linear function - 演算法筆記

eigenbasis - root finding

本章主角是特徵向量與特徵值（真實版本）

矩陣求根

本站文件「root finding」，介紹了「輸入、輸出只有一個變數」的函數，以及根、解、不動點、特徵點。

此處介紹「輸入、輸出有很多個變數」的線性函數，以及根、解、不動向量、特徵向量。

根（核）：Ax = 0，符合條件的x。

線性函數必定有一個根：零向量x = 0，缺乏討論意義。

線性函數通常有許多個根，並且組成空間，空間維度等於A的座標軸數量減去A的rank。白話解釋：座標軸無法觸及之處，當然輸出0囉。另外rank、determinant可以判斷是否有許多個根。

解：Ax = b，符合條件的x。b是已知的、特定的向量。

類似於根。請見本站文件「linear equation」。

不動向量：Ax = x，符合條件的x。

線性函數必定有一個不動向量：零向量x = 0。

特徵向量：Ax = λx，符合條件的x。λ是任意一種倍率。

線性函數通常有許多個特徵向量，後續文章的主角。

eigenvector / eigenvalue

針對一個線性函數，從眾多的輸入向量當中，找到其中一些輸入向量，讓輸出向量恰是輸入向量乘上倍率，此時輸入向量稱作「特徵向量」，倍率稱作「特徵值」。

數學式子是Ax = λx。A是線性函數，x是輸入向量，λx是輸出向量。x是特徵向量，λ是特徵值。

特徵向量變換之後，長度伸縮，方向相同或相反。

特徵值可正、可負、可零。正的伸縮，負的伸縮兼翻轉，零的長度歸零。

特徵向量有無限多個

零向量必定是特徵向量，缺乏討論意義，大家習慣無視它。

特徵向量的任何一種倍率（方向相同或相反），也是特徵向量，也是相同特徵值。大家習慣取長度為一的特徵向量當作代表，方向則隨意。

方向不同的特徵向量，如果恰好擁有相同特徵值，那麼這些特徵向量構成的平面、空間、……當中的任何一個向量，也是特徵向量，也是相同特徵值。大家習慣取互相垂直的特徵向量當作代表，方向則隨意。

characteristic equation / characteristic polynomial

無邊無際、無窮無盡，如何找到特徵向量和特徵值呢？

嘗試解Ax = λx。移項Ax - λx = 0。合併(A - λI)x = 0。

若是唯一解，那麼x顯然是零向量。缺乏討論意義。

若有其他解，那麼令det(A - λI) = 0，稱作「特徵方程式」。展開det(A - λI)，形成λ的多項式，稱作「特徵多項式」。

先求特徵值λ：展開det(A - λI) = 0，特徵多項式求根。

再求特徵向量x：特徵值代入(A - λI)x = 0，線性函數求根。

一種特徵值，得到一種特徵向量。

A = ⎡ 1 -1 ⎤   A - λI = ⎡ (1-λ)  -1   ⎤      det(A - λI) = 0
    ⎣ 2  4 ⎦            ⎣   2   (4-λ) ⎦   => (1-λ)(4-λ) + 2 = 0
                                          => λλ - 5λ + 6 = 0
                                          => λ = 2 or 3

when eigenvalue λ = 2        |   when eigenvalue λ = 3        
                             |                                
    (A - λI)      x  =  0    |       (A - λI)      x  =  0    
                             |                                
⎡ (1-2)  -1   ⎤ ⎡x₁⎤ = ⎡0⎤   |   ⎡ (1-3)  -1   ⎤ ⎡x₁⎤ = ⎡0⎤   
⎣   2   (4-2) ⎦ ⎣x₂⎦   ⎣0⎦   |   ⎣   2   (4-3) ⎦ ⎣x₂⎦   ⎣0⎦   
                             |                                
get ⎰ x₁ =  1k               |   get ⎰ x₁ =  1k               
    ⎱ x₂ = -1k               |       ⎱ x₂ = -2k               
                             |                                
then eigenvector x = ⎡ 1⎤    |   then eigenvector x = ⎡ 1⎤    
                     ⎣-1⎦    |                        ⎣-2⎦

特徵向量最多N種

數學家藉由determinant與characteristic polynomial定義特徵值：N-1次多項式必有N個根，必得N個特徵值。N是方陣邊長。

這種定義方式，儘管湊出了N個特徵值，卻可能湊不出N種特徵向量。詳細分析如下：

一、相異特徵值，各自擁有相異特徵向量，構成相異維度。

二、相同特徵值（重根），不見得都有特徵向量，但至少一種。

縮放矩陣：特徵值和特徵向量都是實數。特徵向量確實存在。

旋轉矩陣：特徵值和特徵向量都是複數。特徵向量確實存在，卻無法作圖。

投影矩陣：特徵值有0。特徵向量確實存在。

歪斜矩陣：特徵向量不足N種。

特徵向量演算法（characteristic polynomial）

一、多項式函數求根，取得特徵值。

二、線性函數求根（求解），取得特徵向量。

然而無法克服：一、根的範圍不明；二、根可能是複數。

沒人用這種方法計算特徵值、特徵向量。

eigenbasis - optimization

本章主角是特徵向量與特徵值（虛擬版本）

矩陣最佳化

本站文件「optimization」，介紹了「輸入有許多個變數、輸出只有一個變數」的函數，以及極值。

本站文件「multivariate optimization」，介紹了「輸入、輸出有許多個變數」的函數，以及極值。

pseudo-eigenvalue【目前稱作Rayleigh quotient】

特徵向量擁有明確的特徵值。其他向量只能擁有虛擬的、近似的特徵值。如何得到近似的特徵值呢？利用最佳化的思維，嘗試讓Ax與λx兩者差異盡量小：對應項的差的平方的總和，越小越好。

（方便起見，x和Ax畫在同一個二維平面。）

minimize ‖Ax - λx‖²

∂
—— ‖Ax - λx‖² = 0                      「一次微分等於零」的地方是極值、鞍點
∂λ                                      二次函數、開口向上，必是最小值
  ∂
[ —— (Ax - λx) ] ∙ [ 2(Ax - λx) ] = 0   微分連鎖律
  ∂λ

[ -x ] ∙ [ 2(Ax - λx) ] = 0             微分


xᵀ A x = xᵀ λ x                         同除以-2、展開、移項

     xᵀ A x
λ = ————————                            Rayleigh quotient
      xᵀ x

虛擬特徵值宛如向量投影公式：Ax投影到x，求得投影量！分子是Ax與x兩者點積，分母是x長度平方。

pseudo-eigenvector【目前稱作Rayleigh quotient導數】

有些矩陣不存在明確的特徵向量，只能擁有虛擬的、近似的特徵向量。如何得到近似的特徵向量呢？利用微分的思維，嘗試求得虛擬特徵值的導數。

微分，就是變化程度。一次微分等於零的地方，就是沒有變化的地方。虛擬特徵值沒有變化的地方，就是虛擬特徵向量。

限制向量長度不為零，避免得到沒有討論意義的虛擬特徵向量。

∂     xᵀ A x
—— [ ———————— ] = 0  and  ‖x‖² = xᵀ x ≠ 0   特徵向量長度不為零
∂x     xᵀ x

 (A + Aᵀ) x     (xᵀ A x) 2x
———————————— - ————————————— = 0    分式微分
    xᵀ x          (xᵀ x)²

                xᵀ A x 
(A + Aᵀ) x = 2 ———————— x           移項，同乘以 xᵀ x，重新整理
                 xᵀ x

(A + Aᵀ) x = 2 λ x                  虛擬特徵值

 A + Aᵀ
——————— x = λ x                     移項，形成特徵向量的格式
   2

虛擬特徵向量即是(A+Aᵀ)/2的特徵向量。

eigenbasis - concatenation

本章主角是特徵分解（宏觀視角）

向量連接【尚無正式名稱】

函數：大量數值併成向量。輸入數值們從上往下併成輸入向量，輸出數值們從上往下併成輸出向量。

⎧ y₁ = f(x₁)        ⎡y₁⎤      ⎡x₁⎤
⎨ y₂ = f(x₂)  ———→  ⎢y₂⎥ = f( ⎢x₂⎥ )  ———→  y⃗ = f(x⃗)
⎩ y₃ = f(x₃)        ⎣y₃⎦      ⎣x₃⎦

線性函數：大量向量併成矩陣。輸入向量從左到右併成輸入矩陣，輸出向量從左到右併成輸出矩陣。

⎧ y₁ = Fx₁
⎨ y₂ = Fx₂  ———→  Y = FX
⎩ y₃ = Fx₃

⎡ |  ⎤   ⎡       ⎤ ⎡ |  ⎤
⎢ y₁ ⎥ = ⎢   F   ⎥ ⎢ x₁ ⎥
⎣ |  ⎦   ⎣       ⎦ ⎣ |  ⎦

⎡ |  ⎤   ⎡       ⎤ ⎡ |  ⎤        ⎡ |  |  |  ⎤   ⎡       ⎤ ⎡ |  |  |  ⎤
⎢ y₂ ⎥ = ⎢   F   ⎥ ⎢ x₂ ⎥  ———→  ⎢ y₁ y₂ y₃ ⎥ = ⎢   F   ⎥ ⎢ x₁ x₂ x₃ ⎥
⎣ |  ⎦   ⎣       ⎦ ⎣ |  ⎦        ⎣ |  |  |  ⎦   ⎣       ⎦ ⎣ |  |  |  ⎦

⎡ |  ⎤   ⎡       ⎤ ⎡ |  ⎤
⎢ y₃ ⎥ = ⎢   F   ⎥ ⎢ x₃ ⎥
⎣ |  ⎦   ⎣       ⎦ ⎣ |  ⎦

重複處理一件事，計算學家弄成迴圈，數學家弄成矩陣。

eigendecomposition（diagonalization）

Ax = λx，λ是數值。亦可寫成Ax = xλ，λ是1×1矩陣。

當特徵向量恰好N種，一口氣考慮所有特徵向量、特徵值。特徵向量們併成矩陣E，特徵值們併成矩陣Λ，得到AE = EΛ。

移項得到A = EΛE⁻¹，稱作「特徵分解」或「對角化」。

縮放矩陣Λ：對角線是特徵值，其餘元素是0。

                                                        -1
⎡ A₁₁ A₁₂ A₁₃ ⎤   ⎡ |  |  |  ⎤ ⎡ λ₁ 0  0  ⎤ ⎡ |  |  |  ⎤
⎢ A₂₁ A₂₂ A₂₃ ⎥ = ⎢ x₁ x₂ x₃ ⎥ ⎢ 0  λ₂ 0  ⎥ ⎢ x₁ x₂ x₃ ⎥
⎣ A₃₁ A₃₂ A₃₃ ⎦   ⎣ |  |  |  ⎦ ⎣ 0  0  λ₃ ⎦ ⎣ |  |  |  ⎦
       A               E            Λ            E⁻¹

J = [3 0 0; 0 2 0; 0 0 1]; E = [1 1 1; 1 2 3; 1 3 6]; A = E * J * inv(E)

當特徵向量恰好N種，就可以特徵分解。舉例來說，恆等矩陣、旋轉矩陣、鏡射矩陣、縮放矩陣、投影矩陣可以特徵分解，歪斜矩陣、移位矩陣不可以特徵分解。

特徵值矩陣Λ，只有唯一一種（重新排列不算在內）。特徵向量矩陣E，通常有許多種（若特徵值相同，則特徵向量方向隨意）。

Jordan–Chevalley decomposition

當特徵向量不足N種，補足成N種特徵向量，對應到N個特徵值。

原始的特徵向量，Ax = λx。補足的特徵向量，Ax ≠ λx，兩邊有差異。數學家採用了「移位shift」：取特徵值相同的、隔壁的特徵向量，填補差異，Axᵢ = λxᵢ + xᵢ₋₁。xᵢ與xᵢ₋₁的特徵值相同。

特徵向量們併成矩陣E，特徵值們併成矩陣Λ，移位們併成分塊移位矩陣S，得到AE = E(Λ + S)。

令J = Λ + S。移項得到A = EJE⁻¹，稱作「喬登分解」。

縮放矩陣Λ：對角線是特徵值。分塊移位矩陣S：次對角線填入1。喬登標準型J：兩者相加。

                                                        -1
⎡ A₁₁ A₁₂ A₁₃ ⎤   ⎡ |  |  |  ⎤ ⎡ λ₁ 1  0  ⎤ ⎡ |  |  |  ⎤
⎢ A₂₁ A₂₂ A₂₃ ⎥ = ⎢ x₁ x₂ x₃ ⎥ ⎢ 0  λ₂ 0  ⎥ ⎢ x₁ x₂ x₃ ⎥
⎣ A₃₁ A₃₂ A₃₃ ⎦   ⎣ |  |  |  ⎦ ⎣ 0  0  λ₃ ⎦ ⎣ |  |  |  ⎦
       A               E            J            E⁻¹

J = [1 1 0; 0 1 0; 0 0 1]; E = [1 1 1; 0 1 2; 0 0 1]; A = E * J * inv(E)

所有方陣都可以喬登分解。此處省略證明。

喬登標準型J，只有唯一一種（人為規定）。特徵向量矩陣E，通常有許多種（若特徵值相同，則特徵向量方向隨意）。

三個問題

可對角化（存在N種特徵向量）：已有演算法。我沒有學會。
特徵分解（找出N種特徵向量）：已有演算法。後續章節陸續介紹。
喬登分解（補足N種特徵向量）：尚無演算法。目前的演算法均不保證數值穩定。

喬登分解採用了「移位shift」，事情變得相當複雜。取哪個特徵向量進行移位，才能填補差異，難以設計演算法。

eigenbasis - dot product🚧

本章主角是方向餘弦矩陣

direction cosine matrix

「方向餘弦矩陣」。兩組向量，每一對向量之間的夾角餘弦值。

    ⎡ |  |  |  ⎤       ⎡ |  |  |  ⎤
A = ⎢ a₁ a₂ a₃ ⎥   B = ⎢ b₁ b₂ b₃ ⎥
    ⎣ |  |  |  ⎦       ⎣ |  |  |  ⎦

             a∙b
