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 |
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 |
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 |
Free forum by Nabble | Edit this page |