2.5.3. Linear System Solvers

  • 稀疏矩阵/特征值问题求解器存在于scipy.sparse.linalg

  • 子模块:
    • dsolve:用于求解线性系统的直接因式分解方法
    • isolve:用于求解线性系统的迭代方法
    • eigen:稀疏特征值问题求解器
  • 所有求解器可从:

    >>> import scipy.sparse.linalg as spla
    
    >>> spla.__all__
    ['LinearOperator', 'Tester', 'arpack', 'aslinearoperator', 'bicg',
    'bicgstab', 'cg', 'cgs', 'csc_matrix', 'csr_matrix', 'dsolve',
    'eigen', 'eigen_symmetric', 'factorized', 'gmres', 'interface',
    'isolve', 'iterative', 'lgmres', 'linsolve', 'lobpcg', 'lsqr',
    'minres', 'np', 'qmr', 'speigs', 'spilu', 'splu', 'spsolve', 'svd',
    'test', 'umfpack', 'use_solver', 'utils', 'warnings']

2.5.3.1. Sparse Direct Solvers

  • 默认解算器:SuperLU 4.0
    • 包括在SciPy
    • 真实和复杂的系统
    • 单精度和双精度
  • 可选:umfpack
    • 真实和复杂的系统
    • 双精度
    • 推荐性能
    • 封装容器现在位于scikits.umfpack
    • 签出Nathaniel Smith的新scikits.suitesparse

