線性方程式求解 在

關於 線性方程式求解

還記得你高中時第一次接觸代數的情景嗎？

通常，我們會拿到一組方程組，然後被要求解出x和y。

4x - 3y = -2
2x + 7y = 16

兩個方程式、兩個未知數：稍微做點代入，很快就能看出x = 1, y = 2。

再回頭看看前面例子的左邊。它看起來就像一個矩陣[4 -3; 2 7]乘上未知向量[x, y]，得到向量[-2, 16]。

以線性代數的記法：A x = b，依循大寫字母代表矩陣、小寫代表向量的慣例。

A x = 0與零空間

一個重要的特例是，當所有方程式的右邊都是零。

4x - 3y = 0
2x + 7y = 0

在上面的例子中，唯一的解是x = y = 0，這通常沒什麼意思。

A x = 0的非平凡解，只有在A的行列式為零時才存在。

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

反之，考慮這些方程式：

4x - 3y = 0
8x + 6y = 0

第二個方程式只是第一個的兩倍，並沒有提供新資訊。凡是滿足x = 0.75y的值都會是解。

這些(x, y)值落在二維空間中的一條直線上。解有無限多個，它們構成矩陣的零空間。

在這種情況下，矩陣中的各列並非線性獨立，而矩陣的秩小於列數或行數。

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

nullspace()函式會回傳單一個點，並盡可能選擇正規化為單位向量，但這個向量的任何純量倍數同樣也在零空間中。

求解A x = b

現在假設你有 1000 個方程式、1000 個未知數？如果這聽起來很蠢，請記得無線數位感測器現在既便宜又多元（量測應變、風速、三軸加速度，什麼都行……）。用 1000 個這樣的感測器來監測像懸索橋這樣的現代結構，再合理不過了。這裡需要有控制系統來解讀資料流，並在情況變得令人擔憂時準備好發出警報。

再次強調，A x = b，其中我們有A（來自工程設計）和b（來自感測器的量測值），但需要求出x。

對不熟悉這件事的人來說，第一個念頭往往是想辦法「把等式兩邊同除以A」，好把A移到右邊。

事實是，A有反矩陣（大多數方陣都有，但並非全部：其行列式必須不為零），而這個計算是可行的（但很慢！）。

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

可惜的是，計算反矩陣對小矩陣來說很慢，對大矩陣則慢如冰川。如果你的橋在計算過程中崩塌了，那可不理想！

幸運的是，其他演算法快得多（在這個例子中是高斯消去法）。

Julia（沿用 Matlab 的做法）只用一個反斜線來當求解器（技術上稱為「左除」）。

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

這是個簡單的例子，但有幾個細節要留意。

• 係數必須按行對齊，缺少的項就補上零。
• 方程式（因此矩陣的各列）必須線性獨立。如果第 3 列只是第 1 列和第 2 列的總和，它就沒有提供新資訊，問題就是欠定的。

rank()指令會給出線性獨立的列數／行數，所以目標是讓它等於未知變數的個數。

julia> rank(A)
2

長方形矩陣

大概在你還小的時候，有人教過你：要解出N個未知數，就需要由N個方程式組成的方程組。

一如數學裡大多數的事情，實際情況要微妙一些（不只是前一節提到的線性獨立這個細節）。

「太少」的列

有N-1個方程式（矩陣A中有N-1列）時，問題是欠定的。我們要絕望地放棄嗎？

這要看情況！

N-1個方程式仍然包含很多資訊。用幾何的角度來說，我們需要N個才能把解確定為 N 維空間中的一個_點_，但pseudoinverse個只會告訴我們解必定落在某條_直線_上。

這麼一來，我們會有無限多個解，但同時也有無限大的參數空間是_沒有_解的。

你解所在的這條直線，是否會穿過一塊有信譽的工程師應該擔心的參數空間區域？或者穿過一塊足以構成理由、讓你把結構對外關閉的區域（懸索橋在強烈側風下可是會晃得很厲害）？也許值得再多做些計算來確定這一點！

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

反斜線求解器給了我們一個解！非官方地說，這顯然是範數最小的解（幾何上最接近原點）。

要得到其他解，我們需要A2的nullspace。

把零空間的任何純量倍數加到原本的解上，都會得到另一個有效的解。

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

