Linear Equation

Linear Function

「一次函數」。多項式函數,變數不相乘,頂多一次方。

f(x, y, z) = 1 x + 2 y - 5 z - 5

數字的一次方呈直線line、二次方呈正方形quadrate、三次方呈正方體cube。因此一次取名linear、二次取名quadratic、三次取名cubic。

Linear Function

「線性函數」。僅由變數的加減、變數的倍率所構成的函數。

f(x, y, z) = 1 x + 2 y - 5 z      線性函數
f(x, y, z) = 1 x + 2 y - 5 z - 5  仿射函數  兩者都是一次函數

線性代數開始流行之後,linear新添了其他意義。當數學家說linear,可能是一次,也可能是線性。

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⃗       y⃗

儘管標題是Linear Equation,但是接下來都是談矩陣求解。

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

Linear Equation: Gaussian Elimination

等量公理

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

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

演算法(Gaussian Elimination)

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

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

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

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

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

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

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

解一次方程組

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

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

UVa 10109 10524 10828 ICPC 3563

LU Decomposition

http://ccjou.wordpress.com/2010/09/01/

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

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

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

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

Cholesky Decomposition

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

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

演算法(Gauss-Jordan Elimination)

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

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

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

求反矩陣

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

求determinant

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

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

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

UVa 684

Linear Equation: Relaxation

移項法則

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

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

演算法(Jacobi Method)

每回合依序計算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₀ = 隨便設定一個初始值
xₖ₊₁ = D⁻¹ [b - (L+U)xₖ]

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

高斯消去法直接得到正解。不動點遞推法逐步逼近正解。

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

x = D⁻¹ [b - (L+U)x]
x = D⁻¹b - D⁻¹(L+U)x

滿足strictly diagonally dominant就保證收斂。不滿足時,可能收斂、也可能不收斂。

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

演算法(Gauss-Seidel Method)

每回合依序計算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 ]

Matrix Splitting

「矩陣分裂」。A = A₁ + A₂ + ...,原本矩陣分裂成多個矩陣相加,使得A₁ A₂ ...清爽暢快,容易計算。

移項法則,換個角度來思考,即是A = D + L + U。

Preconditioner

「預條件」。Ax = b變成MAx = Mb,等號兩邊同乘矩陣M,使得MA清爽暢快,容易計算。

M必須是可逆矩陣(擁有反函數),讓解的數量不變。

移項法則,從某種程度上而言,即是等號兩邊同乘矩陣D⁻¹。

Linear Equation: Optimization

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

等號兩邊同時乘上Aᵀ,得到對稱半正定矩陣AᵀA。

當A有唯一解,則AᵀA從半正定升級成正定。

     solve Ax = b       (A is invertible: unique solution)
---> solve AᵀAx = Aᵀb   (AᵀA is symmetric positive definite)
---> solve A'x = b'     (unique solution)
     solve Ax = b       (overdetermined system: no solution)
---> min ‖Ax - b‖²      (least squares)
---> solve AᵀAx = Aᵀb   (normal equation)
---> solve A'x = b'     (least squares solution)

「一次函數求解」化作「二次函數求極值」

A是對稱正定矩陣。若A是普通矩陣,則下述性質未必成立。

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

xᵀAx是橢圓拋物面函數。等高線是同心橢圓。有唯一最小值。

f(x)是橢圓拋物面函數減去平面函數,仍是橢圓拋物面函數。

f(x)的最小值位於一次微分等於零的地方f′(x) = Ax - b = 0。

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

f(x)的最佳化演算法是「共軛梯度法」。

Linear Equation: Eigendecomposition

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⃗       y⃗

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

                                                 -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⃗       y⃗           x⃗         A⁻¹      y⃗

Inverse

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

solve Ax⃗ = y⃗  ===>  find A⁻¹, then x⃗ = A⁻¹y⃗

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

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

數學公式(Eigendecomposition)

求根、求不動點、求特徵點、求解是等價的,使得矩陣求解有著各式各樣的演算法。

線性函數的精髓:求得在特徵向量上的分量,分別伸縮,伸縮倍率是特徵值。線性反函數的精髓:逆向縮放,縮放倍率變成倒數。

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

矩陣實施特徵分解,求反矩陣,再拿來解一次方程組。

以特徵值來判斷是否存在反矩陣。特徵值皆不為零,則有反矩陣。就這樣子。

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

Linear Equation: Cramer's Rule

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]]]

