[SciPy-User] why is GMRES slow ? [py=3.4.3 numpy=1.10 scipy=0.16]

[zimon toutoune]

Dear,

I am looking for an explanation why GMRes [scipy.sparse.linalg.gmres]
running 'Niter' iterations is much more slower than applying 'Niter'
times 'matvec' (or 'dot').

For sure, there is the Arnoldi orthogonalization, the Givens rotation,
but the gap is more huge than that.

Moreover, for dense matrices [numpy.ndarray], solving with
'scipy.linalg.solve', i.e. performing a LU decomposition, is really,
really faster. (and hopefully same result using 'numpy.linalg.solve')

But LU is about O(N**3) and Krylov about O(Niter*N**2), therefore
fixing 'Niter' small and 'N' big enough should show significant
differences.

I am not talking about convergence, only about the expected performances.
And I know that solving random matrices is not a well-thought idea,
but it appears to me a way to quickly benchmark since I am just
interesting about performances.

Let see the code below.
In short, the matrix is complex and the size is 5000.
There is no restart and maxiter is 50.
And the results are:

> In [1]: run test.py
solve
gmres
mydot
7.203553199768066
26.598630666732788
1.2706046104431152
50

> In [2]: %timeit  scipy.linalg.solve(A, b)
1 loops, best of 3: 7.15 s per loop

> In [3]: %timeit  scipy.sparse.linalg.gmres(A, b, maxiter=m, restart=r)
1 loops, best of 3: 25.7 s per loop

> In [4]: %timeit mydot(A, b, maxiter=info)
1 loops, best of 3: 1.14 s per loop


 - Is it expected ?
 - Is it reproductible ?
(I do with 2 versions of Numpy/Scipy on MacOS and Debian,  and so I am
missing an installation problem ?, and the Debian one runs
python=3.4.3 numpy=1.10 scipy=0.16 installed through conda)

 - Does it come from the constants behind big-O ?
 (size versus time should show confirm or not, right ?)

 - Does it come from a detail about the implementation ?
(all is calling BLAS/LAPACK, right ?)

If I miss a point and my questions are naive, then oups!
Because why ?

Best Regards,
simon


~~python code

#!/usr/env python
# coding: utf8

import numpy
import scipy.linalg
import scipy.sparse.linalg

from time import time

N = 5000
m = 50
r = None

A = numpy.random.rand(N, N) + 1j * numpy.random.rand(N, N)
b = numpy.random.rand(N) + 1j * numpy.random.rand(N)

def mydot(A, x, maxiter=5):
    for i in range(maxiter):
        A.dot(x)
        numpy.sum(x * x)

print('solve')
tt = time()
x = scipy.linalg.solve(A, b)
ts = time() - tt

print('gmres')
tt = time()
y, info = scipy.sparse.linalg.gmres(A, b, maxiter=m, restart=r)
tg = time() - tt

if info <= 0:
    info = m

print('mydot')
tt = time()
mydot(A, b, maxiter=info)
tm = time() - tt

print(ts)
print(tg)
print(tm)
print(info)