巨量計算中的矩陣分解方法 (附Python程式碼)
Matrix Decomposition Method for Enormous Calculation with Python Code
【前言：為甚麼需要矩陣分解】
資料科學中的數據表達方式，大多是以矩陣的形式呈現的，舉個最簡單的例子，我們考慮一個簡單線性迴歸方程式有如下表達方式：
我們建模的目的，在於找出一個最佳的迴歸係數Beta，使得真實的誤差最小，舉一個N=1，也就是輸入的資料只有一個欄位(一個變數)的時候來看，我們的目標是希望找到紅線，能夠達成Perfect_Match。
當然，我們有很多辦法去找到這個Beta，對於機器來說，一開始可以是甚麼都不知道的胡亂嘗試，然後透過如梯度下降等方式，去逐步逼近我們的最佳解，也就是上圖的紅線。然而，回憶一下高中的數學課程，告訴我們要找到最佳Beta，可以由以下的方式來計算：
上述的方法也叫做最小平方法，也就是只要知道了原始資料，我們即有辦法得到最佳迴歸係數，相比於Trial and Error的梯度下降法，好處是我們不需要經過參數調整的過程，也就是不需要花費時間來修正並嘗試改進模型表現，但是壞處是，當我們的資料量非常龐大，也就是m過大的時候，計算時間會變得非常久，因此，我們會需要一些方式來減少這個計算時間，比如我如果能夠將原先計算量為O(n³)的演算法改為O(n²)時，假設O(1)運算佔據0.00001秒的時間，當 n 非常大的時候，我就能夠省下非常多時間，這代表我不會每測試一次程式，就得花費大量的時間等待。
此外，一般的矩陣表示不見得是最有效率的表示方法，因為資料可能會出現過大、雜訊過多等情形，我們會需要更有效率的方式來表示原先的資訊，而這兩個理由就會跟本文的主角－ 矩陣分解(Matrix Decomposition)有關。
舉運算時間為例，為了計算這個係數，需要的時間複雜度為：
在這裡，我必須先帶出另外一個觀念，也就是運算的時間複雜度(Computational Time Complexity)，上圖中用大O符號表示，如果對時間複雜度不熟悉的朋友，我推薦看一下這篇文章：
也就是說我們會依據m, n 的大小不同，產生O(n³)或O(mn²)上界的時間複雜度，而對資料科學家而言同樣一個計算目標，如果能用更短的時間(這裡先不討論空間問題)，就意味著這支演算法更有效率一些，應用開發或實戰時也就能更快一些，那我們要如何知道這支程式的運算效率，以及怎麼樣的運算次序，才能夠達到理論上的效率上限(也就是最有效率)呢?
【演算法：矩陣鏈相乘問題 (Matrix Chain Multiplication)】
對於上述那個問題的答案是，我們可以知道正確的運算次序，但要找到這個次序的過程，往往要花費不匪的代價，這就是一個演算法中很知名的問題，矩陣鏈相乘問題，理論上要找到最佳次序，我們需要花費O(n³)的時間。
先從小的問題開始看，如果我們有三個矩陣相乘，需要花費的時間為：
此時，我們有兩種組合要嘗試，但如果今天有四個矩陣，問題馬上變得複雜許多，需要測試的組合會變成ABCD、(ABC)D、A(BCD)、(AB)(CD)、(AB)CD、AB(CD)、A(BC)D，可以想見當矩陣有N個的時候，找到最優次序本身會太過複雜，如果有n個矩陣的時候，我們要比較的次數最多會達到n³次，也就是O(n³)的時間複雜度。(#這邊就不附程式碼了，因為要吐出結果需要比較麻煩的解釋與編程，牽扯動態規劃，而這不是本文的重點)
因此，我們通常不會為了減少複雜度，先去計算一遍矩陣的相乘最優次序，因為光是做這件事就得花費O(n³)，即使用Strassen算法也必須花費O(n^2.807)的時間複雜度，顯然這個方法是行不通的。於是，就換到我們今天的主角登場，也就是矩陣分解(Matrix Decomposition)，它的概念在於用相似的矩陣相乘去等於原先我們要乘的矩陣，從而達到減少複雜度的方式。
這篇文章裡面，我們介紹三種常見的矩陣分解的方法：
1. QR 分解
2. Cholesky 分解
3. 奇異值分解
當然，還有其他沒提到的分解方式，本文僅介紹相對常見的幾種分解方法的數學原理，並嘗試實作。更完整的推導、性質等，推薦大家看看周老師的線代啟示錄，對於矩陣分解有非常完整的講解。
【QR分解(QR-Decomposition)】
當一個矩陣可以寫成如下形式，我們稱這個矩陣可以被QR分解：
我們針對這三個定義來分別看一下，第一個不用多說，我們需要確保A的行向量彼此不相依。第二個，我們必須找到一個Q矩陣，且其為正交單範矩陣(Orthonormal Matrix)，回憶一下線性代數，我們知道一個矩陣可以透過線性轉換(Linear Transformation)用另外一個矩陣表達，比方旋轉、拉伸、投影等，比如下圖的示意：
數學上來看，這些被轉換的向量或矩陣都經過了另一個矩陣來做線性變換，因此我們需要知道該變換矩陣，來還原開始的資料，比如說Y=AX，我們必須得知道A才能夠做變換。然而，有一種方式可以在只知道X的情況下去做映射，也就是正交投影(Orthogonal Projection)，亦即我們針對一個向量或矩陣，只要讓它與其剩餘量的內積為零就可以了，如下圖示例：
注意雖然上圖都用向量來做示意圖，實際上推廣到矩陣概念亦同，向量可視為行向量或列向量為一維的矩陣。我們先從一維向量的計算開始，根據上圖，我們想要得到 v 向量在 s 上的正交投影，其實等同於我們在做找最小距離的問題，也就是讓 v-cs 這個剩餘量長度最小(嚴謹上來說，我們稱其為 2-範數平方，也就是一個最小平方化問題) 首先，為了不失其一般性，我們將 s 乘上一個常數 c ，因此得到：
推廣到矩陣，假設我們有一個矩陣 A = [a1, a2, …, an] 要找到其正交基底 U = [u1, u2, …, un]，可以透過Gram-Schmidt方法去得到，而再對該矩陣的向量做單位化之後，就可以得到我們的Q矩陣，遵照以下步驟：
此時，我們的R矩陣顯而易見：
至此，我們簡要說明了QR分解的幾個步驟，現在回頭來看原問題，我們就可以透過以下的方式來拆解線性迴歸估計的Normal Equation：
整體的概念上，等於可以將原先的矩陣 XTX，透過矩陣 RTR 來近似。比較一下原先的計算複雜度約是O(2mn²)，透過QR拆解可以將複雜度下降到O(2mn²-2n³)，當m與n的值相對接近的時候，這會是更有效率的方式，但相對地，如果m的值相對n要大得多，下降的幅度就很有限，基本上隨著資料量加大，這個算法就會越來越慢，因此在小樣本時，我們推薦Cholesky算法。
【Cholesky分解(Cholesky-Decomposition)】
第二種常見的分解方式叫做Cholesky分解，在開始介紹以前，首先，我們先理解一個概念較做LU分解。這個分解對原先的資料 X 做出了一些限制，即X必須要是對稱(Symmetric)且正定(XTX >= 1)：
此時X可被LU分解，也就是 X = LU，L代表下三角矩陣，U代表上三角矩陣。
承接這個概念，Cholesky分解在於將 X 分解成一個下三角矩陣及其轉至矩陣的乘積，此時 X 需為正定矩陣：
在這裡就不證明正定矩陣必可Cholesky分解，這種做法的好處是從頭到尾只產生一個 L 矩陣，計算的穩定性相對較高，但壞處是對原先資料做出了一些限制，使得因為輸入資料導致的計算崩潰容易發生，舉例來說，如果元素的第一行第一列的值是負數，計算有可能就會出錯，因此我們使用時最好有一個檢驗其是否為Positive Definite的前置作業。
回到Least Square Estimation，因為要求輸入資料必須正定，我們於是不對 X 矩陣做分解，轉而對 XTX 做分解，原因是一個實對稱矩陣必然可以寫成一個矩陣與其轉置矩陣的乘積，證明如下：
對應地，我們可以將估計式寫成：
當然，直觀來看Cholesky是更有效率一些的，但Cholesky本身比較有可能出現不穩定，比方說矩陣是Rank Deficient的時候，XTX不會存在反矩陣等。
【奇異值分解(Singular Value Decomposition)】
奇異值分解通常簡稱SVD，是非常強大的矩陣分解工具，除了估計迴歸係數以外，一些機器學習的模型也都是基於該分解方式，比如說用於推薦系統的Matrix Factorization等等。
通常覺得SVD好用的原因是在於，本身除了能夠大幅減少資訊儲存量，計算也相對穩定，不用如Cholesky那樣對資料本身有一定要求，對於Rank Deficient (m << n)以及Overdetermined(m >> n)都能有一定表現。
詳細的數學比較冗長，還是推薦看看周老師的線代啟示錄，這裡只針對SVD的性質，以及在線性回歸估計上的應用做探討。
SVD矩陣會將矩陣分解成以下形式：
奇異值是甚麼? 本質上，整個奇異值分解就是對於原資料矩陣的拉伸變換而已，如周老師在文中提到的，將原始資料透過旋轉、伸縮、旋轉的過程，就完成了SVD，而其中的奇異值矩陣，實際上是一種拉伸的方式。
數學意義上，可以想成對原始資料的一種去噪方式，由於：
其中的奇異值是由大到小排列的，等於說，在 n 個奇異值中我們如果只取 k 個(k ≤ n)，我們亦可以近似原始資料 X ，最有名的應用方式如圖片壓縮，比方說我們有一個512*512向素的圖片，可以視為512維的方陣：
如圖，我們可以明顯看到選擇不同數量奇異值時的圖片變化。
回到正題，我們的線性回歸估計可以寫成：
整個推導相對要複雜一些，我將參考連結放在下方，需要注意的是，此時 D 的反矩陣是一個 n * n 的矩陣，也就是我們僅將其中的非零行列取出，稱為Truncated SVD Decomposition，我們可以由結論看到，相較QR分解與Cholesky分解，時間複雜度要來得更高一些，造成整個計算會相對偏慢，因此通常基於SVD的快速計算方法，都會再加上其他的程式加速方式，比如平行計算等來增加效率。雖然如此，幾乎所有資料都適用的情況下，加上有降噪效果，讓這個算法成為機器學習領域中相對熱門的分解方法之一。
其他還有諸多種分解方法，這邊就不一一提及。而以下是這三種方法估計Least Square 問題的 Python 程式碼，Numpy中求解線性代數的Linalg Package非常完整地支援了這幾種方法，當然，你也可以選擇自己手刻一個，整體來說代碼並不複雜：
- 補充：奇異值(Singular Value)與特徵值(Eigenvalue)並沒有特別關聯，但的確有相似之處。本質上，兩者都是給原資料的基變換(Basis Transformation)乘上一個常數，其中奇異值是兩個Orthogonal Basis的變換，而特徵值是基於Eigenvector的基變換。此外，奇異值都是nonnegative，而特徵值沒有該性質。真的要說關聯的話，奇異值會是矩陣 sqrt(T*T) 的特徵值，其中 T 為 Self-adjoint Linear Operator，這方面在線性代數課本裡有詳細的說明，就不再贅述。
順便附上產生資料的程式碼，大家可以搭配上面的分解方式去產生不同的估計值，一樣是用Python 3 做範例：
【結尾：理論 vs. 現實】
在資料量非常龐大的時候，演算法的改進是大幅節省效率的方式，比如我們介紹至今的方法，然而實際上，因為程式本身的運算邏輯如算法運算次序、基本運算的差異，以及為了考慮泛用性而做出的調整等，都會傷害程式的效率，舉例來說，我們對Python的幾種線性回歸模組測試發現：
在這個例子中，簡單的矩陣逆求解反而是最快的方式，大部分這些模型都兼具泛用性以及延展性，比如 Moore-Penrose Inverse 就包含了SVD分解加速，但效能卻不如一般的矩陣逆求解。因此，數學的調整雖然能夠在理論上增進算法效率，考量到實用範圍時，還是要把程式本身的限制等方面考慮進來。
當然，線性回歸只是矩陣分解一個常見也好理解的例子，實際上，矩陣分解在其他更多的演算法中扮演不可或缺的地位，有興趣的朋友也可以去研究!
