73076


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

[code]

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

このコードの自分なりの理解は、`A`を適当なスカラー`u`の列ポインタ`colptr`で指定された部分内(`A.colptr[u]`から`A.colptr[u+1]-1`)に値があった場合の行インデックスを`rowval`で取得する、というものだと考えています。(Juliaのsparse行列に関する[ドキュメンテーション](https://docs.julialang.org/en/v1/stdlib/SparseArrays/)を参考にしました)  
以下に`u=2`とした場合の例があります。  
その例を基に私の考え方と合わせると、`Au`は`A.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`ファイルの中身を記しておきます。

[code]

    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]
    
[/code]

`.mtx`ファイルの中身

[code]

    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
    
[/code]

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

[code]

    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]
    
[/code]

Pythonコード

[code]

    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]
    
[/code]


CC BY-SA 4.0
2021-01-01T07:21:54.207
73007
2021-01-05T02:51:41.837
2021-01-03T05:29:35.803
36750
36750
question
python python3 julia 行列
Juliaの行列オプションrowvalをPythonで表現するには？


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

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

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

* * *

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

[code]

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

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

[code]

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

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

[code]

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

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

* * *

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

[code]

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

Python:

[code]

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

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

[code]

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

Python:

[code]

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

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

[code]

    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    
    
[/code]

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

* * *

根拠情報：

`SparseArrays - rowval`と`scipy.sparse.csc_matrix - indices`について

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

**Julia**  
[SparseArrays](https://docs.julialang.org/en/v1/stdlib/SparseArrays/)

> In Julia, sparse matrices are stored in the [Compressed Sparse Column (CSC)
> format](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29).
> 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:
[code]

>     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
>  
[/code]

上記のリンク先(Sparse matrix - Wikipedia)の該当記事がこちら  
[Compressed sparse column (CSC or
CCS)](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29)

> 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](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html)

>
[code]

>     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
>  
[/code]

つまり英語版記事の回答は(`indices` -> `rowval`は前後が逆ですが)合っていると考えられます。  
[What is “colptr” in Julia and its counterpart in
Python?](https://stackoverflow.com/q/64971750/9014308)

> The typical Python equivalent is in scipy:
> [scipy.sparse.csc_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html);
> simply substitutue `colptr` -> `indptr`, `indices` -> `rowval`, `nzval` ->
> `data` and accommodate the 0-based indexing.

補助的にはこれらの記事も  
[Compressed Sparse Column Format (CSC) - scipy lecture notes](https://scipy-
lectures.org/advanced/scipy_sparse/csc_matrix.html)

>   * 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）で疎行列を生成・変換](https://note.nkmk.me/python-scipy-
sparse-matrix-csr-csc-coo-lil/)

> **CSC: scipy.sparse.csc_matrix**
>
> 属性`data`, `indices`, `indptr`にデータが格納されている。
>
> `data`, `indices`は値、行のインデックスのリスト。`indptr`は列のインデックスを圧縮したリストとなる。
>
> `indptr`のサイズ（要素数）は列数 + 1となる。

* * *

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

**Julia**  
[Arrays with custom indices](https://docs.julialang.org/en/v1/devdocs/offset-
arrays/)

> 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で最低限やっていくための配列操作まとめ](https://qiita.com/A03ki/items/007be353411d19952ef7)

> 1次元配列  
>  一つの要素を取り出す  
>  基本はa[index]で指定する([Linear
> indexing](https://docs.julialang.org/en/v1/manual/arrays/index.html#Linear-
> indexing-1))。indexは0ではなく1から開始することに注意。
>
> 配列の要素にアクセスする  
>  複数の要素を取り出す  
>  スライスが使える。
[code]

>     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
>  
[/code]

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

[code]

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

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

**Python**  
[シーケンス型 --- list, tuple,
range](https://docs.python.org/ja/3/library/stdtypes.html#sequence-types-list-
tuple-range)

> 共通のシーケンス演算
>
> `s[i]` s の **0 から数えて** i 番目の要素 (3)
>
> `s[i:j]` s の i から j までのスライス (3)(4)
>
>   4. s の i から j へのスライスは **`i<= k < j`** となるようなインデックス k
> を持つ要素からなるシーケンスとして定義されます。
>

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


CC BY-SA 4.0
2021-01-04T17:46:42.620
73076
2021-01-05T02:51:41.837
2021-01-05T02:51:41.837
26370
26370
73007
answer