论坛版：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 的心得分享
也希望论坛的功能有吸印到各位来试玩看看