论坛版：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://webptt.com/cn.aspx?n=bbs/Python/M.1639681999.A.96F.html