稀疏正定矩陣的Cholesky分解
本文大部分參考這篇文章。圖片也是從他那里復制的>_<
圖和矩陣的對應
考慮矩陣A，如果A[i][j]=w，那么在i,j之間就有一條長度為w的路徑。由於我們考慮的是無向圖，因此這個矩陣A一定滿足\(A=A^T\)
正定（SPD）矩陣的Cholesky分解
要做的事情是將一個正定矩陣A分解為一個下三角矩陣L和其轉置的乘積，也即\(A=LL^T\)。
考慮這樣一個做法：
for i=1 to n-1
A[i][i]=sqrt(A[i][i])
A[i+1:n][i]=A[i+1:n][i]/A[i][i]
A[i+1:n][i+1:n]=A[i+1:n][i+1:n]-A[i+1:n][i]*(A[i+1:n][i])^T
考慮最后那一步，相當於一個矩陣減去一個稀疏向量和其自己轉置的外積。
舉例說明：
可以發現得到的矩陣一定包含一個稠密子矩陣。
不難知道，如果我們單獨取出這個稠密子矩陣，其對應的圖一定是一個團。
於是這份代碼的最后一個步驟可以理解為：我們將第i個節點刪去，同時向第i個節點所有連接的節點加邊，使他們成為一個團。
這樣子我們就可以得到\(L+L^T\)的所有非零元素的位置。
我們通過這樣一個方法，可以將正定矩陣的Cholesky分解與圖聯系起來。
由於最終我們要解決的是稀疏矩陣上的問題，因此我們希望對消去節點的順序排序，使得最終新增的節點數量最少。但是很可惜，找到最優解已經被證明是一個NP-hard問題。所以我們希望能夠得到一個相對來說比較優秀的順序。
先來看一個3X3分塊矩陣的Cholesky分解。
![[公式]](/image/aHR0cHM6Ly93d3cuemhpaHUuY29tL2VxdWF0aW9uP3RleD0lNUNsZWZ0JTVCKyU1Q2JlZ2luJTdCYXJyYXklN0QrRERfMSslMjYrJTI2U18xJTVFVCU1QyU1QyslMjYrRF8yKyUyNlNfMiU1RVQrJTVDJTVDK1NfJTdCMSU3RCUyNitTXyU3QjIlN0QlMjYrUyU1Q2VuZCU3QmFycmF5JTdEKyU1Q3JpZ2h0JTVEKyUzRCU1Q2xlZnQlNUIrJTVDYmVnaW4lN0JhcnJheSU3RCtMTF8xKyUyNislMjYlNUMlNUMrJTI2K0xfMislMjYrJTVDJTVDK01fJTdCMSU3RCUyNitNXyU3QjIlN0QlMjYrTF9TJTVDZW5kJTdCYXJyYXklN0QrJTVDcmlnaHQlNUQrJTVDdGltZXMrTCU1RVQ=.png)
其中D,S,M都可能是稠密的，而L是下三角矩陣。
如果我們最后才把S對應的行和列處理掉（也就是上面的那個沿着對角線依次處理的算法，由於矩陣可以通過行列變換自由交換，故我們這么做是可行的），那么優先處理的部分就是D1，D2，這兩個矩陣構成了一個對角矩陣。按照圖的說法來討論的話，S對應着圖中的一組節點，我們稱之為separator，即一個分段點，先不考慮S，而如果不通過S的話，D1與D2之間沒有連邊，因此在處理S之前兩部分是獨立的，此時我們只需要遞歸的處理D1和D2兩部分即可。一個好的separator應當使得最終D1與D2的規模盡可能的一致，並且S中包含的點數盡可能的少。這個方法稱之為Nested Dissection。
舉個栗子，還是前面那張圖，我們對其重新進行一個順序的划分。
於是先處理D1，再處理D2，最后再處理S，這樣子最終得到的圖就是，紅色是新加入的邊：
對應到矩陣上，首先通過重新排列行列得到：
而計算完成之后，得到的\(L+L^T\)為：
可見，通過這樣的優化，我們在矩陣中新加入點的個數減少了。
上面的所有過程都是在無向圖中進行考慮的，現在考慮給邊定向，我們強制從較大編號的節點指向到較小編號的節點，不難發現這樣子這樣圖一定是一個DAG，這一點和我們在前面的算法中沿着矩陣的行從小往大依次消去是一致的。
圖中的一條邊表示了一個依賴關系，也就是如果我們要對於矩陣進行分解，那么如果(i,j)之間存在一條i->j的邊，意味着如果我們要計算在矩陣中消去i，我們必須先消去j。我們對其進行拓撲排序，不難發現由於圖聯通，且我們的圖一定是弦圖，因此最終得到的一定是一棵樹。這棵樹被稱作elimination tree，也即消去樹。
接下來我們以一個例子來介紹Multifrontal方法。
這是一個較為簡單的圖，G+是其通過前面介紹的算法添加新邊之后得到的弦圖。A對應着G，其中紅色點對應着G+中新加入的邊，也就是L相比於A多出的元素。T是其消去樹。
我們考慮在消去樹T中找到一些稀疏結構“相同”的節點。如果一個子節點在L的列中除了其本身之外和父親節點有相同的非零元素，那么這兩個節點就可以合並。如上面消去樹T中的節點5和6。
執行完這個合並步驟之后，我們得到這樣一棵樹T'
我們稱T'為Assembly Tree即，組裝樹。由於一個組中的節點具有類似的非零元素，故他們組合在一起代表着一個稠密的子矩陣。
假設在T'中的元素i'包含了T中的點集\(M=\{i,j,k...\}\)構成，則在T'中的節點i'對應的子矩陣是M中所有元素值在L（注意是L，它是一個下三角矩陣）中的所有非零元構成。比如果T'中的2'，其對應的\(M=\{2\}\)，因此其列上的非零元素包含了節點{2,5,6}。因此其對應的子矩陣就是由這幾個點構成的3X3矩陣。
在此之前，我們先考慮一個分解問題，我們先考慮一個2X2的分塊矩陣的Cholesky分解
![[公式]](/image/aHR0cHM6Ly93d3cuemhpaHUuY29tL2VxdWF0aW9uP3RleD0lNUNsZWZ0JTVCKyU1Q2JlZ2luJTdCYXJyYXklN0QrQUFfJTdCTU0lN0QlMjZBXyU3Qk5NJTdEJTVFVCslNUMlNUMrQV8lN0JOTSU3RCUyNitBXyU3Qk5OJTdEJTVDZW5kJTdCYXJyYXklN0QlNUNyaWdodCU1RCslM0QrJTVDbGVmdCU1QislNUNiZWdpbiU3QmFycmF5JTdEK0xMXyU3Qk1NJTdEJTI2KyU1QyU1QytMXyU3Qk5NJTdEJTI2K0xfJTdCTk4lN0QlNUNlbmQlN0JhcnJheSU3RCU1Q3JpZ2h0JTVEKyU1Q3RpbWVzK0wlNUVU.png)
於是有\(L_{MM}L_{MM}^T=A_{MM}\)，是子矩陣\(A_{MM}\)的Cholesky分解，\(L_{NM}=L_{MM}^{-1}A_{NM}\)是矩陣分解之后產生的一些中間元素。
重點考慮\(L_{NN}L_{NN}^T=A_{NN}-L_{NM}L_{NM}^T\)，這里可以類比前面過程中，矩陣減去當前行下方的向量的轉置這個過程。
於是，受到這個過程的啟發，我們按照順序來分解矩陣（順序即從葉子往根節點進行，這與我們消去樹中的依賴關系是一致的）。
上圖為消去1'節點的舉例。首先我們提出其對應的矩陣，然后加上其子節點傳來的矩陣元素（這里子節點，所以沒有這一步），然后直接計算求解這一個小矩陣的Cholesky分解。這一部分也就是前面的\(A_{MM}\)。接着計算\(L_{NM}=L^{-1}_{MM}A_{NM}\),最后計算\(A_{NN}-L_{NM}L_{NM^T}\)傳遞至父親節點。
例如我們給出\(A\)，這個對應的圖就是前面的例子。
顯而易見的，對於不同的分支，其計算在達到LCA之前都是獨立的，故所有的未合並且無依賴關系的分支都可以進行並行計算。