Ошибка в SciPy разреженной конструкции Diags матрицы

Question

При использовании scipy.sparse.spdiags или scipy.sparse.diags я заметил, хочу я считаю, что это ошибка в подпрограммы, например

scipy.sparse.spdiags([1.1,1.2,1.3],1,4,4).toarray() 

возвращается

array([[ 0. , 1.2, 0. , 0. ], 
     [ 0. , 0. , 1.3, 0. ], 
     [ 0. , 0. , 0. , 0. ], 
     [ 0. , 0. , 0. , 0. ]]) 

То есть для положительных диагоналей он отбрасывает первые k данных. Можно было бы утверждать, что для этого есть какая-то грандиозная причина программирования, и мне просто нужно заполнить нулями. ОК раздражает, как это может быть, можно использовать scipy.sparse.diags, который дает правильный результат. Однако эта процедура имеет ошибку, которая не может быть работал вокруг

scipy.sparse.diags([1.1,1.2],0,(4,2)).toarray() 

дает

array([[ 1.1, 0. ], 
     [ 0. , 1.2], 
     [ 0. , 0. ], 
     [ 0. , 0. ]]) 

хорошо, и

scipy.sparse.diags([1.1,1.2],-2,(4,2)).toarray() 

дает

array([[ 0. , 0. ], 
     [ 0. , 0. ], 
     [ 1.1, 0. ], 
     [ 0. , 1.2]]) 

но

scipy.sparse.diags([1.1,1.2],-1,(4,2)).toarray() 

вызывает ошибку ValueError: Диагональная длина (индекс 0: 2 при смещении -1) не согласуется с размером матрицы (4, 2). Очевидно, что ответ

array([[ 0. , 0. ], 
     [ 1.1, 0. ], 
     [ 0. , 1.2], 
     [ 0. , 0. ]]) 

и для дополнительного случайного поведения мы имеем

scipy.sparse.diags([1.1],-1,(4,2)).toarray() 

давая

array([[ 0. , 0. ], 
     [ 1.1, 0. ], 
     [ 0. , 1.1], 
     [ 0. , 0. ]]) 

Каждые знают, если есть функция для построения диагональных разреженных матриц, что на самом деле работает?

Answer 1

Резюме: spdiags работает правильно, даже если вход матрицы не является наиболее интуитивным. diags имеет ошибку, которая влияет на некоторые смещения в прямоугольных матрицах. На scipy github исправлена ​​ошибка.

---

Пример для spdiags является:

>>> data = array([[1,2,3,4],[1,2,3,4],[1,2,3,4]]) 
>>> diags = array([0,-1,2]) 
>>> spdiags(data, diags, 4, 4).todense() 
matrix([[1, 0, 3, 0], 
     [1, 2, 0, 4], 
     [0, 2, 3, 0], 
     [0, 0, 3, 4]]) 

Обратите внимание, что третий столбец data всегда появляется в 3-м столбце разреженный. Остальные столбцы также выстраиваются в линию. Но они опускаются там, где они «падают с края».

Ввод этой функции - это матрица, а входной сигнал diags - это оборванный список. Диагонали разреженной матрицы имеют различное количество значений. Поэтому спецификация должна соответствовать этому в той или иной форме. spdiags делает это, игнорируя некоторые значения, diags, беря ввод списка.

Ошибка sparse.diags([1.1,1.2],-1,(4,2)) является недоумением.

spdiags эквивалент делает работу:

In [421]: sparse.spdiags([[1.1,1.2]],-1,4,2).A 
Out[421]: 
array([[ 0. , 0. ], 
     [ 1.1, 0. ], 
     [ 0. , 1.2], 
     [ 0. , 0. ]]) 

Ошибка возникает в этом блоке кода:

for j, diagonal in enumerate(diagonals): 
    offset = offsets[j] 
    k = max(0, offset) 
    length = min(m + offset, n - offset) 
    if length <= 0: 
     raise ValueError("Offset %d (index %d) out of bounds" % (offset, j)) 
    try: 
     data_arr[j, k:k+length] = diagonal 
    except ValueError: 
     if len(diagonal) != length and len(diagonal) != 1: 
      raise ValueError(
       "Diagonal length (index %d: %d at offset %d) does not " 
       "agree with matrix size (%d, %d)." % (
       j, len(diagonal), offset, m, n)) 
     raise 

Фактическая матрица конструктор в diags является:

dia_matrix((data_arr, offsets), shape=(m, n)) 

Это тот же самый конструктор, который использует spdiags, bu t без каких-либо манипуляций.

In [434]: sparse.dia_matrix(([[1.1,1.2]],-1),shape=(4,2)).A 
Out[434]: 
array([[ 0. , 0. ], 
     [ 1.1, 0. ], 
     [ 0. , 1.2], 
     [ 0. , 0. ]]) 

В формате dia входы хранятся точно так, как дано spdiags (в комплекте с этой матрицей с дополнительными значениями):

In [436]: M.data 
Out[436]: array([[ 1.1, 1.2]]) 
In [437]: M.offsets 
Out[437]: array([-1], dtype=int32) 

---

Как @user2357112 указывает, length = min(m + offset, n - offset неправильна, производство 3 в тестовом примере. Изменение его на length = min(m + k, n - k) делает все случаи для этой (4,2) матрицы работы. Но он терпит неудачу с транспонированной: diags([1.1,1.2], 1, (2, 4))

Коррекция, по состоянию на 5 октября, на этот вопрос:

https://github.com/pv/scipy-work/commit/529cbde47121c8ed87f74fa6445c05d71353eb6c

length = min(m + offset, n - offset, min(m,n)) 

С этим исправлением, diags([1.1,1.2], 1, (2, 4)) работ.