Создать разреженную диагональную матрицу из ряда разреженной матрицы - PullRequest
Новая
       3

Создать разреженную диагональную матрицу из ряда разреженной матрицы

5 голосов
/ 01 декабря 2011

Я обрабатываю довольно большие матрицы в Python / Scipy. Мне нужно извлечь строки из большой матрицы (которая загружается в coo_matrix) и использовать их в качестве диагональных элементов. В настоящее время я делаю это следующим образом: 

import numpy as np
from scipy import sparse

def computation(A):
  for i in range(A.shape[0]):
    diag_elems = np.array(A[i,:].todense())
    ith_diag = sparse.spdiags(diag_elems,0,A.shape[1],A.shape[1], format = "csc")
    #...

#create some random matrix
A = (sparse.rand(1000,100000,0.02,format="csc")*5).astype(np.ubyte)
#get timings
profile.run('computation(A)')

Из вывода profile я вижу, что большую часть времени занимает функция get_csr_submatrix при извлечении diag_elems. Это заставляет меня думать, что я использую либо неэффективное разреженное представление исходных данных, либо неправильный способ извлечения строки из разреженной матрицы. Можете ли вы предложить лучший способ извлечь строку из разреженной матрицы и представить ее в диагональной форме?

EDIT

Следующий вариант удаляет узкое место из извлечения строки (обратите внимание, что простого изменения 'csc' на csr недостаточно, A[i,:] также следует заменить на A.getrow(i)). Однако главный вопрос заключается в том, как исключить материализацию (.todense()) и создать диагональную матрицу из разреженного представления строки. 

import numpy as np
from scipy import sparse

def computation(A):
  for i in range(A.shape[0]):
    diag_elems = np.array(A.getrow(i).todense())
    ith_diag = sparse.spdiags(diag_elems,0,A.shape[1],A.shape[1], format = "csc")
    #...

#create some random matrix
A = (sparse.rand(1000,100000,0.02,format="csr")*5).astype(np.ubyte)
#get timings
profile.run('computation(A)')

Если я создаю матрицу DIAgonal непосредственно из 1-рядной матрицы CSR, как показано ниже: 

diag_elems = A.getrow(i)
ith_diag = sparse.spdiags(diag_elems,0,A.shape[1],A.shape[1])

тогда я не могу ни указать format="csc" аргумент, ни преобразовать ith_diags в формат CSC: 

Traceback (most recent call last):
   File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.6/profile.py", line 70, in run
    prof = prof.run(statement)
  File "/usr/local/lib/python2.6/profile.py", line 456, in run
    return self.runctx(cmd, dict, dict)
  File "/usr/local/lib/python2.6/profile.py", line 462, in runctx
    exec cmd in globals, locals
  File "<string>", line 1, in <module>
  File "<stdin>", line 4, in computation
  File "/usr/local/lib/python2.6/site-packages/scipy/sparse/construct.py", line 56, in spdiags
    return dia_matrix((data, diags), shape=(m,n)).asformat(format)
  File "/usr/local/lib/python2.6/site-packages/scipy/sparse/base.py", line 211, in asformat
    return getattr(self,'to' + format)()
  File "/usr/local/lib/python2.6/site-packages/scipy/sparse/dia.py", line 173, in tocsc
    return self.tocoo().tocsc()
  File "/usr/local/lib/python2.6/site-packages/scipy/sparse/coo.py", line 263, in tocsc
    data    = np.empty(self.nnz, dtype=upcast(self.dtype))
  File "/usr/local/lib/python2.6/site-packages/scipy/sparse/sputils.py", line 47, in upcast
    raise TypeError,'no supported conversion for types: %s' % args
TypeError: no supported conversion for types: object`

1 Ответ

3 голосов
/ 04 декабря 2011

Вот что я придумал: 

def computation(A):
    for i in range(A.shape[0]):
        idx_begin = A.indptr[i]
        idx_end = A.indptr[i+1]
        row_nnz = idx_end - idx_begin
        diag_elems = A.data[idx_begin:idx_end]
        diag_indices = A.indices[idx_begin:idx_end]
        ith_diag = sparse.csc_matrix((diag_elems, (diag_indices, diag_indices)),shape=(A.shape[1], A.shape[1]))
        ith_diag.eliminate_zeros()

Профилировщик Python сказал 1,464 секунды против 5,574 секунды раньше.Он использует преимущества базовых плотных массивов (indptr, indexes, data), которые определяют разреженные матрицы.Вот мой ускоренный курс: A.indptr [i]: A.indptr [i + 1] определяет, какие элементы в плотных массивах соответствуют ненулевым значениям в строке i.A.data - это плотный 1d-массив ненулевых значений A, а A.indptr - столбец, в котором эти значения идут.до.Я проверил только несколько случаев.

...