疎行列の行から疎な対角行列を作成する

Question

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出力から見る何diag_elemsを抽出しながら、ほとんどの時間をget_csr_submatrix機能によって消費されていることです。これは、私が最初のデータの非効率的な疎な表現か、まばらな行列からの行を抽出する間違った方法のどちらかを使用すると考える。疎な行列から行を抽出し、それを斜めの形で表現するより良い方法を提案できますか？

EDIT

は、以下の変異体は、行抽出（csrに'csc'を変更する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)') 

次のように私は、直接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` 

Answer 1

ここに私が思いつきました：

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プロファイラは、5.574秒前と比べて1.464秒と言っています。これは、疎行列を定義する基底の密な配列（indptr、indices、data）を利用します。 A.indptr [i]：A.indptr [i + 1]は、高密度配列のどの要素が行iの非ゼロ値に対応するかを定義します。 A.dataは、AとAの値が0でない密な1次元配列です.indptrは、それらの値が入る列です。

私は以前と同じことをしていることを確かめるために、さらにテストをします。私はいくつかのケースだけをチェックしました。