Wednesday 24 June 2015

Progress Update: Midterm Evaluation


Objective:
Adding functions to the Optim package for Octave using existing back-ends.

Expected deliverables before midterm:
  • 'lsqnonlin' using 'nonlin_residmin'
Done in [1]. Differences in backends, nonlin_residmin uses "lm_svd" algorithm for optimization as currently the only backend. However, lsqnonlin in Matlab can choose between "trust-region-reflective" and "Levenberg-Marquardt" (LM) algorithms.
Another difference is in complex inputs. lm_svd does not support complex valued inputs whereas Matlab's LM algorithm can accept complex input arguments. One way of providing complex inputs to lsqnonlin in Octave is to split the real and imaginary parts into separate variables and running the optimization.  
  • 'lsqcurvefit' using 'nonlin_curvefit', 'nonlin_residmin', or 'lsqnonlin'
Done in [2] using nonlin_curvefit. lsqcurvefit is very similar to lsqnonlin with only a few minor interface differences. lsqcurvefit explicitly takes independent variables and the observations as inputs while these values can be wrapped inside the objective function while using lsqnonlin. Additional bounds for the optimized parameters can be specified. 
  • 'nlinfit' using 'leasqr',
I wrapped nlinfit on nonlin_curvefit and curvefit_stat as leasqr repeats the optimization to compute the additional statistics (Jacobian and Covariance matrix) while curvefit_stat saves this computation overhead. I have partially implemented nlinfit in [3] (It hasn't been thoroughly reviewed yet). Two missing features are: 1) Error models and Error parameters estimation, and 2) Robust Weight function. Meanwhile, no such functionality exists in the Octave's optimization backends for the missing features. My current implementation supports the input of scalar positive array of weights for robust regression.  
Since nlinfit is from the statistics toolbox of Matlab, it uses statset and statget to create and get options respectively. I created additional functions statset, statget and __all_stat_opts__ with minor changes to the code in optimset, optimget and __all_opts__.
  • 'fmincon' using 'nonlin_min',
In progress [4].


Future goals:
    1. Complete fmincon implementation.
    2. Create solver specific options using optimoptions and desirably still be able to use optimget. 
    3. Arranging lambda output for quadprog by wrapping it on __qp__ instead of qp.m
    4. Test cases for all the implemented functions.

[1] https://github.com/AsmaAfzal/octave_workspace/blob/master/lsqnonlin.m
[2] https://github.com/AsmaAfzal/octave_workspace/blob/master/lsqcurvefit.m
[3] https://github.com/AsmaAfzal/octave_workspace/blob/master/nlinfit.m
[4] https://github.com/AsmaAfzal/octave_workspace/blob/master/fmincon.m

Monday 22 June 2015

Week4: fmincon wrapping nonlin_min

Time flies.. A third of the way through already..

fmincon/nonlin_min is the most elaborate function of all that I have previously 
implemented so before actual coding I would like to thoroughly check the the 
mapping of arguments and options. 

[x,fval,exitflag,output,lambda,grad,hessian] = fmincon (fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
[x,fval,exitflag,output] = nonlin_min (fun,x0,settings)

A,b Aeq, beq- Linear inequality and equality constraints
Inequality constraints:
Matlab standard: A * x - b <= 0 ,
Octave standard: m.' * x  + v >= 0. 
This implies: m=-A.' , v=b.
Set in Octave using 
optimset ("inequc", {m,v})
Similar for equality constraints:
optimset ("equc", {m,v})

lb, ub- Lower and upper bounds
Set in Octave using optimset
 optimset ("lbound", ..., "ubound",...)

nonlcon- Nonliner constraint function handle 
In Matlab, nonlinear constraints are given in a function with the following format

function [c,ceq,gradc,gradceq] = mycon(x)
c = ...     % Compute nonlinear inequalities at x.  
ceq = ...   % Compute nonlinear equalities at x.
% Optional output arguments:
gradc = ...   % Gradient of c(x).
gradceq = ...   % Gradient of ceq(x).
options = optimoptions('fmincon','GradConstr','on')

Alternative to nonlcon in nonlin_min
  optimset ("equc", {constraint_function})

Options- Options common to both fmincon and nonlin_min
User-supplied Gradient
In Matlab, the objective function must return the second output when the 
GradObj option is set:
optimset('GradObj','on')
Octave:
optimset("objf_grad",@fun)
User-supplied Hessian
In Matlab, the objective function must return the third output when the Hessian 
option is set:
optimset('Hessian','user-supplied')
Octave:
optimset("objf_hessian",@fun)

Lambda- returned as a structure in Matlab but as a vector field of the "output" 
ouput argument in Octave.

Things to investigate:
1. Algorithm mapping. 
2. Exitflag mapping.
3. Setting gradc (the gradient of general equality/inequality functions). The second 
entry for the equc/inequc setting implements this feature in Octave as stated by 
optim_doc but I was unable to make it work properly.
4. Returning additional statistics (Hessian, Gradient)- I used 
residmin_stat/curvefit_stat previously for lsqnonlin and lsqcurvefit. No such function for nonlin_min.
5. Rearranging values of lambda in the fields of a structure.

P.s. this week I spent quite some time with Mercurial and how to possibly work in the optim package. There was a slight miscommunication/confusion but it is clear now and I will continue to publish my code on github..




Sunday 14 June 2015

Week3: nlinfit, statset, stat_get and __all_stat_opts__

So this week I achieved the following milestones:

Wrapping nlinfit on nonlin_curvefit

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,Y,modelfun,beta0,options,Name,Value)

Implementation [1]:

I chose not to wrap nlinfit on lsqcurvefit because 
  1. We might end up wrapping lsqcurvefit on lsqnonlin eventually so it is undecided.
  2. The default options for lsqnonlin/lsqcurvefit are different from nlinfit. 
Missing features:
  1. RobustWgtFun - The field RobustWgtFun in options can be provided with a function handle which computes the residuals at every iteration. The backend optimization algorithm in Octave currently does not support this functionality.  
  2. Name-Value pairs. Currently the only implemented one is "weights",  which takes an array of weights for the weighted optimization. "ErrorModelInfo" and "ErrorParameters" are not implemented. The possible error models include, "constant", "proportional" and "combined". The error model also translates to a weight function which helps in reducing the effect of outliers.   
  3. ErrorModelInfo- output field which gives information about the error variance, and estimates the parameters of Error models. 
Setting options using statset

In Matlab, options for statistics toolbox are set using statset [2].  The functionality is almost identical to optimset but the different functions in Matlab are because of different toolboxes (statset for statistics and optimset for optimization toolbox).

Added functions:
statset.m, statget.m and __all_stat_opts__.m [3]-[5]
Creating these functions was pretty straight forward.

Still to do:
  1. Have to check if the weighted residual and weighted Jacobian output in octave is consistent with Matlab and further refine the functions with the feedback from my mentors.
  2. Move on to fmincon wrapping nonlin_min.
[1] https://github.com/AsmaAfzal/octave_workspace/blob/master/nlinfit.m
[2] http://uk.mathworks.com/help/stats/statset.html
[3] https://github.com/AsmaAfzal/octave_workspace/blob/master/statset.m
[4] https://github.com/AsmaAfzal/octave_workspace/blob/master/statget.m
[5] https://github.com/AsmaAfzal/octave_workspace/blob/master/__all_stat_opts__.m

Tuesday 9 June 2015

Week 2: lsqnonlin and lsqcurvefit

A bit late blogging about week 2. 

Almost completed functions lsqnonlin and lsqcurvefit. 
Successfully mapped user-specified Jacobian. 
In Matlab, if the Jacobian option is set to "on",  the model function must return 
a second output which is the Jacobian function.
In Octave, the Jacobian function handle is given to the dfdp option using optimset.

Lsqnonlin function description:

[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(fun,x0,lb,ub,options)

This function maps on:

[x, residual, exitflag, output] = nonlin_residmin (fun, x0, options)

Features of lsqnonlin in Octave:
  • Input arguments: Acceptable forms are lsqnonlin(fun,x0), lsqnonlin(fun,x0,lb,ub) and lsqnonlin(fun,x0,lb,ub,options)
  • Outputs
    • x, exitflag, residual and output currently same as nonlin_residmin.
    • resnorm=sum(residual.^2)
    • Lambda is computed using the complementary pivoting in __lm_svd__.
      It's values differ from Matlab's due to the difference in backends.
    • Jacobian is computed using the function residmin_stat ().


Lsqcurvefit function description:

[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

This function maps on

[x, fy, exitflag, output] = nonlin_curvefit (fun, x0, xdata, ydata, options)

Features of lsqcurvefit in Octave:
  • Input arguments: Acceptable forms are lsqcurvefit (fun,x0,xdata,ydata), lsqcurvefit (fun,x0,xdata,ydata,lb,ub) and lsqcurvefit (fun,x0,lb,ub,xdata,ydata,options)
  • Outputs
    • x, exitflag, residual and output currently same as nonlin_curvefit.
    • residual = fy-ydata, resnorm = sum(residual.^2)
    • Lambda and Jacobian same as in lsqnonlin.
There are only minor interface differences between lsqcurvefit and lsqnonlin.

This week's plan:

  • Hopefully, with lsqnonlin and lsqcurvefit wrapped up, I'll move on to nlinfit
  • Three key challenges need to be addressed when wrapping nlinfit using nonlin_curvefit
    and curvefit_stat:
    • Weight functions: Currently, no such functionality exists in nonlin_curvefit,
       where a user can specify weight functions to perform Robust regression
      (weights computed using the specified function in every iteration).
    • Error Models and ErrorModelInfo
    • Setting options using statset instead of optimset or optimoptions.