Сопряженная реализация градиента Python - PullRequest
Новая
       50

Сопряженная реализация градиента Python

0 голосов
/ 07 декабря 2018

Я реализовал Conjugate Gradient в 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

Может кто-нибудь сказать, почему она не работает для разреженной матрицы?

...