linear equation

「一次方程式」。一次函數的等式。

3 x + 4 y - 6 = 2 x + 5 z - 1

移項整理：變數在左邊，常數在右邊。

1 x + 2 y - 5 z = 5

system of linear equations

「一次方程組」。許多道一次方程式同時成立。

⎧ 1 x + 2 y - 5 z = 5
⎨ 2 x + 4 y + 6 z = 1
⎩ 3 x + 1 y + 7 z = 4

一次方程組求解，等同於線性函數求解，或者說矩陣求解。

⎧ 1 x + 2 y - 5 z = 5        ⎡ 1 2 -5 ⎤ ⎡ x ⎤   ⎡ 5 ⎤
⎨ 2 x + 4 y + 6 z = 1  --->  ⎢ 2 4  6 ⎥ ⎢ y ⎥ = ⎢ 1 ⎥
⎩ 3 x + 1 y + 7 z = 4        ⎣ 3 1  7 ⎦ ⎣ z ⎦   ⎣ 4 ⎦
                                 A        x       b

特殊矩陣擁有特殊演算法，稍微省時，詳情請見MATLAB的mldivide說明圖片。簡單起見，以下不介紹特殊矩陣的演算法。

inverse

「反矩陣」。一旦知道反矩陣，即可輕鬆解一次方程組。而且還可以隨意抽換等號右邊數值。

solve Ax = b  --->  find A⁻¹, then x = A⁻¹b

一次方程組求解演算法，稍作修改，即得反矩陣演算法。

等量公理

相信大家已經學過如何解一次方程組：由上往下消除變數，變成階梯狀；由下往上解出變數，變成對角線。

由上往下消除變數的過程，就叫做「高斯消去法」。

演算法（Gaussian elimination）

一個矩陣，化成上三角矩陣，對角線元素皆為一。

由上往下，處理每個橫條，實施三種運算：交換、倍率、相減。

如果橫條的首項係數不是零，以該橫條抵銷下方橫條，令下方橫條的首項係數化成零。如果是零，就往下找到非零橫條，先交換，再抵銷。

pivot row：首項係數不是零的橫條。
pivot：上述橫條的首項係數。
pivoting：找到pivot row，視情況進行交換。

減少誤差的方法：取絕對值最大的pivot row，抵銷其餘row。換句話說，首項係數絕對值最大的橫條，總是交換到上方，再來抵銷下方橫條。倍率小於1，相減誤差少。

矩陣邊長N×M，時間複雜度O(N²M)。大家習慣討論方陣，N = M，時間複雜度O(N³)。

方陣的高斯消去法，程式碼如下所示。矩陣的高斯消去法，留給大家自行練習。

