[SciPy-user] Convergence behaviour of iterative solvers

Nils Wagner nwagner at iam.uni-stuttgart.de
Mon Jan 22 06:54:53 CST 2007


Ed Schofield wrote:
>
>
> On 1/19/07, *Nils Wagner* <nwagner at iam.uni-stuttgart.de
> <mailto:nwagner at iam.uni-stuttgart.de>> wrote:
>
>
>     So how can I use this functionality within the following
>     script
>     from scipy import *
>
>     def func(x):
>          return 0.5*dot(x,dot(A,x))-dot(x,b)
>
>     def callback(x):
>          return linalg.norm(dot(A,x)-b)
>
>     n = 10
>     x0 = zeros(n,float)
>     A = random.rand(n,n)+diag(4*ones(n))
>     A = 0.5*(A+A.T)
>     b = random.rand(n)
>
>     x = optimize.fmin_cg(func,x0)
>
>
> Hi Nils,
>
> The callback function's return value is ignored. Make it do something,
> e.g.:
>
> def callback(x):
>     print "||A.x - b|| = " + str(linalg.norm(dot(A,x)-b))
>
> Then call it with:
>
> x = optimize.fmin_cg(func, x0, callback=callback)
>
> -- Ed
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>   
Hi Ed,

Thank you very much! Where should I place the following two lines

if callback is not None:
     callback(x)

Here is an excerpt of iterative.py (namely cg)

    while 1:
        x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
           revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
        slice1 = slice(ndx1-1, ndx1-1+n)
        slice2 = slice(ndx2-1, ndx2-1+n)
        if (ijob == -1):
            if callback is not None:
                callback(x)
            break
        elif (ijob == 1):
            if matvec is None:
                matvec = get_matvec(A)
            work[slice2] *= sclr2
            work[slice2] += sclr1*matvec(work[slice1])
        elif (ijob == 2):
            if psolve is None:
                psolve = get_psolve(A)
            work[slice1] = psolve(work[slice2])
        elif (ijob == 3):
            if matvec is None:
                matvec = get_matvec(A)
            work[slice2] *= sclr2
            work[slice2] += sclr1*matvec(x)
        elif (ijob == 4):
            if ftflag:
                info = -1
                ftflag = False
            bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
        ijob = 2

    return x, info



More information about the SciPy-user mailing list