3

現在、Juliaで書かれたSparse行列Aの計算に関するコードをPythonに書き直す作業をしています。
また、そのJuliaコードはVer.0.5.2で書かれているものであり、だいぶ古いバージョンです。
そのJuliaで書かれた方のコードにどうしてもわからない部分があり、それをPythonに書き換える事ができずにいます。
以下がそのjuliaコードです。

Au = A.rowval[A.colptr[u]:A.colptr[u+1]-1]

このコードの自分なりの理解は、Aを適当なスカラーuの列ポインタcolptrで指定された部分内(A.colptr[u]からA.colptr[u+1]-1)に値があった場合の行インデックスをrowvalで取得する、というものだと考えています。(Juliaのsparse行列に関するドキュメンテーションを参考にしました)
以下にu=2とした場合の例があります。
その例を基に私の考え方と合わせると、AuA.rowval[15:23]に値がある行のインデックスを返すため、最終的な(重複無しの)returnは[15,16,17,18,19,20,21,22,23,33,34]になるはずですが、実際の結果は[2,3,4,8,14,18,20,22,31]でした。
どこで考え方が間違っているのかわからないため、Pythonコードに変換できません。
どなたか、このrowvalに関する解説と、同じ結果になるようなpythonコードのご教授をよろしくお願いします。

仮にこの問題に取り組んでくれる方がいた場合のために、このコードに該当する部分と実際に使用しているサイズ(34,34)行列の.mtxファイルの中身を記しておきます。

dir  = "ファイルへのディレクトリ"
ln = "行列.mtx"
xxx = MatrixMarketRead(string(dir,strip(ln)));
A = xxx - spdiagm(diag(xxx))
n = size(A,1);
A = speye(n) - A * spdiagm(1./vec(sum(A,1)));
Au = A.rowval[A.colptr[2]:A.colptr[2+1]-1]  # A.colptr[2]=15, A.colptr[2+1]-1 = 23
println(Au)  # [2,3,4,8,14,18,20,22,31]

.mtxファイルの中身

2 1 1
3 1 1
4 1 1
5 1 1
6 1 1
7 1 1
8 1 1
9 1 1
11 1 1
12 1 1
13 1 1
14 1 1
18 1 1
20 1 1
22 1 1
32 1 1
3 2 1
4 2 1
8 2 1
14 2 1
18 2 1
20 2 1
22 2 1
31 2 1
4 3 1
8 3 1
9 3 1
10 3 1
14 3 1
28 3 1
29 3 1
33 3 1
8 4 1
13 4 1
14 4 1
7 5 1
11 5 1
7 6 1
11 6 1
17 6 1
17 7 1
31 9 1
33 9 1
34 9 1
34 10 1
34 14 1
33 15 1
34 15 1
33 16 1
34 16 1
33 19 1
34 19 1
34 20 1
33 21 1
34 21 1
33 23 1
34 23 1
26 24 1
28 24 1
30 24 1
33 24 1
34 24 1
26 25 1
28 25 1
32 25 1
32 26 1
30 27 1
34 27 1
34 28 1
32 29 1
34 29 1
33 30 1
34 30 1
33 31 1
34 31 1
33 32 1
34 32 1

JuliaのrowvalはPythonでindicesに値するみたいなので、シンプルな行列を用いて上の例を基に確認してみたところ、違う結果となりました。以下がその例です。
Juliaコード

a = spzeros(3, 3)
a[1, 1] = 1
a[2, 1] = 4
a[3, 1] = 7
a[1, 2] = 2
a[2, 2] = 5
a[3, 2] = 8
a[1, 3] = 3
a[2, 3] = 6
a[3, 3] = 9
println(a.rowval[a.colptr[1]:a.colptr[3]])  # [1,2,3,1,2,3,1]

Pythonコード

import numpy as np
from scipy.sparse import csc_matrix