cos(a,b) = ———————
           ‖a‖ ‖b‖

           ⎡ cos(a₁,b₁) cos(a₁,b₂) cos(a₁,b₃) ⎤
cos(A,B) = ⎢ cos(a₂,b₁) cos(a₂,b₂) cos(a₂,b₃) ⎥
           ⎣ cos(a₃,b₁) cos(a₃,b₂) cos(a₃,b₃) ⎦

兩個正規正交矩陣的方向餘弦矩陣，可以簡單寫成兩兩交叉點積。此時方向餘弦矩陣恰是正規正交矩陣。

           ⎡ a₁∙b₁ a₁∙b₂ a₁∙b₃ ⎤
cos(A,B) = ⎢ a₂∙b₁ a₂∙b₂ a₂∙b₃ ⎥ = AᵀB   when A and B are both
           ⎣ a₃∙b₁ a₃∙b₂ a₃∙b₃ ⎦         orthonormal matrices

求得變換矩陣，藉由反矩陣。
（find transformation matrix with inverse matrix）

已知多組輸入x輸出y，求變換矩陣F。

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ |  ⎤
⎢ y₁ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ x₁ ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ |  ⎦

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ |  ⎤
⎢ y₂ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ x₂ ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ |  ⎦

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ |  ⎤
⎢ y₃ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ x₃ ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ |  ⎦

方法是向量連接，併成矩陣。

⎡ |  |  |  ⎤   ⎡ |  |  |  ⎤ ⎡ |  |  |  ⎤
⎢ y₁ y₂ y₃ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ x₁ x₂ x₃ ⎥
⎣ |  |  |  ⎦   ⎣ |  |  |  ⎦ ⎣ |  |  |  ⎦
     Y              F            X

然後移項形成反矩陣。

                          -1
⎡ |  |  |  ⎤  ⎡ |  |  |  ⎤    ⎡ |  |  |  ⎤
⎢ y₁ y₂ y₃ ⎥  ⎢ x₁ x₂ x₃ ⎥  = ⎢ v₁ v₂ v₃ ⎥
⎣ |  |  |  ⎦  ⎣ |  |  |  ⎦    ⎣ |  |  |  ⎦
     Y             X⁻¹             F

先算反矩陣，再算矩陣乘法，得到變換矩陣。

求得變換矩陣，藉由標準座標軸。
（find transformation matrix with standard basis）

構造數學式子：輸入是恆等矩陣I，輸出是變換矩陣F。

⎡ |  |  |  ⎤   ⎡ |  |  |  ⎤ ⎡ 1 0 0 ⎤
⎢ v₁ v₂ v₃ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ 0 1 0 ⎥
⎣ |  |  |  ⎦   ⎣ |  |  |  ⎦ ⎣ 0 0 1 ⎦
     F              F           I

向量分離。輸入矩陣、輸出矩陣，拆成向量。

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ 1 ⎤
⎢ v₁ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ 0 ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ 0 ⎦

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ 0 ⎤
⎢ v₂ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ 1 ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ 0 ⎦

⎡ |  ⎤   ⎡ |  |  |  ⎤ ⎡ 0 ⎤
⎢ v₃ ⎥ = ⎢ v₁ v₂ v₃ ⎥ ⎢ 0 ⎥
⎣ |  ⎦   ⎣ |  |  |  ⎦ ⎣ 1 ⎦

如果能夠知道每一個標準座標軸的輸出，那麼輸出併成矩陣，即是變換矩陣。縮放、旋轉、投影、……，皆可如此求得矩陣。

座標變換（coordinate transformation）

座標系固定不變，幾何元件實施變換，座標隨之改變。

例如平移、縮放、旋轉都屬於座標變換。

考慮一個特殊問題：兩組座標A與B，一一對應。從A變成B，請找到變換矩陣F。

座標x轉換成座標y，數學公式y = Fx。

運用向量連接，求得變換矩陣：每個A座標分別旋轉成每個B座標，併成矩陣。數學公式B = FA。移項得到F = BA⁻¹。

當A有反矩陣，例如正規正交矩陣，那麼就存在變換矩陣T。

當A與B恰是正規正交矩陣，那麼T也是正規正交矩陣。正規正交矩陣可以視作旋轉矩陣。問題變成：A轉成B，請找到轉法。

座標系變換（coordinate system transformation）

幾何元件固定不變，重新訂立座標系，座標隨之改變。

例如極座標表示法就屬於座標系變換。

考慮一個特殊問題：兩組直角座標系A與B，原點相同，座標軸長度是1、角度不同。座標系從A改成B，請找到變換矩陣F。

座標x轉換成座標y，數學公式y = Fx。

運用向量分離，解釋變換過程：一、一個A座標拆成多個A分量。二、一個A分量，分別投影到每個B向量，得到多個B分量。每個A分量，以此類推。對應的B分量，實施加總，簡化成一個B分量。三、多個B分量併成一個B座標。

運用標準座標軸，求得變換矩陣：每個標準座標軸的變換結果，併成變換矩陣。

數學公式F = BᵀA，方向餘弦矩陣的轉置。（為了省略轉置二字，許多教材將方向餘弦矩陣直接定義為BᵀA。）

既然A與B是正規正交矩陣，那麼F也是正規正交矩陣。正規正交矩陣可以視作旋轉矩陣。然而F的幾何意義是兩兩交叉投影，而不是旋轉。

【尚無正式名稱】

藉由方向餘弦矩陣，旋轉可以改寫成兩兩交叉投影，反之亦然。

欸，真的有這種東西嗎？

eigenbasis - cross product🚧

本章主角是餘因子矩陣

cofactor matrix / adjugate matrix

「餘因子矩陣」。隨便砍掉一個直條與一個橫條，剩下來的矩陣。根據砍掉的位置，乘上係數+1或-1。

「伴隨矩陣」是其轉置。

Bᵀ變A⁻¹。內積變外積。點積變叉積。投影變容積。

inv(A) = cof(A)ᵀ / det(A) = adj(A) / det(A)

adj(A) = axyᵀ   where Ax = 0 and Aᵀy = 0

transpose與inverse的運算

想要順便講一下基本運算，但是還沒想到幾何證明方式。

(AB)ᵀ = BᵀAᵀ
(AB)⁻¹ = B⁻¹A⁻¹
(A⁻¹)ᵀ = (Aᵀ)⁻¹
https://math.stackexchange.com/questions/340233/

矩陣變換：A變Aᵀ

可以特徵分解：一種答案是特徵基底的外積矩陣的反矩陣(EEᵀ)⁻¹，或者改寫成反矩陣的內積矩陣E⁻¹ᵀE⁻¹。我還不知道幾何意義是什麼。

https://www.physicsforums.com/threads/470322/
https://mathoverflow.net/questions/122345/
https://math.stackexchange.com/questions/1760759/
https://math.stackexchange.com/questions/94599/

A = EΛE⁻¹   --->   Λ = E⁻¹AE

=> Aᵀ = (EΛE⁻¹)ᵀ = E⁻¹ᵀΛᵀEᵀ
=> Aᵀ = E⁻¹ᵀΛEᵀ
=> Aᵀ = E⁻¹ᵀ(E⁻¹AE)Eᵀ
=> Aᵀ = (EEᵀ)⁻¹A(EEᵀ)
=> F = (EEᵀ)⁻¹
=> F = E⁻¹ᵀE⁻¹

只能喬登分解：需要額外考慮移位。喬登標準型J的每個分塊各自反向移位，利用逆排矩陣B。

⎰ A = EJE⁻¹    --->   J = E⁻¹AE
⎱ J = BJᵀB⁻¹   --->   Jᵀ = B⁻¹JB   (Bᵀ = B⁻¹ = B)

=> Aᵀ = (EJE⁻¹)ᵀ = E⁻¹ᵀJᵀEᵀ
=> Aᵀ = E⁻¹ᵀ(B⁻¹JB)Eᵀ
=> Aᵀ = E⁻¹ᵀ(B⁻¹(E⁻¹AE)B)Eᵀ
=> Aᵀ = (EBEᵀ)⁻¹A(EBEᵀ)
=> F = (EBEᵀ)⁻¹
https://math.stackexchange.com/questions/94599/

容積演算法（Laplace expansion）（cofactor expansion）

https://en.wikipedia.org/wiki/Laplace_expansion

容積演算法（Chio's matrix condensation）

Lagrange's identity ‖a‖²‖b‖² = ‖a∙b‖² + ‖a×b‖²
https://en.wikipedia.org/wiki/First_fundamental_form#Calculating_lengths_and_areas

Chio's matrix condensation
http://mathworld.wolfram.com/ChioPivotalCondensation.html
https://ccjou.wordpress.com/2009/11/24/
https://terrytao.wordpress.com/2019/08/13/eigenvectors-from-eigenvalues/
一直叉積  算出容積

eigenbasis - iteration

本章主角是線性動態系統

矩陣次方

藉由特徵分解、喬登分解，矩陣次方化作特徵值次方。

A  = EΛE⁻¹
A² = EΛE⁻¹EΛE⁻¹ = EΛ²E⁻¹
A³ = EΛE⁻¹EΛE⁻¹EΛE⁻¹ = EΛ³E⁻¹

特徵值矩陣Λ的次方，就是對角線元素們的次方。

喬登標準型J的次方，請參考：

UVa 10753

矩陣根號

矩陣平方根：

A = EΛE⁻¹ = E√Λ√ΛE⁻¹ = E√ΛE⁻¹E√ΛE⁻¹ = (E√ΛE⁻¹)²
√A = E√ΛE⁻¹

矩陣立方根：

A = EΛE⁻¹ = E∛Λ∛Λ∛ΛE⁻¹ = E∛ΛE⁻¹E∛ΛE⁻¹E∛ΛE⁻¹ = (E∛ΛE⁻¹)³
∛A = E∛ΛE⁻¹

特徵值矩陣Λ的根號，就是對角線元素們的根號。

喬登標準型J的根號，請參考：

linear function可以畫成箭矢圖

本站文件「function」，介紹了「輸入有很多個變數，輸出只有一個變數」的函數，也畫出了函數圖形。

此處介紹「輸入有很多個變數，輸出有很多個變數」的線性函數，並且畫出函數圖形。

此處以ℝ² ⇨ ℝ²線性函數為例，輸入是兩個變數、輸出是兩個變數，畫在二維平面。為了節省紙張，改成畫在同一個二維平面。

一組輸入輸出仍可辨識，多組輸入輸出就看不清了。因此改成「從輸入往輸出的箭頭：輸出減輸入」。

窮舉所有輸入輸出，繪圖結果一片滿，看不出任何意義。因此改成「輸入等距取樣、保持間隔」。

線性函數是一種特別的函數，函數圖形有著一種難以言喻的整齊。有時是一致朝外，有時是一致螺旋。

這樣的函數圖形，可以用來解釋真實世界的物理現象！例如磁場、氣旋──不過這是後話了。先讓我們專注於線性函數吧！

linear function的本質：朝著目標向前進

整齊的原因是什麼呢？數學家已經初步解決了這個問題！

數學家猜測，所謂的整齊，也許是指：「大家有著共識，方向一致。」再進一步猜測：「這當中有沒有堪稱表率，讓大家群起效尤的輸入輸出呢？是不是有人一路以來始終如一，堅持走在正確的道路上？」於是數學家嘗試尋找：「有哪些輸入向量，套用線性函數之後，方向保持不變。」

套用線性函數。特徵值為實數則縮放，虛數則旋轉，零則消滅，負則翻轉。

反覆套用線性函數。特徵值絕對值小於1則趨近原點，等於1則保持不動，大於1則趨近無限遠，越小則縮短越快，越大則伸長越快。特徵值為正數則連貫移動，負數則來回翻轉、交互跳躍，而整體趨勢仍與正數相同。向量漸漸偏向特徵值絕對值最大的方向。

trace–determinant diagram

以trace和determinant判斷流向。

UVa 720

特徵向量演算法（power iteration）

只能求出「絕對值最大的特徵值的特徵向量」。

反覆套用線性函數，向量朝向「絕對值最大的特徵值的特徵向量」的方向。就這麼簡單。

反覆套用線性函數，向量越來越長或越來越短，超出浮點數範圍。解法是隨時將向量長度縮放為1。這麼做不會影響方向。

時間複雜度O(N²T)，N是矩陣邊長，T是遞推次數。

收斂速度，取決於最大、次大特徵值的絕對值的比值|λ₁|/|λ₂|。比值越大，收斂速度越快，遞推次數越少。

特徵向量演算法（inverse iteration）

有兩種手法可以調整特徵值大小：

一、矩陣對角線加上同一個值，則特徵值們也會加上同一個值：(A + αI)x = Ax + αx = λx + αx = (λ + α)x。

二、矩陣倒數（反矩陣），則特徵值們也會倒數。

猜測特徵值，再將該特徵值調整成無限大，成為最大的特徵值。如此一來power iteration就能得到其特徵向量。

eigenvalues of A: {λ₁, λ₂, λ₃, ...}   ( |λ₁| ≥ |λ₂| ≥ |λ₃| ≥ ... )
let B = (A - αI)⁻¹
eigenvalues of B: {1/(λ₁-α), 1/(λ₂-α), 1/(λ₃-α), ...}

α猜得準，足夠靠近λᵢ，讓1/(λᵢ-α)很大，成為B的絕對值最大的特徵值。
此時B套用power iteration即可得到1/(λᵢ-α)的特徵向量，進而算出λᵢ。

xₖ₊₁ = B xₖ = (A - αI)⁻¹ xₖ

錦上添花，可以使用LU decomposition加速。遞推一次從O(N³)降為O(N²)。

xₖ₊₁ = (A - αI)⁻¹ xₖ
(A - αI) xₖ₊₁ = xₖ
find LU decomposition of (A - αI)

然而特徵值該怎麼猜呢？你還記得虛擬特徵值嗎？

特徵向量演算法（Rayleigh quotient iteration）

每一步以Rayleigh quotient重新估計特徵值。

αₖ = (xₖᵀ A xₖ) / (xₖᵀ xₖ)
xₖ₊₁ = B xₖ = (A - αₖI)⁻¹ xₖ

遺珠之憾，不能使用LU decomposition加速。遞推一次需時O(N³)。

