Python求解T形矩阵(toeplitz矩阵)

文章目录

简介
matmul_toeplitz
solve_toeplitz

简介

T形矩阵，即托普利兹矩阵，其主对角线上的元素相等，平行于主对角线的线上的元素也相等；矩阵中的各元素关于次对角线对称，即T型矩阵为次对称矩阵。

在scipy.linalg中，可通过toeplitz函数来创建T形矩阵

import numpy as np
import scipy.linalg as slm = 10             # 中心项索引
d = np.arange(1, 2*(m+1))
d2 = d[::-1]T = sl.toeplitz(d[m:], d2[m:])
print(T)

结果为

[ 11 10 9 8 7 6 5 4 3 2 1 12 11 10 9 8 7 6 5 4 3 2 13 12 11 10 9 8 7 6 5 4 3 14 13 12 11 10 9 8 7 6 5 4 15 14 13 12 11 10 9 8 7 6 5 16 15 14 13 12 11 10 9 8 7 6 17 16 15 14 13 12 11 10 9 8 7 18 17 16 15 14 13 12 11 10 9 8 19 18 17 16 15 14 13 12 11 10 9 20 19 18 17 16 15 14 13 12 11 10 21 20 19 18 17 16 15 14 13 12 11 ] \begin{bmatrix} 11&10& 9& 8& 7& 6& 5& 4& 3& 2& 1\\ 12&11&10& 9& 8& 7& 6& 5& 4& 3& 2\\ 13&12&11&10& 9& 8& 7& 6& 5& 4& 3\\ 14&13&12&11&10& 9& 8& 7& 6& 5& 4\\ 15&14&13&12&11&10& 9& 8& 7& 6& 5\\ 16&15&14&13&12&11&10& 9& 8& 7& 6\\ 17&16&15&14&13&12&11&10& 9& 8& 7\\ 18&17&16&15&14&13&12&11&10& 9& 8\\ 19&18&17&16&15&14&13&12&11&10& 9\\ 20&19&18&17&16&15&14&13&12&11&10\\ 21&20&19&18&17&16&15&14&13&12&11\\ \end{bmatrix} 11121314151617181920211011121314151617181920910111213141516171819891011121314151617187891011121314151617678910111213141516567891011121314154567891011121314345678910111213234567891011121234567891011

可以将其绘制出来

import matplotlib.pyplot as plt
plt.imshow(t/21)
plt.show()

matmul_toeplitz

由于T型矩阵的特殊形态，所以在解决一些问题时，往往有更加便捷的方法。scipy为T形矩阵提供了两个特殊的函数，分别用于矩阵乘法和 T x = − r Tx=-r Tx=−r问题的求解。

其中，matmul_toeplitz为矩阵乘法函数

import scipy.linalg as sl
import numpy as npc = np.array([1, 3, 6, 10])
r = np.array([1, -1, -2, -3])
b = np.array([1, 2, 2, 5])
T = sl.toeplitz(c,r)np.random.seed(42)
xs = np.random.rand(4)
print(xs)
# array([0.37454012, 0.95071431, 0.73199394, 0.59865848])print(T@xs)
# array([-51.17958147, -41.08476258, -29.59336213,  -9.54874811])
print(matmul_toeplitz((c,r), xs))
# array([-51.17958147, -41.08476258, -29.59336213,  -9.54874811])

solve_toeplitz

对于形如 T x = − r Tx=-r Tx=−r这种问题，其中 T T T为T形矩阵， r r r为已知量，可通过Levinson-Durbin算法对方程进行求解。

在scipy.linalg中，solve_toeplitz函数实现了这个算法，其输入参数主要是c_or_cr和x，其中c_or_r是一个元组，表示T的第一列和第一行。此外还有一个参数是有限校验了。

下面做一下测试

s = T@xs
c, r = np.arange(11,22), np.arange(11,0,-1)
sl.solve_toeplitz((c,r), s)
# array([0.37454012, 0.95071431, 0.73199394, 0.59865848])

这个值与xs完全相同。