x = np.matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
xx = csc_matrix(x)
xxx = xx.indices[xx.indptr[0]:xx.indptr[2]]
print(xxx)  # [0 1 2 0 1 2]
7
  • ご自身の英語版記事 What is “colptr” in Julia and its counterpart in Python? の回答では適応しないのでしょうか? rowvalに関しては逆のようですが。(回答ではindices -> rowvalですが前後に合わせるならrowval -> indicesでしょう) scipy.sparse.csc_matrix
    – kunif
    Commented 2021年1月2日 3:16
  • @kunif その記事で得られた回答を基にJuliaとPythonでそれぞれ同じ行列で確認してみたのですが、結果が異なっていたため、あまり理解できていないrowvalに関して質問してみました。その例を質問に載せておきます
    – S.Ky
    Commented 2021年1月2日 5:02
  • データの数(juliaが7個、pythonが6個)が合ってないのは理由が分かりませんが、数値についてはインデックス値のオリジンがjuliaが1でpythonが0だからでは? pythonのxxをprintして見較べてみてはどうでしょう?
    – kunif
    Commented 2021年1月2日 14:43
  • @kunif アドバイスありがとうございます。その違いを失念していました。しかし、pythonコードの方をxx.indptr[0]:xx.indptr[2]に変更しても結果は変わりませんでした。
    – S.Ky
    Commented 2021年1月3日 5:29
  • 1
    それから、表示されるデータの数が違うのも言語仕様の差であると思われます。Juliaのslice等の範囲指定が、start:endの時にendに指定した位置の値も含むようです。それに対してPythonではendに指定した位置の値は含まれません。質問の例だけに適用されるのかもしれませんが、xxx = xx.indices[xx.indptr[0]:xx.indptr[2]]の部分はxxx = xx.indices[xx.indptr[0]:xx.indptr[2]+1]として、:end部分をxx.indptr[2]+1にすると、表示されるデータ数はJuliaとPythonで同じになるようです。
    – kunif
    Commented 2021年1月3日 15:01

1 件の回答 1

1
+100

JuliaのSparseArrays - rowvalと同等で、Pythonの相当する物は、scipy.sparse.csc_matrix - indicesでしょう。

その上で、JuliaとPythonのインデックス値に関する言語仕様の違いにより、それぞれの言語に従った普通の処理のままでは取得出来るのは同じ数値(のリスト)になり得ないと考えられます。

Julia
インデックス値 = 1から始まる
スライス範囲指定= 開始位置:終了位置
Python
インデックス値 = 0から始まる
スライス範囲指定= 開始位置:終了位置+1


質問の最初にあるJuliaの以下の行は:

Au = A.rowval[A.colptr[u]:A.colptr[u+1]-1]

uの値の考え方がJuliaのまま1オリジンであるならば、Pythonでは以下のようになるでしょう:

Au = A.indices[A.indptr[u-1]:A.indptr[u]]

uの値もPythonの考え方を適用して0オリジンにするなら、以下のようになるでしょう:

Au = A.indices[A.indptr[u]:A.indptr[u+1]]

上記処理により、取得する値のリスト上の範囲は同じになりますが、数値は1オリジン/0オリジンの違いにより、Pythonでは全てのデータがそれぞれ1づつ小さい値になるでしょう。


シンプルな行列の例では、範囲指定は以下のようになります。(Juliaに合わせる)
Julia:

println(a.rowval[a.colptr[1]:a.colptr[3]])  # [1,2,3,1,2,3,1]

Python:

xxx = xx.indices[xx.indptr[0]:xx.indptr[2]+1]
print(xxx)  # [0 1 2 0 1 2 0]

あるいはこちら(Pythonに合わせる)
Julia:

println(a.rowval[a.colptr[1]:a.colptr[3]-1])  # [1,2,3,1,2,3]

Python:

xxx = xx.indices[xx.indptr[0]:xx.indptr[2]]
print(xxx)  # [0 1 2 0 1 2]

JuliaとPythonのシンプルな配列を表示すると以下になります。
実際の表示そのままではなく、対比しやすいように水平に組み合わせています。

julia> a                                           |                      
3×3 sparse matrix with 9 Float64 nonzero entries:  | >>> print(xx)        
        [1, 1]  =  1.0                             |   (0, 0)        1    
        [2, 1]  =  4.0                             |   (1, 0)        4    
        [3, 1]  =  7.0                             |   (2, 0)        7    
        [1, 2]  =  2.0                             |   (0, 1)        2    
        [2, 2]  =  5.0                             |   (1, 1)        5    
        [3, 2]  =  8.0                             |   (2, 1)        8    
        [1, 3]  =  3.0                             |   (0, 2)        3    
        [2, 3]  =  6.0                             |   (1, 2)        6    
        [3, 3]  =  9.0                             |   (2, 2)        9    