然而初始向量該怎麼猜呢？我也不知道，嘻嘻。

eigenbasis - projection

本章主角是線性展開

矩陣內積【尚無正式名稱】

向量內積：元素兩兩相乘，通通相加，形成數值。

【註：向量內積恰是向量點積。】

a∙b = aᵀb = c

⎡ a₁ ⎤   ⎡ b₁ ⎤
⎢ a₂ ⎥ ∙ ⎢ b₂ ⎥ = a₁b₁ + a₂b₂ + a₃b₃
⎣ a₃ ⎦   ⎣ b₃ ⎦

矩陣內積：向量交叉內積，形成矩陣。

【註：此處的定義比較另類。數學家的定義是tr(AᵀB)。】

A∙B = AᵀB = C

⎡ |  |  |  ⎤   ⎡ |  |  |  ⎤   ⎡ a₁∙b₁ a₁∙b₂ a₁∙b₃ ⎤
⎢ a₁ a₂ a₃ ⎥ ∙ ⎢ b₁ b₂ b₃ ⎥ = ⎢ a₂∙b₁ a₂∙b₂ a₂∙b₃ ⎥
⎣ |  |  |  ⎦   ⎣ |  |  |  ⎦   ⎣ a₃∙b₁ a₃∙b₂ a₃∙b₃ ⎦

矩陣外積【尚無正式名稱】

向量外積：元素交叉相乘，形成矩陣。

【註：向量外積恰是克羅內克積。】

a□b = abᵀ = C

⎡ a₁ ⎤   ⎡ b₁ ⎤   ⎡ a₁b₁ a₁b₂ a₁b₃ ⎤
⎢ a₂ ⎥ □ ⎢ b₂ ⎥ = ⎢ a₂b₁ a₂b₂ a₂b₃ ⎥
⎣ a₃ ⎦   ⎣ b₃ ⎦   ⎣ a₃b₁ a₃b₂ a₃b₃ ⎦

矩陣外積：向量兩兩外積，通通相加，形成矩陣。

【註：此處的定義比較另類。數學家的定義採用了分塊矩陣。】

A□B = ABᵀ = C

⎡ |  |  |  ⎤   ⎡ |  |  |  ⎤
⎢ a₁ a₂ a₃ ⎥ □ ⎢ b₁ b₂ b₃ ⎥ = a₁□b₁ + a₂□b₂ + a₃□b₃
⎣ |  |  |  ⎦   ⎣ |  |  |  ⎦

矩陣乘法的内積形式與外積形式

內積形式：橫條配直條，交叉內積。N²個元素合併。

外積形式：直條配橫條，兩兩外積。N個矩陣加總。

     ⎡ —— a₁* —— ⎤ ⎡ |  |  |  ⎤   ⎡ a₁*b₁ a₁*b₂ a₁*b₃ ⎤
AB = ⎢ —— a₂* —— ⎥ ⎢ b₁ b₂ b₃ ⎥ = ⎢ a₂*b₁ a₂*b₂ a₂*b₃ ⎥
     ⎣ —— a₃* —— ⎦ ⎣ |  |  |  ⎦   ⎣ a₃*b₁ a₃*b₂ a₃*b₃ ⎦

                                  informal: view aᵢ*bⱼ as aᵢ∙bⱼ
                                  formal: aᵢ*bⱼ = (aᵢ*ᵀ)∙bⱼ

     ⎡ |  |  |  ⎤ ⎡ —— b₁* —— ⎤
AB = ⎢ a₁ a₂ a₃ ⎥ ⎢ —— b₂* —— ⎥ = a₁b₁* + a₂b₂* + a₃b₃*
     ⎣ |  |  |  ⎦ ⎣ —— b₃* —— ⎦

                                  informal: view aᵢbᵢ* as aᵢ□bᵢ
                                  formal: aᵢbᵢ* = aᵢ□(bᵢ*ᵀ)

*的意義是取橫條。*的文言文是索引轉置。

【註：此處的定義比較另類。數學家的定義是共軛轉置。】

垂直投影（orthogonal projection）

向量a垂直投影到向量b，可以寫成向量內積。

          a∙b     aᵀb
proj_b(a) = ——— b = ——— b
          b∙b     bᵀb

垂直投影是線性函數。改寫成矩陣，恰是外積除以內積。

          (a∙b)b   (aᵀb)b   b(aᵀb)   b(bᵀa)   (bbᵀ)a   (b□b)a
proj_b(a) = —————— = —————— = —————— = —————— = —————— = ——————
            b∙b      b∙b     b∙b      b∙b       b∙b      b∙b

          b□b   bbᵀ
proj_b( ) = ——— = ———
          b∙b   bᵀb

如果b是單位向量，長度為1，可以省略分母。

proj_b( ) = bbᵀ     when b is an unit vector
                  ‖b‖² = b∙b = bᵀb = 1

向量a垂直投影到大量向量B所在空間，分母在中間。

proj_B( ) = B(BᵀB)⁻¹Bᵀ

如果B是正規正交矩陣，長度為I，可以省略分母。

proj_B( ) = BBᵀ     when B is an orthonormal matrix
                  ‖B‖² = B∙B = BᵀB = I

順帶一提，輸入向量減去投影結果，可以進一步消滅空間。

a - proj_B(a) = a - B(BᵀB)⁻¹Bᵀ a = (I - B(BᵀB)⁻¹Bᵀ) a

保存向量b的線性函數：bbᵀ/bᵀb              = bb⁺
消滅向量b的線性函數：I - bbᵀ/bᵀb          = I - bb⁺
保存大量向量B的線性函數：B(BᵀB)⁻¹Bᵀ       = BB⁺
消滅大量向量B的線性函數：I - B(BᵀB)⁻¹Bᵀ   = I - BB⁺

輸入可以是大量向量，大量向量併成矩陣。輸入可以是矩陣。

proj_B(A) = B(BᵀB)⁻¹BᵀA
A - proj_B(A) = (I - B(BᵀB)⁻¹Bᵀ) A

傾斜投影（oblique projection）

兩種做法。

一、反矩陣：向量a傾斜投影到向量x，沿著向量y。細節請見專著《Matrix Analysis for Statistics》。

    ⎡ |      |  ⎤       ⎡ |      |  ⎤
X = ⎢ x₁ ... xᵢ ⎥   Y = ⎢ y₁ ... yⱼ ⎥
    ⎣ |      |  ⎦       ⎣ |      |  ⎦

    ⎡ |      |  |      |  ⎤     V must be invertible
V = ⎢ x₁ ... xᵢ y₁ ... yⱼ ⎥     (size is N×N and rank is N)
    ⎣ |      |  |      |  ⎦     (X and Y are disjoint and complementary)

     ⎡ Iᵢₓᵢ | O ⎤        ⎡ O |   O  ⎤
I_X = ⎢------|---⎥   I_Y = ⎢---|------⎥
     ⎣  O   | O ⎦        ⎣ O | Iⱼₓⱼ ⎦

proj_X,Y( ) = VI_XV⁻¹ = I - VI_YV⁻¹     project to X along Y
proj_Y,X( ) = VI_YV⁻¹ = I - VI_XV⁻¹     project to Y along X

二、外積除以內積：向量a傾斜投影到向量x，取決於向量y的傾斜方向。【尚待確認】

如果是垂直投影，那麼x和y方向相同。

          x□y   xyᵀ
proj_x,y( ) = ——— = ———
          y∙x   yᵀx

proj_X,Y( ) = X(YᵀX)⁻¹Yᵀ

矩陣不動點【目前稱作projection matrix】

「投影矩陣」。P = P²。變換一次等同變換數次。保存的空間一直都在，消滅的空間回不來了。

I - P也是投影矩陣。只需證明I - P = (I - P)²。

如果是垂直投影，那麼P = Pᵀ。

eigendecomposition（spectral decomposition）

特徵分解改寫成矩陣外積形式。再改寫成投影矩陣的加權平均數（權重是特徵值）。

⎡ X X X ⎤   ⎡ |  |  |  ⎤ ⎡ λ₁ 0  0  ⎤ ⎡ —— y₁* —— ⎤
⎢ X X X ⎥ = ⎢ x₁ x₂ x₃ ⎥ ⎢ 0  λ₂ 0  ⎥ ⎢ —— y₂* —— ⎥
⎣ X X X ⎦   ⎣ |  |  |  ⎦ ⎣ 0  0  λ₃ ⎦ ⎣ —— y₃* —— ⎦
    A            E            Λ            E⁻¹

1. A = λ₁x₁y₁* + ... + λₙxₙyₙ* = sum λᵢxᵢyᵢ*
     = λ₁P₁    + ... + λₙPₙ    = sum λᵢPᵢ
2. I =   x₁y₁* + ... +   xₙyₙ* = sum xᵢyᵢ*
     =   P₁    + ... +   Pₙ    = sum Pᵢ
3. E⁻¹E = I    --->  yᵢ*xᵢ = 1 , yᵢ*xⱼ = 0
4. xᵢyᵢ* = Pᵢ  --->  PᵢPᵢ = Pᵢ , PᵢPⱼ = 0

Jordan–Chevalley decomposition

喬登分解如法炮製。基本單位從向量變成了正規正交矩陣（不是方陣），矩陣數量符合Jordan block數量。另外多了餘數。細節請見專著《Eigenvalues of Matrices》。

               X₁   X₂
              ____  _     ______
⎡ X X X ⎤   ⎡|X  X||X|⎤ ⎡|λ₁ 1 |0  ⎤ ⎡|X̅‾‾X̅‾‾X̅|⎤ Y₁*
⎢ X X X ⎥ = ⎢|X  X||X|⎥ ⎢|0  λ₂|0̲_ ⎥ ⎢|X̲__X̲__X̲|⎥
⎣ X X X ⎦   ⎣|X  X||X|⎦ ⎣ ‾‾‾‾‾|λ₃|⎦ ⎣|X̲__X̲__X̲|⎦ Y₂*
    A         ‾‾‾‾E ‾        J  ‾‾       E⁻¹

1. A = sum λᵢXᵢYᵢ* + Dᵢ = sum λᵢPᵢ + Dᵢ     where Dᵢ = (A - λᵢI)P
2. I = sum XᵢYᵢ*        = sum Pᵢ
3. E⁻¹E = I    --->  Yᵢ*Xᵢ = I , Yᵢ*Xⱼ = O
4. XᵢYᵢ* = Pᵢ  --->  PᵢPᵢ = Pᵢ , PᵢPⱼ = O
5. Xᵢ is orthonormal
6. sum over k Jordan block (each Jordan block has same eigenvalue)

特徵值歸零演算法（matrix deflation）

deflation在英文當中是指洩氣消掉。已知任何一個特徵值暨特徵向量，藉由投影與縮放，讓特徵值歸零。

消掉λ₁x₁y₁*，讓特徵值λ₁歸零。然而難以求得y₁*，只好換成其他向量v。令v∙x₁ = 1，減去λ₁x₁vᵀ，讓特徵值λ₁歸零。細節請見專著《A Friendly Introduction to Numerical Analysis》。

eigenvalues  of A: {λ₁, λ₂, λ₃, ...}     |λ₁| ≥ |λ₂| ≥ |λ₃| ≥ ...
eigenvectors of A: {x₁, x₂, x₃, ...}     ‖x₁‖ = ‖x₂‖ = ‖x₃‖ = ... = 1
let v∙x₁ = vᵀx₁ = 1
let B = A - λ₁x₁vᵀ
eigenvalues  of B: { 0, λ₂, λ₃, ...}
eigenvectors of B: {x₁, x̕₂, x̕₃, ...}     xᵢ = (λᵢ-λ₁)x̕ᵢ + λ₁(vᵀx̕ᵢ)x₁

該特徵值歸零，該特徵向量不變。其餘特徵值皆不變，其餘特徵向量皆改變。

N次power iteration，逐次消滅絕對值最大的特徵值。不必猜測特徵值了！

特徵值歸零演算法（Wielandt deflation）

取A的橫條當作v，調整一下倍率。

let v = [ith row of A]ᵀ / (λ₁x₁ᵢ)

B的對應橫條將歸零，減少誤差。

特徵值歸零演算法（Hotelling deflation）

僅適用於對稱矩陣：特徵向量正規正交。x₁ᵀ = y₁*。

取特徵向量當作v，調整一下長度。

let B = A - λ₁x₁x₁ᵀ and ‖x₁‖ = 1
eigenvalues  of B: { 0, λ₂, λ₃, ...}
eigenvectors of B: {x₁, x₂, x₃, ...}

垂直投影，特徵向量皆不變。

eigenbasis - transformation

本章主角是相似矩陣：特徵值相同（嚴格來說是喬登標準型相同）

矩陣變換【目前稱作similar transformation】

本站文件「transformation」，介紹了數據的變換。

一、數據：此處是向量x。

二、數據變換：此處採用線性函數F，向量x變換成向量Fx。

三、為了還原數據，於是令變換函數F擁有反函數F⁻¹。

x ———F———> Fx
  <——F⁻¹——

此處介紹線性函數的變換。事情稍微複雜一點。

四、函數：此處是線性函數A。

五、為了對應運算，於是令函數A擁有對應運算B。

六、函數變換：函數A變換成函數B，推導得到FAF⁻¹ = B。

  ———F———>        ⎰ y = A(x)
x <——F⁻¹—— Fx      ⎱  Fy = B(Fx)
|           |
A ———?———>  B     => FAx = BFx
↓           ↓     => FA = BF
y ———F———> Fy     => FAF⁻¹ = B
  <——F⁻¹——        => A = F⁻¹BF

也有人寫成A = F⁻¹BF。矩陣連乘要從右往左讀，意義是繞一圈回來仍是A。

|          |
| ———F———> | 
A          B      A = F⁻¹BF
| <——F⁻¹—— |
↓          ↓

特徵分解、喬登分解，都是矩陣變換。

|          |                    |          |
| ———E⁻¹—> |                    | ———E⁻¹—> |
A          Λ      A = EΛE⁻¹     A          J      A = EJE⁻¹
| <——E———— |                    | <——E———— |
↓          ↓                    ↓          ↓

矩陣變換，特徵值不變。

特徵向量隨之變換，特徵值保持不變。

