Eigenbasis

Eigenbasis

「特徵基底」來龍去脈盤根錯節,林林總總各有千秋。以下將從各種面向,描述「特徵基底」及演算法,分成數章介紹。

這些章節的標題都是我自己定的。這種章節架構僅供參考。這種章節架構著重數感直觀,但是不依循公理系統。數學系的線性代數教科書絕不會採用這種章節架構。

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 (Separation)

本章主角是特徵分解(微觀視角)

向量分離【尚無正式名稱】

矩陣拆成向量。

[ |  |  |  ]         [ |  ]   [ |  ]   [ |  ]
[ x₁ x₂ x₃ ]  ―――→ { [ x₁ ]   [ x₂ ]   [ x₃ ] }
[ |  |  |  ]         [ |  ] , [ |  ] , [ |  ]

Eigendecomposition

以座標軸的觀點來重新看待特徵分解。

當輸入向量與特徵向量的方向相同,那麼輸出向量就是輸入向量乘上特徵值,長度伸縮,方向不變。

當輸入向量與特徵向量的方向不同,那麼可以運用分量的概念:Ax = A(x₁+x₂) = Ax₁ + Ax₂。輸入向量,在特徵向量上面的分量,分別乘上特徵值,最後合而為一,得到輸出向量。

拆成三個步驟,嘗試改成矩陣:

一、「輸入向量,在特徵向量上面的分量」:套用反矩陣,座標軸是特徵向量。
二、「各個分量分別按照特徵值來縮放」:套用縮放矩陣,縮放比例是特徵值。
三、「各個分量合而為一,得到輸出向量」:套用矩陣,座標軸是特徵向量。

矩陣連乘要從右往左乘:

                                                     -1
[  4 -1  0 ]   [ |  |  |  ] [ λ₁ 0  0  ] [ |  |  |  ]
[  0  5 -2 ] = [ x₁ x₂ x₃ ] [ 0  λ₂ 0  ] [ x₁ x₂ x₃ ]
[ -3  9 -3 ]   [ |  |  |  ] [ 0  0  λ₃ ] [ |  |  |  ]
      A             E            Λ            E⁻¹

特徵向量們,宛如座標軸,稱作「特徵基底eigenbasis」。

特徵分解A = EΛE⁻¹:分量、座標、合量。

喬登分解A = EJE⁻¹:分量、座標與移位、合量。

三個變換

一、主對角線:
 甲、實數:縮放
 乙、虛數:旋轉
二、次對角線:
 甲、元素0:無作用
 乙、元素1:移位

Eigenbasis (Iteration)

本章主角是線性動態系統

矩陣次方

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

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

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

喬登標準型J的次方,請參考:

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

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的根號,請參考:

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

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 + aI)x = Ax + ax = λx + ax = (λ + a)x。

二、矩陣倒數(反矩陣),則特徵值們也會倒數。

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

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

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

xk+1 = B xk = (A - αI)⁻¹ xk

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

xk+1 = (A - αI)⁻¹ xk
(A - αI) xk+1 = xk
find LU decomposition of (A - αI)

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

特徵向量演算法(Rayleigh Quotient Iteration)

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

αk = (xkᵀ A xk) / (xkᵀ xk)
xk+1 = B xk = (A - αkI)⁻¹ xk

遺珠之憾,不能使用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
projb(a) = ——— b = ——— b
           b∙b     bᵀb

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

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

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

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

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

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

projB( ) = B(BᵀB)⁻¹Bᵀ

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

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

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

