Crear una matriz diagonal dispersa a partir de la fila de una matriz dispersa

Question

Procedo matrices bastante grandes en Python/Scipy. Necesito extraer filas de la matriz grande (que está cargada en coo_matrix) y usarlas como elementos diagonales. Actualmente lo hago de la siguiente manera:

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)') 

Lo que veo desde la salida profile es que la mayoría de las veces es consumida por get_csr_submatrix función mientras que la extracción diag_elems. Eso me hace pensar que utilizo una representación dispersa ineficiente de datos iniciales o una forma incorrecta de extraer una fila de una matriz dispersa. ¿Puede sugerir una mejor manera de extraer una fila de una matriz dispersa y representarla en forma diagonal?

EDITAR

La siguiente variante elimina cuello de botella de la extracción fila (aviso de sencilla cambiando 'csc'-csr no es suficiente, A[i,:] debe ser reemplazado con A.getrow(i) también). Sin embargo, la pregunta principal es cómo omitir la materialización (.todense()) y crear la matriz diagonal a partir de la representación dispersa de la fila.

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)') 

Si creo matriz diagonal de 1-fila de la matriz RSE directamente, de la siguiente manera:

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

entonces puedo ni especificar format="csc" argumento, ni convertir ith_diags a formato 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` 

Answer 1

Esto es lo que se me ocurrió:

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() 

El generador de perfiles de Python dijo 1,464 segundos contra 5.574 segundos antes. Aprovecha las matrices densas subyacentes (indptr, índices, datos) que definen matrices dispersas. Aquí está mi curso intensivo: A.indptr [i]: A.indptr [i + 1] define qué elementos en las matrices densas corresponden a los valores distintos de cero en la fila i. A.data es una matriz 1d densa de valores distintos de cero de A y A.indptr es la columna a la que van esos valores.

Haría algunas pruebas más para asegurarme de que esto hace lo mismo que antes. Solo revisé algunos casos.