Блочная трехдиагональная матрица питона

Question

Я хотел бы создать блочную трехдиагональную матрицу, начиная с трех numpy.ndarray. Есть ли (прямой) способ сделать это в Python?

Заранее спасибо!

Приветствия

Answer 1

С "обычными" массивами 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 2

Вы также можете сделать это с «обычными» массивами с помощью необычного индексирования:

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.arange(3)+4, np.arange(3)] используйте 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.(Обратите внимание, что вам нужно добавить фальшивые данные в значения строк, если вы помещаете данные в диагональную позицию с положительным значением (например, 3 в позиции 4 в примере))

В качестве быстрого примера:

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

Используйте функцию 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.]])

Answer 4

@ 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 5

Мой ответ строит из ответа @ 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)

Answer 6

Что бы там ни было, все остальные ответы, похоже, отвечают о трехдиагональных матрицах, а не block трехдиагональных матрицах.

Я не думаю, что есть встроенная поддержка трехдиагональных матриц, поэтому я написал свой собственный код. У меня были нули на главной диагонали, и моя матрица была симметричной.

Вот мой код.

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

Надеюсь, это поможет будущим людям, которые ищут block tridiagonal matrices.

Answer 7

Поскольку трехдиагональная матрица является разреженной матрицей, использование разреженного пакета может быть хорошим вариантом, см. http://pysparse.sourceforge.net/spmatrix.html#matlab-implementation,, есть некоторые примеры и сравнения с MATLAB даже ...