使用 scipy 稀疏矩阵求解方程组

Question

这是对How to set up and solve simultaneous equations in python 的跟进，但我觉得任何回答都值得拥有自己的声誉积分。

对于一个固定整数n，我有一组2(n-1) 联立方程如下。

M(p) = 1+((n-p-1)/n)*M(n-1) + (2/n)*N(p-1) + ((p-1)/n)*M(p-1)

N(p) = 1+((n-p-1)/n)*M(n-1) + (p/n)*N(p-1)

M(1) = 1+((n-2)/n)*M(n-1) + (2/n)*N(0)

N(0) = 1+((n-1)/n)*M(n-1)

M(p) 是为1 <= p <= n-1 定义的。 N(p) 是为 0 <= p <= n-2 定义的。另请注意，p 在每个方程中只是​​一个常数整数，因此整个系统是线性的。

对于如何在 python 中建立方程组给出了一些非常好的答案。但是，该系统是稀疏的，我想为大 n 解决它。例如，如何使用 scipy 的稀疏矩阵表示和http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html？

Answer 1

这里有一些我以前看过的代码：http://jkwiens.com/heat-equation-using-finite-difference/ 他的函数实现了一种有限差分法，使用 scipy 稀疏矩阵包求解热方程。

Answer 2

这是一个使用 scipy.sparse 的解决方案。不幸的是，这里没有说明问题。所以为了理解这个解决方案，未来的访问者必须首先在问题中提供的链接下查找问题。

使用 scipy.sparse 的解决方案：

from scipy.sparse import spdiags, lil_matrix, vstack, hstack
from scipy.sparse.linalg import spsolve
import numpy as np


def solve(n):
    nrange = np.arange(n)
    diag = np.ones(n-1)

    # upper left block
    n_to_M = spdiags(-2. * diag, 0, n-1, n-1)

    # lower left block
    n_to_N = spdiags([n * diag, -nrange[-1:0:-1]], [0, 1], n-1, n-1)

    # upper right block
    m_to_M = lil_matrix(n_to_N)
    m_to_M[1:, 0] = -nrange[1:-1].reshape((n-2, 1))

    # lower right block
    m_to_N = lil_matrix((n-1, n-1))
    m_to_N[:, 0] = -nrange[1:].reshape((n-1, 1))

    # build A, combine all blocks
    coeff_mat = hstack(
                       (vstack((n_to_M, n_to_N)),
                        vstack((m_to_M, m_to_N))))

    # const vector, right side of eq.
    const = n * np.ones((2 * (n-1),1))

    return spsolve(coeff_mat.tocsr(), const).reshape((-1,1))

Answer 3

我通常不会一直打死马，但碰巧我解决您的其他问题的非矢量化方法在事情变得大时有一些优点。因为我实际上是一次填充系数矩阵，所以很容易转换成COO稀疏矩阵格式，可以高效地转换为CSC并求解。以下是这样做的：

import scipy.sparse

def sps_solve(n) :
    # Solution vector is [N[0], N[1], ..., N[n - 2], M[1], M[2], ..., M[n - 1]]
    n_pos = lambda p : p
    m_pos = lambda p : p + n - 2
    data = []
    row = []
    col = []
    # p = 0
    # n * N[0] + (1 - n) * M[n-1] = n
    row += [n_pos(0), n_pos(0)]
    col += [n_pos(0), m_pos(n - 1)]
    data += [n, 1 - n]
    for p in xrange(1, n - 1) :
        # n * M[p] + (1 + p - n) * M[n - 1] - 2 * N[p - 1] +
        #  (1 - p) * M[p - 1] = n
        row += [m_pos(p)] * (4 if p > 1 else 3)
        col += ([m_pos(p), m_pos(n - 1), n_pos(p - 1)] +
                ([m_pos(p - 1)] if p > 1 else []))
        data += [n, 1 + p - n , -2] + ([1 - p] if p > 1 else [])
        # n * N[p] + (1 + p -n) * M[n - 1] - p * N[p - 1] = n
        row += [n_pos(p)] * 3
        col += [n_pos(p), m_pos(n - 1), n_pos(p - 1)]
        data += [n, 1 + p - n, -p]
    if n > 2 :
        # p = n - 1
        # n * M[n - 1] - 2 * N[n - 2] + (2 - n) * M[n - 2] = n
        row += [m_pos(n-1)] * 3
        col += [m_pos(n - 1), n_pos(n - 2), m_pos(n - 2)]
        data += [n, -2, 2 - n]
    else :
        # p = 1 
        # n * M[1] - 2 * N[0] = n
        row += [m_pos(n - 1)] * 2
        col += [m_pos(n - 1), n_pos(n - 2)]
        data += [n, -2]
    coeff_mat = scipy.sparse.coo_matrix((data, (row, col))).tocsc()
    return scipy.sparse.linalg.spsolve(coeff_mat,
                                       np.ones(2 * (n - 1)) * n)

当然，它比从矢量化块构建它要冗长得多，就像 TheodorosZelleke 所做的那样，但是当你对这两种方法都计时时会发生一件有趣的事情：

首先，这（非常）好，时间在两种解决方案中都是线性缩放的，正如人们对使用稀疏方法所期望的那样。但是我在这个答案中给出的解决方案总是更快，对于更大的ns 更是如此。只是为了好玩，我还从另一个问题中计时了 TheodorosZelleke 的密集方法，它给出了这个很好的图表，显示了两种类型的解决方案的不同比例，以及多早，在 n = 75 附近的某个地方，这里的解决方案应该是你的选择：

我对@9​​87654326@ 的了解还不够，无法真正弄清楚为什么这两种稀疏方法之间存在差异，尽管我严重怀疑使用 LIL 格式的稀疏矩阵。通过将 TheodorosZelleke 的答案转换为 COO 格式，尽管代码非常紧凑，但可能会有一些非常小的性能提升。但这留给 OP 做练习！