[心得] Interpolative Decomposition 分享

[mikemike1021 (mike)]

論壇版：https://forum.community.tw/t/topic/200
這裡是我自架的論壇，可以用 markdown, latex，並在旁邊建立列表，程式碼的部分也
可以自動上色
除此之外，如果有自己建立部落格也可以將此論壇來當做留言系統（可顯示在原頁中
目前以電腦相關討論開箱評測為主
希望這些功能能夠讓大家討論上更加方便，也歡迎大家註冊分享自己的心得、問問題或
幫忙解答XD

以下正文，會盡量將 latex 以方便 ptt 描述的方式打出來，也可以到前面的論壇觀看
不太確定發在這個版合不合適，但有用 python 來做驗證

直到前陣子才知道有這樣一種矩陣分解方式，之前知道的大概就是 QR, SVD, NMF 等這種
分解，但這些都沒有直接使用原始矩陣作為基底向量；然而最近暸解到的
Interpolative Decomposition - ID，他就是利用原始矩陣的向量直接來表示。

Interpolative Decomposition
給定矩陣 A(mxn)，在 rank r (r < rank(A)) 下，他的 Interpolative Decomposition
可表示為
A ~ A[:, J]Z
J 為 r 個索引，也就是從 A 挑出 r 個向量，Z 則為 rxn 矩陣，並含有一個單位子矩陣

原理概述
其實 ID 是可以從 QR 分解的結果再稍微組合一下組出 ID 的，因為在 QR 分解中，我們
是逐個向量去把 Q, R 組出來，Qk = sum(i=0->k)R{ik}A{:,k}，那把前面的部分組回去是不是就得到了原先的 A 向量了呢？但為了能得到更準確的估計，所以會加
上 Pivoting 或者說重新排列。這邊我們先假設 A 是不需要 pivoting 的因為我們只打
算用前面 A 的幾條向量來表示，所以 Q 只取 r 個，R 相對應的也把沒用掉的部分砍掉
A = Q(mxr)R(rxn) ()內為矩陣大小
就如同剛剛所說，前面得部分組回去就會得到原始的 A 向量，因此我們這裡在把 R 分
為 R1 前面 rxr 部分，以及 R2 後面的 rx(n-r)
A = Q[R1 R2] = QR1[I, R1^{-1}R2] = A1[I, R1^{-1}R2]
把 R1 題到外面去，跟 Q 結合後產生 A1，而由於 R1 是一個上三角矩陣，所以一定有他
的反矩陣。因此我們就有一個由 A 向量表示的分解
A = A1[I, R1^{-1}R2] = A1 Z

