插值矩阵分解 (scipy.linalg.interpolative )¶
0.13 新版功能.
矩阵的插值分解(ID) \(A \in \mathbb{{C}}^{{m \times n}}\) 级别的 \(k \leq \min \{{ m, n \}}\) 是因式分解
哪里 \(\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}}\) 都是关于 \((k + 1)\) -第3个最大奇异值 \(A\) 。请注意， \(\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计算矩阵的奇异值分解。 |
|
用随机化幂法估计矩阵的谱范数。 |
|
用随机化幂法估计两个矩阵之差的谱范数。 |
|
使用随机化方法将矩阵秩估计到指定的相对精度。 |
支持功能：
|
为此ID包中使用的内部随机数生成器设定种子。 |
|
通过一种非常有效的滞后Fibonacci方法产生标准的均匀伪随机数。 |
参考文献¶
本模块使用ID软件包 [R5a82238cdab4-1] 由Martinsson、Rokhlin、Shkolnisky和Tygert编写，Tygert是一个Fortran库，用于使用各种算法计算ID，包括 [R5a82238cdab4-2] 中描述的较新的随机化方法 [R5a82238cdab4-3], [R5a82238cdab4-4], 和 [R5a82238cdab4-5]. 此模块以便于Python用户使用的方式公开其功能。请注意，除了组织更简单、更一致的界面之外，此模块不会添加任何功能。
我们建议用户同时查阅 documentation for the ID package 。
- R5a82238cdab4-1
首页--期刊主要分类--期刊细介绍--期刊题录与文摘--期刊详细文摘内容ID：用于通过插值分解进行矩阵低阶近似的软件包，版本0.2。http://tygert.com/id_doc.4.pdf.
- R5a82238cdab4-2
首页--期刊主要分类--期刊细介绍--期刊题录与文摘--期刊详细文摘内容“关于低秩矩阵的压缩。” 暹罗·J·科学。电脑。 26(4)：1389--1404,2005。 DOI:10.1137/030602678 。
- R5a82238cdab4-3
首页--期刊主要分类--期刊细介绍--期刊题录与文摘--期刊详细文摘内容“矩阵低阶近似的随机化算法。” 程序纳特尔。阿卡德。SCI。美国。 第104(51)：20167--20172,2007年。 DOI:10.1073/pnas.0709640104 。
- R5a82238cdab4-4
首页--期刊主要分类--期刊细介绍--期刊题录与文摘--期刊详细文摘内容“矩阵分解的随机化算法。” 应用程序。电脑。哈蒙。肛门。 30(1)：47--68,2011。 DOI:10.1016/j.acha.2010.02.003 。
- R5a82238cdab4-5
题名/责任者：Required/J.“一种矩阵近似的快速随机化算法。” 应用程序。电脑。哈蒙。肛门。 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(m):
>>> A[i,j] = 1. / (i + j + 1)
注意标志的用法 order='F' 在……里面 numpy.empty 。这以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计算为固定精度，请键入以下内容：
>>> 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)
计算奇异值分解¶
ID可以通过以下命令转换为SVD：
>>> U, S, V = sli.id_to_svd(B, idx, proj)
因此，SVD近似值为：
>>> C = np.dot(U, np.dot(np.diag(S), np.dot(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 具体如下：
>>> 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')
要指定种子值，请使用：
>>> sli.seed(s)
哪里 s 必须是整数或55个浮点数的数组。如果是整数，则通过使用 numpy.random.rand 使用给定的整数种子。
要简单地生成一些随机数，请键入以下内容：
>>> sli.rand(n)
哪里 n 要生成的随机数的数量。
备注¶
上述函数都会自动检测适当的接口，并处理实数和复数数据类型，将输入参数传递给适当的后端例程。