PyCall.jl with SparseArrays

General Usage
,
1

I am trying to diagonalize a huge sparse Hermitian matrix (original thread here: Suggestions needed: diagonalizing large Hermitian sparse matrix) and am exploring the possibility of utilizing the PRIMME package. Since the Primme.jl (GitHub - andreasnoack/Primme.jl: Julia wrapper for Primme: Preconditioned Iterative Multimethod Eigensolver) is outdated and does not contain eigenvalue solver interface, I am trying to call PRIMME by PyCall.

It seems PyCall does not know how to translate SparseArrays to scipy sparse arrays:

julia> using PyCall, SparseArrays

julia> sp = pyimport("scipy")
PyObject <module 'scipy' from '/path/to/myhome/.julia/conda/3/lib/python3.9/site-packages/scipy/__init__.py'>

julia> A = sprandn(1000, 1000, 0.0025);

julia> B = sprandn(1000, 1000, 0.0025);

julia> C = A + im * B;

julia> D = C' + C;

julia> sp.sparse.linalg.eigs(D)
ERROR: PyError ($(Expr(:escape, :(ccall(#= /path/to/myhome/.julia/packages/PyCall/twYvK/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'AttributeError'>
AttributeError("'list' object has no attribute 'shape'")
  File "/path/to/myhome/.julia/conda/3/lib/python3.9/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py", line 1256, in eigs
    if A.shape[0] != A.shape[1]:

Stacktrace:
  [1] pyerr_check
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/exception.jl:96
  [4] macro expansion
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:473 [inlined]
  [7] __pycall!
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyObject, o::PyObject, args::Tuple{SparseMatrixCSC{ComplexF64, Int64}}, nargs::Int64, kw::Ptr{Nothing})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:29
  [9] _pycall!
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:11 [inlined]
 [10] #_#114
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:86 [inlined]
 [11] (::PyObject)(args::SparseMatrixCSC{ComplexF64, Int64})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:86
 [12] top-level scope
    @ REPL[34]:1

Any ideas about how I can utilize sparse eigenvalue solvers from the Python ecosystem?

2

one of these step converted type, check the type of im*B and A + im*B and C'

3

The types are exactly what I expected. What do you mean?

julia> typeof(A)
SparseMatrixCSC{Float64, Int64}

julia> typeof(B)
SparseMatrixCSC{Float64, Int64}

julia> typeof(C)
SparseMatrixCSC{ComplexF64, Int64}

julia> typeof(D)
SparseMatrixCSC{ComplexF64, Int64}
4

See Convert to sparse · Issue #204 · JuliaPy/PyCall.jl · GitHub for conversion in the opposite direction. You should just be able to call the scipy.sparse.csc_array constructor directly with the corresponding data from a Julia sparse arrray, e.g. something like:

pysparse = pyimport("scipy.sparse")
Apy = pysparse.csc_array((A.nzval, A.rowval .- 1, A.colptr .- 1), shape=size(A))

where A is a Julia SparseMatrixCSC. (The .- 1 is to convert from 1-based to 0-based indices.)

1 Like
5

I would call that something I didn’t expect, because I called Scipy function, and got back Julia matrix. This is not what you want if the intention is to call Scipy function again.

Btw maybe PythonCall.jl with its no-conversion default would serve you better

6

I think you probably misunderstood. sprandn is a Julia function from SparseArrays.

1 Like
7

Yep I totally got that backwards sorry