上記の中で、どの部分が通知されているかというと、以下の図の赤い枠がJuliaに合わせた場合、緑の枠がPythonに合わせた場合のデータです。
青の枠は、Pythonの当初のxxx = xx.indices[xx.indptr[1]:xx.indptr[3]]で取得していた時のデータです。
JuliaPythonSparse配列


根拠情報:

SparseArrays - rowvalscipy.sparse.csc_matrix - indicesについて

説明としては全く同じ文言という訳ではないですが、以下それぞれの仕様記述と関連を見れば同等なものと考えられます。

Julia
SparseArrays

In Julia, sparse matrices are stored in the Compressed Sparse Column (CSC) format. Julia sparse matrices have the type SparseMatrixCSC{Tv,Ti}, where Tv is the type of the stored values, and Ti is the integer type for storing column pointers and row indices. The internal representation of SparseMatrixCSC is as follows:

struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int                  # Number of rows
    n::Int                  # Number of columns
    colptr::Vector{Ti}      # Column j is in colptr[j]:(colptr[j+1]-1)
    rowval::Vector{Ti}      # Row indices of stored values
    nzval::Vector{Tv}       # Stored values, typically nonzeros
end

上記のリンク先(Sparse matrix - Wikipedia)の該当記事がこちら
Compressed sparse column (CSC or CCS)

For example, CSC is (val, row_ind, col_ptr), where val is an array of the (top-to-bottom, then left-to-right) non-zero values of the matrix; row_ind is the row indices corresponding to the values; and, col_ptr is the list of val indexes where each column starts.

上記の段落の最後の方にSee scipy.sparse.csc_matrixとリンクされている先がこちら
Python
scipy.sparse.csc_matrix

Attributes
    shape : 2-tuple
        Get shape of a matrix.

    data
        Data array of the matrix

    indices
        CSC format index array

    indptr
        CSC format index pointer array

つまり英語版記事の回答は(indices -> rowvalは前後が逆ですが)合っていると考えられます。
What is “colptr” in Julia and its counterpart in Python?

The typical Python equivalent is in scipy: scipy.sparse.csc_matrix; simply substitutue colptr -> indptr, indices -> rowval, nzval -> data and accommodate the 0-based indexing.

補助的にはこれらの記事も
Compressed Sparse Column Format (CSC) - scipy lecture notes

  • column oriented
    • three NumPy arrays: indices, indptr, data
      • indices is array of row indices
      • data is array of corresponding nonzero values
      • indptr points to column starts in indices and data
      • length is n_col + 1, last item = number of values = length of both indices and data

Python, SciPy(scipy.sparse)で疎行列を生成・変換

CSC: scipy.sparse.csc_matrix

属性data, indices, indptrにデータが格納されている。

data, indicesは値、行のインデックスのリスト。indptrは列のインデックスを圧縮したリストとなる。

indptrのサイズ(要素数)は列数 + 1となる。


JuliaとPythonのインデックス値に関する言語仕様の違い

Julia
Arrays with custom indices

Conventionally, Julia's arrays are indexed starting at 1, whereas some other languages start numbering at 0, and yet others (e.g., Fortran) allow you to specify arbitrary starting indices. While there is much merit in picking a standard (i.e., 1 for Julia), there are some algorithms which simplify considerably if you can index outside the range 1:size(A,d) (and not just 0:size(A,d)-1, either). To facilitate such computations, Julia supports arrays with arbitrary indices.

Juliaで最低限やっていくための配列操作まとめ

1次元配列
一つの要素を取り出す
基本はa[index]で指定する(Linear indexing)。indexは0ではなく1から開始することに注意。

配列の要素にアクセスする
複数の要素を取り出す
スライスが使える。

julia> a = Vector(1:5)
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> a[1:end]
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

例えば続けてa[2:4]とすると以下になる

julia> a[2:4]
3-element Array{Int64,1}:
 2
 3
 4

インデックス値 = 1から始まる
スライス範囲指定= 開始位置:終了位置

Python
シーケンス型 --- list, tuple, range

共通のシーケンス演算

s[i]  s の 0 から数えて i 番目の要素 (3)

s[i:j] s の i から j までのスライス  (3)(4)

  1. s の i から j へのスライスは i <= k < j となるようなインデックス k を持つ要素からなるシーケンスとして定義されます。

インデックス値 = 0から始まる
スライス範囲指定= 開始位置:終了位置+1

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。