用 scipy 求解稀疏和复杂的线性系统

Question

下午好，

我正在开发一些数字代码，在关键点上我需要求解一个稀疏线性系统。我已经用真正的浮点数设置并测试了这一切。当我需要求解系统 Ax = b 时会出现问题，其中 A 是稀疏实数矩阵（以 scipy 的 CSR 格式存储）和复杂的 RHS b，存储为 numpy 数组。当我解决系统时，我只得到解决方案的实际值，以及警告

 /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)

这是一个最小的工作示例：

 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)

这是一个 scipy 问题还是我做错了什么？有可能的解决方法吗？任何帮助将不胜感激。

编辑：我发现了一些解决方法，我们可以只用实数右侧求解两个线性系统（做明显的事情，将 RHS 分成实部和虚部）。有没有更好的方法来做到这一点？这会使代码不太清晰。

Answer 1

求解器正在尝试将矩阵和右侧转换为相同的数据类型以便进行计算。由于涉及复数，因此应将它们转换为复数，但由于某种原因，它试图转换为真正的浮点数（错误？或者只是假设矩阵的数据类型是要使用的数据类型）。制作复杂数据类型的A解决了问题：

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

现在解算器returns正确的解，

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

用 scipy 求解稀疏和复杂的线性系统

Solving sparse and complex linear systems with scipy

numpy

scipy

sparse-matrix