內插矩陣分解 (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}\) 與 \((k + 1)\)-th 個最大奇異值 \(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\) 的單位子矩陣，計算效率更高。

常式

主要功能

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 函數中使用的隨機演算法產生均勻分佈的隨機數。

參考文獻

此模組使用 ID 軟體套件 [R5a82238cdab4-1] 中的演算法，該套件由 Martinsson、Rokhlin、Shkolnisky 和 Tygert 開發，它是一個 Fortran 程式庫，用於使用各種演算法計算 ID，包括 [R5a82238cdab4-2] 的秩揭示 QR 方法，以及 [R5a82238cdab4-3]、[R5a82238cdab4-4] 和 [R5a82238cdab4-5] 中描述的較新的隨機方法。

我們建議使用者也查閱 ID 套件 的文件。

[R5a82238cdab4-1]

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.

[R5a82238cdab4-2]

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.

[R5a82238cdab4-3]

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.

[R5a82238cdab4-4]

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.

[R5a82238cdab4-5]

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。這些演算法主要分為兩種二分法

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

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。

備註

上述函數都會自動偵測適當的介面，並適用於實數和複數資料類型，將輸入引數傳遞到正確的後端常式。