[SciPy-dev] SciPy improvements

Michael McNeil Forbes mforbes@physics.ubc...
Thu Apr 12 16:31:54 CDT 2007

```On 11 Apr 2007, at 1:56 PM, Matthieu Brucher wrote:

> 2) Nonlinear solvers
>
> I have written these nonlinear solvers for the problem R(x) = 0, where
> x and R has a dimension "n":
>
>     broyden1 - Broyden's first method - is a quasi-Newton-Raphson
> method for
>         updating an approximate Jacobian and then inverting it
>     broyden2 - Broyden's second method - the same as broyden1, but
>         inverse Jacobian directly
>     broyden3 - Broyden's second method - the same as broyden2, but
>         directly computing the inverse Jacobian, it remembers how
> to construct
>         it using vectors, and when computing inv(J)*F, it uses
> those vectors to
>         compute this product, thus avoding the expensive NxN matrix
>         multiplication.
>     broyden_generalized - Generalized Broyden's method, the same as
> broyden2,
>         but instead of approximating the full NxN Jacobian, it
> construct it at
>         every iteration in a way that avoids the NxN matrix
> multiplication.
>         This is not as precise as broyden3.
>     anderson - extended Anderson method, the same as the
> broyden_generalized,
>         but added w_0^2*I to before taking inversion to improve the
> stability
>     anderson2 - the Anderson method, the same as anderson, but
> formulated
>         differently
>     linear_mixing
>     exciting_mixing
>
> I use them in the self-consistent cycle of the Density Functional
> Theory (so I use a terminology of DFT literature in the names of the
> methods).
>
> Could the part that computes the step be separated from the
> function iself and the optimizer ? I'm trying to "modularize" non
> linear solvers so as to select more efficiently what is needed -
> kind of optimizer, kind of step, kind of stopping criterion, ... -
>
> Matthieu

It should be possible to modularize these with a step class that
maintains state
(the Jacobian, or its inverse etc.).

to look at it yet, but have been meaning to and would like to look at
Perhaps we should make a page on the Wiki to collect suggestions and
code samples.)

I have been meaning to get a good Broyden based algorithm coded for
python for a
while.  I have a MATLAB version of a globally convergent Broyden
implementation
using a linesearch as a base with a couple of unusual features that
might be
useful (not specific to Broyden based methods).  (The code is based
on the
presentation of the algorithm given in Numerical Recipies with some
modifications
suggested by Dennis and Schnabel's book and is partially documented
using noweb.)

http://www.phys.washington.edu/users/mforbes/projects/broyden/

1) Variable tolerances.  The idea is to quickly estimate the starting
Jacobian
with low tolerance calculations and then improve the tolerances
as the code
converges to the solution.  This is useful if the function R(x)
is computed
with numerical integration or some similar technique that is
quick for low
tolerances but expensive for high tolerance functions.  (I have
rarely seen
this technique in generic optimization algorithms, but found it
invaluable for
several problems.)
2) Real-time bailout.  This allows you to compute for a specified
length of time
and then return the best solution found within that time frame.
Most algorithms
simply count function calls.
3) Selective refinement of the Jacobian as the iteration proceeds.
This amounts
to monitoring the condition number of the Jacobian and
recomputing parts of it
selectively if it becomes ill-conditioned.  (For this I update
the QR
factorization of J rather than maintaining inv(J)).

These things all slow down the fundamental algorithm, but are very
useful when the
function R(x) is extremely expensive to evaluate.

Michael.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-dev/attachments/20070412/06caa228/attachment.html
```