Как я могу запомнить разреженную матрицу с scipy.sparse.dia_matrix, давая ей массивы (диагонали) разных размеров?

Question

Я хочу запомнить пятиугольную квадратную матрицу с scipy.sparse.dia_matrix. Первая и последняя диагонали находятся далеко от главной диагонали, поэтому они содержат меньше элементов, чем главная диагональ. Кажется, dia_matrix принимает только диагонали с такой же длиной главной диагонали. Я хочу сохранить память и дать ей массив разной длины. Является ли это возможным? Ниже приведен пример (трехдиагональный):

diag1 = [88,99]
main_diag = [1,2,3,4]
diag2 = [ 101,202]



A = [ 1,   0, 88, 0
      0,   2,  0, 99
      101, 0,  3, 0
      0,  202, 0, 4 ] 
A = scipy.sparse.dia_matrix(([[diag1],[main_diag],[diag2]],offsets), shape= (4,4))

Но это выдаст мне ошибку, если я не запомню диагонали как:

diag1 = [0,0,88,99]
main_diag = [1,2,3,4]
diag2 = [ 101,202,0,0]

Но это не экономия памяти.

Answer 1

Ваша плотная матрица (почему вы не сделали эту совместимую вырезку-н-пасту?):

In [589]: A = np.array([[ 1,   0, 88, 0],
     ...:       [0,   2,  0, 99],
     ...:       [101, 0,  3, 0],
     ...:       [0,  202, 0, 4] ])
     ...:       
In [590]: A
Out[590]: 
array([[  1,   0,  88,   0],
       [  0,   2,   0,  99],
       [101,   0,   3,   0],
       [  0, 202,   0,   4]])

Мы можем сделать матрицу dia непосредственно из A:

In [591]: M = sparse.dia_matrix(A)
In [592]: M
Out[592]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 8 stored elements (3 diagonals) in DIAgonal format>

dia хранилища формата - это значения в двух массивах: 2d data и 1d offsets:

In [593]: M.data
Out[593]: 
array([[101, 202,   0,   0],
       [  1,   2,   3,   4],
       [  0,   0,  88,  99]])
In [594]: M.offsets
Out[594]: array([-2,  0,  2], dtype=int32)

Мы можем использовать функцию sparse.diags для построения матрицы dia из ваших 3 диагоналей:

In [605]: M1 = sparse.diags([diag1, main_diag, diag2],[2,0,-2], dtype=int)
In [606]: M1
Out[606]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 8 stored elements (3 diagonals) in DIAgonal format>
In [607]: M1.A
Out[607]: 
array([[  1,   0,  88,   0],
       [  0,   2,   0,  99],
       [101,   0,   3,   0],
       [  0, 202,   0,   4]])
In [608]: M1.data
Out[608]: 
array([[  0,   0,  88,  99],
       [  1,   2,   3,   4],
       [101, 202,   0,   0]])

Но основное хранилище то же самое.

Я не знаю деталей того, почему dia_matrix использует этот формат. Я подозреваю, что это связано с балансом затрат на хранение и обработку. Чтобы хранить data как двумерный массив, он должен иметь своего рода заполнение. В противном случае вы должны использовать список списков (или массив объектов из массивов / списков), как в [diag1, main_diag, diag2]. Использует ли это больше или меньше памяти, может зависеть от размера матрицы и того, как далеко находятся offsets.

В формате lil та же матрица хранится как 2 массива списков:

In [611]: M.tolil().data
Out[611]: 
array([list([1, 88]), list([2, 99]), list([101, 3]), list([202, 4])],
      dtype=object)
In [612]: M.tolil().rows
Out[612]: 
array([list([0, 2]), list([1, 3]), list([0, 2]), list([1, 3])],
      dtype=object)

и coo:

In [615]: M.tocoo().data
Out[615]: array([101, 202,   1,   2,   3,   4,  88,  99])
In [616]: M.tocoo().row
Out[616]: array([2, 3, 0, 1, 2, 3, 0, 1], dtype=int32)
In [617]: M.tocoo().col
Out[617]: array([0, 1, 0, 1, 2, 3, 2, 3], dtype=int32)

csr, широко используемый для расчетов, сжимает массив row.

Если вы выполняете все вычисления в плотной форме, вы можете построить A из чего-то вроде:

np.diag(main_diag,0)+np.diag(diag1,2)+np.diag(diag2,-2)

Сравнивая это с sparse опосредованным подходом:

In [623]: timeit np.diag(main_diag,0)+np.diag(diag1,2)+np.diag(diag2,-2)
27.5 µs ± 1.28 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [624]: timeit sparse.diags([diag1, main_diag, diag2],[2,0,-2], dtype=int).A
288 µs ± 13 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)