Search code examples
pythonscipylinear-algebra

Print current residual from callback in scipy.sparse.linalg.cg


I am using scipy.sparse.linalg.cg to solve a large, sparse linear system, and it works fine, except that I would like to add a progress report, so that I can monitor the residual as the solver works. I've managed to set up a callback, but I can't figure out how to access the current residual from inside the callback. Calculating the residual myself is possible, of course, but that is a rather heavy operation, which I'd like to avoid. Have I missed something, or is there no efficient way of getting at the residual?


Solution

  • The callback is only sent xk, the current solution vector. So you don't have direct access to the residual. However, the source code shows resid is a local variable in the cg function.

    So, with CPython, it is possible to use the inspect module to peek at the local variables in the caller's frame:

    import inspect
    import numpy as np
    import scipy as sp
    import scipy.sparse as sparse
    import scipy.sparse.linalg as splinalg
    import random
    
    def report(xk):
        frame = inspect.currentframe().f_back
        print(frame.f_locals['resid'])
    
    N = 200
    A = sparse.lil_matrix( (N, N) )
    for _ in xrange(N):
        A[random.randint(0, N-1), random.randint(0, N-1)] = random.randint(1, 100)
    
    b = np.random.randint(0, N-1, size = N)
    x, info = splinalg.cg(A, b, callback = report)