插值矩阵分解 (scipy.linalg.interpolative
)#
在 0.13 版中新增。
在 1.15.0 版中更改:底层算法已从原始 Fortran77 代码移植到 Python。更多详情请参阅下面的参考文献。
矩阵 \(A \in \mathbb{C}^{m \times n}\) 的秩为 \(k \leq \min \{ m, 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}\) 是带有非负项的对角矩阵。使用 ID 而非 SVD 的主要优点是:
构建成本更低;
它保留了 \(A\) 的结构;并且
鉴于 \(P\) 的单位子矩阵,用它进行计算更有效率。
例程#
主要功能
|
计算矩阵的 ID。 |
|
从其 ID 重构矩阵。 |
|
从 ID 重构插值矩阵。 |
|
从 ID 重构骨架矩阵。 |
|
将 ID 转换为 SVD。 |
|
通过 ID 计算矩阵的 SVD。 |
|
通过随机幂法估计矩阵的谱范数。 |
|
通过随机幂法估计两个矩阵差的谱范数。 |
|
使用随机方法,以指定的相对精度估计矩阵秩。 |
以下支持函数已弃用,并将在 SciPy 1.17.0 中移除
|
此函数,在历史上,用于设置 Fortran77 编写的 |
|
此函数,在历史上,用于为 Fortran77 编写的 |
参考文献#
本模块使用 Martinsson、Rokhlin、Shkolnisky 和 Tygert 的 ID 软件包 [R5a82238cdab4-1] 中的算法,这是一个用于使用各种算法计算 ID 的 Fortran 库,包括 [R5a82238cdab4-2] 的秩显式 QR 方法和 [R5a82238cdab4-3]、[R5a82238cdab4-4] 和 [R5a82238cdab4-5] 中描述的最新随机方法。
我们建议用户也查阅 ID 软件包的文档。
P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert. “ID: a software package for low-rank approximation of matrices via interpolative decompositions, version 0.2.” http://tygert.com/id_doc.4.pdf.
H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. “On the compression of low rank matrices.” SIAM J. Sci. Comput. 26 (4): 1389–1404, 2005. DOI:10.1137/030602678.
E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, M. Tygert. “Randomized algorithms for the low-rank approximation of matrices.” Proc. Natl. Acad. Sci. U.S.A. 104 (51): 20167–20172, 2007. DOI:10.1073/pnas.0709640104.
P.G. Martinsson, V. Rokhlin, M. Tygert. “A randomized algorithm for the decomposition of matrices.” 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. “A fast randomized algorithm for the approximation of matrices.” 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
SVD 也可以通过将 ID 和转换步骤合并到一个命令中来“重新”计算。按照上面各种 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
控制数值秩的定义。
最后,所有随机例程所需的随机数生成都可以通过提供带有固定种子的 NumPy 伪随机生成器来控制。有关详细信息,请参阅 numpy.random.Generator
和 numpy.random.default_rng
。
备注#
上述所有函数都会自动检测合适的接口,并支持实数和复数数据类型,将输入参数传递给适当的后端例程。