[SciPy-user] fsolve with sparse matrices

[Jan Wicijowski]

Hi all!

I would like to ask you, if anyone was ever confronted with solving
nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
say 10000.

The usual approach fails, because the size of jacobian matrix is
proportional to the square of equations number. The paging kills my OS.

In my problem, which is solving a nonlinear electrical circuit, the matrix
is mostly zeros, which has immediately led me to trying sparse matrices.

The first obstacle was numpy.atleast_1d contained in
scipy.optimize.minpack.check_func - instead of preserving sparse matrix as
is, it has a suprising behaviour of boxing the matrix in a list. This may
be a bug, but I may post it elsewhere.

Injecting custom function instead of scipy.optimize.minpack.check_func
didn't work as well with fsolve. This wasn't surprising, as I guessed,
that FORTRAN hybrj won't be able to deal with interfacing with scipy
matrices.

So, am I left with rewriting the original FORTRAN hybrj source to python,
or is there somebody, who dealt with such problem?

Regards,
Jan Wicijowski

[Gergely Imreh]

Post by Jan Wicijowski
Hi all!
I would like to ask you, if anyone was ever confronted with solving
nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
say 10000.
The usual approach fails, because the size of jacobian matrix is
proportional to the square of equations number. The paging kills my OS.
In my problem, which is solving a nonlinear electrical circuit, the matrix
is mostly zeros, which has immediately led me to trying sparse matrices.
The first obstacle was numpy.atleast_1d contained in
scipy.optimize.minpack.check_func - instead of preserving sparse matrix as
is, it has a suprising behaviour of boxing the matrix in a list. This may
be a bug, but I may post it elsewhere.
Injecting custom function instead of scipy.optimize.minpack.check_func
didn't work as well with fsolve. This wasn't surprising, as I guessed,
that FORTRAN hybrj won't be able to deal with interfacing with scipy
matrices.
So, am I left with rewriting the original FORTRAN hybrj source to python,
or is there somebody, who dealt with such problem?
Regards,
Jan Wicijowski

Hi,

I was working on a similar problem (large sparse matrix, solving
optical Bloch-equations), and been planning to check out interfacing
with Fortran. The idea was, that I know the NAG library [1] works very
well for these kinds of calculations, just have to get the interface
between the sparse matrices and Fortran working. I checked it now, and
they have some documents showing how to use the NAG with F2PY [2].
Though they only use numpy, so no mention of sparse matrices, but NAG
does sparse.
This supposed to be just a heads up, need to check it out. Though in
my experience even the NAG people can be contacted, and they helped to
resolve issues (or extend functionality) before.

Cheers,
Greg

[1] http://www.nag.co.uk/
[2] http://www.nag.co.uk/doc/techrep/index.asp#np3665

[Jan Wicijowski]

Hi all!

I would like to ask you, if anyone was ever confronted with solving
nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
say 10000.

The usual approach fails, because the size of jacobian matrix is
proportional to the square of equations number. The paging kills my OS.

In my problem, which is solving a nonlinear electrical circuit, the matrix
is mostly zeros, which has immediately led me to trying sparse matrices.

The first obstacle was numpy.atleast_1d contained in
scipy.optimize.minpack.check_func - instead of preserving sparse matrix as
is, it has a suprising behaviour of boxing the matrix in a list. This may
be a bug, but I may post it elsewhere.

Injecting custom function instead of scipy.optimize.minpack.check_func
didn't work as well with fsolve. This wasn't surprising, as I guessed,
that FORTRAN hybrj won't be able to deal with interfacing with scipy
matrices.

So, am I left with rewriting the original FORTRAN hybrj source to python,
or is there somebody, who dealt with such problem?

Regards,
Jan Wicijowski

[Pauli Virtanen]

Fri, 17 Apr 2009 17:29:29 +0200, Jan Wicijowski kirjoitti:
[clip]

Post by Jan Wicijowski
I would like to ask you, if anyone was ever confronted with solving
nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
say 10000.

Sure, but as you noticed, fsolve won't cut it, as it assumes dense
matrices.

[clip: fsolve & sparse jacobian]

Post by Jan Wicijowski
Injecting custom function instead of scipy.optimize.minpack.check_func
didn't work as well with fsolve. This wasn't surprising, as I guessed,
that FORTRAN hybrj won't be able to deal with interfacing with scipy
matrices.
So, am I left with rewriting the original FORTRAN hybrj source to
python, or is there somebody, who dealt with such problem?

Translating hybrj directly to sparse matrices is a bit of work; it's a
trust-region Newton method, so it doesn't simply invert the Jacobian, but
solves a restricted linear programming problem involving it. (I think
some of the tricks it does work well only with sparse matrices.) In
principle, writing something like this with sparse matrices should
nevertheless be possible with the tools in Scipy (though Scipy does not
have a sparse QR decomposition).

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

***

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.

[David Huard]

Have you looked at pyamg ? It has a number of sparse matrix solvers based on
multigrid.

David

Post by Pauli Virtanen
[clip]

Post by Jan Wicijowski
I would like to ask you, if anyone was ever confronted with solving
nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
say 10000.

Sure, but as you noticed, fsolve won't cut it, as it assumes dense
matrices.
[clip: fsolve & sparse jacobian]

Post by Jan Wicijowski
Injecting custom function instead of scipy.optimize.minpack.check_func
didn't work as well with fsolve. This wasn't surprising, as I guessed,
that FORTRAN hybrj won't be able to deal with interfacing with scipy
matrices.
So, am I left with rewriting the original FORTRAN hybrj source to
python, or is there somebody, who dealt with such problem?

Translating hybrj directly to sparse matrices is a bit of work; it's a
trust-region Newton method, so it doesn't simply invert the Jacobian, but
solves a restricted linear programming problem involving it. (I think
some of the tricks it does work well only with sparse matrices.) In
principle, writing something like this with sparse matrices should
nevertheless be possible with the tools in Scipy (though Scipy does not
have a sparse QR decomposition).
If you write a sparse trust-region Newton algorithm in Python, I'm sure
there's a place in Scipy for it :)
***
The easier way would be just to implement a Newton method combined with
x = x0
maxiter = 100
abs_tolerance = 1e-8
J = jacobian(x)
r = residual(x)
d = -sparse.linalg.spsolve(J, r) # or, iterative
r_norm = norm(r)
# check convergence
break
# naive line search, could consider using
# scipy.optimize.line_search instead
x2 = x + alpha*d
r2 = residual(x2)
break
x = x2
raise RuntimeError("Didn't converge")
It's very naive, but if your problem is easy enough, it works.
--
Pauli Virtanen
_______________________________________________
SciPy-user mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user

[Jan Wicijowski]

Post by David Huard
Have you looked at pyamg ? It has a number of sparse matrix solvers based on
multigrid.
David

Just have taken a quick peek. Scipy also has got them, so when I find the
proper method in which I _could_ use sparse matrix solver I may consider
it.

Thanks

Jan Wicijowski

[Jan Wicijowski]

Post by Pauli Virtanen
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.

Post by Pauli Virtanen
The easier way would be just to implement a Newton method combined with
x = x0
maxiter = 100
abs_tolerance = 1e-8
J = jacobian(x)
r = residual(x)
d = -sparse.linalg.spsolve(J, r) # or, iterative
r_norm = norm(r)
# check convergence
break
# naive line search, could consider using
# scipy.optimize.line_search instead
x2 = x + alpha*d
r2 = residual(x2)
break
x = x2
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