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

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[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)
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: why is GMRES slow ? [py=3.4.3 numpy=1.10 scipy=0.16]

Pauli Virtanen-3
Tue, 01 Dec 2015 21:48:27 -0300, zimon toutoune kirjoitti:
> 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').

Each gmres outer iteration requires ~ `restart` matrix-vector products.
restart=None means the default value, restart=20.
The returned `info` contains the number of outer iterations.

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: why is GMRES slow ? [py=3.4.3 numpy=1.10 scipy=0.16]

zimon toutoune
On 2 December 2015 at 12:26, Pauli Virtanen <[hidden email]> wrote:
> Tue, 01 Dec 2015 21:48:27 -0300, zimon toutoune kirjoitti:
>> 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').
>
> Each gmres outer iteration requires ~ `restart` matrix-vector products.
> restart=None means the default value, restart=20.
> The returned `info` contains the number of outer iterations.

Sorry for the naive question.
Now completely clear.

This means that by default gmres is not gmres(full).
And the docstring of 'info' is a bit disturbing.
So clearly the total of number of iterations is always: Niter = maxiter*restart.

However, the only way to use gmres(full) is to set: 'maxiter=1' and
'restart=val'.
Which is also a bit confusing.
The matlab behavior seems clearer, no ?
 - if restart==None then gmres does not restart.
 - if maxiter!=None and restart==None then the maximum number of total
iterations is maxiter (instead of restart*maxiter).

Last, why not collect more information in info, as the residual
history, the number of outer iterations, the number of inner
iteration, etc. ?
Is there a reason ? or just because nobody has added it ?

Thanks for your times and your quick and clear answer.

All the best,
simon
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user