I had chosen the classical half-step approach because it is robust and easy to implement (I thought I used Richardson extrapolation? Then it should also be second order). Anyway, I guess that the difference is negligible and we could go for odebwenk here without any loss.

I propose odeset('nlssolver',@function) as an option for odebwe with fallback to a private newton scheme. Then we don't need a dependency on optim because it's the user's responsibility.

Checking for the package is also fine (I personall check for functions because I want to keep matlab compatibility).

> <snip>
> Yes it is meant for problems with large, sparse jacobians
> <snip>
> indeed it is a zero-finder but many ptimization functions could made
> to optionally use it when dealing with large problems.

Principally yes, but this probably has some issues:

- It would be an algorithm without constraints (but this could be
o.k. due to its specialization).

- The user could not, as usual, specify a Hessian function of the
parameters directly, but would have to do it indirectly for
this algorithm by specifying a jacobian (for inexact_newton)
of the gradient function (of the outer algorithm).

- A user-specified (in which way whatever) Hessian would not
necessarily be positive definite. This can spoil the outer
algorithm which uses inexact_newton.

- A user-specified Hessian isn't even necessarily non-singular.
I'm not sure that inexact_newton would cope with that.

> <snip>
> so you suggest adding inexact_newton both in odepkg and optim? I
> think it could be done if we make it private in odepkg to avoid
> conflicts ...

Yes, something like that. I don't suggest it, but it would be a way if you value independence. Which Sebastian seems to do. I think this way would also be compatible with his suggestion to make inexact_newton rather a configurable sub-algorithm of an existing solver.

Anyway I think it could be good if this function were also in 'optim'.

Re-posted for Carlo de Falco <kingcrimson@tiscali.it>
(Actually the order of this and the previous post was reversed, sorry.):
---------------------------------------------------------------

On 13 May 2013, at 22:10, Olaf Till <INVALID.NOREPLY@gnu.org> wrote:

> Follow-up Comment #1, patch #8047 (project octave):
>
> As I understand, the inexact Newton method, as well as the ODE solver which
> uses it, is meant for large scale problems --- is that right?

Yes it is meant for problems with large, sparse jacobians

> If specially suitable for large scale problems, I think 'inexact_newton' could
> be a valuable alternative to 'fsolve', and so would fit well into the 'optim'
> package (I believe it would be the first zero-finder there, except a
> vectorized clone of 'fzero').

indeed it is a zero-finder but many ptimization functions could made to optionally
use it when dealing with large problems. I'd like the function to be given more testing
before such changes are done, though.

> The dependency on other packages can sometimes be a problem, though. 'optim'
> sometimes makes problems on some systems due to its use of LAPACK, and it in
> turn depends on further packages, partially due to some older code. Maybe
> eventually some reorganization will happen; until that, duplicating the code
> of 'inexact_newton' in the 'odepkg' package could be an (ugly) solution for
> being independent …

so you suggest adding inexact_newton both in odepkg and optim?
I think it could be done if we make it private in odepkg to avoid conflicts ...

> Or 'inexact_newton' might actually go into Octave core, if it is really good
> for large scale problems.

I don't think it fits in Octave core, its scope does not seem general enough to me
and there is not an equivalent function in Matlab

> But that is not mine to decide. Its interface would
> probably have to be changed for that in favor of 'optimset' (and maybe this
> would even be better anyway).

The backward Euler that I have written (odebwe.m) uses a private Newton scheme to avoid dependencies. I would encourage to merge the two integrators.

Instead there should be a switch in odeset to choose the nonlinear solver (fixed point, Newton, simplified Newton, Krylov) as suggested by Olaf. I would prefer
+to keep the classical Newton build-in and the others could be made optional with something like "if !exist(...,'file') warning('Install Optim')".

Best regards,
Sebastian

Am 14.05.2013 um 08:25 schrieb "c." <kingcrimson@tiscali.it>:

>
> On 13 May 2013, at 22:10, Olaf Till <INVALID.NOREPLY@gnu.org> wrote:
>
>> Follow-up Comment #1, patch #8047 (project octave):
>>
>> As I understand, the inexact Newton method, as well as the ODE solver which
>> uses it, is meant for large scale problems --- is that right?
>
> Yes it is meant for problems with large, sparse jacobians
>
>> If specially suitable for large scale problems, I think 'inexact_newton' could
>> be a valuable alternative to 'fsolve', and so would fit well into the 'optim'
>> package (I believe it would be the first zero-finder there, except a
>> vectorized clone of 'fzero').
>
> indeed it is a zero-finder but many ptimization functions could made to optionally
> use it when dealing with large problems. I'd like the function to be given more testing
> before such changes are done, though.
>
>> The dependency on other packages can sometimes be a problem, though. 'optim'
>> sometimes makes problems on some systems due to its use of LAPACK, and it in
>> turn depends on further packages, partially due to some older code. Maybe
>> eventually some reorganization will happen; until that, duplicating the code
>> of 'inexact_newton' in the 'odepkg' package could be an (ugly) solution for
>> being independent …
>
> so you suggest adding inexact_newton both in odepkg and optim?
> I think it could be done if we make it private in odepkg to avoid conflicts ...
>
>> Or 'inexact_newton' might actually go into Octave core, if it is really good
>> for large scale problems.
>
> I don't think it fits in Octave core, its scope does not seem general enough to me
> and there is not an equivalent function in Matlab
>
>> But that is not mine to decide. Its interface would
>> probably have to be changed for that in favor of 'optimset' (and maybe this
>> would even be better anyway).
>
> that is a good idea.
>
>> Olaf
> c.

As I understand, the inexact Newton method, as well as the ODE solver which uses it, is meant for large scale problems --- is that right?

If specially suitable for large scale problems, I think 'inexact_newton' could be a valuable alternative to 'fsolve', and so would fit well into the 'optim' package (I believe it would be the first zero-finder there, except a vectorized clone of 'fzero').

The dependency on other packages can sometimes be a problem, though. 'optim' sometimes makes problems on some systems due to its use of LAPACK, and it in turn depends on further packages, partially due to some older code. Maybe eventually some reorganization will happen; until that, duplicating the code of 'inexact_newton' in the 'odepkg' package could be an (ugly) solution for being independent ...

Or 'inexact_newton' might actually go into Octave core, if it is really good for large scale problems. But that is not mine to decide. Its interface would probably have to be changed for that in favor of 'optimset' (and maybe this would even be better anyway).

Please find attached a new submission I'd like to be included
in ODEPKG it implements an adaptive timestep Backward-Euler solver that uses a Newton-Krylov solver for non-linear system solution, implemented by Simone Panza and Fabio Cisaria.

The main reason for posting this contribution here rather than commiting to the subversion repository directly is that I need help deciding what should be done with the "inexact_newton" function.

Personally I believe the inexact_newton function has use beyond the scope of ODEPKG so maybe it should rather go into the optim package.

This though would impose a dependence of odepkg on optim, which is something that the previous maintainer of odepkg has always tried to avoid (odepkg currently has no additional dependencies). Therefore I'd like you to comment on where this function fits better before committing.