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)