⎰ Ax = λx       eigenvector and eigenvalue (if exist)
⎱ A = F⁻¹BF     similarity transformation

=> (F⁻¹BF)x = λx
=> BFx = Fλx
=> B(Fx) = λ(Fx)

特徵向量隨之變換，喬登標準型保持不變。

⎰ A = EJE⁻¹     Jordan–Chevalley decomposition
⎱ B = FAF⁻¹     similarity transformation

=> B = F(EJE⁻¹)F⁻¹
=> B = (FE)J(E⁻¹F⁻¹)
=> B = (FE)J(FE)⁻¹

similar transformation / similar matrices

FAF⁻¹ = B稱作「相似變換」。A與B稱作「相似矩陣」。

相似變換，特徵值不變。運用相似變換調整矩陣，讓矩陣更容易求得特徵值，甚至一口氣求得大量特徵值。

Schur decomposition

「上三角矩陣」。對角線之下的元素必須是零。

    ⎡ T₁₁ T₁₂ T₁₃ T₁₄ T₁₅ ⎤
    ⎢  0  T₂₂ T₂₃ T₂₄ T₂₅ ⎥
T = ⎢  0   0  T₃₃ T₃₄ T₃₅ ⎥     upper triangular matrix
    ⎢  0   0   0  T₄₄ T₄₅ ⎥
    ⎣  0   0   0   0  T₅₅ ⎦

「上三角分解」。實施相似變換，變成上三角矩陣。變換矩陣必須是（實數版本）正規正交矩陣、（複數版本）么正矩陣。

實數版本不一定可以上三角分解。複數版本一定可以。

Q⁻¹AQ = T     Schur decomposition

上三角分解可以一口氣求得大量特徵值：上三角矩陣的對角線元素們，即是特徵值們。上三角分解無法求得特徵向量。

【註：以特徵向量求特徵值，檢查長度縮放倍率即可。以特徵值求特徵向量，目前沒有簡單的演算法。】

順帶一提，如果A是對稱矩陣，那麼T也是對稱矩陣。T是對角線矩陣，即是特徵值。Q是正規正交矩陣，即是特徵向量。

上三角分解演算法（QR iteration）

矩陣A實施QR分解。

A = QR

A實施相似變換，相似變換矩陣是Q⁻¹，得到相似矩陣RQ。

Q⁻¹AQ = RQ

RQ擁有相同特徵值。然而RQ不是什麼特別的矩陣，依然難以求得特徵值。那就多做幾遍吧。不斷實施相似變換，碰碰運氣，看看會不會變成上三角矩陣。

遞推一次O(N³)。時間複雜度O(N³T)。

收斂條件：特徵值皆相異。導致特徵向量方向明確。

收斂速度：取決於特徵值的絕對值的比值|λᵢ|/|λⱼ|。比值越遠離1，收斂速度越快，遞推次數越少。

對稱矩陣：收斂至對角線矩陣，即是特徵值。每回合的Q依序連乘，Q₁Q₂Q₃...Qₜ，即是特徵向量。

上三角分解演算法（Hessenberg QR iteration）

幾乎上三角矩陣的情況下，幾乎都是零了，遞推次數T大幅減少，但是收斂速度依然相同。

幾乎上三角矩陣的情況下，遞推一次降為O(N²)，時間複雜度降為O(N²T)。

一、幾乎上三角矩陣的QR分解，採用「軸面旋轉矩陣Givens matrix」，時間複雜度從O(N³)降為O(N²)。

二、幾乎上三角矩陣的Q，可以由Q的第一個直條遞推而得，以此快速算出RQ，時間複雜度從O(N³)降為O(N²)。

三、RQ仍然是幾乎上三角矩陣。

Hessenberg decomposition

「幾乎上三角矩陣」。次對角線之下的元素必須是零。

    ⎡ H₁₁ H₁₂ H₁₃ H₁₄ H₁₅ ⎤
    ⎢ H₂₁ H₂₂ H₂₃ H₂₄ H₂₅ ⎥
H = ⎢  0  H₃₂ H₃₃ H₃₄ H₃₅ ⎥     Hessenberg matrix
    ⎢  0   0  H₄₃ H₄₄ H₄₅ ⎥
    ⎣  0   0   0  H₅₄ H₅₅ ⎦

「幾乎上三角分解」。實施相似變換，變成幾乎上三角矩陣。變換矩陣必須是（實數版本）正規正交矩陣、（複數版本）么正矩陣。

Q⁻¹AQ = H     Hessenberg decomposition

任意矩陣一定可以幾乎上三角分解。

幾乎上三角分解演算法，其實就是QR分解那三個演算法的改良版，QR分解改成了相似變換。

大量特徵值演算法，可以分成兩步驟：第一步先變成幾乎上三角矩陣，第二步再變成上三角矩陣。時間複雜度較低。

Hessenberg decomposition: Q⁻¹AQ = H
Schur decomposition:      Q⁻¹AQ = T

幾乎上三角分解演算法（Arnoldi iteration）

矩陣AQ實施QR分解，將Q移項即是相似變換。強制使用同一個Q的緣故，分解之後不是上三角矩陣，而是幾乎上三角矩陣。

AQ = QH
Q⁻¹AQ = H

矩陣A實施QR分解。

⎡ A₁₁ A₁₂ A₁₃ A₁₄ ⎤   ⎡  |  |  |  | ⎤ ⎡ R₁₁ R₁₂ R₁₃ R₁₄ ⎤
⎢ A₂₁ A₂₂ A₂₃ A₂₄ ⎥ = ⎢ q₁ q₂ q₃ q₄ ⎥ ⎢  0  R₂₂ R₂₃ R₂₄ ⎥
⎢ A₃₁ A₃₂ A₃₃ A₃₄ ⎥   ⎣  |  |  |  | ⎦ ⎢  0   0  R₃₃ R₃₄ ⎥
⎣ A₄₁ A₄₂ A₄₃ A₄₄ ⎦                   ⎣  0   0   0  R₄₄ ⎦
         A                   Q                 R

矩陣AQ實施QR分解。

                                         ⎡ H₁₁ H₁₂ H₁₃ H₁₄ ⎤
⎡  |   |   |   |  ⎤   ⎡  |  |  |  |  | ⎤ ⎢ H₂₁ H₂₂ H₂₃ H₂₄ ⎥
⎢ Aq₁ Aq₂ Aq₃ Aq₄ ⎥ = ⎢ q₁ q₂ q₃ q₄ q₅ ⎥ ⎢  0  H₃₂ H₃₃ H₃₄ ⎥
⎣  |   |   |   |  ⎦   ⎣  |  |  |  |  | ⎦ ⎢  0   0  H₄₃ H₄₄ ⎥
                                         ⎣  0   0   0  H₅₄ ⎦
        AQ                     Q̃                  H̃

數學式子。H̃切成直條，視作座標。

Aq₁ = H₁₁q₁ + H₂₁q₂                           (1st column of H)
Aq₂ = H₁₂q₁ + H₂₂q₂ + H₃₂q₃                   (2nd column of H)
Aq₃ = H₁₃q₁ + H₂₃q₂ + H₃₃q₃ + H₄₃q₄           (3rd column of H)
Aq₄ = H₁₄q₁ + H₂₄q₂ + H₃₄q₃ + H₄₄q₄ + H₅₄q₅   (4th column of H)
where Hᵢⱼ = (Aqⱼ ∙ qᵢ) = qᵢᵀAqⱼ

q₁ q₂ ... q₄已經形成N維空間，q₅是多餘的，H₅₄ = 0也是多餘的。刪除之後即是方陣。

⎡  |   |   |   |  ⎤   ⎡  |  |  |  | ⎤ ⎡ H₁₁ H₁₂ H₁₃ H₁₄ ⎤
⎢ Aq₁ Aq₂ Aq₃ Aq₄ ⎥ = ⎢ q₁ q₂ q₃ q₄ ⎥ ⎢ H₂₁ H₂₂ H₂₃ H₂₄ ⎥
⎣  |   |   |   |  ⎦   ⎣  |  |  |  | ⎦ ⎢  0  H₃₂ H₃₃ H₃₄ ⎥
                                      ⎣  0   0  H₄₃ H₄₄ ⎦
        AQ                   Q                 H

遞推法得到q₁ q₂ ... q₄。隨機選擇初始向量x₀，扳正為q₁。Aq₁作為下一個向量，扳正為q₂。以此類推。

q₁ = x₀/‖x₀‖     normalization
for t = 1⋯N
   x = Aq_t
   x = x - H_1tq₁ - H_2tq₂ - ... - H_ttq_t     orthonormalization
   q_t+1 = x/‖x‖     normalization
   H_t+1,k = ‖x‖     length

遞推一次O(N²)，時間複雜度O(N³)。

幾乎上三角分解演算法（Lanczos iteration）

如果A是對稱矩陣，那麼H也是對稱矩陣。

導致H變成「三對角線矩陣tridiagonal matrix」。

⎡  |   |   |   |  ⎤   ⎡  |  |  |  | ⎤ ⎡ H₁₁ H₁₂  0   0  ⎤
⎢ Aq₁ Aq₂ Aq₃ Aq₄ ⎥ = ⎢ q₁ q₂ q₃ q₄ ⎥ ⎢ H₂₁ H₂₂ H₂₃  0  ⎥
⎣  |   |   |   |  ⎦   ⎣  |  |  |  | ⎦ ⎢  0  H₃₂ H₃₃ H₃₄ ⎥
                                      ⎣  0   0  H₄₃ H₄₄ ⎦
        AQ                   Q                 H

矩陣求值的部分，時間複雜度仍是O(N³)；向量扳正的部分，計算公式精簡成三項，從O(N²)降為O(N)，時間複雜度從O(N³)降為O(N²)。

對稱矩陣改用Wilkinson shift。

幾乎上三角分解演算法（Householder reflection）

實施N-1次相似變換，變換矩陣是鏡射矩陣。

遞推一次O(N²)，時間複雜度O(N³)。

幾乎上三角分解演算法（Jacobi rotation）

實施O(N²)次相似變換，變換矩陣是軸面旋轉矩陣。

僅適用於對稱矩陣。遞推一次O(N)，時間複雜度O(N³)。

矩陣相同的成立條件

一、特徵值相同、特徵向量相同。

二、特徵值與特徵向量的對應方式相同。（特徵對相同。）

三、特徵向量湊足N種。（一種方式是令特徵值皆異。）

同時滿足三點，那麼矩陣必定相同。只滿足其中幾點，那麼矩陣可能相同、可能不相同！舉一個範例：

⎡ 1 0 ⎤ and ⎡   1 0 ⎤ are not identical.
⎣ 1 1 ⎦     ⎣ 200 1 ⎦

same eigenvalues:
1 and 1

same eigenvector:
⎡ 0 ⎤
⎣ 1 ⎦   only one eigenvector

https://math.stackexchange.com/questions/2535288/
https://www.wolframalpha.com/input/?i={{1,0},{1,1}}
https://www.wolframalpha.com/input/?i={{1,0},{200,1}}

矩陣相似的成立條件

一、特徵值相同。

二、特徵向量湊足N種。（一種方式是令特徵值皆異。）

同時滿足兩點，那麼矩陣必定相似。只滿足其中幾點，那麼矩陣可能相似、可能不相似！舉一個範例：

⎡ 1 0 ⎤ and ⎡ 1 0 ⎤ are not similar.
⎣ 1 1 ⎦     ⎣ 0 1 ⎦

same eigenvalues:
1 and 1

https://math.stackexchange.com/questions/3330203/

eigenbasis - congruence

本章主角是相合矩陣：特徵值正負號數量相同

congruence transformation / congruent matrices

FAFᵀ = B稱作「相合變換」。A與B稱作「相合矩陣」。

也有人寫成PᵀAP = B，其中P = Fᵀ。

相似變換：FAF⁻¹ = B
相似矩陣：A與B特徵值相同
相合變換：FAFᵀ = B
相合矩陣：A與B特徵值正負號數量相同

所知甚少，只知道是二次型。

估計大量特徵向量（subspace iteration）

power iteration，一口氣處理大量向量。

大量向量併成矩陣X。

「向量長度縮放為1」改成了「矩陣扳正為Q」。

loop t times
   X := AX                  power iteration
   X := Q  (let X = QR)     orthonormalization

QR分解代價太大，改成只做一次，無傷大雅。

X := AᵗX                    power iteration
X := Q     (let X = QR)     orthonormalization

特徵向量們最後被扳正了，嚴格來說是估計特徵空間。

估計大量特徵值（Rayleigh–Ritz projection）

Rayleigh quotient，一口氣處理大量向量。

大量向量併成矩陣X。大量虛擬特徵值併成矩陣Λ。

對整個矩陣微分，令矩陣是I，以利微分。

minimize ‖AX - ΛX‖²   subject to XᵀX = I (X is orthonormal)

最佳解位於一次微分等於零的地方。

    XᵀAX
Λ = ———— = XᵀAX
    XᵀX

X是正規正交矩陣，寫成Q。Λ通常不是對角線矩陣，寫成B。由於Qᵀ = Q⁻¹，形成相似變換，特徵值相同。

Q⁻¹AQ = B

虛擬特徵值們最初被整併了，嚴格來說是估計虛擬特徵值矩陣。

大量特徵值演算法（subspace iteration）

事先實施subspace iteration和Rayleigh–Ritz projection，調整矩陣，讓矩陣接近正解，減少遞推次數。集大成。

X := AᵗX                      power iteration
X := Q       (let X = QR)     orthonormalization
B := Q₁⁻¹AQ₁ (let Q₁ = Q)     Rayleigh–Ritz projection
H := Q₂⁻¹BQ₂                  Hessenberg decomposition
T := Q₃⁻¹HQ₃                  Schur decomposition

順便整理一下

            QR iteration：A的QR分解拿來當作Q。
       Arnoldi iteration：AQ的QR分解移項當作Q。
Rayleigh–Ritz projection：任意正規正交矩陣拿來當作Q。

eigenbasis - representation

本章主角是可交換矩陣：特徵向量相同（嚴格來說是特徵空間相同）

commuting matrices

「可交換矩陣」。AB = BA。複合順序不影響結果。

