2017年9月24日 星期日

閱讀筆記 - matrix computation (3e)

讀完之後打下來以免哪天又要重新看一次
不是什麼教學性質的,只是單純的筆記而已


Chapter 2

  • 說明關於 float point rounding error 的問題
  • 點出 SVD 的存在
  • condition number(SVD 最大最小的 singular value 比值)
    一個矩陣數值穩定度好不好可以由這個數字判斷,越大越爛
其他太數學了,不理他

Chapter 3

基本教學,不過其實不是很實用
點出了 A=LU 分解的存在(L 下三角,U 上三角)
其實就是我們熟悉的高斯消去法的代數表現
L 就是我們要做 reduced row echelon form 的係數
U 就是 reduced row echelon form

這個方法其實很爛,簡單的例如說

[0 1]
[1 0]

就不會動了
所以就有 PA=LU 這種方法,P 是 permuate matrix
就是說 PA 是 A 的 row permutation
概念是在要進行高斯消去的時候(AKA 找 L 的係數的時候)
先找 column 剩下的數字中絕對值最大的 swap 上來
這個方法相比之下好了不少

另外因為 LU 分解是唯一的,所以如果說 A 對稱的話
A = LU, LT = U
而且如果是對稱正定,這個還算穩定
不是的話還是要 permute
此外,對稱矩陣的 LDL 分解就是把 L 的對角項目 normalize
把該係數放到中間而已
同理 LDU 分解也是存在的

Chapter 4~5

非常實用的兩個章節,屌打 chapter 3 的方法
如果不考慮其實這些方法都可以用 numpy 或是 matlab 一行搞定的話 XD

Householder

假設我們有個向量是 v = [4 1 8]
用一個鏡射,想把他反射到 [9 0 0] (PS 保持長度,剩下第一個分量)
矩陣要怎麼寫
我們知道對於一個「單位法向量是 a 經過原點的平面」的鏡射
可以寫作 I-2aaT
所以 v-2a(aTv) = [|v| 0 0 ...0]
得到 a = normalize(v-[|v| 0 0 ...0])
(也可以是 a = normalize(v+[|v| 0 0 ...0]),實際上建議取v第一分量異號)

可以試著用 a = [-5 1 8]/sqrt(90) 建構出反射矩陣
乘到這個上(注意第一個 row)

[4 2]     [* *]
[1 0] --> [0 *]
[8 0]     [0 *]

Given

就是兩個 row 的旋轉,像這個形式

[ cos sin]
[-sin cos]

[ 0.8 0.6]  [4 1] = [* *]
[-0.6 0.8]  [3 1]   [0 *]

Householder / Given 的特性
都是 orthonormal 轉換
都有把矩陣變成 0 的功能
Householder O(n2) 把一個 column 變成 0,比較適合平行處理
Given 可以在 O(n) 把一個值變成 0
而且因為可以部份保持原本 0 的還是 0
(像是 Householder 的 [2 0 0]轉完之後全部都不是0了)

因此
Householder 比較適合 general matrix
Given 適合處理特殊結構

QR 分解

A = QR (Q is orthonormal, R is upper triangular)
根據前面寫的方法我們可以建造一系列的 orthonormal matrix 乘到 A 上面
最終得到 R = MMMMMA 的形式
我們就能很容易求出 Q
順序如下:

Householder | Given
[* * *]     | [* * *]
[0 * *]     | [2 * *]
[0 1 *]     | [1 4 *]
[0 1 2]     | [0 3 5]

另外,Householder 可以累積多個 reflection 一起轉
使得 GEMM blas3 級的運算佔大多數

不過實際上在解 least square 問題 Ax=b 的時候我們也不需要 Q
只要把 b 做一樣的轉換即可
R = MMMMA = MMMMb
剩下的問題是上三角矩陣的 inverse
可以在 O(n2) 完成
這樣可以大幅增加數值穩定度
基本上只要矩陣好,該方法的數值非常穩定

另外還有 gram schmidt 跟 modified gram schmidt
基本上前者就是大一課本裡面的 gram schmidt 的代數表示,基本上沒用
後者還行,不過運算複雜度較大
而且得到的 Q 的穩定度差了 condition number 的 order,也不實用

如果矩陣性質太差的時候就要用 SVD 解
A=USV -> A-1 = VTS-1UT
如果遇到 S 裡面有 0 或是太小的值直接捨棄
這個會使 I - AAT 的平方和最小

另外一般常見的 Ax=b 解法
ATA = ATb
這個方法其實超級爛
因為前面那個 (ATA) condition number 是 A 的平方

一個很合理的說法是,明確的建構出不需要的矩陣都會「喪失資訊」
導致穩定度下降

Chapter 8

像是 Householder 這種東西其實可以清理出更多 0
只要乘到左邊也可以清理出一個 row 的 0
順序如下:

[* * 1 1]
[0 * * 3]
[0 2 * *]
[0 2 4 *]
[0 2 4 6]

要注意的是如果我們想清出更多 0
像是 step 1 如果我們也想把紅色星號變成 0
我們就會不小心破壞 step 0 清理出來的 0
所以這個演算法就變成 A = UBV
其中 B 只有對角線跟上面有值
可作為 SVD 的前置步驟

BTW 對於長的矩陣的分解可以先做 QR 分解
A    = Q       R
100x10 100x10  10x10
再做 R = UBV → A = (QU)BV
基本上只要高大於 5/3 寬的話就有節省
基於這個方法的 SVD 叫 R-SVD

要把 bidiagnoal 矩陣變成對角線會用收斂法
這時候就是 given rotation 發揮的時間了
只要 col12 -> row12 -> col23 -> row23... 交互做 given rotation
就會得到每個時間點內,在 ABCDE 這些位置有一個不是 0

 [* * B 0]
 [A * * D]
 [0 C * *]
 [0 0 E *]
 [0 0 0 0]

聽起來做完後還是 bidiagonal,一點用也沒有
不過這個方法下「非對角線的平方總和保證越來越小」
所以這個方法遲早會收斂
雖然是 linear time 收斂,很慢

很多情形下其實也不需要建構出完整的 USV
例如我們只想找最大的 singular value 對應的 vector

沒有留言:

張貼留言