插值矩阵分解 (scipy.linalg.interpolative)#
版本 0.13 新增。
秩为 \(k \leq \min \{ m, n \}\) 的矩阵 \(A \in \mathbb{C}^{m \times n}\) 的插值分解 (ID) 是一个因子分解
其中 \(\Pi = [\Pi_{1}, \Pi_{2}]\) 是一个置换矩阵，其中 \(\Pi_{1} \in \{ 0, 1 \}^{n \times k}\)，即 \(A \Pi_{2} = A \Pi_{1} T\)。这等效于写成 \(A = BP\)，其中 \(B = A \Pi_{1}\) 和 \(P = [I, T] \Pi^{\mathsf{T}}\) 分别是骨架和插值矩阵。
如果 \(A\) 没有精确的秩 \(k\)，那么存在一个以 ID 形式的近似值，使得 \(A = BP + E\)，其中 \(\| E \| \sim \sigma_{k + 1}\) 与 \(A\) 的第 \((k + 1)\) 个最大奇异值的数量级相同。请注意，\(\sigma_{k + 1}\) 是秩为 \(k\) 的近似的最佳可能误差，并且实际上，由奇异值分解 (SVD) \(A \approx U S V^{*}\) 实现，其中 \(U \in \mathbb{C}^{m \times k}\) 和 \(V \in \mathbb{C}^{n \times k}\) 具有正交列，并且 \(S = \mathop{\mathrm{diag}} (\sigma_{i}) \in \mathbb{C}^{k \times k}\) 是对角线，其对角线元素非负。与 SVD 相比，使用 ID 的主要优点是
构造成本更低；
它保留了 \(A\) 的结构；并且
鉴于 \(P\) 的单位子矩阵，它在计算方面更有效。
例程#
主要功能
|
计算矩阵的 ID。 |
|
从其 ID 重构矩阵。 |
|
从 ID 重构插值矩阵。 |
|
从 ID 重构骨架矩阵。 |
|
将 ID 转换为 SVD。 |
|
通过 ID 计算矩阵的 SVD。 |
|
通过随机幂法估计矩阵的谱范数。 |
|
通过随机幂法估计两个矩阵差的谱范数。 |
|
使用随机方法将矩阵秩估计到指定的相对精度。 |
支持函数
参考文献#
此模块使用 Martinsson、Rokhlin、Shkolnisky 和 Tygert 编写的 ID 软件包 [R5a82238cdab4-1]，这是一个用于使用各种算法计算 ID 的 Fortran 库，包括 [R5a82238cdab4-2] 中的秩揭示 QR 方法以及 [R5a82238cdab4-3]、[R5a82238cdab4-4] 和 [R5a82238cdab4-5] 中描述的更近期的随机方法。此模块以适合 Python 用户使用的方式公开其功能。请注意，此模块除了组织更简单和更一致的接口之外，没有添加任何其他功能。
我们建议用户也参考 ID 包的文档。
P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert。 “ID：通过插值分解对矩阵进行低秩近似的软件包，版本 0.2。” http://tygert.com/id_doc.4.pdf。
H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin。 “关于低秩矩阵的压缩。” SIAM J. Sci. Comput. 26 (4): 1389–1404, 2005. DOI:10.1137/030602678。
E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, M. Tygert。 “用于矩阵的低秩近似的随机算法。” Proc. Natl. Acad. Sci. U.S.A. 104 (51): 20167–20172, 2007. DOI:10.1073/pnas.0709640104。
P.G. Martinsson, V. Rokhlin, M. Tygert。 “用于矩阵分解的随机算法。” Appl. Comput. Harmon. Anal. 30 (1): 47–68, 2011. DOI:10.1016/j.acha.2010.02.003。
F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert。 “用于矩阵近似的快速随机算法。” Appl. Comput. Harmon. Anal. 25 (3): 335–366, 2008. DOI:10.1016/j.acha.2007.12.002。
教程#
初始化#
第一步是通过发出以下命令导入 scipy.linalg.interpolative
>>> import scipy.linalg.interpolative as sli
现在让我们构建一个矩阵。为此，我们考虑一个希尔伯特矩阵，该矩阵众所周知具有低秩
>>> from scipy.linalg import hilbert
>>> n = 1000
>>> A = hilbert(n)
我们也可以通过以下方式显式地执行此操作
>>> import numpy as np
>>> n = 1000
>>> A = np.empty((n, n), order='F')
>>> for j in range(n):
... for i in range(n):
... A[i,j] = 1. / (i + j + 1)
请注意在 numpy.empty 中使用标志 order='F'。这会以 Fortran 连续顺序实例化矩阵，并且对于在传递到后端时避免数据复制非常重要。
然后，我们通过将其视为一个 scipy.sparse.linalg.LinearOperator 来定义矩阵的乘法例程
>>> from scipy.sparse.linalg import aslinearoperator
>>> L = aslinearoperator(A)
这会自动设置描述矩阵及其伴随矩阵对向量的作用的方法。
计算 ID#
我们有几种算法选择来计算 ID。这些算法主要根据两个二分法进行划分
矩阵是如何表示的，即通过其条目还是通过其对向量的作用；以及
是将其近似到固定的相对精度还是近似到固定的秩。
我们将在下面依次介绍每种选择。
在所有情况下，ID 由三个参数表示
一个秩
k；一个索引数组
idx；以及插值系数
proj。
ID 由关系 np.dot(A[:,idx[:k]], proj) == A[:,idx[k:]] 指定。
从矩阵条目#
我们首先考虑一个以其条目形式给出的矩阵。
要计算到固定精度的 ID，请键入
>>> eps = 1e-3
>>> k, idx, proj = sli.interp_decomp(A, eps)
其中 eps < 1 是所需的精度。
要计算到固定秩的 ID，请使用
>>> idx, proj = sli.interp_decomp(A, k)
其中 k >= 1 是所需的秩。
两种算法都使用随机抽样，通常比相应的旧的确定性算法更快，可以通过以下命令访问
>>> k, idx, proj = sli.interp_decomp(A, eps, rand=False)
和
>>> idx, proj = sli.interp_decomp(A, k, rand=False)
分别。
从矩阵操作#
现在考虑一个矩阵，它根据其对向量的操作给出为一个scipy.sparse.linalg.LinearOperator。
要计算到固定精度的 ID，请键入
>>> k, idx, proj = sli.interp_decomp(L, eps)
要计算到固定秩的 ID，请使用
>>> idx, proj = sli.interp_decomp(L, k)
这些算法是随机的。
重建 ID#
上面的 ID 例程不会显式输出骨架和插值矩阵，而是以更紧凑（有时更有用）的形式返回相关信息。要构建这些矩阵，请编写
>>> B = sli.reconstruct_skel_matrix(A, k, idx)
用于骨架矩阵，以及
>>> P = sli.reconstruct_interp_matrix(idx, proj)
用于插值矩阵。然后可以计算 ID 近似值，如下所示
>>> C = np.dot(B, P)
这也可以使用以下方法直接构建
>>> C = sli.reconstruct_matrix_from_id(B, idx, proj)
无需先计算 P。
或者，也可以使用以下方法显式完成
>>> B = A[:,idx[:k]]
>>> P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
>>> C = np.dot(B, P)
计算 SVD#
可以通过以下命令将 ID 转换为 SVD
>>> U, S, V = sli.id_to_svd(B, idx, proj)
然后 SVD 近似值为
>>> approx = U @ np.diag(S) @ V.conj().T
还可以通过将 ID 和转换步骤组合到一个命令中来“全新”计算 SVD。在上面的各种 ID 算法之后，相应地存在可以使用各种 SVD 算法。
从矩阵条目#
我们首先考虑用于根据其条目给出的矩阵的 SVD 算法。
要计算到固定精度的 SVD，请键入
>>> U, S, V = sli.svd(A, eps)
要计算到固定秩的 SVD，请使用
>>> U, S, V = sli.svd(A, k)
两种算法都使用随机抽样；对于确定性版本，请如上所述发出关键字 rand=False。
从矩阵操作#
现在考虑一个矩阵，它根据其对向量的操作给出。
要计算到固定精度的 SVD，请键入
>>> U, S, V = sli.svd(L, eps)
要计算到固定秩的 SVD，请使用
>>> U, S, V = sli.svd(L, k)
实用例程#
还提供了一些实用例程。
要估计矩阵的谱范数，请使用
>>> snorm = sli.estimate_spectral_norm(A)
该算法基于随机幂方法，因此只需要矩阵向量乘积。可以使用关键字 its 设置要执行的迭代次数（默认值：its=20）。该矩阵被解释为一个scipy.sparse.linalg.LinearOperator，但也可以将其提供为numpy.ndarray，在这种情况下，它使用 scipy.sparse.linalg.aslinearoperator 微不足道地转换。
相同的算法还可以估计两个矩阵 A1 和 A2 的差的谱范数，如下所示
>>> A1, A2 = A**2, A
>>> diff = sli.estimate_spectral_norm_diff(A1, A2)
这通常用于检查矩阵近似的准确性。
中的一些例程scipy.linalg.interpolative 还需要估计矩阵的秩。这可以使用以下两种方法完成
>>> k = sli.estimate_rank(A, eps)
或
>>> k = sli.estimate_rank(L, eps)
取决于表示形式。参数 eps 控制数值秩的定义。
最后，可以通过 scipy.linalg.interpolative.seed 控制所有随机例程所需的随机数生成。要将种子值重置为其原始值，请使用
>>> sli.seed('default')
要指定种子值，请使用
>>> s = 42
>>> sli.seed(s)
其中 s 必须是整数或 55 个浮点数的数组。如果为整数，则浮点数数组是通过使用 numpy.random.rand 使用给定的整数种子获得的。
要简单地生成一些随机数，请键入
>>> arr = sli.rand(n)
其中 n 是要生成的随机数的个数。
备注#
以上所有函数都会自动检测合适的接口，并使用实数和复数数据类型，将输入参数传递给适当的后端例程。