Ctrl+K

插值矩阵分解 (scipy.linalg.interpolative)#

在 0.13 版本中添加。

在 1.15.0 版本中更改: 底层算法已从原始 Fortran77 代码移植到 Python。 详情请参阅下面的参考文献。

秩为 \(k \leq \min \{ m, n \}\) 的矩阵 \(A \in \mathbb{C}^{m \times n}\) 的插值分解 (ID) 是一个因式分解

\[A \Pi = \begin{bmatrix} A \Pi_{1} & A \Pi_{2} \end{bmatrix} = A \Pi_{1} \begin{bmatrix} I & T \end{bmatrix},\]

其中 \(\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\) 的单位子矩阵，计算起来效率更高。

例程#

主要功能

interp_decomp(A, eps_or_k[, rand, rng])

计算矩阵的 ID。

reconstruct_matrix_from_id(B, idx, proj)

从其 ID 重构矩阵。

reconstruct_interp_matrix(idx, proj)

从 ID 重构插值矩阵。

reconstruct_skel_matrix(A, k, idx)

从 ID 重构骨架矩阵。

id_to_svd(B, idx, proj)

将 ID 转换为 SVD。

svd(A, eps_or_k[, rand, rng])

通过 ID 计算矩阵的 SVD。

estimate_spectral_norm(A[, its, rng])

通过随机幂方法估计矩阵的谱范数。

estimate_spectral_norm_diff(A, B[, its, rng])

通过随机幂方法估计两个矩阵之差的谱范数。

estimate_rank(A, eps[, rng])

使用随机方法估计矩阵的秩，达到指定的相对精度。

以下支持函数已弃用，将在 SciPy 1.17.0 中删除

seed([seed])

此函数在历史上曾用于设置 Fortran77 中编写的 scipy.linalg.interpolative 函数中使用的随机化算法的种子。

rand(*shape)

此函数在历史上曾用于为 Fortran77 中编写的 scipy.linalg.interpolative 函数中使用的随机化算法生成均匀分布的随机数。

参考文献#

此模块使用 Martinsson、Rokhlin、Shkolnisky 和 Tygert 的 ID 软件包 [R5a82238cdab4-1] 中找到的算法，该软件包是一个 Fortran 库，用于使用各种算法计算 ID，包括 [R5a82238cdab4-2] 的秩揭示 QR 方法和 [R5a82238cdab4-3]、[R5a82238cdab4-4] 和 [R5a82238cdab4-5] 中描述的更新的随机方法。

我们建议用户查阅 ID 软件包的文档。

[R5a82238cdab4-1]

P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert. “ID：用于通过插值分解对矩阵进行低秩近似的软件包，版本 0.2。” http://tygert.com/id_doc.4.pdf。

[R5a82238cdab4-2]

H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. “关于低秩矩阵的压缩。” SIAM J. Sci. Comput. 26 (4): 1389–1404, 2005. DOI:10.1137/030602678。

[R5a82238cdab4-3]

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。

[R5a82238cdab4-4]

P.G. Martinsson, V. Rokhlin, M. Tygert. “用于矩阵分解的随机算法。” Appl. Comput. Harmon. Anal. 30 (1): 47–68, 2011. DOI:10.1016/j.acha.2010.02.003。

[R5a82238cdab4-5]

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。这些算法大致可以根据两个二分法进行分类

  1. 矩阵的表示方式，即通过其条目还是通过其对向量的作用；以及

  2. 是将其近似到固定的相对精度还是固定的秩。

下面我们将依次介绍每个选项。

在所有情况下，ID 由三个参数表示

  1. 一个秩 k；

  2. 一个索引数组 idx；以及

  3. 插值系数 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#

>>> U, S, V = sli.id_to_svd(B, idx, proj)

可以通过以下命令将 ID 转换为 SVD

>>> approx = U @ np.diag(S) @ V.conj().T

然后 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。

备注#

上述函数都会自动检测相应的接口，并使用实数和复数数据类型，将输入参数传递给正确的后端例程。

在本页