數學公式(Cramer's Rule)

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

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

linear equation:
  Ax⃗ = y⃗

solution:
  { x = det(Aˣ) / det(A)
  { y = det(Aʸ) / det(A)
  { z = det(Aᶻ) / det(A)

unisolvence:
   Ax⃗ = y⃗ 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⃗       y⃗
         _ y⃗                  _ y⃗                  _ y⃗
       [|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是容積。

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

http://mathworld.wolfram.com/DeterminantExpansionbyMinors.html

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

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

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

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

Linear Least Squares

Linear Least Squares(Linear Least Squares Equation)

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

solve Ax = b

 overdetermined system: 等式太多 -> 無解 -> 改求‖Ax - b‖²最小的解
underdetermined system: 等式太少 -> 多解 -> 改求‖x‖²最小的解

嚴格來說,必須預先消去所有等價的、多餘的等式,以rank大小、矩陣邊長來區分這兩種情況。等式太多太少,無解多解唯一解,其實沒有直接關係。詳情請參考「Rouché-Capelli Theorem」。

等式太多、等式太少(中學數學的講法是:變數少於等號、變數多於等號)所導致的無解、多解,兩種情況分別處理。

一、無解:方程組每一道等式,求得等號左右兩邊的差的平方;累計所有等式,總和越小越好。

二、多解:解的每一項的平方,總和越小越好。

Least意指「盡量小」,Squares意指「平方和」。

平方誤差的優勢是:循規蹈矩,成為一個參考指標,誤差高低可以拿來比較,科學多了。缺陷是:計算速度慢。

平方誤差比起絕對值誤差,有三個好處:一、讓每道等式的誤差保持均勻,不會有某道等式誤差特別高。二、「一次微分等於零」容易推導數學公式。三、誤差函數是開口向上的拋物線函數,沒有鞍點,容易最佳化。

Linear Least Squares: Decomposition

三種數學公式

函式庫兜一兜,答案就出來了,大可不必深究細節。

solve overdetermined system Ax = b   minimize ‖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   minimize ‖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

Timus 1668

數學公式(Normal Equation)

http://people.csail.mit.edu/bkph/articles/Pseudo_Inverse.pdf

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

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

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

solve Ax = b

minimize ‖Ax - b‖²

∂/∂x ‖Ax - b‖² = 0                    「一次微分等於零」的地方是最小值

[ ∂/∂x (Ax - b) ] [ 2(Ax - b) ] = 0   微分連鎖律

Aᵀ [ 2(Ax - b) ] = 0                  微分

Aᵀ A x = 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)

http://www.cs.cornell.edu/courses/cs322/2007sp/notes/qr.pdf

A = QR。將矩陣拆開成正規正交矩陣Q、零餘部分R。正規正交矩陣Q不影響最小值,最小值取決於零餘部分R。

以下只證明無解的情況。多解的情況,改為分解A的轉置矩陣Aᵀ = QR。

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

solve Ax = b

min ‖ Ax - b ‖²

min ‖ Qᵀ (Ax - b) ‖²          正規正交矩陣,變換後長度不變

min ‖ Qᵀ A x - Qᵀ b ‖²        展開

min ‖ R x - Qᵀ b ‖²           Qᵀ A = Qᵀ Q R = R

min ‖ [ R₁ x - Q₁ᵀ b ] ‖²     R = [ R₁ ]  區分出零,讓R₁是方陣
    ‖ [    0 - Q₂ᵀ b ] ‖          [ 0  ]  區分上段和下段

令 R₁ x - Q₁ᵀ b = 0           此式有唯一解,可為零

令最小值是 ‖ Q₂ᵀ b ‖²

R₁ x = Q₁ᵀ b                  移項

x = R₁⁻¹ Q₁ᵀ b                移項

數學公式(Singular Value Decompostion)

和QR分解的手法如出一轍。

以下只證明無解的情況。多解的情況,改為分解A的轉置矩陣Aᵀ = UΣVᵀ。

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

solve Ax = b

min ‖ Ax - b ‖²

min ‖ Uᵀ (Ax - b) ‖²       正規正交矩陣,變換後長度不變

min ‖ Uᵀ A x - Uᵀ b ‖²     展開

min ‖ Σ Vᵀ x - Uᵀ b ‖²     Uᵀ A = Uᵀ U Σ Vᵀ = Σ Vᵀ

令 Σ Vᵀ x - Uᵀ b = 0       此式有唯一解,可為零

Σ Vᵀ x = Uᵀ b              移項

x = V Σ⁻¹ Uᵀ b             移項

Linear Least Squares: Optimization

Tikhonov Regularization

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

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

甚至利用Regularization,追加其他最佳化目標。

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

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

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

solve Ax = b   --->   minimize ‖Ax - b‖² + α ‖x‖²   (α ≥ 0)
∂/∂x [ ‖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的情況。當b = 0,則x = 0,缺乏討論意義。於是添加限制「x長度(的平方)為1」,增進討論意義。

solve Ax = 0   --->   minimize ‖Ax‖² subject to ‖x‖² = 1
minimize ‖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最小的且大於零的特徵值的特徵向量。

欲求最小的特徵值,可以採用QR Iteration或Lanczos Iteration演算法求得所有特徵值,再挑出最小的,時間複雜度O(N³ + N²K)。亦可採用Singular Value Decomposition的演算法,不必計算AᵀA,節省一點時間。

Basis Pursuit(Lasso)

改成L¹ norm,討論多解的情況。NP-hard。

solve Ax = b   --->   minimize ‖x‖₁ subject to Ax = b
                      [underdetermined system]

解法是改寫成一次規劃:

http://www4.ncsu.edu/~kksivara/masters-thesis/kristen-thesis.pdf

Basis Pursuit Denoising

無解採用L² norm、多解採用L¹ norm。走火入魔。

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

Linear Inequality

System of Linear Inequalities

「一次不等式組」。許多道一次不等式同時成立。

正是計算幾何「Half-plane Intersection」推廣到高維度,所有解形成一個凸多胞形,也可能形成開放區間、退化、空集合。

目前沒有演算法。有人適度乘上負號,調整成Ax > b的格式,套用高斯消去法,但是不知道正不正確。

大家習慣採用一次規劃,將凸多胞形硬是位移至第一象限(各個變數加上一個足夠大的數值,代換成新變數),以符合一次規劃的格式。