4

I am attempting to use Scipy's gmres command. My matrix is dense and large, so my plan was to use the LinearOperator command to return only the matrix vector product. Please see my previous question here LinearOperator with Two Inputs. With help from that question, I was able to build a LinearOperator object which successfully does the computation A*x where A is the matrix and x is the vector.

The problem is when I call gmres. I run the command:

x, info = scipy.sparse.linalg.gmres(A, b)

and it returns the error that the operator A has no dtype. This is true, as A.dtype returns an error. My problem is that I have no idea what to specify for dtype. When I construct my linear operator, there is an optional parameter to pass for dtype, but I do not know what to give it.
I tried passing dtype='float64' and that froze my IDE, so I suspect I'm wrong there.
Trying dtype = None just yields the default, where dtype is not defined.

I also tried simply defining A leaving dtype blank, and then typing A.dtype = None. This actually gives the attribute A.dtype and yields another error when calling A.

This seems to be linked to another problem, which is that the gmres seems to want a preconditioner given to it. I don't actually have a pre-conditioner I want to give it, so it tries to construct one and it tries to use the same dtype as A but, since A has no dtype attribute it errors out. Any suggestion would be greatly appreciated.

A = sparse.scipy.linalg.linearoperator(shape = (n,n), matvec = mv, dtype = 'float64')

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/isolve/iterative.py", line 393, in gmres
    A,M,x,b,postprocess = make_system(A,M,x0,b,xtype)

  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/isolve/utils.py", line 119, in make_system

M = LinearOperator(A.shape, matvec=psolve, rmatvec=rpsolve, dtype=A.dtype)
AttributeError: LinearOperator instance has no attribute 'dtype'
Community
  • 1
  • 1
user35959
  • 359
  • 2
  • 13
  • Is there a particular reason you're using the sparse module for dense matrix operations? – askewchan Nov 28 '13 at 02:57
  • No, I would prefer not to, but I could only find it in the sparse library. There was a scipy.linalg.gmres, but it seems to either be deprecated or not showing up in my version of scipy. The documentation claims it works for sparse, dense or LinearOperator type matrices. – user35959 Nov 28 '13 at 04:08
  • 1
    Note that there is no need to use LinearOperator, you can just directly pass the dense matrix to `gmres(a, b)`. – pv. Nov 28 '13 at 09:52
  • 1
    I cannot load the matrix directly because it takes up too much RAM for the computer. – user35959 Nov 28 '13 at 15:47

1 Answers1

4

(This is more of an extended comment than an answer.)

What version of scipy are you using? I'm using scipy 0.13.0. If I don't specify the dtype when I create the LinearOperator, I get the dtype error that you get. But specifying dtype='float64' works for me:

In [1]: import numpy as np

In [2]: from scipy.sparse.linalg import LinearOperator, gmres

In [3]: def mymatvec(v):
   ...:     a = np.array([[4,2,1],[2,2,1],[1,1,1]])
   ...:     return a.dot(v)
   ...: 

In [4]: A = LinearOperator((3,3), mymatvec, dtype='float64')

In [5]: b = np.array([1,2,3])

In [6]: gmres(A, b)
Out[6]: (array([-0.5, -0.5,  4. ]), 0)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • I am using Scipy version 0.9.0. I unfortunately think I'm stuck with this version, as it is the version on the computer I plan on running this code on which I do not have the ability to install new versions on. When I add dtype='float64' and try calling A from the interpreter, my Python IDE completely freezes. I can add code if that would be helpful. – user35959 Nov 28 '13 at 04:45
  • 1
    Interestingly enough, when I run this in my Python IDE (Spyder), it freezes if I ever try calling A, but if I just execute from the command line python my_file.py it runs without a problem after setting it to dtype='float64'. Maybe the problem is now just with my IDE. I haven't confirmed yet its giving me the "right" result, but otherwise it looks like your suggestion works for me. Thanks! – user35959 Nov 28 '13 at 04:53