[SciPy-user] fsolve with sparse matrices

[Jan Wicijowski]

Dnia 17-04-2009 o 18:28:14 Pauli Virtanen <pav at iki.fi> napisał(a):

> If you write a sparse trust-region Newton algorithm in Python, I'm sure
> there's a place in Scipy for it :)

In my work I immediately picked ready solution of scipy.optimize.fsolve,  
and didn't look at simpler methods. Fortran hybrj algorithm works  
blazingly fast for my problem, and I expect simpler methods to converge  
much slower. Although it may be a reasonable trade-off with sparse  
matrices.

> The easier way would be just to implement a Newton method combined with
> line search:
>
> 	x = x0
> 	maxiter = 100
> 	abs_tolerance = 1e-8
> 	for k in xrange(maxiter):
> 	    J = jacobian(x)
> 	    r = residual(x)
> 	    d = -sparse.linalg.spsolve(J, r)  # or, iterative
>
> 	    r_norm = norm(r)
>
> 	    # check convergence
> 	    if r_norm < abs_tolerance:
> 	        break
>
> 	    # naive line search, could consider using
> 	    # scipy.optimize.line_search instead
> 	    for alpha in [1, 0.5, 0.1, 0.01]:
> 		x2 = x + alpha*d
> 	        r2 = residual(x2)
> 	        if norm(r2) < r_norm:
> 	            break
> 	    x = x2
> 	else:
> 	    raise RuntimeError("Didn't converge")
>
> It's very naive, but if your problem is easy enough, it works.

Cool, it looks nice and self-explanatory. I'll try it ASAP.

If the results for my problem won't satisfy me, expect to see hybrj  
algorithm ported to python in two weeks or so ;)

Thanks for your help - you've pointed out many important facts for my work.

Regards,
Jan Wicijowski

-- 
Jan Wicijowski
a29jaGFtIEp1bGnqIEcu