那如果 A 需要 pivoting 的狀況呢？
AP = A' = QR = A'1 Z' = AJ Z'
A = AJ (Z' P^*) = AJ Z
其實差不多，所用的 A 向量就變成被 pivoting 調到前面的那些向量
A'1 = AJ = A[:, J] 是 A' 的前幾個向量，也就是 A 某些向量。而最後 Z' 把
permutation matrix (置換矩陣) 給涵蓋進去，得到 Z 就可以了，那整體誤差其實就
跟 QR 分解一樣，這邊就不多談了

python 試試
那 scipy 其實有提供 ID  - scipy.linalg.interpolative 以及 QR 分解 -
scipy.linalg.qr 那我們就可以來用 python 來試試看上述的內容
python 程式碼
- 引入相關函式庫
import scipy.linalg.interpolative as sli
from scipy.linalg import hilbert
from scipy.linalg import qr
from numpy.linalg import inv
import numpy as np

- 初始化矩陣跟設定
# Set the matrix and setting
n = 5
# use hilbert to generate matrix
A = hilbert(n)
r = 3
print("A:")
print(A)
https://imgur.com/9cqfxiw.png

- 使用 scipy 的 Interpolative Decomposition
# Use scipy to generate the interpolative decomposition
idx, proj = sli.interp_decomp(A, r)
# A[:, idx[:r]] x proj = A[:, idx[r:]]
print("idx:")
print(idx)
print("proj:")
print(proj)
print("A[:, idx[:r]]:")
print(A[:, idx[:r]])
print("reconstruct A from ID:")
print((A[:, idx[:r]] @ np.hstack([np.eye(r), proj]))[:, np.argsort(idx)])
print("original A:")
print(A)
https://imgur.com/ZjBIGzA.png
可以看到 0, 2, 4 就是所用基底，所以他們 ID 所組出來的結果就會跟原始一樣，那
1, 3 就會是由這些 0, 2, 4 所組成，而係數就由 proj 來表述。

- 使用 scipy QR 來求 Interpolative Decomposition
# Use scipy QR to generate interpolative decomposition
Q, R, P = qr(A, pivoting=True)
print("P:")
print(P)
A_J = Q[:, :r] @ R[:r, :r]
print("A_J:")
print(A_J)
T = inv(R[:r, :r]) @ R[:r, r:n]
print("T:")
print(T)
# np.argsort(P) transpose the projection
Z = np.concatenate((np.eye(r), T), axis=1)[:, np.argsort(P)]
print("Z:")
print(Z)
print("A_J x Z:")
print(A_J @ Z)
print("A:")
print(A)
https://imgur.com/nlRSfDN.png

跟 ID 一樣可以看到 0, 2, 4 就是所用基底，所以他們 ID 所組出來的結果就會跟原始
一樣，那 1, 3 就會是由這些 0, 2, 4 所組成，而係數就由 T 來表述。

這邊 QR 跟 ID 給出的 pivoting 不一樣，應該是 ID 只要求算到 rank 3 所以後面就隨
便排了，如果把 ID 的 rank 改成 4 就會得到一樣的 pivoting，由於是 rank 外的所以
順序不會影響結果。 那由於 pivoting 的順序不一樣，那 T 其實也就跟 proj 不一樣，
但也只差在順序而已，係數部分是一樣的。

因為基底就是Ａ本身的某幾條向量，係數就是一個單位矩陣配 T，最後再把
permutation 倒回去，因此 Z 很明顯地含有一個單位子矩陣。

結語
一開始有人問到說要怎麼只利用 A 的向量做分解，且要近似，一開始只有想到 QR/SVD
可以做到近似，QR 前面的基底的確是用 A 的特定向量做成，但沒仔細想原來 QR 組一組
，就可以做到這件事 (ID)，且離 QR 並不遠。事後想一想這個過程也都很合理，能夠組
出來 A 的向量，那 A 的向量也可以組回去那些基底。過程並不複雜，即邊其他函式庫沒
有直接提供 ID，也可以從 QR 推導過去。另外，也有用 row 表示而非使用 column 的
ID，或者同時使用，詳情可以再閱讀參考資料暸解。那 ID 比起 QR 多的好處，就是直接
用 A 的向量表示分解後的結果，可能在解釋上會更加直觀。

參考資料 (連結也可直接從論壇版點選)
wiki - Interpolative decomposition
course note - The Interpolative Decomposition
On the Compression of Low Rank Matries
scipy Interpolative matrix decomposition
scipy QR

以上就是這次 Interpolative Decomposition 的心得分享
也希望論壇的功能有吸印到各位來試玩看看

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 112.78.81.99 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Python/M.1639681999.A.96F.html
推 allexj: 論壇可以使用 markdown / latex 是怎麼做的? 很好奇
→ mikemike1021: 我是用現有的把他架起來，markdown 本身他就可以用
→ mikemike1021: 了，另外裝主題的元件可以讓使用者插入自動生成的 m
→ mikemike1021: arkdown ，而另外裝 plugin 讓他可以用 mathjax 或
→ mikemike1021: kajax   將 latex 弄出來
→ allexj: 了解，我也來試試看