A and B are commute: AB = BA

以矩陣變換的觀點來看，A是保B變換，B是保A變換。

|          |      |          |
| ———B———> |      | ———B⁻¹—> |     A = B⁻¹AB
A          A      A          A        and
| <——B⁻¹—— |      | <——B———— |     A = BAB⁻¹
↓          ↓      ↓          ↓

A commute with B => A is B-preserving transformation

順帶一提，A和B是對稱矩陣的情況下，AB是對稱矩陣、A和B是可交換矩陣，兩者等價。

A and B are symmetric. AB = BA iff AB is symmetric.

simultaneous diagonalization

「同步特徵分解」、「同步對角化」。A和B擁有相同的N種特徵向量。

A and B are simultaneous diagonalizable: A = EΛE⁻¹ and B = EΓE⁻¹

可同步特徵分解，可交換矩陣，兩者等價。伸縮方向（特徵向量）一致，伸縮倍率（特徵值）擁有交換律。

可同步喬登分解，似乎沒有特別意義。

AB = BA iff A = EΛE⁻¹ and B = EΓE⁻¹

特殊案例

兩個二維旋轉矩陣，可交換。三維以上則不可交換。
兩個縮放矩陣，可交換。
均勻縮放矩陣、旋轉矩陣，可交換。

縮放矩陣、旋轉矩陣，通常不可交換，除非遇到特例。
縮放矩陣、其他矩陣，通常不可交換，除非遇到特例。
旋轉矩陣、其他矩陣，通常不可交換，除非遇到特例。
投影矩陣、其他矩陣，通常不可交換，除非遇到特例。

convolution theorem

矩陣宛如數列，矩陣乘法宛如數列卷積，特徵值相乘宛如頻譜點積，同步特徵基底宛如傅立葉轉換。

  |  forward   |            |  forward   |
  | ————E⁻¹——→ |            | ————E⁻¹——→ |
  A            Λ            B            Γ     AB = E(ΛΓ)E⁻¹
  | ←———E————— |            | ←———E————— |
  ↓  backward  ↓            ↓  backward  ↓
 time       frequency      time       frequency
domain       domain       domain       domain

commuting matrices性質證明

https://math.stackexchange.com/questions/1736639/

Eigenbasis - orthogonal axes

本章主角是正規矩陣：特徵向量互相垂直

normal matrix

「（實數版本）正規矩陣」。A與Aᵀ可交換。AᵀA = AAᵀ。

「（複數版本）正規矩陣」。A與Aᴴ可交換。AᴴA = AAᴴ。

特徵向量互相垂直。真．正交矩陣。

註：normal matrix的normal應該是指法線。orthonormal matrix的normal應該是指長度為1。

orthonormal eigendecomposition

「正規正交特徵分解」。正規矩陣的特徵分解。

正規矩陣（包含了對稱矩陣、反對稱矩陣）一定可以特徵分解。其特徵向量互相垂直，再令特徵向量長度為一，如此一來特徵基底正規正交：反矩陣等於轉置矩陣。

A = EΛE⁻¹     when A is diagonalizable
A = EΛEᵀ      when A is real normal

normal matrix性質證明【尚待確認】

(⇒)正規矩陣，則特徵向量互相垂直。

一、轉置矩陣，特徵值相同。共軛轉置矩陣，特徵值共軛。

二、可交換矩陣，特徵向量相同。

1. eigval(A) = eigval(Aᵀ)
2. eigval(A) = eigval(Aᴴ)
3. AB = BA  ⟺  eigvec(A) = eigvec(B)

（實數版本）正規矩陣、及其轉置，特徵值相同，特徵向量相同。

（複數版本）正規矩陣、及其共軛轉置，特徵值共軛，特徵向量相同。

   AᵀA = AAᵀ                    A is real-version normal
=> ⎰ eigval(A) = eigval(Aᵀ)
   ⎱ eigvec(A) = eigvec(Aᵀ)
≠> A = Aᵀ                       A is real symmetric ???

   AᴴA = AAᴴ                    A is normal
=> ⎰ eigval(A) = eigval(Aᴴ)
   ⎱ eigvec(A) = eigvec(Aᴴ)

但是這樣沒有證明特徵值、特徵向量一一對應！

（實數版本）正規矩陣不一定是對稱矩陣。【尚待確認】

normal matrix性質證明

(⇒)正規矩陣，則特徵向量互相垂直。

正規矩陣、及其共軛轉置，向量經過變換，向量長度一致。

   AᴴA = AAᴴ
=> ‖Ax‖² = (Ax)ᴴ(Ax) = xᴴAᴴAx = xᴴAAᴴx = (Aᴴx)ᴴ(Aᴴx) = ‖Aᴴx‖²
=> ‖Ax‖ = ‖Aᴴx‖

正規矩陣、及其共軛轉置，特徵值共軛，特徵向量相同，一一對應。

   Ax = λx
=> Ax - λI = 0
=> (A - λI)x = 0
=> ‖(A - λI)x‖ = 0
=> ‖(A - λI)ᴴx‖ = 0
=> ‖(Aᴴ - λI)x‖ = 0
=> Aᴴx - λI = 0
=> Aᴴx = λx

https://math.stackexchange.com/questions/461958/

正規矩陣，若特徵值相異，則特徵向量互相垂直。

正規矩陣，若特徵值相同，大家習慣取互相垂直的特徵向量。

   ⎰ Ax₁ = λ₁x₁
   ⎱ Ax₂ = λ₂x₂

=> ⎰ Aᴴx₁ = λ₁x₁
   ⎱ Aᴴx₂ = λ₂x₂

=> ⎰ ❮Ax₁,x₂❯ = ❮λ₁x₁,x₂❯ = λ₁❮x₁,x₂❯
   ⎱ ❮Ax₁,x₂❯ = ❮x₁,Aᴴx₂❯ = ❮x₁,λ₂x₂❯ = λ₂❮x₁,x₂❯

=> since λ₁ ≠ λ₂ thus ❮x₁,x₂❯ = 0

https://math.stackexchange.com/questions/778946/

(⇐)特徵向量互相垂直，則是正規矩陣。

   A = EΛEᵀ
=> Aᵀ = EΛᵀEᵀ = EΛEᵀ
=> A = Aᵀ
=> AAᵀ = AᵀA

   A = UΛUᴴ
=> Aᴴ = UΛᴴUᴴ
=> ⎰ AᴴA = UΛᴴUᴴUΛUᴴ = UΛᴴΛUᴴ     where UᴴU = I
   ⎱ AAᴴ = UΛUᴴUΛᴴUᴴ = UΛΛᴴUᴴ
=> since ΛΛᴴ = ΛᴴΛ thus AᴴA = AAᴴ

normal matrix相關定理

正規矩陣、及其轉置，若特徵值相同，則特徵向量必定垂直。

⎰ A xᵢ = λᵢxᵢ  ⟹  xᵢᴴyⱼ = 0
⎱ Aᴴyᵢ = λᵢyᵢ
http://macs.citadel.edu/chenm/344.dir/08.dir/lect4_2.pdf

實數正規矩陣、實數特徵值，必是實數對稱矩陣。

https://math.stackexchange.com/questions/2991679/
https://math.stackexchange.com/questions/459207/

備忘

想要順便講一下特殊觀念，但是我還沒想到如何整合思路。

A and Aᵀ are similar (have same eigenvalue)
AᵀA and AAᵀ are similar (have same eigenvalue)
https://math.stackexchange.com/questions/1726932/

eigenbasis - symmetry

本章主角是對稱矩陣：特徵向量互相垂直、特徵值為實數

矩陣轉置

實數的反面：變號。

虛數的反面：共軛。

矩陣的反面：轉置。

symmetric matrix / Hermitian matrix

「對稱矩陣」是實數版本。A = Aᵀ。

「共軛對稱矩陣」是複數版本。A = Aᴴ。

特徵向量互相垂直，特徵值是實數。真．縮放矩陣。

skew-symmetric matrix / skew-Hermitian matrix

「反對稱矩陣」是實數版本。A = -Aᵀ。

「反共軛對稱矩陣」是複數版本。A = -Aᴴ。

特徵向量互相垂直，特徵值是虛數。真．旋轉矩陣。

變號、共軛、轉置，三個通通用上。

symmetric and skew-symmetric decomposition
（even–odd decomposition）

「對稱與反對稱分解」。一個矩陣拆成兩個矩陣相加：對稱矩陣（真．縮放矩陣）、反對稱矩陣（真．旋轉矩陣）。也有複數版本。

A = S + R
S = (A + Aᵀ) / 2   symmetric part (even part) (scale)
R = (A - Aᵀ) / 2   skew-symmetric part (odd part) (rotate)

即是「奇偶分解」。變號是主角，共軛和轉置是配角。

複數的實部與虛部：原本複數加減共軛複數再除以二。
函數的偶函數與奇函數：原本函數加減鏡射函數再除以二。
線性函數的的偶部與奇部：原本矩陣加減共軛轉置矩陣再除以二。

眼尖的讀者應該注意到了，虛擬特徵向量就是偶部(A+Aᵀ)/2。

Hermitian matrix性質證明

我只證明複數版本。實數版本如法炮製。

(⇒)共軛對稱矩陣，則特徵值是實數。

⎰ A = Aᴴ      Hermitian matrix
⎱ Ax = λx     eigenvector and eigenvalue

=> ‖Ax‖² = (Ax)ᴴ(Ax) = xᴴAᴴAx = xᴴA²x = xᴴ(λ²x) = λ²‖x‖²
=> λ² = ‖Ax‖² / ‖x‖²
=> λ² is real and λ² ≥ 0
=> the square root λ is real 

https://math.stackexchange.com/questions/354115/

⎰ A = Aᴴ
⎱ eigval(A) = eigval(Aᴴ)

=> eigval(A) = eigval(Aᴴ) = eigval(Aᴴ)
=> imaginary part of eigval(Aᴴ) is 0
=> imaginary part of eigval(A) is 0

https://math.stackexchange.com/questions/2592295/

(⇐)特徵向量互相垂直（正規矩陣），特徵值是實數，則是共軛對稱矩陣。

   A = UΛUᴴ
=> Aᴴ = UΛᴴUᴴ
=> Aᴴ = UΛᴴUᴴ = UΛUᴴ     since Λ is real
=> A = Aᴴ

備忘

想要順便講一下特殊觀念，但是我還沒想到如何整合思路。

1. An odd sized skew-symmetric matrix cannot be invertible.
2. skew-symmetric matrix A exhibits xᵀAx = 0.

eigenbasis - conjugate

本章主角是正定矩陣：特徵值為正數

矩陣共軛分解【目前稱作conjugate decomposition】

複數的共軛分解：分解成複數乘以共軛複數。

線性函數的共軛分解：分解成矩陣乘以轉置矩陣。

注意到，矩陣共軛分解⟌A = √ΛEᵀ、矩陣平方根√A = E√ΛE⁻¹，兩者不同。

矩陣共軛分解目前沒有正式數學符號。

矩陣共軛分解【目前稱作conjugate decomposition】

「長度平方」按照慣例必須是正值或零，不可為負值。我們姑且將「矩陣長度平方是正值或零」定義成「特徵值全是正值或零」。

        def
A = MᵀM === ‖M‖²   A cannot have negative eigenvalues.

我們希望矩陣A是「矩陣長度平方」、M與Mᵀ是「矩陣共軛對」。然而A的特徵值可能為負數、M的特徵值可能為複數。我們必須追加限制條件，讓事情變得合理。

normal matrix

「正規矩陣」。特徵向量互相垂直。

正規矩陣擁有正規正交特徵基底。

A = EΛE⁻¹ = EΛEᵀ

正規矩陣可以分解成「矩陣自身內積」或「矩陣自身外積」。

A = EΛE⁻¹ = EΛEᵀ = E√Λ√ΛEᵀ = E√Λᵀ√ΛEᵀ = (√ΛEᵀ)ᵀ(√ΛEᵀ) = MᵀM
A = EΛE⁻¹ = EΛEᵀ = E√Λ√ΛEᵀ = E√Λ√ΛᵀEᵀ = (E√Λ)(E√Λ)ᵀ = NNᵀ

symmetric matrix

「對稱矩陣」。特徵向量互相垂直、特徵值全是實數。

對稱矩陣是正規矩陣的特例，也可以分解成「矩陣自身內積」。

positive definite matrix / positive semidefinite matrix

「正定矩陣」。虛擬特徵值全是正數。特徵值當然也是。

「半正定矩陣」。虛擬特徵值全是正數或零。特徵值當然也是。

    xᵀAx
λ = ————— ≥ 0   when x ≠ 0     虛擬特徵值全是正數或零
     xᵀx

數學式子簡化成xᵀAx > 0和xᵀAx ≥ 0。

xᵀAx ≥ 0        when x ≠ 0     同乘xᵀx。元素平方和，必定為正數，不必變號。

數學符號標記成A ≻ 0和A ≽ 0。

A ≽ 0                          positive semidefinite

symmetric positive semidefinite matrix

「對稱半正定矩陣」。同時具備對稱和半正定兩種性質。

可以分解成「矩陣自身內積」，而且是實數矩陣。

可以改寫成「矩陣長度平方」，而且是實數矩陣。

        def
A = MᵀM === ‖M‖²     when A is symmetric positive semidefinite

A = (√ΛEᵀ)ᵀ(√ΛEᵀ) = MᵀM   where M is real
A = ‖√ΛEᵀ‖² = ‖M‖²        where M is real

特徵值是非負實數，特徵值開根號仍是非負實數，特徵向量是實數，因而得到實數矩陣M = √ΛEᵀ。

conjugate decomposition【尚無正式名稱】

對稱半正定矩陣的共軛分解，答案無限多個。常見解法：特徵分解、Cholesky分解。

eigendecomposition:             Cholesky decomposition:
A = EΛEᵀ                        A = LLᵀ
A = E√Λ√ΛEᵀ = (√ΛEᵀ)ᵀ(√ΛEᵀ)     A = (Lᵀ)ᵀ(Lᵀ)

eigenbasis - absolute value🚧

本章主角是奇異值分解

singular value decomposition（SVD）

