Неожиданная ошибка от sparse.spdiags()

В Python 3 я пытаюсь запустить следующую строку кода, чтобы получить конкретную разреженную матрицу.

sparse.spdiags(np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1), [0, 1], 9, 10)

Это дает следующее сообщение об ошибке:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/scipy/sparse/construct.py", line 61, in spdiags
    return dia_matrix((data, diags), shape=(m,n)).asformat(format)
  File "/usr/lib/python3/dist-packages/scipy/sparse/dia.py", line 138, in __init__
    % (self.data.shape[0], len(self.offsets)))
ValueError: number of diagonals (9) does not match the number of offsets (2)

Выполнение того, что я понимаю как эквивалентный код в Octave, похоже, дает мне разреженную матрицу.

spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)

Compressed Column Sparse (rows = 9, cols = 10, nnz = 18 [20%])

(1, 1) -> -1
(1, 2) ->  1
(2, 2) -> -1
(2, 3) ->  1
(3, 3) -> -1
(3, 4) ->  1
(4, 4) -> -1
(4, 5) ->  1
(5, 5) -> -1
(5, 6) ->  1
(6, 6) -> -1
(6, 7) ->  1
(7, 7) -> -1
(7, 8) ->  1
(8, 8) -> -1
(8, 9) ->  1
(9, 9) -> -1
(9, 10) ->  1

Любые идеи о том, почему они ведут себя по-разному, и как это исправить?

ДОПОЛНЕНИЕ

У меня есть дополнительная проблема с выводом Scipy.sparse по сравнению с Octave.

ПИТОН

>>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A.shape                                                      
(9, 10)

>>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A                         
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A.shape
(9, 10) 

>>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A.shape
(9, 10)

ОКТАВА

>full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10))
ans =

  -1   1   0   0   0   0   0   0   0   0
   0  -1   1   0   0   0   0   0   0   0
   0   0  -1   1   0   0   0   0   0   0
   0   0   0  -1   1   0   0   0   0   0
   0   0   0   0  -1   1   0   0   0   0
   0   0   0   0   0  -1   1   0   0   0
   0   0   0   0   0   0  -1   1   0   0
   0   0   0   0   0   0   0  -1   1   0
   0   0   0   0   0   0   0   0  -1   1

>size(full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)))
ans =

9   10

Почему scipy и Octave не дают одинакового значения в последнем столбце последней строки?


python scipy sparse-matrix octave diagonal
person Galen    schedule 09.08.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Ваш concatenate создает матрицу (9,2):

In [310]: np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1)
Out[310]: 
array([[-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.]])

In [311]: _.shape
Out[311]: (9, 2)

spdiags doc описывает этот data параметр как matrix diagonals stored row-wise. То есть каждой строке матрицы соответствует диагональ. 9 строк, но только 2 значения в [0,1].

Это важное отличие, о котором я упоминал в своем предыдущем ответе, хотя, возможно, я недостаточно подчеркнул это.

Если вам нужны 2 диагонали, вам нужно указать массив (2,9), например транспонирование этой матрицы:

In [317]: sparse.spdiags(np.concatenate((-np.ones((9,1)),
    np.ones((9,1))), axis=1).T, [0,1],9,10)
Out[317]: 
<9x10 sparse matrix of type '<class 'numpy.float64'>'
    with 18 stored elements (2 diagonals) in DIAgonal format>

Вы также можете построить диагонали с помощью:

In [321]: np.concatenate([-np.ones((1,9)), np.ones((1,9))],axis=0)
Out[321]: 
array([[-1., -1., -1., -1., -1., -1., -1., -1., -1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]])

Or np.vstack([-np.ones(9),np.ones(9)]) or np.ones(9)*[[-1],[1]].

Посмотрите на мой предыдущий ответ, но измените окончательную форму, чтобы столбцов было больше, чем строк):

octave:17> reshape (1:12, 4, 3)
ans =
    1    5    9
    2    6   10
    3    7   11
    4    8   12

octave:18> full(spdiags (reshape (1:12, 4, 3), [-1 0 1], 4,5))
ans =
    5    9    0    0    0
    1    6   10    0    0
    0    2    7   11    0
    0    0    3    8   12

In [327]: np.arange(1,13).reshape(3,4)
Out[327]: 
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])

In [328]: sparse.spdiags(np.arange(1,13).reshape(3,4), [-1, 0, 1], 4,5).A
Out[328]: 
array([[ 5, 10,  0,  0,  0],
       [ 1,  6, 11,  0,  0],
       [ 0,  2,  7, 12,  0],
       [ 0,  0,  3,  8,  0]])

В Octave (и предположительно MATLAB) диагональ +1 начинается с 9, заканчивается 12, полным последним столбцом входной матрицы. Посмотрите на [2,6,10] - в расположении под прямым углом.

Диагональ scipy начинается с 10, заканчивается добавленным 0. 9 невидим в несуществующей строке выше. Посмотрите на [2,6,10] - в один столбец.

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

Другая функция scipy устраняет неоднозначность, ожидая правильного количества элементов для каждой диагонали (в виде списка списков):

In [337]: sparse.diags([[1,2,3],[5,6,7,8],[9,10,11,12]],[-1,0,1],(4,5),dtype=int).A
Out[337]: 
array([[ 5,  9,  0,  0,  0],
       [ 1,  6, 10,  0,  0],
       [ 0,  2,  7, 11,  0],
       [ 0,  0,  3,  8, 12]])

Обратите внимание, что мне пришлось опустить 4.

Посмотрите на метод tocoo в scipy/sparse/dia.py, чтобы узнать больше о том, как dia_matrix сопоставляет данные диагоналей с разреженными координатами (формат coo).

person hpaulj    schedule 09.08.2015