TW Community

Interpolative Decomposition - ID

直到前陣子才知道有這樣一種矩陣分解方式,之前知道的大概就是 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 不一樣,但也只差在順序而已,係數部分是一樣的。

因為基底就是A本身的某幾條向量,係數就是一個單位矩陣配 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