共轭梯度实现 Python答案

【问题标题】：Conjugate Gradient implementation Python共轭梯度实现 Python
【发布时间】：2018-12-07 07:59:57
【问题描述】：

我通过查看维基百科参考在 python 中实现了共轭梯度 - https://en.wikipedia.org/wiki/Conjugate_gradient_method

实现应该解决

ax = b

我的应用程序输入如下，

a = <400x400 sparse matrix of type '<class 'numpy.float64'>'
    with 1920 stored elements in Compressed Sparse Row format>
b = vector of shape (400, ) and dtype = float64
x = vector of random numbers of shape (400, )

这是我的实现 -

def ConjGrad(a, b, x):
       r = (b - np.dot(np.array(a), x));
       p = r;
       rsold = np.dot(r.T, r);

       for i in range(len(b)):
           a_p = np.dot(a, p);
           alpha = rsold / np.dot(p.T, a_p);
           x = x + (alpha * p);
           r = r - (alpha * a_p);
           rsnew = np.dot(r.T, r);
           if (np.sqrt(rsnew) < (10 ** -5)):
               break;
           p = r + ((rsnew / rsold) * p);
           rsold = rsnew;        
       return p

当我调用上述 CG 函数时，我在该行的函数中得到一个错误 -

r = (b - np.dot(np.array(a), x));

错误是这样的 -

NotImplementedError: subtracting a sparse matrix from a nonzero scalar is 
not supported

在运行时，以下是 CG 函数中变量的属性 -

np.dot(np.array(a), x).shape
(400,)
b.shape
(400,)

我想知道为什么减法没有发生？？？

我用下面的示例输入参数测试了相同的函数，它运行良好。

a = np.array([[3, 2, -1], [2, -1, 1], [-1, 1, -1]]) # 3X3 symmetric matrix
b = (np.array([1, -2, 0])[np.newaxis]).T  # 3X1 matrix
x = (np.array([0, 1, 2])[np.newaxis]).T

有人可以告诉我为什么它不适用于稀疏矩阵吗？

【问题讨论】：

标签： python numpy optimization sparse-matrix matrix-multiplication

【解决方案1】：

当将稀疏矩阵乘以数组时，您不应使用：np.dot(np.array(a), x)) 但 a.dot(x)。请参阅以下文档：

https://docs.scipy.org/doc/scipy/reference/sparse.html

遵循正确的惯例：

    def conjGrad(A,x,b,tol,N):

    r = b - A.dot(x)
    p = r.copy()
    for i in range(N):
        Ap = A.dot(p)
        alpha = np.dot(p,r)/np.dot(p,Ap)
        x = x + alpha*p
        r = b - A.dot(x)
        if np.sqrt(np.sum((r**2))) < tol:
            print('Itr:', i)
            break
        else:
            beta = -np.dot(r,Ap)/np.dot(p,Ap)
            p = r + beta*p
    return x

【讨论】：