無法特徵分解的矩陣，可以先平方再開根號，硬是分解。

大家不使用平方A²，而是使用內積AᵀA、外積AAᵀ：

一、兩者都是對稱矩陣：特徵向量正規正交、特徵值是實數。

二、兩者都是半正定矩陣：特徵值是正數或零。

AᵀA = VΛVᵀ = VΣᵀΣVᵀ = VΣᵀUᵀUΣVᵀ = (UΣVᵀ)ᵀ(UΣVᵀ)
AAᵀ = UΛUᵀ = UΣΣᵀUᵀ = UΣVᵀVΣᵀUᵀ = (UΣVᵀ)(UΣVᵀ)ᵀ
A = UΣVᵀ          (Λ = ΣᵀΣ = ΣΣᵀ = Σ² , √Λ = Σ)

硬湊出A = UΣVᵀ，稱作「奇異值分解」。

一、U是AAᵀ特徵向量，V是AᵀA特徵向量。U每個直條稱作「左奇異向量」，V每個直條稱作「右奇異向量」。正規正交。

二、Σ是AAᵀ特徵值開根號，也是AᵀA特徵值開根號。Σ每個對角線元素稱作「奇異值」。正數或零。

「奇異值分解」是「特徵分解取絕對值」。SVD = |EVD|。

複數的長度平方（共軛相乘）的平方根，複數變成了非負實數。

複數矩陣的長度平方（內積與外積）的平方根，複數特徵值變成了非負實數奇異值，複數特徵向量變成了正規正交實數奇異向量。

旋轉矩陣（虛數特徵向量）、歪斜矩陣（缺少特徵向量）、不是方陣（沒有特徵向量），無論哪種矩陣，通通可以奇異值分解！

singular vector與singular value

預條件矩陣是轉置矩陣，實施投影，調整維度。

左預條件，輸出降維。右預條件，輸入升維。

left preconditioner Aᵀ and right singular vector v.
Aᵀ(Av) = Aᵀ(σv)     where Aᵀ(σv) = σ(Aᵀv) = σ(σv) = σ²v

right preconditioner Aᵀ and left singular vector u.
A(Aᵀu) = σ(Aᵀu)     where σ(Aᵀu) = σ(σu) = σ²u

QR transformation / RQ transformation【尚無正式名稱】

相似變換保留特徵值，QR變換、RQ變換保留奇異值。

P⁻¹AP = B                         Similar transformation
Qᵀ(AAᵀ)Q = (RRᵀ)  --->  QᵀA = R   QR transformation
Q(AᵀA)Qᵀ = (RᵀR)  --->  AQᵀ = R   RQ transformation

Golub–Kahan bidiagonalization

「雙對角線矩陣」。主對角線、次對角線（通常是上方）以外必須為零。

    ⎡ B₁₁ B₁₂  0   0   0  ⎤
    ⎢  0  B₂₂ B₂₃  0   0  ⎥
B = ⎢  0   0  B₃₃ B₃₄  0  ⎥     bidiagonal matrix
    ⎢  0   0   0  B₄₄ B₄₅ ⎥
    ⎣  0   0   0   0  B₅₅ ⎦

「雙對角線分解」。實施QR變換，變成雙對角線矩陣。左右矩陣必須是（實數版本）正規正交矩陣、（複數版本）么正矩陣。

UᵀAV = B     Golub–Kahan bidiagonalization

任意矩陣一定可以雙對角線分解。

雙對角線分解演算法，其實就是QR分解那三個演算法的改良版，QR分解改成了交互QR變換、RQ變換。

奇異值分解演算法，可以分成兩步驟：第一步先變成雙對角線矩陣，第二步再變成對角線矩陣。時間複雜度較低。

Golub–Kahan bidiagonalization: UᵀAV = B
Schur decomposition:           Qᵀ(BᵀB)Q = Σ²

雙對角線分解演算法（Golub–Kahan–Lanczos method）

雙對角線分解演算法（Householder reflection）

⎡ X X X X ⎤     ⎡ X X X X ⎤     ⎡ X X     ⎤     ⎡ X X     ⎤     ⎡ X X     ⎤
⎢ X X X X ⎥ U₁ᵀ ⎢   X X X ⎥ V₁  ⎢   X X X ⎥ U₂ᵀ ⎢   X X X ⎥ V₂  ⎢   X X   ⎥
⎢ X X X X ⎥ ——→ ⎢   X X X ⎥ ——→ ⎢   X X X ⎥ ——→ ⎢     X X ⎥ ——→ ⎢     X X ⎥ ...
⎢ X X X X ⎥     ⎢   X X X ⎥     ⎢   X X X ⎥     ⎢     X X ⎥     ⎢     X X ⎥
⎢ X X X X ⎥     ⎢   X X X ⎥     ⎢   X X X ⎥     ⎢     X X ⎥     ⎢     X X ⎥
⎣ X X X X ⎦     ⎣   X X X ⎦     ⎣   X X X ⎦     ⎣     X X ⎦     ⎣     X X ⎦
     A             U₁ᵀA            U₁ᵀAV₁        U₂ᵀU₁ᵀAV₁      U₂ᵀU₁ᵀAV₁V₂

雙對角線分解演算法（Jacobi rotation）

嗯哼。喔嗬。哎呀。【待補文字】

https://people.eecs.berkeley.edu/~demmel/
https://www.fernuni-hagen.de/MATHPHYS/veselic/
Jacobi's method is more accurate than QR   (1989)
New fast and accurate Jacobi SVD algorithm (2007)

one-side

AV = W      make V orthonormal and W orthogonal by Jacobi rotation
AV = DU     normalization of each vector in W

奇異值分解演算法（Golub–Kahan–Reinsch algorithm）

兩階段。

⎡ X X X X ⎤     ⎡ X X     ⎤     ⎡ X       ⎤
⎢ X X X X ⎥     ⎢   X X   ⎥     ⎢   X     ⎥
⎢ X X X X ⎥ ——→ ⎢     X X ⎥ ——→ ⎢     X   ⎥
⎢ X X X X ⎥     ⎢       X ⎥     ⎢       X ⎥
⎢ X X X X ⎥     ⎢         ⎥     ⎢         ⎥
⎣ X X X X ⎦     ⎣         ⎦     ⎣         ⎦
     A               B               Σ

奇異值分解演算法（Lawson–Hanson–Chan algorithm）

QR上三角分解，時間複雜度高，但是可以快速清空下方。

三階段。事先清空下方，只保留方陣，進行後續計算。

⎡ X X X X ⎤     ⎡ X X X X ⎤     ⎡|X̅‾X̅‾‾‾‾|⎤       ⎡ X̅‾‾‾‾‾‾|⎤
⎢ X X X X ⎥     ⎢   X X X ⎥     ⎢|  X X  |⎥       ⎢ |  X    |⎥
⎢ X X X X ⎥ ——→ ⎢     X X ⎥ ——→ ⎢|    X X|⎥  ——→  ⎢ |    X  |⎥
⎢ X X X X ⎥     ⎢       X ⎥     ⎢|      X|⎥       ⎢ |      X|⎥
⎢ X X X X ⎥     ⎢         ⎥     ⎢ ‾‾‾‾‾‾‾ ⎥       ⎢  ‾‾‾‾‾‾‾ ⎥
⎣ X X X X ⎦     ⎣         ⎦     ⎣         ⎦       ⎣          ⎦
     A              QᵀA            UᵀQᵀAV            Σ

四階段。根據矩陣長寬比例，決定QR上三角分解的使用範圍。

⎡ X X X X ⎤   ⎡ X X     ⎤   ⎡ X X     ⎤   ⎡ X X     ⎤   ⎡ X       ⎤
⎢ X X X X ⎥   ⎢   X X X ⎥   ⎢  |X̅‾X̅‾X̅|⎥   ⎢  |X̅‾X̅‾‾|⎥   ⎢   X     ⎥
⎢ X X X X ⎥ → ⎢   X X X ⎥ → ⎢  |  X X|⎥ → ⎢  |  X X|⎥ → ⎢     X   ⎥
⎢ X X X X ⎥   ⎢   X X X ⎥   ⎢  |    X|⎥   ⎢  |    X|⎥   ⎢       X ⎥
⎢ X X X X ⎥   ⎢   X X X ⎥   ⎢  |     |⎥   ⎢   ‾‾‾‾‾ ⎥   ⎢         ⎥
⎣ X X X X ⎦   ⎣   X X X ⎦   ⎣  |_____|⎦   ⎣         ⎦   ⎣         ⎦
     A           U₁ᵀAV₁       QᵀU₁ᵀAV₁    U₂ᵀQᵀU₁ᵀAV₁V₂      Σ

eigenbasis - product🚧

本章主角是共相關矩陣：左右奇異向量

correlation matrix

「共相關矩陣」。矩陣外積。XYᵀ。

AᵀA瘦奇異值分解演算法（NIPALS）

註：AᵀA是對稱半正定矩陣，特徵分解恰好等於奇異值分解。

power iteration求特徵向量，Hotelling deflation消滅該維度。

每回合先代入A、再代入Aᵀ，不必事先求得AᵀA。

時間複雜度O(RCT)，R是矩陣高度，C是矩陣寬度，T是遞推次數。

若事先求得AᵀA，矩陣乘法O(CRC)，奇異值分解O(CCT)。當T遠小於C，當R遠小於C，則適合採用此演算法。換句話說，扁平矩陣適合採用此演算法。

加速技巧：檢查v是否收斂，改為檢查t是否收斂。v是矩陣寬度、t是矩陣高度，選擇較小者，以減少檢查時間。

加速技巧的缺點：t沒有正規化，數值範圍較大。檢查t是否收斂，浮點數精確度必須調低，以容忍誤差，以避免無窮迴圈。

XᵀY瘦奇異值分解演算法（NIPALS）

每回合輪流套用XᵀY和YᵀX，得到左右奇異向量。

matrix normalize ⧚A⧛ =  UVᵀ

repeat rank(XᵀY) times
    v := the longest row vector of Y
    repeat until y converges
        u := ⧚Xᵀ(Yv)⧛               // power iteration of XᵀY
        v := ⧚Yᵀ(Xu)⧛               // power iteration of YᵀX
    end
    U := U ∪ u                      // left singular vector
    V := V ∪ v                      // right singular vector
    s := s ∪ ‖Xᵀ(Yv)‖               // singular value
    X := X - (Xu)uᵀ                 // Hotelling deflation of X
    Y := Y - (Yv)vᵀ                 // Hotelling deflation of Y
end

u = ⧚Xᵀ(Yv)⧛ = ⧚Xᵀ(Y⧚Yᵀ(Xu)⧛)⧛ = ⧚XᵀYYᵀXu⧛ = ⧚(XᵀYYᵀX)u⧛
u is eigenvector of XᵀYYᵀX = (XᵀY)(XᵀY)ᵀ = AAᵀ
u is left singular vector of XᵀY = A

v = ⧚Yᵀ(Xu)⧛ = ⧚Yᵀ(X⧚Xᵀ(Yv)⧛)⧛ = ⧚YᵀXXᵀYv⧛ = ⧚(YᵀXXᵀY)v⧛
v is eigenvector of YᵀXXᵀY = (XᵀY)ᵀ(XᵀY) = AᵀA
v is right singular vector of XᵀY = A

void NIPALS()
{
	MatrixXf X(3, 5);
	X << 1, 4, 5, 3, 7,
		 2, 5, 0, 3, 5,
		 3, 6, 4, 3, 9;

MatrixXf Y(3, 7);
	Y << 4, 16,  9,  3,  8, 17, 23,
		 5,  4,  8, 16, 14, 11,  9,
		 8, 16, 14, 11,  9, 19,  1;

int rank = min(X.rows(), X.cols(), Y.cols());
	MatrixXf U(X.cols(), rank);
	MatrixXf V(Y.cols(), rank);
	VectorXf s(rank);

for (int i=0; i<rank; ++i)
	{
		const float ε = 1e-15;
		MatrixXf u(X.cols(), 1);
		MatrixXf v(Y.cols(), 1);
		MatrixXf x(X.rows(), 1);
		MatrixXf y(Y.rows(), 1);
		MatrixXf v0(Y.cols(), 1);

// 找到長度最長的向量，當作初始向量。以避免零向量。
		float max = -1e9, pos = 0;
		for (int j=0; j<Y.rows(); ++j)
		{
			float length2 = Y.row(j).squaredNorm();
			if (length2 > max) {max = length2; pos = j;}
		}
		if (max <= ε) {rank = i; break;}
		v = Y.row(pos).transpose().normalized();

// power iteration：尋找奇異向量暨奇異值。
		do
		{
			v0 = v;

y = Y * v;
			u = X.transpose() * y;
			u.normalize();
			x = X * u;
			v = Y.transpose() * x;
			v.normalize();
		} while ((v - v0).squaredNorm() > ε);

// 紀錄奇異向量暨奇異值
		U.col(i) = u;
		V.col(i) = v;
		s(i) = (X.transpose() * y).norm();

// Hotelling deflation：消滅該維度。
		X -= x * u.transpose();
		Y -= y * v.transpose();
	}

U.conservativeResize(Eigen::NoChange, rank);
	V.conservativeResize(Eigen::NoChange, rank);
	s.conservativeResize(rank);
	cout << U << V << s;
}

使用加速技巧的版本。NIPALS的相關文獻通常記載此版本。

repeat rank(XᵀY) times
    y := the longest column vector of Y
    repeat until y converges
        x := X ⧚Xᵀy⧛
        y := Y ⧚Yᵀx⧛
    end
    U := U ∪ ⧚Xᵀy⧛                  // left singular vector
    V := V ∪ ⧚Yᵀx⧛                  // right singular vector
    s := s ∪ ‖Xᵀy‖                  // singular value
    X := X - x ⧚Xᵀy⧛ᵀ               // Hotelling deflation of X
    Y := Y - y ⧚Yᵀx⧛ᵀ               // Hotelling deflation of Y
end

eigenbasis - reciprocal length🚧

本章主角是虛擬反矩陣：奇異值倒數

inverse

「反矩陣」。可以特徵分解：特徵值通通倒數。只能喬登分解：喬登矩陣求反矩陣。

