塊三對角矩陣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.]])