块三对角矩阵python

Question

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

预先感谢您！

干杯

Answer 1

由于三角矩阵是用一个稀疏的包应该是一个不错的选择，看http://pysparse.sourceforge.net/spmatrix.html#matlab-implementation，也有一些实例和比较与MATLAB甚至...

Answer 2

你也可以做到这一点稀疏矩阵通过花哨的索引“常规” 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.r_如果你想更简洁替换那些调用np.arange例如，而不是。，使用data[np.r_[:3]+4, np.r_[: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位），则需要将前缀假数据拖到您的行值上）

作为快速例如：

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]] 

Answer 3

用 “常规” 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) 

Answer 4

我的答案建立@ TheCorwoodRep的答案。我只是发布它，因为我做了一些更改，使它更加模块化，以便它可以用于不同的矩阵顺序，并且还可以更改​​,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) 

Answer 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.]]) 

Answer 6

使用功能scipy.sparse.diags。

例子：

from scipy.sparse import diags 
import numpy as np 
# 
n = 10 
k = np.array([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.]])