Умножение разреженных матриц , При использовании функции dot для умножения разреженной матрицы выдается ошибка как обойти избежать

[Urung]

V etoj programmme modul sg1 v kotoroj funkcia dot

__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
from scipy.sparse import linalg
from sg1 import sg1
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
def f(x, y, r):
    return 2*c*(y-y**2+r*(x**2-x))
def poisson(n, h, r):
    V = []
    I = []
    J = []
    b = []
    for j in range(n):
        for i in range(n):
            k = i+n*j
            b.append(h**2*f(i*h+h, j*h+h, r))
            V.append(2.-2*r)
            I.append(k)
            J.append(k)
            if i>0:
                V.append(-1.)
                I.append(k)
                J.append(k-1)
            if i<n-1:
                V.append(-1.)
                I.append(k)
                J.append(k+1)
            if j>0:
                V.append(r)
                I.append(k)
                J.append(k-n)
            if j<n-1:
                V.append(r)
                I.append(k)
                J.append(k+n)
    return sparse.coo_matrix((V, (I,J))), b
N =100
h = 1./N
c = 4.
r = 0.4
tst = time.clock()
A, b = poisson(N-1, h, r)
B = A.tocsr()
v, info = sg1(B, b,tol=1.e-08)
dt = time.clock() - tst
 
if iter == 0:
    print "N = %i, time solution = %1.3e" % (N, dt)
 
fig = plt.figure()
ax = Axes3D(fig)
x = np.linspace(h/2, 1.-h/2, N-1)
y = np.linspace(h/2, 1.-h/2, N-1)
X,Y = np.meshgrid(x, y)
Z = np.reshape(v, (N-1, N-1))
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral')
plt.show()

Vot modul

__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
 
def sg1(B,b,tol=1.0e-8):
    """
    Solve the linear system Ax=b by conjugate gradient method
    """
    n=len(b)
    x=np.zeros((n),'float')
    r=np.copy(b)
    for i in range(n):
        r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
    s=np.copy(r)
    As=np.zeros((n),'float')
    for k in range(n):
        for i in range(n):
            As[i]=np.dot(B[i,0:n],s[0:n])
        alpha=np.dot(r,r)/np.dot(s,As)
        x=x-alpha*s
        for i in range(n):
            r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
        if np.dot(r,r)<tol**2:
            break
        else:
            beta=-np.dot(r,As)/np.dot(s,As)
            s=r+beta*s
    return x,k

A vot shto soobchaet

/usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py
Traceback (most recent call last):
  File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module>
    v, info = sg1(B, b,tol=1.e-08)
  File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 193, in __sub__
    raise NotImplementedError('adding a scalar to a sparse '
NotImplementedError: adding a scalar to a sparse matrix is not supported
 
Process finished with exit code 1

[Славян]

Навряд ли это от умножения выдаётся. Действительно вычитание скаляра не поддерживаемо разреженными матрицами, а вы b вычитаете. Надо или как-то указать, что результат будет обычная матрица, или менять сам метод.
dot тут не при чём.

[Славян]

Что-то я прогоняю или не понимаю: там вообще строка из B на вектор множится, а потом число вычитается. В итоге всё в вектор уходит. Да, он не разрежен, но ошибку то пишет про матрицу, а не про вектор!

[Urung]

Я не знаю насколько правильно но я думаю так: В основной программе матрица V формируется в формате return sparse.coo_matrix((V, (I,J))), b то есть хранятся лишь координаты и значения ненулевых элементов матрицы. А dot ждет обычную запись. Я искал в этом направлении но не нашел

[Славян]

А просто попробуйте убрать вычитание b и увидите в чём проблема.

[Urung]

Я пробовал убирать это вычитание не получалось

Добавлено 21.07.14, 23:54

/usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py
Traceback (most recent call last):
  File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module>
    v, info = sg1(B, b,tol=1.e-08)
  File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=np.dot(B[i,0:n],x[0:n])
ValueError: setting an array element with a sequence.
 
Process finished with exit code 1

Добавлено 22.07.14, 00:18
А может подскажете как переформатировать с sparce cool в обычный формат чтобы dot заработал мне бы не хотелось менять основную программу

[Славян]

У вас одинаковая строка в двух местах, а ошибку пишет про одно место:

r[i]=np.dot(B[i,0:n],x[0:n])-b[i]

Попробуйте в одном месте написать так и посмотрим, в чём проблема:

r[i]=-b[i] + np.dot(B[i,0:n],x[0:n])

[Urung]

Если так

__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
 
def sg1(B,b,tol=1.0e-8):
    """
    Solve the linear system Ax=b by conjugate gradient method
    """
    n=len(b)
    x=np.zeros((n),'float')
    r=np.copy(b)
    for i in range(n):
        r[i]=-b[i]+np.dot(B[i,0:n],x[0:n])
    s=np.copy(r)
    As=np.zeros((n),'float')
    for k in range(n):
        for i in range(n):
            As[i]=np.dot(B[i,0:n],s[0:n])
        alpha=np.dot(r,r)/np.dot(s,As)
        x=x-alpha*s
        for i in range(n):
            r[i]=np.dot(B[i,0:n],x[0:n])
        if np.dot(r,r)<tol**2:
            break
        else:
            beta=-np.dot(r,As)/np.dot(s,As)
            s=r+beta*s
    return x,k

то так

Traceback (most recent call last):
  File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module>
    v, info = sg1(B, b,tol=1.e-08)
  File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=-b[i]+np.dot(B[i,0:n],x[0:n])
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 187, in __radd__
    return self.__add__(other)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 173, in __add__
    raise NotImplementedError('adding a scalar to a CSC or CSR '
NotImplementedError: adding a scalar to a CSC or CSR matrix is not supported

если так то так

/usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py
Traceback (most recent call last):
  File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module>
    v, info = sg1(B, b,tol=1.e-08)
  File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 193, in __sub__
    raise NotImplementedError('adding a scalar to a sparse '
NotImplementedError: adding a scalar to a sparse matrix is not supported
 
Process finished with exit code 1

Добавлено 22.07.14, 05:14
если убрать совсем строку и т д то второй дот все равно находит

[Urung]

Это программа прекрасно работает с библиотечным модулем linalg.cg который одинаково вычисляет с модулем с которым мы мучаемся в другой тестовой программе. То есть библиотечный модуль одинаково хорошо работает с любыми матрицами. Наш с вами модуль хорошо вычисляет в другой программе. Исходный код библиотечного модуля я не сумел найти. И код нашего модуля для меня более понятен и я хотел используя этот код попробовать похожие методы модули так вот приведу код который работает.

__author__ = 'kardashevskiy'
#Hyperbolic equation
import numpy as np
from scipy import sparse
from scipy.sparse import linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
def f(x, y, r):
    return 2*c*(y-y**2+r*(x**2-x))
def poisson(n, h, r):
    V = []
    I = []
    J = []
    b = []
    for j in range(n):
        for i in range(n):
            k = i+n*j
            b.append(h**2*f(i*h+h, j*h+h, r))
            V.append(2.-2*r)
            I.append(k)
            J.append(k)
            if i>0:
                V.append(-1.)
                I.append(k)
                J.append(k-1)
            if i<n-1:
                V.append(-1.)
                I.append(k)
                J.append(k+1)
            if j>0:
                V.append(r)
                I.append(k)
                J.append(k-n)
            if j<n-1:
                V.append(r)
                I.append(k)
                J.append(k+n)
    return sparse.coo_matrix((V, (I,J))), b
 
N =100
h = 1./N
c = 4.
r = 0.4
tst = time.clock()
A, b = poisson(N-1, h, r)
B = A.tocsr()
v, info = linalg.cg(B, b, tol=1.e-08)
dt = time.clock() - tst
 
if iter == 0:
    print "N = %i, time solution = %1.3e" % (N, dt)
 
fig = plt.figure()
ax = Axes3D(fig)
x = np.linspace(h/2, 1.-h/2, N-1)
y = np.linspace(h/2, 1.-h/2, N-1)
X,Y = np.meshgrid(x, y)
Z = np.reshape(v, (N-1, N-1))
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral')
plt.show()

[Славян]

Цитата Urung @ 22.07.14, 05:11

если убрать совсем строку и т д то второй дот все равно находит

А давайте так попробуем: вместо одной итерации:

for i in range(n):
   r[i]=np.dot(B[i,0:n],x[0:n])-b[i]

Запишем пару:

for i in range(n):
    r[i] = np.dot(B[i,0:n],x[0:n])
    r[i] = r[i]-b[i]

Я не знаю как у вас begin/end, потому напишите как надо.
Вот тогда и посмотрим на что будет ругаться, на вычитание или на умножение.

Добавлено 22.07.14, 16:27
Ай, тьфу, вы ж в №6 написали, что не работает так. Виноват. Думаю.

Добавлено 22.07.14, 16:36
Хм... смотрите, интересный момент:
1.Тут пишет ошибку:

Цитата Urung @ 21.07.14, 23:48

File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=np.dot(B[i,0:n],x[0:n])
ValueError: setting an array element with a sequence.

2.а тут норма:

As[i]=np.dot(B[i,0:n],s[0:n])

хотя и x и s - одинаковые вектора, как и As и r. В чём тонкость?

[Славян]

Забавно, но по этому коду:

Цитата Urung @ 20.07.14, 23:37

x=np.zeros((n),'float')
r=np.copy(b)
for i in range(n):
    r[i]=np.dot(B[i,0:n],x[0:n])-b[i]

можно сделать вывод, что r = -b, т.к. x - нулевой.

[Urung]

СлавянСпасибо за участие но тут кажется тупик. Но вот может быть посоветуете найти исходный код библиотечного ьодуля linalg.cg я пытался искал но не смог найти я попробовал бы преобразовать тот модуль наверняка для вас это не проблема наверное вы вовсю используете модифицируете и новые делаете. А мы чайники даже не можем найти.

[Славян]

Навряд ли "посоветую". Я ни в питоне ни в Руби "не шарю", а вам советую создать тему в разделе С/С++ Алгоритмы про linalg, а там много умного народа бывает, они явно укажут. Даже если найдёте не питоновский, то переделать на него будет не шибко сложно. Так что "рвите когти" туда!

[Urung]

Спасибо я вам очень признателен за то что приняли очень живое участие. Спасибо попробую надеюсь не побьют.

[amk]

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

Если так, то надо учитывать, что например в MathCAD'е или MATLAB'е матрица размером 1х1 может быть использована там, где нужен скаляр (отдельное число). В Python'овских библиотеках numpy и scipy не так - матрица остаётся матрицей, даже если она вообще не содержит элементов или содержит единственный элемент. В частности чтобы получить из такой одноэлементной матрицы её единственный элемент в пакете numpy есть функция asscalar.
Можешь попробовать написать

r[i] = np.asscalar(np.dot(B[i,0:n],x[0:n])) - b[i]

может быть заработает. Я бы сам проверил, но у меня стоит Python 3, а программа написана на каком-то из 2-х