Python SciPy linalg.matmul_toeplitz用法及代码示例

本文简要介绍 python 语言中 scipy.linalg.matmul_toeplitz 的用法。

用法:

scipy.linalg.matmul_toeplitz(c_or_cr, x, check_finite=False, workers=None)#

使用 FFT 的高效 Toeplitz Matrix-Matrix 乘法

此函数返回 Toeplitz 矩阵和密集矩阵之间的矩阵乘法。

Toeplitz 矩阵有恒定的对角线，c 为第一列，r 为第一行。如果没有给出 r，则假定为 r == conjugate(c)。

参数 ：：

c_or_cr： 数组 或 (数组, 数组) 的元组

向量 c 或数组元组( c ， r )。无论 c 的实际形状如何，它将被转换为一维数组。如果未提供，则假定为 r = conjugate(c)；在这种情况下，如果 c[0] 是实数，则 Toeplitz 矩阵是 Hermitian。 r[0] 被忽略； Toeplitz 矩阵的第一行是 [c[0], r[1:]] 。无论 r 的实际形状如何，它将被转换为一维数组。

x： (M,) 或 (M, K) 数组

要与之相乘的矩阵。

check_finite： 布尔型，可选

是否检查输入矩阵是否仅包含有限数。禁用可能会提高性能，但如果输入确实包含无穷大或 NaN，则可能会导致问题(结果完全是 NaN)。

workers： 整数，可选

传递给scipy.fft.fft 和ifft。用于并行计算的最大工作线程数。如果为负，则该值从 os.cpu_count() 环绕。有关详细信息，请参阅scipy.fft.fft。

返回 ：：

T @ x： (M,) 或 (M, K) ndarray

矩阵乘法的结果T @ x.返回的形状匹配的形状x.

注意：

Toeplitz 矩阵嵌入循环矩阵中，FFT 用于有效计算matrix-matrix 乘积。

由于计算基于 FFT，整数输入将导致浮点输出。这与 NumPy 的 matmul 不同，后者保留了输入的数据类型。

这部分基于可在 [1] 中找到的实现，该实现在 MIT 许可下获得许可。有关该方法的更多信息，请参见引用 [2]。参考文献 [3] 和 [4] 在 Python 中有更多参考实现。

参考：

[1]

Jacob R Gardner、Geoff Pleiss、David Bindel、Kilian Q Weinberger、Andrew Gordon Wilson，“GPyTorch：Blackbox Matrix-Matrix GPU 加速的高斯过程推理”，Max Balandat 和 Ruihan Wu 贡献。在线提供：https://github.com/cornellius-gp/gpytorch

[2]

J. Demmel、P. Koev 和 X. Li，“直接线性求解器的简要调查”。在 Z. Bai、J. Demmel、J. Dongarra、A. Ruhe 和 H. van der Vorst，编辑。代数特征值问题求解的模板：实用指南。 SIAM，费城，2000 年。可在：http://www.netlib.org/utk/people/JackDongarra/etemplates/node384.html

[3]

R. Scheibler、E. Bezzam、I. Dokmanic、Pyroomacoustics：用于音频室模拟和阵列处理算法的 Python 包，Proc。 IEEE ICASSP，加利福尼亚州卡尔加里，2018 年。https://github.com/LCAV/pyroomacoustics/blob/pypi-release/pyroomacoustics/adaptive/util.py

[4]

Marano S、Edwards B、Ferrari G 和 Fah D(2017 年)，“拟合地震频谱：有色噪声和不完整数据”，美国地震学会公报，2017 年 1 月。 107(1)，第 276-291 页。

例子：

将 Toeplitz 矩阵 T 与矩阵 x 相乘：

[ 1 -1 -2 -3]       [1 10]
T = [ 3  1 -1 -2]   x = [2 11]
    [ 6  3  1 -1]       [2 11]
    [10  6  3  1]       [5 19]

要指定 Toeplitz 矩阵，只需要第一列和第一行。

>>> import numpy as np
>>> c = np.array([1, 3, 6, 10])    # First column of T
>>> r = np.array([1, -1, -2, -3])  # First row of T
>>> x = np.array([[1, 10], [2, 11], [2, 11], [5, 19]])

>>> from scipy.linalg import toeplitz, matmul_toeplitz
>>> matmul_toeplitz((c, r), x)
array([[-20., -80.],
       [ -7.,  -8.],
       [  9.,  85.],
       [ 33., 218.]])

通过创建完整的 Toeplitz 矩阵并将其乘以 x 来检查结果。

>>> toeplitz(c, r) @ x
array([[-20, -80],
       [ -7,  -8],
       [  9,  85],
       [ 33, 218]])

完整的矩阵从未明确形成，因此该例程适用于非常大的 Toeplitz 矩阵。

>>> n = 1000000
>>> matmul_toeplitz([1] + [0]*(n-1), np.ones(n))
array([1., 1., 1., ..., 1., 1., 1.])