块三对角矩阵python

Posted 2023-03-12

【中文标题】块三对角矩阵python

【英文标题】：Block tridiagonal matrix python

【发布时间】：2011-08-16 02:53:33

【问题描述】：

我想从三个 numpy.ndarray 开始创建一个块三对角矩阵。 在 python 中是否有任何（直接）方法可以做到这一点？

提前谢谢你！

干杯

【问题讨论】：

您希望结果是另一个 ndarray，还是愿意为结果使用稀疏数组？

这个问题应该更准确，应该提供一个框架示例。

【参考方案1】：

为了从三个单独的块构建一个块状三对角矩阵（并将这些块重复 N 次），一种解决方案可以是：

import numpy as np
from scipy.linalg import block_diag

def tridiag(c, u, d, N): 
    # c, u, d are center, upper and lower blocks, repeat N times
    cc = block_diag(*([c]*N)) 
    shift = c.shape[1]
    uu = block_diag(*([u]*N)) 
    uu = np.hstack((np.zeros((uu.shape[0], shift)), uu[:,:-shift]))
    dd = block_diag(*([d]*N)) 
    dd = np.hstack((dd[:,shift:],np.zeros((uu.shape[0], shift))))
    return cc+uu+dd

例如：

c = np.matrix([[1,1],[1,1]])
u = np.matrix([[2,2],[2,2]])
d = -1*u
N =4
H = tridiag(c,u,d,N)
print(H)

给出答案

[[ 1.  1.  2.  2.  0.  0.  0.  0.]
 [ 1.  1.  2.  2.  0.  0.  0.  0.]
 [-2. -2.  1.  1.  2.  2.  0.  0.]
 [-2. -2.  1.  1.  2.  2.  0.  0.]
 [ 0.  0. -2. -2.  1.  1.  2.  2.]
 [ 0.  0. -2. -2.  1.  1.  2.  2.]
 [ 0.  0.  0.  0. -2. -2.  1.  1.]
 [ 0.  0.  0.  0. -2. -2.  1.  1.]]

【参考方案2】：

使用函数scipy.sparse.diags。

例子：

from scipy.sparse import diags
import numpy as np

n = 10
k = [np.ones(n-1),-2*np.ones(n),np.ones(n-1)]
offset = [-1,0,1]
A = diags(k,offset).toarray()

这会返回：

array([[-2.,  1.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  1., -2.]])

【参考方案3】：

我只是需要类似的东西。我也分享了我的解决方法。

from  functools import reduce
import numpy as np
size = 5
vals = [1,3,4]
res = reduce(lambda a, b: a+np.eye(size,k=b[0])*b[1], enumerate(vals), np.zeros((size, size)))
print(res)
# [[1. 3. 4. 0. 0.]
#  [0. 1. 3. 4. 0.]
#  [0. 0. 1. 3. 4.]
#  [0. 0. 0. 1. 3.]
#  [0. 0. 0. 0. 1.]]

size = 7
vals = [1,3,4,7,6]
offset = -2
res = reduce(lambda a, b: a+np.eye(size,k=b[0]+offset)*b[1], enumerate(vals), np.zeros((size, size)))
print(res)
# [[4. 7. 6. 0. 0. 0. 0.]
#  [3. 4. 7. 6. 0. 0. 0.]
#  [1. 3. 4. 7. 6. 0. 0.]
#  [0. 1. 3. 4. 7. 6. 0.]
#  [0. 0. 1. 3. 4. 7. 6.]
#  [0. 0. 0. 1. 3. 4. 7.]
#  [0. 0. 0. 0. 1. 3. 4.]]

【参考方案4】：

无论好坏，所有其他答案似乎都是关于三对角矩阵而不是块三对角矩阵。

我认为没有对三对角矩阵的原生支持，所以我编写了自己的代码。我在主对角线上有零，我的矩阵是对称的。

这是我的代码。

n1 = 784
n2 = 256
n3 = 128
n4 = 10
M1 = np.ones((n1,n2))
M2 = np.ones((n2,n3))
M3 = np.ones((n3, n4))

