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

Interpolative Decomposition

給定矩陣 A (\( m \times n \)) ，在 rank r (\( r \leq rank(A) \)) 下，他的 Interpolative Decomposition 可以表示為

\[ A \approx A[:, J] Z \]

其中 J 為 r 個索引，也就是 A[:, J] 從 A 當中挑出 r 個向量出來，Z 則為 \( r \times n \) 的矩陣，且含有一個單位子矩陣

原理概述

其實 ID 是可以從 QR 分解的結果再稍微組合一下組出 ID 的，因為在 QR 分解中，我們是逐個向量去把 Q, R 組出來，\( Q_0 = R_{00} A_{:,0} \), \( Q_k = \sum_{i = 0}^k(R_{ik}A_{:,k}) \) ，那把前面的部分組回去是不是就得到了原先的 A 向量了呢？但為了能得到更準確的估計，所以會加上 Pivoting 或者說重新排列。這邊我們先假設 A 是不需要 pivoting 的

因為我們只打算用前面 A 的幾條向量來表示，所以 Q 只取 r 個，R 相對應的也把沒用掉的部分砍掉

\[ A = Q_{m\times r}R_{r\times n} \]

就如同剛剛所說，前面得部分組回去就會得到原始的 A 向量，因此我們這裡在把 R 分成 \(R_1\) \(R_2\), \(R_1\) 是前面 \( r \times r \) 的部分，那 \(R_2\) 就是後面 \( r \times (n-r) \) 了

\[ A = Q [R_1, R_2] = QR_1[I, R_1^{-1}R_2] = A_1 [I, R_1^{-1}R_2] \]

我們把 \(R_1\) 提到外面去，跟 Q 結合後就會產生 \(A_1\)，而由於 \(R_1\)是一個上三角矩陣，所以一定有他的反矩陣。因此我們就得到了一個由 A 向量表示的分解

\[ A = A_1 [ I, R_1^{-1}R_2] = A_1 Z \]

那如果 A 是需要 pivoting 的狀況呢？

\[ AP = A' = QR = A'_1 Z' = A_J Z' \]

\[ A = A_J (Z' P^*) = A_J Z \]

其實差不多，所用的 A 向量就變成被 pivoting 調到前面的那些向量，\(A'_1=A_J = 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)

• 使用 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)

可以看到 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)

跟 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