CODEWITHSUNDEEP
numpyscipysparse-matrix
Matt
Matt

Reputation: 217

Solving sparse and complex linear systems with scipy

Good afternoon,

I'm working on developing some numerical code, and at a key point I need to solve a sparse linear system. I have this all set up and tested with real floating point numbers. The issue occurs when I need to solve a system Ax = b where A is a sparse and real matrix (stored in CSR format with scipy) and a complex RHS b, stored as a numpy array. When I solve the system, I only get the real values of the solution, along with the warning

 /System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/core/numeric.py:460: 
 ComplexWarning: Casting complex values to real discards the imaginary part 
 return array(a, dtype, copy=False, order=order)

Here's a minimal working example:

 import numpy as np
 import scipy.sparse.linalg
 import scipy.sparse

 A = scipy.sparse.csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])

 b = np.array([1,1,1]) +1j*np.array([1,1,1])

 # throws warning and returns erroneous result
 x_complex = scipy.sparse.linalg.spsolve(A,b)

Is this a scipy issue or am I doing something wrong? Is there a possible workaround? Any help would be greatly appreciated.

Edit: I've found a bit of a workaround where we can just solve two linear systems with real right hand sides (doing the obvious thing of splitting the RHS into real and imaginary parts). Is there a better way to do this though? It will make the code somewhat less clear.

Upvotes: 0

Views: 2790

Answers (1)

user6655984
user6655984

Reputation:

The solver is trying to cast the matrix and the right-hand side to the same data type in order to do the computations. Since you have complex numbers involved, they should be cast to complex, but for some reason it's trying to cast to real floating point numbers (a bug? or just the assumption that the data type of the matrix is the one to use). Making A of complex data type solves the problem:

A = scipy.sparse.csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]], dtype=np.cfloat)

Now the solver returns the correct solution,

[-0.16666667-0.16666667j  0.58333333+0.58333333j  0.33333333+0.33333333j]

Upvotes: 3

Related Questions