def blockTri(Ms):
    #Takes in a list of matrices (not square) and returns a tridiagonal block matrix with zeros on the diagonal
    count = 0
    idx = []
    for M in Ms:
        #print(M.shape)
        count += M.shape[0]
        idx.append(count)
    count += Ms[-1].shape[-1]
    mat = np.zeros((count,count))
    count = 0
    for i, M in enumerate(Ms):
        mat[count:count+M.shape[0],idx[i]:idx[i]+M.shape[1]] = M
        count = count + M.shape[0]
    mat = mat + mat.T    
    return mat

M = blockTri([M1, M2, M3])

希望这可以帮助未来寻找块三对角矩阵的人。

【参考方案5】：

@TheCorwoodRep 的回答实际上可以在一行中完成。不需要单独的函数。

np.eye(3,3,k=-1) + np.eye(3,3)*2 + np.eye(3,3,k=1)*3

这会产生：

array([[ 2.,  3.,  0.],
       [ 1.,  2.,  3.],
       [ 0.,  1.,  2.]])

【参考方案6】：

我的答案建立在@TheCorwoodRep 的答案之上。我只是发布它，因为我做了一些更改以使其更加模块化，以便它适用于不同的矩阵顺序，并更改 k1,k2,k3 的值，即决定对角线出现的位置, 会自动处理溢出。在调用该函数时，您可以指定哪些值应出现在对角线上。

import numpy as np
def tridiag(T,x,y,z,k1=-1, k2=0, k3=1):
    a = [x]*(T-abs(k1)); b = [y]*(T-abs(k2)); c = [z]*(T-abs(k3))
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

D=tridiag(10,-1,2,-1)

【参考方案7】：

使用“常规”numpy 数组，使用numpy.diag：

def tridiag(a, b, c, k1=-1, k2=0, k3=1):
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

a = [1, 1]; b = [2, 2, 2]; c = [3, 3]
A = tridiag(a, b, c)

【参考方案8】：

您也可以通过花哨的索引对“常规”numpy 数组执行此操作：

import numpy as np
data = np.zeros((10,10))
data[np.arange(5), np.arange(5)+2] = [5, 6, 7, 8, 9]
data[np.arange(3)+4, np.arange(3)] = [1, 2, 3]
print data

（如果您想更简洁，可以将那些对np.arange 的调用替换为np.r_。例如，使用data[np.r_[:3]+4, np.r_[:3]] 代替data[np.arange(3)+4, np.arange(3)]）

这会产生：

[[0 0 5 0 0 0 0 0 0 0]
 [0 0 0 6 0 0 0 0 0 0]
 [0 0 0 0 7 0 0 0 0 0]
 [0 0 0 0 0 8 0 0 0 0]
 [1 0 0 0 0 0 9 0 0 0]
 [0 2 0 0 0 0 0 0 0 0]
 [0 0 3 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]]

但是，如果您仍然要使用稀疏矩阵，请查看scipy.sparse.spdiags。 （请注意，如果您将数据放置到具有正值的对角线位置（例如示例中位置 4 中的 3），则需要在行值上 预先添加假数据） p>

举个简单的例子：

import numpy as np
import scipy as sp
import scipy.sparse

diag_rows = np.array([[1, 1, 1, 1, 1, 1, 1],
                      [2, 2, 2, 2, 2, 2, 2],
                      [0, 0, 0, 0, 3, 3, 3]])
positions = [-3, 0, 4]
print sp.sparse.spdiags(diag_rows, positions, 10, 10).todense()

这会产生：

[[2 0 0 0 3 0 0 0 0 0]
 [0 2 0 0 0 3 0 0 0 0]
 [0 0 2 0 0 0 3 0 0 0]
 [1 0 0 2 0 0 0 0 0 0]
 [0 1 0 0 2 0 0 0 0 0]
 [0 0 1 0 0 2 0 0 0 0]
 [0 0 0 1 0 0 2 0 0 0]
 [0 0 0 0 1 0 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 1 0 0 0]]

【参考方案9】：

由于三对角矩阵是一个稀疏矩阵，使用稀疏包可能是一个不错的选择，请参阅http://pysparse.sourceforge.net/spmatrix.html#matlab-implementation，甚至还有一些示例和与 MATLAB 的比较...