Создайте разреженную циркулянтную матрицу в python

Question

Я хочу создать в Python большую (скажем, 10^5 x 10^5) разреженную циркулянтную матрицу. Он имеет 4 элемента в строке в позициях [i,i+1], [i,i+2], [i,i+N-2], [i,i+N-1], где я принял периодические граничные условия для индексов (т. Е. [10^5,10^5]=[0,0], [10^5+1,10^5+1]=[1,1] и т. Д.). Я посмотрел документацию на scipy редких матрицах, но я очень смущен (я новичок в Python).

я могу создать матрицу с NumPy

import numpy as np 

def Bc(i, boundary): 
    """(int, int) -> int 

    Checks boundary conditions on index 
    """ 
    if i > boundary - 1: 
     return i - boundary 
    elif i < 0: 
     return boundary + i 
    else: 
     return i 

N = 100 
diffMat = np.zeros([N, N]) 
for i in np.arange(0, N, 1): 
    diffMat[i, [Bc(i+1, N), Bc(i+2, N), Bc(i+2+(N-5)+1, N), Bc(i+2+(N-5)+2, N)]] = [2.0/3, -1.0/12, 1.0/12, -2.0/3] 

Однако, это довольно медленно и для больших N использует много памяти, поэтому я хочу, чтобы избежать создания с NumPy и преобразование в разреженные матрицы и перейти непосредственно к последнему.

Я знаю, как это сделать в Mathematica, где можно использовать SparseArray и шаблоны индексов - здесь что-то похожее?

Answer 1

Чтобы создать плотный циркулянтная матрица, вы можете использовать scipy.linalg.circulant. Например,

In [210]: from scipy.linalg import circulant 

In [211]: N = 7 

In [212]: vals = np.array([2.0/3, -1.0/12, 1.0/12, -2.0/3]) 

In [213]: offsets = np.array([1, 2, N-2, N-1]) 

In [214]: col0 = np.zeros(N) 

In [215]: col0[offsets] = -vals 

In [216]: c = circulant(col0) 

In [217]: c 
Out[217]: 
array([[ 0. , 0.6667, -0.0833, 0. , 0. , 0.0833, -0.6667], 
     [-0.6667, 0. , 0.6667, -0.0833, 0. , 0. , 0.0833], 
     [ 0.0833, -0.6667, 0. , 0.6667, -0.0833, 0. , 0. ], 
     [ 0. , 0.0833, -0.6667, 0. , 0.6667, -0.0833, 0. ], 
     [ 0. , 0. , 0.0833, -0.6667, 0. , 0.6667, -0.0833], 
     [-0.0833, 0. , 0. , 0.0833, -0.6667, 0. , 0.6667], 
     [ 0.6667, -0.0833, 0. , 0. , 0.0833, -0.6667, 0. ]]) 

Как вы отмечаете, для больших N, что требует много памяти и большинство значений равно нуль. Чтобы создать скудную разреженную матрицу, вы можете использовать scipy.sparse.diags. Мы должны создать смещения (и соответствующие значения) для диагоналей выше и ниже главной диагонали:

In [218]: from scipy import sparse 

In [219]: N = 7 

In [220]: vals = np.array([2.0/3, -1.0/12, 1.0/12, -2.0/3]) 

In [221]: offsets = np.array([1, 2, N-2, N-1]) 

In [222]: dupvals = np.concatenate((vals, vals[::-1])) 

In [223]: dupoffsets = np.concatenate((offsets, -offsets)) 

In [224]: a = sparse.diags(dupvals, dupoffsets, shape=(N, N)) 

In [225]: a.toarray() 
Out[225]: 
array([[ 0. , 0.6667, -0.0833, 0. , 0. , 0.0833, -0.6667], 
     [-0.6667, 0. , 0.6667, -0.0833, 0. , 0. , 0.0833], 
     [ 0.0833, -0.6667, 0. , 0.6667, -0.0833, 0. , 0. ], 
     [ 0. , 0.0833, -0.6667, 0. , 0.6667, -0.0833, 0. ], 
     [ 0. , 0. , 0.0833, -0.6667, 0. , 0.6667, -0.0833], 
     [-0.0833, 0. , 0. , 0.0833, -0.6667, 0. , 0.6667], 
     [ 0.6667, -0.0833, 0. , 0. , 0.0833, -0.6667, 0. ]]) 

матрица хранится в «диагонального» формат:

In [226]: a 
Out[226]: 
<7x7 sparse matrix of type '<class 'numpy.float64'>' 
    with 28 stored elements (8 diagonals) in DIAgonal format> 

Вы можете использовать преобразование методы разреженной матрицы, чтобы преобразовать ее в другой разреженный формат. Например, следующие результаты в матрице в формате CSR:

In [227]: a.tocsr() 
Out[227]: 
<7x7 sparse matrix of type '<class 'numpy.float64'>' 
    with 28 stored elements in Compressed Sparse Row format>