使用稀疏稀疏矩阵求解方程组 [英] Using scipy sparse matrices to solve system of equations

问题描述

这是对的后续操作并在python中解决联立方程，但我觉得任何答案都应该拥有自己的声誉点.

This is a follow up to How to set up and solve simultaneous equations in python but I feel deserves its own reputation points for any answer.

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

For a fixed integer n, I have a set of 2(n-1) simultaneous equations as follows.

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)

为1 <= p <= n-1定义了

M(p).为0 <= p <= n-2定义了N(p).还请注意，p在每个方程式中都是一个常数整数，因此整个系统是线性的.

M(p) is defined for 1 <= p <= n-1. N(p) is defined for 0 <= p <= n-2. Notice also that p is just a constant integer in every equation so the whole system is linear.

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

Some very nice answers were given for how to set up a system of equations in python. However, the system is sparse and I would like to solve it for large n. How can I use scipy's sparse matrix representation and http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html for example instead?

推荐答案

我通常不会继续击败一匹死马，但是碰巧我的非矢量化方法可以解决您的其他问题，当事情变大时，它有一些优点. .因为实际上我一次只填充一个系数矩阵，所以将其转换为COO稀疏矩阵格式非常容易，可以将其有效地转换为CSC并求解.做到这一点:

I wouldn't normally keep beating a dead horse, but it happens that my non-vectorized approach to solving your other question, has some merit when things get big. Because I was actually filling the coefficient matrix one item at a time, it is very easy to translate into COO sparse matrix format, which can efficiently be transformed to CSC and solved. The following does it:

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一样，它当然比从矢量化块中构建更为冗长，但是当您同时使用这两种方法时，会发生一件有趣的事情:

It is of course much more verbose than building it from vectorized blocks, as TheodorosZelleke does, but an interesting thing happens when you time both approaches:

首先，这(非常)不错，两种解决方案中的时间都是线性增长的，就像使用稀疏方法所期望的那样.但是我在这个答案中给出的解决方案总是更快，对于更大的n来说更是如此.只是为了好玩，我还从另一个问题中选择了TheodorosZelleke的密集方法，它给出了这张漂亮的图，显示了两种类型的解决方案的不同缩放比例，并且在n = 75左右的某个位置，该解决方案应该是很早的选择:

First, and this is (very) nice, time is scaling linearly in both solutions, as one would expect from using the sparse approach. But the solution I gave in this answer is always faster, more so for larger ns. Just for the fun of it, I also timed TheodorosZelleke's dense approach from the other question, which gives this nice graph showing the different scaling of both types of solutions, and how very early, somewhere around n = 75, the solution here should be your choice:

我对scipy.sparse的了解还不足以真正弄清楚两种稀疏方法之间的区别的原因，尽管我怀疑大量使用LIL格式的稀疏矩阵.通过将TheodorosZelleke的答案转换为COO格式，尽管代码中有很多紧凑性，但是可能会获得非常微不足道的性能提升.但这是OP的练习！

I don't know enough about scipy.sparse to really figure out why the differences between the two sparse approaches, although I suspect heavily of the use of LIL format sparse matrices. There may be some very marginal performance gain, although a lot of compactness in the code, by turning TheodorosZelleke's answer into COO format. But that is left as an exercise for the OP!

这篇关于使用稀疏稀疏矩阵求解方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！