Re: fmincon wrapping nonlin_min

[Olaf Till]

On Fri, Jun 26, 2015 at 08:56:39PM +0100, Asma Afzal wrote:
> Hi Olaf,
> 
> I need some help in wrapping fmincon to nonlin_min:
> 
> - User supplied gradient of nonlinear constraints.
> Consider the example:
> 
> objective_function = @ (p) p(1)^2 + p(2)^2;
>  pin = [-2; 5];
>  constraint_function = @ (p) p(1)^2 + 1 - p(2);
>  gradc = @(p) [2*p(1);-1];
>  [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset
> ("equc", {constraint_function, gradc}))
> 
> It gives the following error:
> error: function handle type invalid as index value

A bug in nonlin_min.m (and __nonlin_residmin__.m). Should be fixed now
in central repository.

> What is the right way pass gradc to nonlin_min?

It was the right way, with the important exception that the gradient
must be a row vector. (This is consistant with the Jacobians in
residual estimation.)

> - Another simple example:
> 
> f = @(x) -x(1) * x(2) * x(3);
> S = [1  -1;   2  -2;   2  -2]
> b = [0;72];
> x0 = [10;10;10];
> [x,fval] = nonlin_min( f, x0, optimset ("inequc",{S,b}) )
> 
> Gives an error with lm_feasible
> error: __lm_feasible__: operator -: nonconformant arguments (op1 is
> 3x1, op2 is 3x0)
> __lm_feasible__.m at line 208, column 11

A bug in __lm_feasible__.m. Should be fixed now in central repository.

> but it works fine with octave_sqp.
> 
> - Output Gradient and Hessian at the solution as in Matlab.

The situation is a bit different here than for
nonlin_residmin/residmin_stat, since the Hessian approximation is (by
default) refined over the course of the whole optimization and can not
by default be computed without repeating the optimization. And if we
compute the Hessian for the final parameters, by default we also have
to compute the gradient for the final parameters, so returning the
gradient would not be an additional effort. And, of course, we
currently have no equivalent of residmin_stat for scalar optimization.

On the other hand, the algorithm doesn't need the Hessian refined for
the final parameters, so doing this final refinement _is_ an
additional effort.

What I currently think is we should introduce the options
'ret_objf_grad' and 'ret_hessian' to nonlin_min, and compute a
gradient and Hessian for the final parameters if set. They could be
returned as .objf_grad and .hessian fields of the 'outp' structure
(fourth output argument of nonlin_min). And I'd say both or nothing
should be returned, regardless whether ret_hessian or ret_objf_grad or
both had been set. As long as I've not implemented this, you could
just set these options and check if the output fields are present.

> - Output Lambda as a structure: changing backend to octave_sqp when
> this output is requested? but the ordering issue again comes down to
> the lambda returned from qp.

I probably should make __lm_feasible__ return the lambda the same way
as __lm_svd__ does it now. It's a bit less nonsensical (?) than in
__lm_svd__, since the Hessian approximation does not inherently wildly
deviate from the real Hessian. Still if you should now what a user can
do with the returned lambda, please inform me. For now, you could test
for the field lambda.

I probably can't do the same for __octave_sqp_wrapper__, do to the
already discussed order problem. Maybe I should disable outputting
lambda in __octave_sqp_wrapper__.

> My code for fmincon:
> 
> https://github.com/AsmaAfzal/octave_workspace/blob/master/fmincon.m

My impression is that you have a good understanding of the general
principals of these wrappings now, so that my advice is rather needed
in small details, or if something strategical has to be decided. This
is good, because it would get rather much for me if I had to proofread
every line. But it would be good if you would provide some examples
for easy testing of your code.

Some of these small details:

- You have a mistake in some functions (lsqnonlin, fmincon, ...?)
  which always sets the default for TypicalX, never a user-defined
  value:

  e.g lsqnonlin, that'd be a correction:

147a148
>       TypicalX = optimget (settings, "TypicalX", TypicalX_default);
152c153
<                            "TypicalX", TypicalX_default,
---
>                            "TypicalX", TypicalX,

- In my recent remark on the line-continuation marks, I seemingly had
  some 'thinko'. You don't need them in
  'optimset ("FinDiffRelStep", [],...'

- I'd use:

  [~, ~, GC] = nonlcon (p);

  and not

  [~, ~, GC, ~] = nonlcon (p);

  because nonlcon might test for 'nargout () == 4' instead of
  'isargout (4)', possibly leading to some needless computations.

Olaf

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc
Description: Digital signature