同樣地，N-2個方程式會把（無限多的）解限制在一個平面上，同樣的原則也適用。

「太多」的列

相反的情況是方程式_比_未知數還多，但它們卻仍然線性獨立。

在真實世界的工程裡，這完全是正常的，而且被視為一件好事！

感測器的精度有限，b向量的精度有限，計算中還有雜訊。這時的解就變成對含雜訊資料的最小平方法擬合。

最簡單的技術是用矩陣的偽逆，Julia 把它實作成pinv()函式。

它可以像（非奇異）方陣的反矩陣那樣使用，給出變數的最小平方估計，而不是精確解。

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

特徵值與特徵向量

前一節談的是A x = b的解，等同於我們熟悉的代數。

本節談的是A x = λ x的解，其中 λ 是純量。

這是線性代數特有的概念，但我們可以從幾何上理解它：

對於適當的x和 λ 值，方陣A會把非零向量x的長度縮放 λ 倍，但不改變它的方向（λ 為負時則反向）。

這聽起來很小眾，但結果卻_荒謬地_有用！

術語： λ 的有效值稱為A的_特徵值_，而對應的x值稱為A的_特徵向量_。很遺憾，我們只能接受有些詞會中途換語言這件事。

學生通常會學到如何用手算 2×2 矩陣的特徵值／特徵向量，但用電腦要輕鬆得多（雖然仍舊相當慢，一般情況下對一個n×n矩陣的複雜度是O(n^3)）。

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

一般來說，一個n×n矩陣會有n個特徵值，不過這些值未必_互異_。可以把它們想成一個n次多項式（稱為特徵多項式）的根，這些根可能重複，而且即使矩陣是實數值，它們也常常是複數。

每個特徵向量代表一個方向，而任何純量倍數同樣也是有效的特徵向量。為了方便後續計算，Julia 會回傳範數為 1 的單位向量。

應用

要解釋特徵向量為什麼重要，理想上需要一本 500 頁的教科書，而不是短短幾段文字。這個主題滲透到現代應用數學的_極多_層面。

從最高的層次來看，特徵向量代表資料集中「最重要」的軸（方向）。特徵值則指出每個軸的「相對重要性」（前提是對輸入的正規化做了一些假設）。

這實際上_代表_什麼，取決於應用。

主成分分析

在典型的資料科學問題裡，我們可能會有 100 個以上的「特徵」，以資料行的形式儲存。幾乎無可避免地，會有雜訊、冗餘，以及不想要的相關性。

我們需要做降維，而 PCA 是把混亂整理出秩序的一種方法：

• 計算資料的共變異數矩陣。
• 現在我們有一個方陣了，接著計算特徵值與特徵向量。
• 依特徵值（更精確地說，是絕對值|λ|）由大到小排序特徵向量。
• 前k個特徵向量現在就是你的principal components，其中k明顯小於原本的特徵數量。
• 把整個資料集投影到這k個軸上，然後開始尋找有趣的模式。

以這種形式，PCA 被各種不同的團體使用，從試圖在藥物設計中改善分子性質的醫學科學家，到試圖從民調資料理解選民偏好的政治競選團隊。大概也有想賣你東西的行銷團隊，但任何技術都可以被用在好的或壞的用途上。

PCA 不是最小安裝的 Julia 的一部分，但（在 Exercism 之外）MultivariateStats 套件裡有你需要的東西。

影像處理

把 PCA 的想法再往前推一步：數位影像不過是像素值的矩陣，我們可以為它們計算主成分。

這在影像壓縮上已經用了幾十年，其中 PCA 是判斷縮減檔案大小時最該保留什麼的指引。

越來越多的情況是，PCA 是影像分類器（例如人臉辨識）的關鍵一環。下次你走過機場時，老大哥不只是在看著你，他還用線性代數來理解他所看到的！

機械工程

比較不具爭議的是，機械元件的主軸是它的慣性矩張量I的特徵向量（這是矩陣換了個稍微不同的名字）。

你車上的每個輪子大概都有一或多個配重塊，用來把這個張量I的非對角元素歸零。修車廠的技師無疑會避開這些數學（與飛機設計師和火箭工程師相反），但「抖動」只是用來描述張量中交叉項的日常用語；在這個例子裡，這些交叉項會讓你的旅程不那麼舒適，並增加軸承的機械磨損。