a - projB(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⁺

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

projB(A) = B(BᵀB)⁻¹BᵀA
A - projB(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  ]
IX = [------|---]   IY = [---|------]
     [  O   | O ]        [ O | Iⱼₓⱼ ]

projX,Y( ) = VIXV⁻¹ = I - VIYV⁻¹     project to X along Y
projY,X( ) = VIYV⁻¹ = I - VIXV⁻¹     project to Y along X

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

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

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

projX,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―――― |
↓          ↓                    ↓          ↓

矩陣變換,特徵值不變。

特徵向量們隨之變換,特徵值們、喬登標準型都保持不變。

{ A = EJE⁻¹     Jordan-Chevalley decomposition
{ B = FAF⁻¹     similarity transformation

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

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

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

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

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)

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

幾乎上三角矩陣的情況下,幾乎都是零了,遞推次數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分解。

[ X X X X X ]   [                ] [ R₁₁ R₁₂ R₁₃ R₁₄ R₁₅ ]
[ X X X X X ]   [  |  |  |  |  | ] [  0  R₂₂ R₂₃ R₂₄ R₂₅ ]
[ X X X X X ] = [ q₁ q₂ q₃ q₄ q₅ ] [  0   0  R₃₃ R₃₄ R₃₅ ]
[ X X X X X ]   [  |  |  |  |  | ] [  0   0   0  R₄₄ R₄₅ ]
[ X X X X X ]   [                ] [  0   0   0   0  R₅₅ ]
      A                  Q                    R

矩陣AQ實施QR分解。

[                     ]   [                   ] [ H₁₁ H₁₂ H₁₃ H₁₄ H₁₅ ]
[  |   |   |   |   |  ]   [  |  |  |  |  |  | ] [ H₂₁ H₂₂ H₂₃ H₂₄ H₂₅ ]
[ Aq₁ Aq₂ Aq₃ Aq₄ Aq₅ ] = [ q₁ q₂ q₃ q₄ q₅ q₆ ] [  0  H₃₂ H₃₃ H₃₄ H₃₅ ]
[  |   |   |   |   |  ]   [  |  |  |  |  |  | ] [  0   0  H₄₃ H₄₄ H₄₅ ]
[                     ]   [                   ] [  0   0   0  H₅₄ H₅₅ ]
                                                [  0   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)
Aq₅ = H₁₅q₁ + H₂₅q₂ + H₃₅q₃ + H₄₅q₄ + H₅₅q₅ + H₆₅q₆ (5th column of H)
where Hᵢⱼ = (Aqⱼ ∙ qᵢ) = qᵢᵀAqⱼ

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

[                     ]   [                ] [ H₁₁ H₁₂ H₁₃ H₁₄ H₁₅ ]
[  |   |   |   |   |  ]   [  |  |  |  |  | ] [ H₂₁ H₂₂ H₂₃ H₂₄ H₂₅ ]
[ Aq₁ Aq₂ Aq₃ Aq₄ Aq₅ ] = [ q₁ q₂ q₃ q₄ q₅ ] [  0  H₃₂ H₃₃ H₃₄ H₃₅ ]
[  |   |   |   |   |  ]   [  |  |  |  |  | ] [  0   0  H₄₃ H₄₄ H₄₅ ]
[                     ]   [                ] [  0   0   0  H₅₄ H₅₅ ]
          AQ                       Q                    H

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

q₁ = x/‖x‖     normalization
for t = 1~N
   x = Aqt
   x = x - H1tq₁ - H2tq₂ - ... - Httqt     orthonormalization
   qt+1 = x/‖x‖     normalization
   Ht+1,k = ‖x‖     length

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

幾乎上三角分解演算法(Lanczos Iteration)

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

導致H變成「三對角線矩陣Tridiagonal Matrix」。

[                     ]   [                ] [ H₁₁ H₁₂  0   0   0  ]
[  |   |   |   |   |  ]   [  |  |  |  |  | ] [ H₂₁ H₂₂ H₂₃  0   0  ]
[ Aq₁ Aq₂ Aq₃ Aq₄ Aq₅ ] = [ q₁ q₂ q₃ q₄ q₅ ] [  0  H₃₂ H₃₃ H₃₄  0  ]
[  |   |   |   |   |  ]   [  |  |  |  |  | ] [  0   0  H₄₃ H₄₄ H₄₅ ]
[                     ]   [                ] [  0   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)

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

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

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

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

Eigenbasis (Axis)(Under Construction!)

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

Normal Matrix

「正規矩陣」。A與Aᵀ可交換。AᵀA = AAᵀ。

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

命名錯誤。本是正交(垂直),卻稱正規(一單位量)。

Orthonormal Eigendecomposition

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

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

A = EΛE⁻¹     when A is diagonalizable
A = EΛEᵀ      when A is normal : AᵀA = AAᵀ

備忘

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

1. A and Aᵀ are similar.
2. An odd sized skew-symmetric matrix cannot be invertible.
3. skew-symmetric matrix A exhibits xᵀAx = 0.
AᵀA has full row rank iff det(AᵀA) != 0
AAᵀ has column row rank iff det(AAᵀ) != 0
https://math.stackexchange.com/questions/903028/

row rank = column rank
https://en.wikipedia.org/wiki/Rank_(linear_algebra)

【尚待確認】

可交換矩陣保留特徵向量。因此A與Aᵀ特徵向量相同。

轉置矩陣保留特徵值。因此A與Aᵀ特徵值相同。

因此,特徵向量暨特徵值相同,是相同矩陣?我懵逼了。

所以說,應該需要考慮特徵向量與特徵值的對應關係吧。

Ax = λx
❮x,Aᴴx❯ = ❮Ax,x❯ = ❮λx,x❯ = ❮x,λx❯
Aᴴx = λx ???
https://math.stackexchange.com/questions/461958/
Ax₁ = λ₁x₁
Ax₂ = λ₂x₂
❮Ax₁,x₂❯ = ❮λ₁x₁,x₂❯ = λ₁❮x₁,x₂❯
❮Ax₁,x₂❯ = ❮x₁,Aᴴx₂❯ = ❮x₁,λ₂x₂❯ = λ₂❮x₁,x₂❯
λ₁ ≠ λ₂ thus ❮x₁,x₂❯ = 0
https://math.stackexchange.com/questions/778946/
{ Axᵢ = λᵢxᵢ   => xᵢᴴyⱼ = 0
{ Aᴴyᵢ = λᵢyᵢ
http://macs.citadel.edu/chenm/344.dir/08.dir/lect4_2.pdf

矩陣變換:A變Aᵀ

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

可以特徵分解:一種答案是特徵基底的外積矩陣的反矩陣(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/

Transpose與Inverse的運算

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

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

Eigenbasis (Symmetry)

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

矩陣轉置

實數的反面:變號。

虛數的反面:共軛。

矩陣的反面:轉置。

Symmetric Matrix / Hermitian Matrix

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

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

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

{ A = Aᵀ      real symmetric 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 

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。

Eigenbasis (Absolute Value)(Under Construction!)

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

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                          半正定矩陣

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ᵀ。

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)

http://www.netlib.org/utk/people/JackDongarra/etemplates/node198.html

雙對角線分解演算法(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)

http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf

兩階段。

[ 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)

https://books.google.com.tw/books?id=4Mou5YpRD_kC&pg=PA238

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 (Reciprocal)(Under Construction!)

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

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。

Adjoint Matrix

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

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 (Length)(Under Construction!)

本章主角是內積矩陣:奇異值平方

矩陣長度平方【尚無正式名稱】

複數的長度平方:複數乘以共軛複數。

線性函數的長度平方:矩陣乘以轉置矩陣。

注意到,矩陣長度平方‖A‖² = AᵀA、矩陣平方A² = AA,兩者不同。

矩陣長度平方【尚無正式名稱】

「向量自身點積」定義成「向量長度平方」。相信大家都學過。

    def
xᵀx === ‖x‖²

「矩陣自身點積」定義成「向量長度平方」。如法炮製。

    def
MᵀM === ‖M‖²   

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

複數的內積分解:分解成複數乘以共軛複數。

線性函數的內積分解:分解成矩陣乘以轉置矩陣。

注意到,矩陣內積分解⟌A = √ΛEᵀ、矩陣平方根√A = E√ΛEᵀ,兩者不同。

只有對稱矩陣允許內積分解。答案無限多個。常見解法:特徵分解、Cholesky分解。

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

內積分解目前沒有數學符號。內積分解類似於開根號。

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

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

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

我們希望對稱矩陣A是「矩陣長度平方」,然而特徵值Λ可能為負數,√Λ可能為複數,M可能為複數。長度平方很怪,開根號更怪。我們必須追加限制條件,讓事情變得合理。

Symmetric Matrix

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

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

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ᵀ

內積矩陣相關定理

rank(AᵀA) = rank(A) = rank(Aᵀ) = rank(AAᵀ)
https://math.stackexchange.com/questions/349738/

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,得到左右奇異向量。

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

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

http://users.cecs.anu.edu.au/~kee/pls.pdf

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 (Unit)(Under Construction!)

本章主角是單元矩陣:奇異值為一

Unit Matrix【尚無正式名稱】

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

unit vector譯做單位向量、identity matrix譯做單位矩陣,顯然有一個是錯的。identity matrix按照字面翻譯,應該譯做恆等矩陣。啊就古人數學水平較低,無法辨別unit和identity的差別,翻譯成相同詞彙。

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

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

‖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⧛ = Ib×b when a ≥ b。

扁平矩陣:當直條多於等於橫條,自身外積是單位矩陣,使得橫條正規正交。

⧚A⧛⧚A⧛ᵀ = Ia×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ᵀ
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 (Metric)

本章主角是距離

Symmetric Matrix

「對稱矩陣」。yᵀAx = xᵀAy。

這是仿照dot(y, Ax) = dot(x, Ay)的概念,幾何意義是兩個向量x與y,無論取哪一個實施線性變換,點積仍相同。

Positive Definite Matrix / Positive Semidefinite Matrix

「正定矩陣」。xᵀAx > 0,不討論x = 0。

「半正定矩陣」。xᵀAx ≥ 0,不討論x = 0。

這是仿照dot(x, Ax) > 0的概念,幾何意義是Ax投影到x上面,投影量總是大於0,無論哪種x。換句話說,向量經過線性變換,轉彎幅度少於±90˚,未達半個面。

眼尖的讀者應該注意到了,dot(x, Ax)再除以x的長度平方,就是虛擬特徵值。正定矩陣的虛擬特徵值恆正,於是特徵值恆正。

Symmetric Positive Definite Matrix

「對稱正定矩陣」。既對稱,又正定。同時具備兩者性質。

原本定義可以改寫成xᵀAx = xᵀMᵀMx = (Mx)ᵀMx = ‖Mx‖² > 0。代數意義是線性變換之後的長度的平方、恆正。

A-orthogonal

「矩陣正交」。yᵀAx = 0。xᵀAy = 0。兩式等價。

Ax與y互相垂直;x經過A變換之後,與y互相垂直。

知名範例:切線速度與加速度互相垂直;一次微分與二次微分互相垂直。

定義y與Ax的內積空間,A必須是對稱正定矩陣。內積空間開根號可以定義距離。

Length / Distance

l(x) = ‖x‖ = √xᵀx   |   d(x,y) = ‖x - y‖ = √(x-y)ᵀ(x-y)
lᴍ(x) = ‖Mx‖        |   dᴍ(x,y) = ‖Mx - My‖ = ‖M(x - y)‖
      = √(Mx)ᵀ(Mx)  |           = √(M(x-y))ᵀ(M(x-y))
      = √xᵀMᵀMx     |           = √(x-y)ᵀMᵀM(x-y)
      = √xᵀAx       |           = √(x-y)ᵀA(x-y)

長度:元素平方和開根號。自己與自己的點積、再開根號。

線性變換之後的長度:恰是二次型開根號,是對稱半正定矩陣。

距離:相減之後的長度。

線性變換之後的距離:宛如二次型開根號,是對稱半正定矩陣。

Eigenbasis (Polynomial)(Under Construction!)

本章主角是特徵多項式

矩陣多項式

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

多項式,變數x原本是實數(複數),現在改成實數矩陣(複數矩陣)。

Cayley-Hamilton Theorem

https://en.wikipedia.org/wiki/Cayley–Hamilton_theorem

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

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

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

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

Krylov Matrix

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

    [  |   |        |   ]
K = [ A⁰v A¹v ... Aⁿ⁻¹v ]
    [  |   |        |   ]

靈感源自Power Iteration。用來求特徵值。

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

Krylov Subspace

Krylov Matrix的N個向量所構成的空間。

如果初始向量是K個特徵向量的線性組合,那麼Krylov Matrix將生成這K個特徵向量的線性組合。特徵值皆異,那麼Krylov Subspace是K維。特徵值相同,那麼維度短少。

當矩陣A擁有N個特徵向量,特徵值皆異,那麼Krylov Subspace才是N維。如果不同特徵向量擁有相同特徵值,那麼維度短少。

當矩陣A不足N個特徵向量,那麼Krylov Subspace的維度,需要考慮Jordan Block。初始向量所屬的K個特徵向量,落在哪些Jordan Block裡面,該Jordan Block擁有的特徵向量都會藉由移位出現。【尚待確認】

隨機生成初始向量,通常可以涵蓋所有特徵向量,運氣不會差到只涵蓋K個特徵向量。

矩陣A擁有反矩陣,那麼Krylov Subspace是N維。

Cayley-Hamilton theorem:
(A - λ₁I) (A - λ₂I) ... (A - λₙI) = 0 
c₀A⁰ + c₁A¹ + ... + cₙAⁿ = 0            展開
c₀A⁻¹ + c₁A⁰ + ... + cₙAⁿ⁻¹ = 0         同乘A⁻¹
A⁻¹ = (-c₁/c₀)A⁰ + ... + (-cₙ/c₀)Aⁿ⁻¹   移項
也就是說,A⁻¹是A⁰...Aⁿ⁻¹的線性組合,只有唯一解,A⁻¹的維度是n。

Krylov subspace:
Ax = b  =>  x = A⁻¹b  =>  x = (-c₁/c₀)A⁰b + ... + (-cₙ/c₀)Aⁿ⁻¹b
Krylov(A, b) = span(A⁰b, A¹b, ... , Aⁿ⁻¹b)

Krylov Basis

QR分解,扳正矩陣K,仍然得到相同子空間。

QR分解,遞推法,生成向量、扳正向量,仍然得到相同子空間。

1. shift: A + λI
2. scale: kA
3. change of basis:

Power Iteration、Arnoldi Iteration、Conjugate Gradient Descent都可以使用Krylov Subspace進行演算法分析。

有些數學教科書將Krylov Subspace當作主角,反客為主,本末倒置。

Companion Matrix

    [  |   |        |   ]
K = [ A⁰b A¹b ... Aⁿ⁻¹b ]
    [  |   |        |   ]
     [  |   |       |  ]     [  |  |        |   ]
AK = [ A¹b A²b ... Aⁿb ] = K [ e₂ e₃ ... K⁻¹Aⁿb ] = KC
     [  |   |       |  ]     [  |  |        |   ]
K⁻¹AK = C   similarity transformation
          [ 0 0 ... 0 c₀   ]
          [ 1 0 ... 0 c₁   ]
where C = [ 0 1 ... 0 c₂   ] , c = K⁻¹Aⁿb
          [ : :     : :    ]
          [ 0 0 ... 1 cₙ₋₁ ]

任意矩陣實施相似變換,Krylov Matrix做為變換矩陣,得到Companion Matrix。

Cayley-Hamilton Theorem:矩陣多項式的根,就是特徵值!

Companion Matrix:多項式的係數與變數,化作矩陣與輸入向量。多項式的根,化作矩陣的特徵值。

特徵多項式就是原本多項式。

特徵向量是(λ⁰, λ¹, ... , λⁿ⁻¹),可以用特徵值λ求得,非常方便。

Minimal Polynomial

內插矩陣,特徵值的各個次方。特徵值相異時,內插矩陣可逆,可以用來說明Krylov Subspace維度大小。

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

延伸閱讀: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

Eigenbasis (Field)(Under Construction!)

本章主角是矩陣代數

矩陣代數

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

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步之後座標軸必盡數消滅。

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

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

Idempotent Matrix

「冪等矩陣」。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 : Ai,j = 0 when i ≠ j

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

upper triangular matrix A : Ai,j = 0 when i > j
lower triangular matrix A : Ai,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

幾何意義

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

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

實數矩陣=正規正交.喬登標準型.正規正交       喬登分解
實數矩陣=正規正交.對角線.正規正交         特徵分解
實數矩陣=正規正交.實數對稱半正定          極分解
實數矩陣=正規正交.上三角              QR分解(N次鏡射)
實數矩陣=下三角.上三角               LU分解(N次歪斜)

複數矩陣=特徵正交基底縮放+特徵正交基底旋轉     對稱分解
實數矩陣=特徵正交基底剛體旋轉.特徵正交基底正向縮放 極分解

Eigenbasis應用(Under Construction!)

Linear Approximation【尚無正式名稱】

特徵值絕對值最小者歸零,得到最小平方誤差矩陣。

消滅任一項:該特徵值歸零,該特徵向量不變。維度減一。

消滅特徵值最小項:矩陣元素誤差平方和最小的降維方式。

‖A‖F² = trace(AᵀA) = trace((sum λᵢPᵢ)ᵀ(sum λᵢPᵢ))
=========== trace(sum λᵢ²Pᵢ) = sum λᵢ² trace(Pᵢ) = sum λᵢ²
when Pᵀ = P aka A is normal
pseudoinverse = projection = least squares = linear regression
= dot product of arbitrary basis

linear approximation
= dot product of othronormal basis
galerkin method: orthonormal basis dot
conjugated gradient method: A-orthonormal residual dot
linear feedback system: time series dot (fixedpoint)
companion matrix: x x^2 x^3 dot (differential equation eigenfunction)
taylor expansion: order-1 to order-n gradient dot
kalman filter: residual dot
https://www.encyclopediaofmath.org/index.php/Galerkin_method
正規正交基底的線性組合的線性變換,分別對各個基底點積。
Ax = A(c1q1 + c2q2 + ...)
q*Ax = c q*Aq   when dot(qi,qj) = 0

點積結果是二次型。正規化之後得到虛擬特徵值,乘上c倍。
c q*Aq/q*q

https://en.wikipedia.org/wiki/QR_decomposition
QR分解之GS正交化 Q當作正規正交基底 R是點積
dot(q,a) = q*a
[galerkin method solve linear equation]
xk = x0 + subspace = x0 + (α0 b + α1 Ab + α2 A²b)
such that dot(b - Axk, q) = 0   餘數通通通過基底q的測試
如果是 linear least squares  應該要改成 min dot(b - Axk, q)
xₜ₊₁ = xₜ - dot(J, F) / dot(J, J)
視作投影公式:F(x)投影到J。
向量投影到梯度,再求得向量到梯度的差距。
can we use them? when use the linearization
1. fixpoint iteration (linear equation)
2. power iteration (linear eigenvalue)
n->n求根   其實就是求特徵值? power iteration?
           companion matrix? Cayley-Hamilton Theorem?