2.5.3.1.1. Examples

  • 导入整个模块,并查看其docstring:

    >>> from scipy.sparse.linalg import dsolve
    
    >>> help(dsolve)
  • 可以使用superlu和umfpack(如果后者安装)如下:

    • 准备线性系统:

      >>> import numpy as np
      
      >>> from scipy import sparse
      >>> mtx = sparse.spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
      >>> mtx.todense()
      matrix([[ 1, 5, 0, 0, 0],
      [ 0, 2, 8, 0, 0],
      [ 0, 0, 3, 9, 0],
      [ 0, 0, 0, 4, 10],
      [ 0, 0, 0, 0, 5]])
      >>> rhs = np.array([1, 2, 3, 4, 5], dtype=np.float32)
    • 求解为单精度实数:

      >>> mtx1 = mtx.astype(np.float32)
      
      >>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
      >>> print(x)
      [ 106. -21. 5.5 -1.5 1. ]
      >>> print("Error: %s" % (mtx1 * x - rhs))
      Error: [ 0. 0. 0. 0. 0.]
    • 求解为双精度实数:

      >>> mtx2 = mtx.astype(np.float64)
      
      >>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
      >>> print(x)
      [ 106. -21. 5.5 -1.5 1. ]
      >>> print("Error: %s" % (mtx2 * x - rhs))
      Error: [ 0. 0. 0. 0. 0.]
    • 求解为单精度复数:

      >>> mtx1 = mtx.astype(np.complex64)
      
      >>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
      >>> print(x)
      [ 106.0+0.j -21.0+0.j 5.5+0.j -1.5+0.j 1.0+0.j]
      >>> print("Error: %s" % (mtx1 * x - rhs))
      Error: [ 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
    • 求解为双精度复数:

      >>> mtx2 = mtx.astype(np.complex128)
      
      >>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
      >>> print(x)
      [ 106.0+0.j -21.0+0.j 5.5+0.j -1.5+0.j 1.0+0.j]
      >>> print("Error: %s" % (mtx2 * x - rhs))
      Error: [ 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
"""
Solve a linear system
=======================
Construct a 1000x1000 lil_matrix and add some values to it, convert it
to CSR format and solve A x = b for x:and solve a linear system with a
direct solver.
"""
import numpy as np
import scipy.sparse as sps
from matplotlib import pyplot as plt
from scipy.sparse.linalg.dsolve import linsolve
rand = np.random.rand
mtx = sps.lil_matrix((1000, 1000), dtype=np.float64)
mtx[0, :100] = rand(100)
mtx[1, 100:200] = mtx[0, :100]
mtx.setdiag(rand(1000))
plt.clf()
plt.spy(mtx, marker='.', markersize=2)
plt.show()
mtx = mtx.tocsr()
rhs = rand(1000)
x = linsolve.spsolve(mtx, rhs)
print('rezidual: %r' % np.linalg.norm(mtx * x - rhs))

2.5.3.2. Iterative Solvers

  • isolve模块包含以下解算器:
    • bicg(BIConjugate Gradient)
    • bicgstab(BIC交联梯度稳定)
    • cg(共轭梯度) - 仅对称正定矩阵
    • cgs(Conjugate Gradient Squared)
    • gmres(广义最小残次)
    • minres(MINimum RESidual)
    • qmr(准最小残差)

2.5.3.2.1. Common Parameters

  • 强制:

    一个
    {稀疏矩阵,密集矩阵,LinearOperator}

    线性系统的N×N矩阵。

    b
    {array,matrix}

    线性系统的右手侧。具有形状(N,)或(N,1)。

  • 可选的:

    x0
    {array,matrix}

    开始猜测解决方案。

    托尔
    浮点数

    相对容差在终止前实现。

    maxiter
    整数

    最大迭代次数。迭代将在maxiter步骤后停止,即使尚未达到指定的容差。

    M
    {稀疏矩阵,密集矩阵,LinearOperator}

    A的预处理程序预处理器应近似于A的逆。有效预处理显着地提高了收敛速率,这意味着需要更少的迭代来达到给定的容差。

    回电话
    功能

    用户提供的函数在每次迭代后调用。它被称为回调(xk),其中xk是当前的解向量。

2.5.3.2.2. LinearOperator Class

from scipy.sparse.linalg.interface import LinearOperator
  • 公共接口用于执行矩阵向量乘积
  • 有用的抽象使得在求解器中使用密集和稀疏的矩阵,以及无矩阵的
  • 具有形状matvec()(+一些可选参数)
  • 例:
>>> import numpy as np
>>> from scipy.sparse.linalg import LinearOperator
>>> def mv(v):
... return np.array([2*v[0], 3*v[1]])
...
>>> A = LinearOperator((2, 2), matvec=mv)
>>> A
<2x2 LinearOperator with unspecified dtype>
>>> A.matvec(np.ones(2))
array([ 2., 3.])
>>> A * np.ones(2)
array([ 2., 3.])

2.5.3.2.3. A Few Notes on Preconditioning

  • 问题具体
  • 往往很难发展
  • 如果不确定,请尝试ILU
    • 可在dsolve中作为spilu()

2.5.3.3. Eigenvalue Problem Solvers

2.5.3.3.1. The eigen module

  • arpack *一组Fortran77子程序,用于解决大规模特征值问题

  • lobpcg(局部最优块预处理的共轭梯度法)*与Nathan Bell的PyAMG

    """
    
    Compute eigenvectors and eigenvalues using a preconditioned eigensolver
    ========================================================================
    In this example Smoothed Aggregation (SA) is used to precondition
    the LOBPCG eigensolver on a two-dimensional Poisson problem with
    Dirichlet boundary conditions.
    """
    import scipy
    from scipy.sparse.linalg import lobpcg
    from pyamg import smoothed_aggregation_solver
    from pyamg.gallery import poisson
    N = 100
    K = 9
    A = poisson((N,N), format='csr')
    # create the AMG hierarchy
    ml = smoothed_aggregation_solver(A)
    # initial approximation to the K eigenvectors
    X = scipy.rand(A.shape[0], K)
    # preconditioner based on ml
    M = ml.aspreconditioner()
    # compute eigenvalues and eigenvectors with LOBPCG
    W,V = lobpcg(A, X, M=M, tol=1e-8, largest=False)
    #plot the eigenvectors
    import pylab
    pylab.figure(figsize=(9,9))
    for i in range(K):
    pylab.subplot(3, 3, i+1)
    pylab.title('Eigenvector %d' % i)
    pylab.pcolor(V[:,i].reshape(N,N))
    pylab.axis('equal')
    pylab.axis('off')
    pylab.show()
  • 示例由Nils Wagner:

  • 输出:

    $ python examples/lobpcg_sakurai.py
    
    Results by LOBPCG for n=2500
    [ 0.06250083 0.06250028 0.06250007]
    Exact eigenvalues
    [ 0.06250005 0.0625002 0.06250044]
    Elapsed time 7.01
../../_images/lobpcg_eigenvalues.png