必須避免1/0 = ∞。特徵值不為0，才能倒數，才有反矩陣。

A⁻¹ = (EΛE⁻¹)⁻¹ = E⁻¹⁻¹Λ⁻¹E⁻¹ = EΛ⁻¹E⁻¹
A⁻¹ = (EJE⁻¹)⁻¹ = E⁻¹⁻¹J⁻¹E⁻¹ = EJ⁻¹E⁻¹

奇異值分解：左右奇異向量交換並轉置，奇異值通通倒數。

A⁻¹ = (UΣVᵀ)⁻¹ = Vᵀ⁻¹Σ⁻¹U⁻¹ = VΣ⁻¹Uᵀ

pseudoinverse

「虛擬反矩陣」、「偽反矩陣」。

奇異值分解：擱置奇異值0，不算倒數。有多少rank，做多少事。

瘦奇異值分解：截斷奇異值0暨奇異向量。

A⁺ = (UΣVᵀ)⁺ = Vᵀ⁺Σ⁺U⁺ = VΣ⁺Uᵀ

可以使用MATLAB的pinv。

condition number

「條件數」。κ(A) = ‖A‖ ‖A⁻¹‖ = σₘₐₓ / σₘᵢₙ。

相似變換的大致的縮放倍率。用來判斷演算法收斂速度，例如QR iteration。

truncated singular value decomposition

「截斷奇異值分解」。刪除奇異值0，刪除對應的奇異向量，縮小矩陣。

原始版本暱稱「滿奇異值分解full SVD」，截斷版本暱稱「瘦奇異值分解thin SVD」。

https://ccjou.wordpress.com/2009/09/01/

外積形式，刪除奇異值0的項。

svd/evd cross version  sum sUiVi*/sEiEi*

瘦奇異值分解演算法（eigendecomposition）

問題化作特徵分解。大可不必計算AᵀA和AAᵀ。

⎡ O Aᵀ⎤ = 1/_ ⎡ V  V ⎤ ⎡ Σ  O ⎤ 1/_ ⎡ V V ⎤ᵀ
⎣ A O ⎦   /√2 ⎣ U -U ⎦ ⎣ O -Σ ⎦ /√2 ⎣ U U ⎦

eigenbasis - normalization🚧

本章主角是單元矩陣：長度為一

unit matrix【尚無正式名稱】

「單位矩陣」。長度平方等於恆等矩陣I。

矩陣長度倒數【尚無正式名稱】

矩陣長度倒數：矩陣自身內積、開根號、反矩陣。

‖A‖²  =  AᵀA   = VΣ²Vᵀ
‖A‖   = √AᵀA   = VΣVᵀ
‖A‖⁻¹ = √AᵀA⁻¹ = VΣ⁻¹Vᵀ

矩陣長度倒數的存在條件【尚無正式名稱】

長度不可為零（奇異值均不為零），才擁有倒數。

det(‖A‖) ≠ 0 ⇒ ‖A‖⁻¹。

如果A是方陣，結論比較簡單。A必須有反矩陣。

A⁻¹ ⇒ ‖A‖⁻¹ when a = b。

如果A是矩陣，需要考慮維度。A直條線性獨立，且A橫條多於等於直條。

rank(A) = b ⇒ ‖A‖⁻¹。

  rank(A) = b       A:a×b
⇒ rank(AᵀA) = b     AᵀA:b×b
⇒ (AᵀA)⁻¹
⇒ √AᵀA⁻¹
⇒ ‖A‖⁻¹

改成虛擬反矩陣，允許奇異值有零，不必斤斤計較維度。

‖A‖⁺ = √AᵀA⁺ = VΣ⁺Vᵀ

矩陣正規化【尚無正式名稱】

矩陣正規化：矩陣除以自身長度。矩陣乘以自身長度倒數。

兩種定義方式：先分子後分母、先分母後分子。因為後者答案比較漂亮、用途比較多元，所以我採用後者。

 A      A
——— = ————— = √AᵀA⁻¹A = VΣ⁻¹VᵀUΣVᵀ         ✘
‖A‖   √AᵀA

 A      A
——— = ————— = A√AᵀA⁻¹ = UΣVᵀVΣ⁻¹Vᵀ = UVᵀ   ✔
‖A‖   √AᵀA

兩種寫作風格：內積風格、外積風格。採用瘦奇異值分解，使得矩陣們順利相消。

⧚A⧛ = A√AᵀA⁻¹ = UΣVᵀVΣ⁻¹Vᵀ = UVᵀ   ✔
⧚A⧛ = √AAᵀ⁻¹A = UΣ⁻¹UᵀUΣVᵀ = UVᵀ   ✔

因為大家從未討論矩陣正規化，所以我只好自創數學符號。

⧚A⧛ = A/‖A‖ = A/√AᵀA
⧚A⧛ = A‖A‖⁻¹ = A√AᵀA⁻¹

矩陣正規化，目前沒有高速演算法。等你幫忙翻身！

矩陣正規化的性質

方陣：自身內積暨自身外積皆是恆等矩陣。

當維度不足，1可能短少，數量只有rank(A)。

⧚A⧛ᵀ⧚A⧛ = ⧚A⧛⧚A⧛ᵀ = I when a = b。

高瘦矩陣：當直條少於等於橫條，自身內積是恆等矩陣，使得直條正規正交。

⧚A⧛ᵀ⧚A⧛ = I_b×b when a ≥ b。

扁平矩陣：當直條多於等於橫條，自身外積是恆等矩陣，使得橫條正規正交。

⧚A⧛⧚A⧛ᵀ = I_a×a when a ≤ b。

奇異值變1（倍率），奇異向量不變（旋轉暨鏡射）。

A = UΣVᵀ。⧚A⧛ = UVᵀ。

先轉置再正規化＝先正規化再轉置。

⧚Aᵀ⧛ = ⧚A⧛ᵀ。

矩陣正規化的用途

白化。主成分分析。離差矩陣正規化。⧚X̃⧛。

⧚X̃⧛ = √X̃X̃ᵀ⁻¹X̃ = UΣ⁻¹UᵀX     let X̃ = X(I - 𝟏/n) = XC

X到Y的最近似仿射保距變換。⧚YXᵀ⧛。

可以想成：除以X，乘以Y。消去倍率，保留旋轉暨鏡射。

data is row
C = YᵀX = UΣVᵀ
Q = ⧚YᵀX⧛ = UVᵀ

data is column
C = YXᵀ = UΣVᵀ
Q = ⧚YXᵀ⧛ = UVᵀ

順帶一提，當X與Y已經正規化，變換矩陣即是YXᵀ = UVᵀ。

Q = YXᵀ = UVᵀ  when X and Y are normalized

矩陣版本的cosine。先正規化再內積。【尚待確認】

  cos(X,Y)
= ⧚X⧛ᵀ⧚Y⧛
= (X/‖X‖)ᵀ(Y/‖Y‖)
= (X‖X‖⁻¹)ᵀ(Y‖Y‖⁻¹)
= ‖X‖⁻¹ XᵀY ‖Y‖⁻¹
:= XᵀY / ‖X‖‖Y‖

一次迴歸的預測結果是y⧚Xᵀ⧛⧚Xᵀ⧛ᵀ。【尚待確認】

一次迴歸也許是矩陣版本的cosine。

data is row  (X:nxd, y:nx1, w:nx1), column is linearly independent
w = (XᵀX)⁻¹Xᵀy
Xw = X(XᵀX)⁻¹Xᵀy = X√XᵀX⁻¹√XᵀX⁻¹Xᵀy = ⧚X⧛⧚X⧛ᵀy
error = y - Xw = y - ⧚X⧛⧚X⧛ᵀy = (I - ⧚X⧛⧚X⧛ᵀ)y

data is column, coefficent is row vector (X:dxn, y:1xn, w:1xn)
wᵀ = (XXᵀ)⁻¹Xyᵀ
Xᵀwᵀ = Xᵀ(XXᵀ)⁻¹Xyᵀ
wX = yXᵀ(XXᵀ)⁻¹X = yXᵀ√XXᵀ⁻¹√XXᵀ⁻¹X = y⧚Xᵀ⧛⧚Xᵀ⧛ᵀ
      ^^^^^^^^^^
      quadratic form: metric ‖Xᵀ‖⁻² ?
error = y - wX = y - y⧚Xᵀ⧛⧚Xᵀ⧛ᵀ = y(I - ⧚Xᵀ⧛⧚Xᵀ⧛ᵀ)

y⧚Xᵀ⧛⧚Xᵀ⧛ᵀ = y⧚X⧛ᵀ⧚X⧛ = y‖⧚X⧛‖² = y   when X is squared matrix

polar decomposition

A = QP。Q是正規正交矩陣。P是對稱半正定矩陣。

Q是旋轉暨鏡射，P是伸展。P設定為半正定，使得鏡射（縮放倍率為負值）放在Q。

Q與P宛如角度與長度，宛如極座標，因而得名「極分解」。

極分解可以用奇異值分解硬湊而得。

A = UΣVᵀ= U(VᵀV)ΣVᵀ = (UVᵀ)(VΣVᵀ) = QP
let Q = UVᵀ, P = VΣVᵀ

極分解可以用正規化硬湊而得。

A = UΣVᵀ = QP
let Q = ⧚A⧛ = UVᵀ, P = ⧚A⧛⁻¹A

X = USV*
U是主成分  移項U^-1=U*是投影矩陣  U*X是投影量
SV*是投影量:資料與主成分的夾角cosine並且乘上長度
S是X不相關化之後(即是U*X)，維度(橫條)每項平方和
SVD其實就是cosine

let Q0 = M, Qi+1 = (Qi + Qi⁻ᵀ)/2 until Qi+1 - Qi ≈ 0
Newton algorithm for the square root of I
converges quadratically when Qi is nearly orthogonal.

eigenbasis - polynomial🚧

本章主角是同伴矩陣：特徵向量是特徵值各種次方

Krylov matrix

唸作krill-ˋ off而非ˋ cry-love。

    ⎡  |   |        |   ⎤
K = ⎢ A⁰x A¹x ... Aⁿ⁻¹x ⎥
    ⎣  |   |        |   ⎦

矩陣A。初始向量x。每次遞推結果，併成一個方陣。

初始向量x又稱作種子seed。亂數演算法也有種子這個詞彙。

靈感源自power iteration。用來求特徵值。

Krylov subspace

Krylov matrix的所有向量，構成的空間。

span(K) = span(A⁰x, A¹x, ..., Aⁿ⁻¹x)

一、考慮A的特徵基底，推導Krylov subspace的維度。

初始向量是K個特徵向量的線性組合，那麼Krylov subspace也是這K個特徵向量的線性組合。

特徵值是零，那麼Krylov subspace維度短少。

特徵值相同，那麼Krylov subspace維度短少。

二、考慮A的喬登標準型，推導Krylov subspace的維度。

初始向量涵蓋的特徵向量，落在哪些Jordan block裡面。每個Jordan block藉由移位，生成其所有維度。【尚待確認】

Krylov subspace維度短少

Krylov subspace維度短少原因：

一、初始向量的特徵向量短少。

二、特徵值是零、特徵值相同，其特徵向量無法分離出來。

Krylov subspace維度短少其實不太常見：

一、隨機初始向量，通常涵蓋所有特徵向量。

二、隨機矩陣，通常特徵值皆異、特徵值不是零。

Krylov basis / orthonormal Krylov basis

Krylov matrix的所有向量，當作基底，甚至扳正。

Krylov matrix經過下述操作，Krylov subspace仍相同。

一、以QR分解扳正向量。

K = QR  ⇒  span(K) = span(Q)

二、以扳正後的向量，即時生成下一個向量。

     ⎡  |    |          |     ⎤
K' = ⎢ A⁰q₀ A¹q₁ ... Aⁿ⁻¹qₙ₋₁ ⎥  ⇒  span(K) = span(K')
     ⎣  |    |          |     ⎦

q₀ = normalize(x) = x/‖x‖
q₁ = normalize(Aq₀ - proj_q₀(Aq₀))
q₂ = normalize(Aq₁ - proj_q₀,q₁(Aq₁))
 :       :

三、特徵值同加一數、矩陣乘上倍率。

eigenvalue shift: A + αI
matrix scaling: kA

power iteration、Arnoldi iteration、conjugate gradient descent、generalized minimum residual method，儘管細節不同，但是均可使用Krylov subspace進行演算法分析。

Krylov subspace進階主題

enlarged Krylov subspace
https://who.rocq.inria.fr/Laura.Grigori/
https://people.eecs.berkeley.edu/~demmel/

rational Krylov subspace
http://guettel.com/rktoolbox/examples/html/index.html

【查無正式名稱】

任意矩陣實施相似變換，Krylov matrix作為變換矩陣，恰是「同伴矩陣」。

    ⎡  |   |        |   ⎤
K = ⎢ A⁰x A¹x ... Aⁿ⁻¹x ⎥
    ⎣  |   |        |   ⎦

     ⎡  |   |       |  ⎤
AK = ⎢ A¹x A²x ... Aⁿx ⎥
     ⎣  |   |       |  ⎦

       ⎡  |  |        |   ⎤
AK = K ⎢ e₂ e₃ ... K⁻¹Aⁿx ⎥
       ⎣  |  |        |   ⎦

       ⎡ 0 0 ... 0 c₀   ⎤
       ⎢ 1 0 ... 0 c₁   ⎥
AK = K ⎢ 0 1 ... 0 c₂   ⎥     where c = K⁻¹Aⁿx
       ⎢ : :     : :    ⎥
       ⎣ 0 0 ... 1 cₙ₋₁ ⎦

AK = KC
K⁻¹AK = C   similarity transformation

companion matrix

同伴矩陣：特徵多項式化作同伴矩陣，使用Horner's rule。特徵多項式的係數，化作矩陣最右直條。特徵多項式的根，化作矩陣的特徵值。

det(λI - C) = p(λ)
det(C - λI) = (-1)ⁿ p(λ)

    ⎡  0  0  0  0 -c₀ ⎤      ⎡ -c₄ -c₃ -c₂ -c₁ -c₀ ⎤
    ⎢  1  0  0  0 -c₁ ⎥      ⎢  1   0   0   0   0  ⎥