double matrix[10][10]; void gaussian_elimination() { for (int i=0; i<10; ++i) { // 如果橫條的首項係數為零，嘗試與下方橫條交換。 if (matrix[i][i] == 0) for (int j=i+1; j<10; ++j) if (matrix[j][i] != 0) { // 交換上方橫條與下方橫條。 for (int k=i; k<10; ++k) swap(matrix[i][k], matrix[j][k]); break; } // 如果首項係數都是零，那就略過。 if (matrix[i][i] == 0) continue; // 橫條的首項係數調整成一。為了讓對角線皆為一。 double t = matrix[i][i]; // 首項係數 for (int k=i; k<10; ++k) matrix[i][k] /= t; // 抵銷下方橫條，令下方橫條的首項係數化成零。 for (int j=i+1; j<10; ++j) if (matrix[j][i] != 0) { double t = matrix[j][i]; for (int k=i; k<10; ++k) matrix[j][k] -= matrix[i][k] * t; } } }

解一次方程組

首先實施高斯消去法，求得上三角矩陣。

由下往上，處理每個橫條。把先前解出的變數，代入到目前橫條，解出變數。時間複雜度O(N²)。

double matrix[10][10+1]; // 參數化的矩陣 double x[10]; // 一次方程組的解 void back_substitution() { for (int i=10-1; i>=0; --i) { if (matrix[i][i] == 0) return; // 無解、多解 double t = 0; for (int k=i+1; k<10; ++k) t += matrix[i][k] * x[k]; x[i] = (matrix[i][10] - t) / matrix[i][i]; } }

UVa 10109 10524 10828 ICPC 3563

演算法（Gauss–Jordan elimination）

「高斯喬登消去法」是延伸版本。對角線化成一、其餘化成零。

時間複雜度仍是O(N³)。

高斯喬登消去法也可以解一次方程組，但是步驟數量比較多。主要用途是求反矩陣。

求反矩陣

利用「高斯喬登消去法」。

double matrix[10][10*2]; // 參數化的矩陣 bool gauss_jordan_elimination() { // 填好參數化的部分 for (int i=0; i<10; ++i) for (int j=0; j<10; ++j) matrix[i][10+j] = 0; for (int i=0; i<10; ++i) matrix[i][10+i] = 1; // 開始進行高斯喬登消去法 // 內容幾乎與高斯消去法相同 for (int i=0; i<10; ++i) { if (matrix[i][i] == 0) for (int j=i+1; j<10; ++j) if (matrix[j][i] != 0) { for (int k=i; k<10*2; ++k) swap(matrix[i][k], matrix[j][k]); break; } // 反矩陣不存在。 if (matrix[i][i] == 0) return false; double t = matrix[i][i]; for (int k=i; k<10*2; ++k) matrix[i][k] /= t; // 消去時，所有橫條通通消去。 for (int j=0; j<10; ++j) if (i != j && matrix[j][i] != 0) { double t = matrix[j][i]; for (int k=i; k<10*2; ++k) matrix[j][k] -= matrix[i][k] * t; } } return true; }

求determinant

利用「高斯消去法」，對角線不化成一，保留原有數字。

上三角矩陣，對角線元素的乘積，便是determinant。

如果矩陣裡都是整數，那麼determinant也是整數。想避免浮點數誤差，可以使用輾轉相除法進行消去。時間複雜度O(N³logC)，C是絕對值最大的首項係數。

int matrix[10][10]; void rowgcd(int* a, int* b, int i) { if (abs(a[i]) < abs(b[i])) swap(a, b); while (b[i] != 0) { int t = a[i] / b[i]; for (int k=i; k<10; ++k) a[k] -= b[k] * t; swap(a, b); } } int determinant() { int det = 1; for (int i=0; i<10; ++i) { for (int j=i+1; j<10; ++j) { // 輾轉相除法進行消去。 rowgcd(matrix[i], matrix[j], i); // 把非零橫條，挪到當前橫條。 if (matrix[i][i] == 0 && matrix[j][i] != 0) { for (int k=i; k<10; ++k) swap(matrix[i][k], matrix[j][k]); // 兩行交換，determinant會相差一負號。 det *= -1; } } // determinant等於上三角方陣的對角線乘積。 det *= matrix[i][i]; if (det == 0) break; } return det; }

UVa 684

LU decomposition

高斯消去法可以改寫成LUP分解！時間複雜度O(N³)。

交換橫條、抵銷橫條，可以改寫成矩陣乘法。一連串交換橫條、抵銷橫條，可以整併成三個矩陣連乘：下三角矩陣L、上三角矩陣U、列交換矩陣P。稱作「LUP分解」。

如果恰好都沒有交換橫條，則可忽略P。稱作「LU分解」。

LU分解的用途是解大量一次方程組Ax = b，A固定，b有許多組。這種情況下，LU分解，每次求解需時O(N²)；高斯消去法，每次求解需時O(N³)。

Cholesky decomposition

「對稱正定矩陣symmetric positive definite matrix」：特徵值均為正數。

對稱正定矩陣的LU分解。L與U將互相對稱。

時間複雜度仍是O(N³)，但是步驟數量較少。

移項法則

所有等式一齊求反矩陣，繁文縟節，慢慢吞吞。

每道等式一一求反函數，化整為零，簡單明快。

演算法（Jacobi iteration）

每回合依序計算x y z，一次處理一種變數。可以視作不動點遞推法的進化版本。

二維：

⎡ 4  3 ⎤ ⎡ x ⎤ = ⎡ 1 ⎤
⎣ 2  5 ⎦ ⎣ y ⎦   ⎣ 2 ⎦

⎰ 4x + 3y = 1  => ⎰ x = (1 - 3y) / 4
⎱ 2x + 5y = 2     ⎱ y = (2 - 2x) / 5

⎡ x₀ ⎤ = ⎡ 0 ⎤ 隨便設定一解，作為初始值。
⎣ y₀ ⎦   ⎣ 0 ⎦
    
⎡ x₁ ⎤ = ⎡ (1 - 3y₀) / 4 ⎤
⎣ y₁ ⎦   ⎣ (2 - 2x₀) / 5 ⎦
    
⎡ x₂ ⎤ = ⎡ (1 - 3y₁) / 4 ⎤
⎣ y₂ ⎦   ⎣ (2 - 2x₁) / 5 ⎦

三維：

⎡  4  3 -1 ⎤ ⎡ x ⎤   ⎡ 1 ⎤
⎢  2  5  1 ⎥ ⎢ y ⎥ = ⎢ 2 ⎥
⎣ -2 -2  6 ⎦ ⎣ z ⎦   ⎣ 3 ⎦

⎧  4x + 3y -  z = 1     ⎧ x = (1 - 3y +  z) / 4
⎨  2x + 5y +  z = 2  => ⎨ y = (2 - 2x -  z) / 5
⎩ -2x - 2y + 6z = 3     ⎩ z = (3 + 2x + 2y) / 6

⎡ x₀ ⎤   ⎡ 0 ⎤
⎢ y₀ ⎥ = ⎢ 0 ⎥ 隨便設定一解，作為初始值。
⎣ z₀ ⎦   ⎣ 0 ⎦

⎡ x₁ ⎤   ⎡ (1 - 3y₀ +  z₀) / 4 ⎤
⎢ y₁ ⎥ = ⎢ (2 - 2x₀ -  z₀) / 5 ⎥
⎣ z₁ ⎦   ⎣ (3 + 2x₀ + 2y₀) / 6 ⎦

任意維度：

      Ax = b
(D+L+U)x = b                  D是對角線、L是嚴格下三角、U是嚴格上三角
      Dx = b - (L+U)x
       x = D⁻¹ [b - (L+U)x]
       x = D⁻¹b - D⁻¹(L+U)x

時間複雜度O(N²T)，N是方陣維度，T是遞推次數。

高斯消去法直接得到正解。鬆弛法逐步逼近正解。

當b = 0，則完全等價於D⁻¹(L+U)實施不動點遞推法。

判斷收斂：檢查D⁻¹(L+U)的特徵值的絕對值是否都小於1。

「嚴格對角優勢矩陣strictly diagonally dominant matrix」：每個橫條，對角線元素大於其餘元素總和。

SDD矩陣，保證收斂；不滿足時，可能收斂、也可能不收斂。

for each row, |Aii| > ∑ |Aij|
                     j≠i

演算法（Gauss–Seidel iteration）

每回合依序計算x y z，剛出爐的數字，馬上拿來使用，加快收斂速度。

⎡  4  3 -1 ⎤ ⎡ x ⎤   ⎡ 1 ⎤
⎢  2  5  1 ⎥ ⎢ y ⎥ = ⎢ 2 ⎥
⎣ -2 -2  6 ⎦ ⎣ z ⎦   ⎣ 3 ⎦

⎧  4x + 3y -  z = 1     ⎧ x = (1 - 3y +  z) / 4
⎨  2x + 5y +  z = 2  => ⎨ y = (2 - 2x -  z) / 5
⎩ -2x - 2y + 6z = 3     ⎩ z = (3 + 2x + 2y) / 6

⎡ x₀ ⎤   ⎡ 0 ⎤
⎢ y₀ ⎥ = ⎢ 0 ⎥ 隨便設定一解，作為初始值。
⎣ z₀ ⎦   ⎣ 0 ⎦

⎡ x₁ ⎤   ⎡ (1 - 3y₀ +  z₀) / 4 ⎤ 依序計算x₁、y₁、z₁，
⎢ y₁ ⎥ = ⎢ (2 - 2x₁ -  z₀) / 5 ⎥ 剛出爐的數字，馬上拿來使用，
⎣ z₁ ⎦   ⎣ (3 + 2x₁ + 2y₁) / 6 ⎦ 加快收斂速度。

xₖ₊₁ = D⁻¹ (b - Uxₖ - Lxₖ₊₁)

演算法（successive over-relaxation）

原數值、新數值，以固定比例混合。

可令原數值權重小於零、新數值權重大於一，硬是催出收斂速度。

⎡ x₁ ⎤   ⎡ (1-w) ⋅ x₀ + w ⋅ (1 - 3y₀ +  z₀) / 4 ⎤
⎢ y₁ ⎥ = ⎢ (1-w) ⋅ y₀ + w ⋅ (2 - 2x₁ -  z₀) / 5 ⎥
⎣ z₁ ⎦   ⎣ (1-w) ⋅ z₀ + w ⋅ (3 + 2x₁ + 2y₁) / 6 ⎦

linear equation與geometry

一次方程式得視作幾何元件：點、線、面、……。

ContourPlot3D[1 x + 2 y - 5 z == 5, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, Boxed -> False, Axes -> False, Mesh -> None]

一次方程組的解得視作一堆幾何元件的交集。

f := 1 x + 2 y - 5 z - 5; g := 2 x + 4 y + 6 z - 1; h := 3 x + 1 y + 7 z - 4; ContourPlot3D[{f == 0, g == 0, h == 0}, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, Boxed -> False, Axes -> False, Mesh -> None, ContourStyle -> Directive[Opacity[0.5]]]

演算法（Kaczmarz's method）

隨便設定一解，依序投影到第一道等式、第二道等式、……，不斷循環，逐步逼近正解。

時間複雜度O(N²T)，N是方陣維度，T是遞推次數。

無解時，最後形成無窮迴圈，稱作「極限環limit cycle」。

數學公式（Cramer's rule）

兩線交點的演算法：求得平行四邊形的面積，以面積比例求得交點位置。請見本站文件「intersection」。

「克拉瑪公式」則是此演算法的高維度版本，形成了非常漂亮的公式解！求得超平行體的容積，以容積比例求得解。

linear equation:
  Ax = b

solution:
  ⎧ x = det(Aₓ) / det(A)
  ⎨ y = det(A₝) / det(A)
  ⎩ z = det(A₞) / det(A)

unisolvence:
   Ax = b has a unique solution   iff   det(A) ≠ 0

example:
  ⎧ 1 x + 2 y - 5 z = 5        ⎡ 1 2 -5 ⎤ ⎡ x ⎤   ⎡ 5 ⎤
  ⎨ 2 x + 4 y + 6 z = 1  --->  ⎢ 2 4  6 ⎥ ⎢ y ⎥ = ⎢ 1 ⎥
  ⎩ 3 x + 1 y + 7 z = 4        ⎣ 3 1  7 ⎦ ⎣ z ⎦   ⎣ 4 ⎦
                                   A        x       b
         _ b                  _ b                  _ b
       ⎡|5| 2 -5 ⎤       ⎡ 1 |5|-5 ⎤       ⎡ 1  2 |5|⎤
  Aₓ = ⎢|1| 4  6 ⎥  A₝ = ⎢ 2 |1| 6 ⎥  A₞ = ⎢ 2  4 |1|⎥
       ⎣|4| 1  7 ⎦       ⎣ 3 |4| 7 ⎦       ⎣ 3  1 |4|⎦
         ‾                    ‾                    ‾

determinant是矩陣當中所有向量所構成的超平行體的容積。時間複雜度等於N+1次determinant的時間複雜度，O(N⁴)。

determinant

determinant起初用來判定一個一次方程組是否有解、解是多少，因而稱作「決定因子」。古人沒有意識到determinant是容積。

字面意義是「決定因子」，中文教科書卻譯作「行列式」。真是異想天開的翻譯啊！

行列式的計算過程是：先刪除一橫行，接著分別刪除每一直行，形成N-1個(N-1)×(N-1)子矩陣，添上正負號。原矩陣的行列式，等於這些子矩陣的行列式總和。每個子矩陣各自遞迴下去，直到N=1。1×1矩陣的行列式，等於矩陣元素。時間複雜度O(N!)。

N = 2 or 3的時候比較特別，可以直接累加所有「左上右下斜線」的乘積、累減「右上左下斜線」的乘積。中學數學課程有教。

計算行列式，也可以使用高斯消去法，時間複雜度O(N³)。

與其採用高斯消去法求行列式、再用行列式解一次方程組，不如直接採用高斯消去法解一次方程組。就當作是學個想法吧。

linear equation與algebra

一次方程組可以寫成矩陣乘法的形式。

⎧ 1 x + 2 y - 5 z = 5        ⎡ 1 2 -5 ⎤ ⎡ x ⎤   ⎡ 5 ⎤
⎨ 2 x + 4 y + 6 z = 1  --->  ⎢ 2 4  6 ⎥ ⎢ y ⎥ = ⎢ 1 ⎥
⎩ 3 x + 1 y + 7 z = 4        ⎣ 3 1  7 ⎦ ⎣ z ⎦   ⎣ 4 ⎦
                                 A        x       b

一次方程組得視作線性函數，解一次方程組得視作線性反函數。

                                                 -1
⎡ 1 2 -5 ⎤ ⎡ x ⎤   ⎡ 5 ⎤        ⎡ x ⎤   ⎡ 1 2 -5 ⎤ ⎡ 5 ⎤
⎢ 2 4  6 ⎥ ⎢ y ⎥ = ⎢ 1 ⎥  --->  ⎢ y ⎥ = ⎢ 2 4  6 ⎥ ⎢ 1 ⎥
⎣ 3 1  7 ⎦ ⎣ z ⎦   ⎣ 4 ⎦        ⎣ z ⎦   ⎣ 3 1  7 ⎦ ⎣ 4 ⎦
    A        x       b            x         A⁻¹      b

數學公式（inverse matrix）

一旦求得反矩陣，即可輕鬆解一次方程組。

solve Ax = b  --->  find A⁻¹, then x = A⁻¹b

計算反矩陣，使用高斯喬登消去法，時間複雜度O(N³)。

與其採用高斯喬登消去法求反矩陣、再用反矩陣解一次方程組，不如直接採用高斯消去法解一次方程組。就當作是學個想法吧。

數學公式（eigendecomposition）

一旦求得特徵分解，即可輕鬆解一次方程組。

A   = EΛE⁻¹
A⁻¹ = EΛ⁻¹E⁻¹
x = A⁻¹b = EΛ⁻¹E⁻¹b

與其採用特徵分解求反矩陣、再用反矩陣解一次方程組，不如直接採用高斯消去法解一次方程組。就當作是學個想法吧。

演算法（Newton's method）

Ax = b的解，即是f(x) = Ax - b的根。

牛頓法。x = x - Aᵀ⁻¹(Ax - b)。

光是求反矩陣就需時O(N³)。無人使用。

與其求反矩陣、不斷遞推，不如直接採用高斯消去法。就當作是學個想法吧。

演算法（Richardson iteration）

Ax = b的解，即是g(x) = Ax - b + x的不動點。

不動點遞推法。可以添加權重α。x = x + α(Ax - b)。

時間複雜度O(N²T)，N是方陣維度，T是遞推次數。

與其求餘數、不斷遞推，不如直接採用Jacobi iteration。就當作是學個想法吧。

「一次函數求解」等於「二次函數求極值」

A必須是對稱正定矩陣，使得二次函數有唯一最小值。

     solve Ax = b               (if A has unique solution)
---> solve Ax - b = 0           (transposition)
---> min 1/2 xᵀAx - bᵀx         (quadratic optimization)

let f(x)  = 1/2 xᵀAx - bᵀx
let f′(x) = Ax - b

     solve Ax = b               (if A has no solution)
---> min ‖Ax - b‖²              (least squares)
---> min xᵀAᵀAx - 2bᵀAx + bᵀb   (expansion)
---> min xᵀAᵀAx - 2bᵀAx         (quadratic optimization)

let f(x)  = xᵀAᵀAx - 2bᵀAx
let f′(x) = Aᵀ(Ax - b)

f(x)的最小值位置，即是f′(x)的根，即是Ax = b的解。

f(x)的最佳化演算法是「共軛梯度法」。優於高斯消去法。

時間複雜度O(N²T)，N是方陣維度，T是遞推次數。

「矩陣求解」化作「對稱正定矩陣求解」

當A不是對稱正定矩陣，那麼等號兩邊同時乘上Aᵀ，得到對稱半正定矩陣AᵀA。當A有唯一解，則得到對稱正定矩陣AᵀA。

     solve Ax = b       (if A has unique solution)
---> solve AᵀAx = Aᵀb   (AᵀA is symmetric positive definite)
---> solve A'x = b'     (unique solution)

     solve Ax = b       (if A has no solution)
---> min ‖Ax - b‖²      (least squares)
---> solve AᵀAx = Aᵀb   (AᵀA is symmetric positive semidefinite)
                        (normal equation)
---> solve A'x = b'     (least squares solution)

與其變成對稱正定矩陣、再用共軛梯度法解一次方程組，不如直接採用高斯消去法解一次方程組。就當作是學個想法吧。

linear least squares

唯一解是稀奇的，無解、多解是普遍的。無解、多解時，改為找到平方誤差最小的解。

no solution：無解。改求‖Ax - b‖²最小的解。
infinitely many solutions：多解。改求‖x‖²最小的解。
unique solution：唯一解。求Ax = b的解。

min ‖Ax - b‖²：方程組每一道等式，求得等號左右兩邊的差的平方；累計所有等式，總和越小越好。

min ‖x‖²：解的每一項的平方，總和越小越好。

no solution / infinitely many solutions

請見「Rouché-Capelli theorem」。

矩陣的rank、column，得以區分無解、多解。

rank：消去冗餘的等式。column：直條數量、變數數量。

no solution              : rank(A) ≠ rank([A|b])
infinitely many solutions: rank(A) = rank([A|b]) ≠ column(A)
unique solution          : rank(A) = rank([A|b]) = column(A)

overdetermined system / underdetermined system

⎡ X X ⎤               ⎡ X X X ⎤
⎢ X X ⎥   ⎡ X X X ⎤   ⎢ X X X ⎥
⎣ X X ⎦   ⎣ X X X ⎦   ⎣ X X X ⎦
 thin        fat        square

overdetermined system  (thin matrix)  : row(A) > column(A)
underdetermined system (fat matrix)   : row(A) < column(A)
welldetermined system  (square matrix): row(A) = column(A)

超定方程組（瘦矩陣）：等式太多。等式多於變數。
欠定方程組（胖矩陣）：等式太少。等式少於變數。
良定方程組（方陣）　：等式剛好。等式等於變數。

等式太多太少剛好，答案無解多解唯一解，兩邊沒有直接關係。

方便起見，許多人修改定義，兩邊直接一一對應。包括下文。

一般而言，等式太多，很可能無解；等式太少，很可能多解。

pseudoinverse（Moore–Penrose inverse）

「虛擬反矩陣」、「偽反矩陣」。反矩陣從方陣推廣成矩陣。

solve Ax = b  --->  find A⁺, then x = A⁺b

三種數學公式

solve overdetermined system Ax = b  --->  min ‖Ax - b‖²

x = (Aᵀ A)⁻¹ Aᵀ b     ( Aᵀ A x = Aᵀ b )   normal equation

x = R⁻¹ Q b           ( A = Q R )         QR decomposition

x = V Σ⁻¹ Uᵀ b        ( A = U Σ Vᵀ )      singular value decomposition

solve underdetermined system Ax = b  --->  min ‖x‖²

x = Aᵀ (A Aᵀ)⁻¹ b                         normal equation

x = Q (Rᵀ)⁻¹ b        ( Aᵀ = Q R )        QR decomposition

x = U Σ⁻¹ Vᵀ b        ( Aᵀ = U Σ Vᵀ )     singular value decomposition

https://eigen.tuxfamily.org/dox-devel/group__LeastSquares.html

數學公式（normal equation）

線性代數經典公式！視作最佳化問題，以微分求極值。

「一次微分等於零」的地方是極值、鞍點。因為平方誤差是開口向上的拋物面，所以「一次微分等於零」的地方必是最小值。

以下只證明無解的情況。時間複雜度O(N³)。

solve Ax = b
min ‖Ax - b‖²
d/dx ‖Ax - b‖² = 0                    「一次微分等於零」的地方是極值、鞍點
[ d/dx (Ax - b) ] [ 2(Ax - b) ] = 0   微分連鎖律
Aᵀ [ 2(Ax - b) ] = 0                  微分
AᵀAx = Aᵀb                            同除以2、展開、移項
x = (AᵀA)⁻¹Aᵀb                        移項。注意到A的向量們必須線性獨立！

注意到最後一步，A的向量們必須線性獨立（事先清除冗餘的、無意義的變數），AᵀA才有反矩陣。

結果宛如向量投影公式：b投影到A。

                 Aᵀb   dot(A,b)
x = (AᵀA)⁻¹Aᵀb = ——— = ———————— = projAb
                 AᵀA   dot(A,A)

平方誤差盡量小＝垂直投影！分析與幾何兩大領域打通了！

     solve Ax = b
---> min ‖Ax - b‖²
---> 2Aᵀ(Ax - b) = 0
---> x = projAb

數學公式（QR decompostion）

以下只證明無解的情況。A = QR，正規正交矩陣Q不影響最小值，最小值取決於零餘部分R。

至於多解的情況，改為分解A的轉置矩陣Aᵀ = QR。

時間複雜度O(N³)。但是計算量比normal equation少。

solve Ax = b
min ‖Ax - b‖²
min ‖Qᵀ(Ax - b)‖²           正規正交矩陣，變換後長度不變
min ‖QᵀAx - Qᵀb‖²           展開
min ‖Rx - Qᵀb‖²             QᵀA = QᵀQR = R
min ‖ ⎡ R₁x - Q₁ᵀb ⎤ ‖²     R = ⎡ R₁ ⎤  區分出零，讓R₁是方陣
    ‖ ⎣   0 - Q₂ᵀb ⎦ ‖          ⎣ 0  ⎦  區分上段和下段
solve R₁x - Q₁ᵀb = 0        此式有唯一解，可為零（因此最小值是 ‖Q₂ᵀb‖²）
R₁x = Q₁ᵀb                  移項
x = R₁⁻¹Q₁ᵀb                移項

數學公式（singular value decompostion）

以下只證明無解的情況。證明手法如出一轍。A = UΣVᵀ。

至於多解的情況，改為分解A的轉置矩陣Aᵀ = UΣVᵀ。

時間複雜度O(N³)。我不確定實務上是否比較快。

solve Ax = b
min ‖Ax - b‖²
min ‖Uᵀ(Ax - b)‖²        正規正交矩陣，變換後長度不變
min ‖UᵀAx - Uᵀb‖²        展開
min ‖ΣVᵀx - Uᵀb‖²        UᵀA = UᵀUΣVᵀ = ΣVᵀ
solve ΣVᵀx - Uᵀb = 0     此式有唯一解，可為零
ΣVᵀx = Uᵀb               移項
x = VΣ⁻¹Uᵀb              移項

演算法（gradient descent）（Richardson iteration）

梯度下降法＝不動點遞推法。

f(x)  = ‖Ax - b‖²
f′(x) = Aᵀ(Ax - b)

xnext = x + α f′(x)     α is step size

自訂步伐大小α。有人設定成α = 2 / σ_max²，其中σ_max是最大的奇異值。

演算法（steepest descent）

最陡下降法＝梯度下降法＆直線搜尋。

xnext = x + α f′(x)
α = argmin f(x + α f′(x)) = ‖f′(x)‖² / ‖A f′(x)‖²

步伐大小是‖f′(x)‖² / ‖A f′(x)‖²。證明省略。

Tikhonov regularization

一次方程組，無解、多解時，改為最佳化問題，以得到唯一解。

solve Ax = b   --->   min ‖Ax - b‖²   [overdetermined system]
solve Ax = b   --->   min ‖x‖²        [underdetermined system]

甚至利用regularization，追加其他最佳化目標。

solve Ax = b   --->   min ‖Ax - b‖² + α f(x)   (α ≥ 0)
solve Ax = b   --->   min ‖x‖² + α f(x)   (α ≥ 0)

一次方程組，有許多式子和變數。可能有其中一群變數與式子構成無解、另一群構成唯一解、剩下一群構成多解。更有甚者，切割一些群結果無解變多解、整併某些群結果多解變無解。

無法釐清是無解、多解的時候，那就兩個一起上吧。

solve Ax = b   --->   min ‖Ax - b‖² + α ‖x‖²   (α ≥ 0)

d/dx [ ‖Ax - b‖² + α ‖x‖² ] = 0   「一次微分等於零」的地方是極值、鞍點
                                  二次函數、恆正，必得最小值
2 AᵀA x - 2 Aᵀ b + 2 α x = 0      展開

( AᵀA + α I ) x = Aᵀ b            移項，左式即是 AᵀA 的對角線加上 α

左式是實數對稱正定矩陣，有唯一解。時間複雜度O(N³)。

homogeneous linear equation

討論特例b = 0的情況。顯然答案包含x = 0。

可能多解、唯一解。不可能無解。至少有一解x = 0。

solve Ax = 0

homogeneous linear least squares

討論特例b = 0的情況。當b = 0，則x = 0，缺乏討論意義。於是添加限制「x長度（的平方）為1」，增進討論意義。

如果原本是多解，那麼添加限制可以找到x = 0以外的解。如果原本是唯一解x = 0，那麼添加限制可以找到原本不存在的解。

solve Ax = 0   --->   min ‖Ax‖² subject to ‖x‖² = 1

min ‖Ax‖² - λ ( ‖x‖² - 1 )             Lagrange multiplier

∂/∂x [ ‖Ax‖² - λ ( ‖x‖² - 1 ) ] = 0   「一次微分等於零」的地方是極值、鞍點
                                       二次函數，必得極值
2 AᵀA x - 2 λ x = 0                    展開

AᵀA x = λ x                            移項，此即特徵向量的格式

答案是AᵀA的最小的特徵值的特徵向量！

AᵀA是對稱半正定矩陣，特徵值即是奇異值，特徵向量即是奇異向量。

AᵀA的特徵值即是A的奇異值，AᵀA的特徵向量即是A的右奇異向量。