C = ⎢  0  1  0  0 -c₂ ⎥  or  ⎢  0   1   0   0   0  ⎥
    ⎢  0  0  1  0 -c₃ ⎥      ⎢  0   0   1   0   0  ⎥
    ⎣  0  0  0  1 -c₄ ⎦      ⎣  0   0   0   1   0  ⎦

  det(λI - C)

      ⎡  λ  0  0  0   c₀ ⎤
      ⎢ -1  λ  0  0   c₁ ⎥
= det ⎢  0 -1  λ  0   c₂ ⎥
      ⎢  0  0 -1  λ   c₃ ⎥
      ⎣  0  0  0 -1 λ+c₄ ⎦

              ⎡ *  *  *  *   *  ⎤
              ⎢ *  λ  0  0   c₁ ⎥
= λ (-1)⁰ det ⎢ * -1  λ  0   c₂ ⎥
              ⎢ *  0 -1  λ   c₃ ⎥
              ⎣ *  0  0 -1 λ+c₄ ⎦
                ⎡  *  *  *  *  * ⎤
                ⎢ -1  λ  0  0  * ⎥
 + c₀ (-1)⁴ det ⎢  0 -1  λ  0  * ⎥
                ⎢  0  0 -1  λ  * ⎥
                ⎣  0  0  0 -1  * ⎦

= det(λI - Cꜜ) ⋅ λ + c₀

= ((((λ + c₄) ⋅ λ + c₃) ⋅ λ + c₂) ⋅ λ + c₁) ⋅ λ + c₀

= c₀λ⁰ + c₁λ¹ + c₂λ² + c₃λ³ + c₄λ⁴ + λ⁵

= p(λ)

Vandermonde matrix

同伴矩陣的特徵基底，恰是「內插矩陣」。

內插矩陣：同伴矩陣的特徵向量是(λ⁰, λ¹, ... , λⁿ⁻¹)。以特徵值λ求得特徵向量(λ⁰, λ¹, ... , λⁿ⁻¹)，非常方便。

C = VΛV⁻¹     eigendecomposition

    ⎡ λ₁⁰   λ₂⁰   ... λₙ⁰   ⎤
    ⎢ λ₁¹   λ₂¹   ... λₙ¹   ⎥
V = ⎢ λ₁²   λ₂²   ... λₙ²   ⎥     det(V) = prod {λⱼ - λᵢ}
    ⎢ :     :         :     ⎥              i < j
    ⎣ λ₁ⁿ⁻¹ λ₂ⁿ⁻¹ ... λₙⁿ⁻¹ ⎦

eigenbasis - factorization🚧

本章主角是矩陣多項式

矩陣多項式

多項式，變數x原本是實數（複數），現在改成實數矩陣（複數矩陣）。

Cayley–Hamilton theorem

矩陣A的特徵多項式，變數原本是特徵值λ，現在改成矩陣A，結果仍是零。

det(A - λI) = p(λ) = c₀λ⁰ + c₁λ¹ + ... + cₙλⁿ = 0
              p(A) = c₀A⁰ + c₁A¹ + ... + cₙAⁿ = 0

一、考慮特徵分解Aᵏ = EΛᵏE⁻¹。多項式每一項使用了相同的特徵基底E，可以提項出來，乘在括號外面。

一口氣考慮所有特徵值。特徵值們併成對角線矩陣，特徵多項式變成矩陣多項式，最後補上特徵基底，即得此定理。

二、考慮喬登分解Aᵏ = EJᵏE⁻¹。Jᵏ每個元素分別處理即可。

Cayley–Hamilton theorem

矩陣多項式的因式分解。矩陣多項式的根，就是λI！

p(A) = c₀A⁰ + c₁A¹ + ... + cₙAⁿ = 0
p(A) = (A - λ₁I) (A - λ₂I) ... (A - λₙI) = 0

power decomposition【尚無正式名稱】

Aⁿ是A⁰ ... Aⁿ⁻¹的線性組合。

c₀A⁰ + c₁A¹ + ... + cₙAⁿ = 0
(c₀/cₙ)A⁰ + (c₁/cₙ)A¹ + ... = Aⁿ        線性遞迴公式

A⁻¹是A⁰ ... Aⁿ⁻¹的線性組合。

c₀A⁰ + c₁A¹ + ... + cₙAⁿ = 0
c₀A⁻¹ + c₁A⁰ + ... + cₙAⁿ⁻¹ = 0         同乘A⁻¹
A⁻¹ = (-c₁/c₀)A⁰ + ... + (-cₙ/c₀)Aⁿ⁻¹   反矩陣公式

Ax = b的唯一解是A⁰b ... Aⁿ⁻¹b的線性組合。

Ax = b  ⟹  x = A⁻¹b  ⟹  x = (-c₁/c₀)A⁰b + ... + (-cₙ/c₀)Aⁿ⁻¹b
x ∈ span(A⁰b, A¹b, ... , Aⁿ⁻¹b)

A的任意次方是A⁰ ... Aⁿ⁻¹的線性組合。同乘、相減。

invariants of tensor

特徵多項式的各項係數。

c₀   = λ₁λ₂ ... λₙ               = det(A) 
 :
cₙ₋₂ = λ₁λ₂ + λ₁λ₃ + ... + λₙ₋₁λₙ
cₙ₋₁ = λ₁ + λ₂ + ... + λₙ         = tr(A)
cₙ   = 1

minimal polynomial

monic polynomial and m(A) = 0

A的最小多項式：首項係數為一、代入A等於0。

m(λ) | p(λ) and m(λ) = (x-λ₁)^n₁ (x-λ₂)^n₂ ... where n₁ n₂ ... are positive

特徵多項式的因式、至少包含每一種特徵值。

https://math.stackexchange.com/questions/2667415/

size of maximal Jordan block = degree of minimal polynomial
https://math.stackexchange.com/questions/2382023/

annihilator

細節請見專著《Basic Matrix Algebra with Algorithms and Applications》。

krylov matrix: K = [x Ax AAx ...]
T-Annihilator of x: c(A)x = 0
minimal polynomial: c(A) = 0 for all x
1. Kc = 0      (basic solution)
2. K⁻¹AK = C

eigenbasis - group🚧

本章主角是矩陣代數

矩陣代數

我還沒搞懂，大家自己揣摩吧。

Jordan–Chevalley decomposition

A = S + N where S and N commute
S is semisimple (complex diagonalizable)
N is nilpotent (vanish via shift) (fixpoint iteration has stable point)

特徵基底：整數x
主對角線矩陣（縮放暨旋轉）：倍率c
上對角線矩陣（移位）：餘數r
喬登分解：n = cx + r
可交換：互質

( (A - λI) + λI )
  ^^^^^^^^   ^^
  nilpotent  diagonal

如何換基底  換到左邊剛好變成 shift matrix   目前沒有演算法

nilpotent matrix / semisimple matrix

「冪零矩陣」。特徵值皆零。對角線皆零、上三角非零。

「半單矩陣」。特徵值皆異。最小多項式不含重根。

A matrix is diagonalizable over C
iff its minimal polynomial has distinct roots.

雖然特徵空間以外的空間，可以隨便找基底，然後對應的倍率設成零。
但是冪零矩陣，特徵值矩陣的對角線都是零，怎麼乘都乘不出上三角元素。

nilpotent matrix / unipotent matrix

「冪零矩陣」。Aᵏ = 0。特徵值都是0。不可逆。

「冪么矩陣」。(A - I)ᵏ = 0。特徵值都是1。可逆。

兩者皆不可對角化（除非是零矩陣、恆等矩陣）。

nilpotent matrix / permutation matrix

「冪零矩陣」。Aᵏ = 0。

「排列矩陣」。Aᵏ = I。

冪零矩陣不含循環節，於是消滅座標軸。N步之後座標軸必盡數消滅。

排列矩陣都是循環節，於是不消滅座標軸。大家各自自轉，總有一天（最小公倍數）回到原點。

也有同時擁有兩者特性的矩陣。例如設計一種分塊矩陣，一部分冪零、一部分排列。

nilpotent matrix / idempotent matrix

「冪零矩陣」。Aᵏ = 0。特徵值都是0。可對角化。

「冪等矩陣」。A² = A。特徵值都是0或1。可對角化。

鏡射矩陣相乘，等同向量角度相加，等同旋轉。

3 shear matrix = 1 rotation matrix.
http://www.matrix67.com/blog/archives/5453

AB = A
https://math.stackexchange.com/questions/973013/
https://math.stackexchange.com/questions/929741/

A = sum λᵢXᵢYᵢᵀ + Dᵢ = sum λᵢPᵢ + Dᵢ     where Dᵢ = (A - λᵢI)P
P is idempotent
D is nilpotent
jordan form = shift + nonshift
不動點就是不斷shift把nilpotent變不見
https://books.google.com.tw/books?id=xt136Bvd4zgC&pg=PA28

eigenbasis - 總結

特殊的矩陣

（複數版本）么正矩陣unitary matrix：特徵向量形成么正矩陣，特徵值的絕對值都是1。

（實數版本）正規正交矩陣orthonormal matrix：特徵向量形成正規正交矩陣，特徵值的絕對值都是1。注意到特徵值可為複數，例如旋轉矩陣。

    unitary matrix A : Aᴴ = A⁻¹
orthonormal matrix A : Aᵀ = A⁻¹

（複數版本）共軛對稱矩陣Hermitian matrix：特徵向量形成么正矩陣，特徵值都是實數。

（實數版本）對稱矩陣symmetric matrix：特徵向量形成正規正交矩陣，特徵值都是實數。

Hermitian matrix A : Aᴴ = A
symmetric matrix A : Aᵀ = A

（複數版本）反共軛對稱矩陣skew-Hermitian matrix：特徵向量形成么正矩陣，特徵值都是虛數。

（實數版本）反對稱矩陣skew-symmetric matrix：特徵向量形成正規正交矩陣，特徵值都是虛數。

skew-Hermitian matrix A : Aᴴ = -A
skew-symmetric matrix A : Aᵀ = -A

正規矩陣normal matrix：特徵向量形成么正矩陣。

normal matrix A : Aᴴ A = A Aᴴ

半正定矩陣positive semidefinite matrix：特徵值都是正數或零。

正定矩陣positive definite matrix：特徵值都是正數。

positive semidefinite matrix A : xᴴAx ≥ 0
    positive definite matrix A : xᴴAx > 0 when x ≠ 0

symmetric positive semidefinite matrix A : A = Mᴴ M
    symmetric positive definite matrix A : A = Mᴴ M and M⁻¹ exists

可交換矩陣commuting matrices：特徵向量相同。如果有特徵向量的話。

commuting matrices A and B : AB = BA

相似矩陣similar matrices：特徵值相同。這是單向推論。例外是恆等矩陣和歪斜矩陣，特徵值都是兩個1，但是不相似。

similar matrices A and B : A = M⁻¹ B M

對角線矩陣diagonal matrix：特徵向量即是矩陣本身，特徵值即是對角線。這是單向推論。

diagonal matrix A : A_i,j = 0 when i ≠ j

三角矩陣triangular matrix：特徵值即是對角線。這是單向推論。

upper triangular matrix A : A_i,j = 0 when i > j
lower triangular matrix A : A_i,j = 0 when i < j

特殊的分解

特徵分解：特徵向量矩陣（合量）．特徵值矩陣（座標）．特徵向量矩陣的反矩陣（分量）。

A = E Λ E⁻¹

喬登分解：半單矩陣（倍數）＋冪零矩陣（餘數）。

A = M + N     where MN = NM    (commuting)
                    M  = EΛE⁻¹ (semisimple) (diagonalizable)
                    Nᵏ = 0     (nilpotent)

奇異值分解：內積與外積的特徵分解，融合之後開根號。

A = U Σ Vᵀ    where AᵀA = VΛVᵀ
                    AAᵀ = UΛUᵀ
                    Λ = ΣΣᵀ = ΣᵀΣ = Σ²

極分解：正規正交矩陣（旋轉翻轉）．正定矩陣（縮放）。

A = Q P

對稱分解：對稱矩陣（真．縮放）＋反對稱矩陣（真．旋轉）。

A = S + R     where S = (A + Aᵀ) / 2
                    R = (A - Aᵀ) / 2

幾何意義

對角線　　　＝縮放/鏡射
正規正交　　＝旋轉/鏡射（兩鏡射＝一旋轉）
對稱正規正交＝鏡射
次對角線　　＝移位
可逆　　　　＝Ｎ維
特徵值無零　＝可逆
特徵值實數　＝特徵基底縮放/鏡射
特徵值虛數　＝特徵基底旋轉/鏡射
特徵值負號　＝特徵基底鏡射
可對角化　　＝特徵基底Ｎ維
可交換　　　＝特徵基底相同
正規　　　　＝特徵基底垂直
對稱　　　　＝特徵基底垂直、特徵值實數
反對稱　　　＝特徵基底垂直、特徵值虛數
正規正交　　＝特徵基底垂直、特徵值長度一
對稱正規正交＝特徵基底垂直、特徵值正負一
轉置　　　　＝投影（測量向量長度？）
三角/冪零 　＝歪斜（改變基底夾角？）

複數矩陣　　＝特徵基底Ｎ維縮放旋轉歪斜　　　　　　　　喬登分解
複數可對角化＝特徵基底Ｎ維縮放旋轉　　　　　　　　　　特徵分解
實數可對角化＝特徵基底Ｎ維縮放　　　　　　　　　　　　特徵分解

實數矩陣＝正規正交．喬登標準型．正規正交　　　　　　　喬登分解
實數矩陣＝正規正交．對角線．正規正交　　　　　　　　　特徵分解
實數矩陣＝正規正交．實數對稱半正定　　　　　　　　　　極分解
實數矩陣＝正規正交．上三角　　　　　　　　　　　　　　ＱＲ分解（Ｎ次鏡射）
實數矩陣＝下三角．上三角　　　　　　　　　　　　　　　ＬＵ分解（Ｎ次歪斜）

複數矩陣＝特徵正交基底縮放＋特徵正交基底旋轉　　　　　對稱分解
實數矩陣＝特徵正交基底剛體旋轉．特徵正交基底正向縮放　極分解