TOMLAB Solving Least Squares and Parameter Estimation Problems: Difference between revisions
No edit summary |
|||
(10 intermediate revisions by 2 users not shown) | |||
Line 1: | Line 1: | ||
{{Part Of Manual|title=the TOMLAB Manual|link=[[TOMLAB|TOMLAB Manual]]}} | |||
[[Category:Manuals]] | [[Category:Manuals]] | ||
This section describes how to define and solve different types of linear and nonlinear least squares and parameter estimation problems. Several examples are given on how to proceed, depending on if a quick solution is wanted, or more advanced tests are needed. TOMLAB is also compatible with MathWorks Optimization TB. See [[TOMLAB Appendix E]] for more information and test examples. | |||
All demonstration examples that are using the TOMLAB format are collected in the directory ''examples''. The examples relevant to this section are ''lsDemo ''and ''llsDemo''. The full path to these files are always given in the text. | |||
[[TOMLAB_Solving_Least_Squares_and_Parameter_Estimation_Problems#Large_Scale_LS_problems_with_Tlsqr|Large Scale LS problems with Tlsqr]] contains information on solving extreme large-scale '''ls '''problems with ''Tlsqr''. | |||
== | ==Linear Least Squares Problems== | ||
This section | This section shows examples how to define and solve linear least squares problems using the TOMLAB format. As a first illustration, the example ''lls1Demo ''in file ''llsDemo ''shows how to fit a linear least squares model with linear constraints to given data. This test problem is taken from the User's Guide of ''LSSOL''. | ||
<source lang="matlab"> | |||
Name='LSSOL test example'; | |||
% In TOMLAB it is best to use Inf and -Inf, not big numbers. | |||
n = 9; % Number of unknown parameters | |||
% In TOMLAB it is best to use Inf and -Inf, not | |||
n = 9; % Number | |||
x_L = [-2 -2 -Inf, -2*ones(1,6)]'; | x_L = [-2 -2 -Inf, -2*ones(1,6)]'; | ||
x_U = 2*ones(n,1); | x_U = 2*ones(n,1); | ||
Line 31: | Line 24: | ||
y = ones(10,1); | y = ones(10,1); | ||
C = | C = [ ones(1,n); 1 2 1 1 1 1 2 0 0; 1 1 3 1 1 1 -1 -1 -3; ... | ||
1 1 1 4 1 1 1 1 1;1 1 1 3 1 1 1 1 1;1 1 2 1 1 0 0 0 -1; ... | 1 1 1 4 1 1 1 1 1;1 1 1 3 1 1 1 1 1; 1 1 2 1 1 0 0 0 -1; ... | ||
1 1 1 1 0 1 1 1 1;1 1 1 0 1 1 1 1 1;1 1 0 1 1 1 2 2 3; ... | 1 1 1 1 0 1 1 1 1;1 1 1 0 1 1 1 1 1; 1 1 0 1 1 1 2 2 3; ... | ||
1 0 1 1 1 1 0 2 2 | 1 0 1 1 1 1 0 2 2]; | ||
x_0 = 1./[1:n]'; | x_0 = 1./[1:n]'; | ||
t = []; % No time set for y(t) (used | t = []; % No time set for y(t) (used for plotting) | ||
weightY = []; % No weighting | weightY = []; % No weighting | ||
weightType = []; % No weighting type | weightType = []; % No weighting type set | ||
x_min = []; % No lower bound for plotting | x_min = []; % No lower bound for plotting | ||
Line 47: | Line 40: | ||
Prob = llsAssign(C, y, x_L, x_U, Name, x_0, t, weightType, weightY, ... | Prob = llsAssign(C, y, x_L, x_U, Name, x_0, t, weightType, weightY, ... | ||
A, b_L, b_U, x_min, x_max); Result = tomRun('lsei',Prob,2); | A, b_L, b_U, x_min, x_max); Result = tomRun('lsei',Prob,2); | ||
</ | </source> | ||
It is trivial to change the solver in the call to ''tomRun '' to a nonlinear least squares solver, e.g. ''clsSolve'', or a general nonlinear programming solver. | It is trivial to change the solver in the call to ''tomRun'' to a nonlinear least squares solver, e.g. ''clsSolve'', or a general nonlinear programming solver. | ||
==Linear Least Squares Problems using the SOL Solver LSSOL== | |||
The example ''lls2Demo ''in file ''llsDemo ''shows how to fit a linear least squares model with linear constraints to given data using a direct call to the SOL solver ''LSSOL''. The test problem is taken from the | The example ''lls2Demo'' in file ''llsDemo'' shows how to fit a linear least squares model with linear constraints to given data using a direct call to the SOL solver ''LSSOL''. The test problem is taken from the User's Guide of ''LSSOL''. | ||
< | <source lang="matlab"> | ||
% Note that when calling the | % Note that when calling the LSSOL MEX interface directly, avoid using | ||
% Inf and -Inf. Instead use big | % Inf and -Inf. Instead use big numbers that indicate Inf. | ||
% The standard | % The standard for the MEX interfaces is 1E20 and -1E20, respectively. | ||
n = 9; % There are | n = 9; % There are nine unknown parameters, and 10 equations | ||
x_L = [-2 -2 -1E20, -2*ones(1,6)]'; | x_L = [-2 -2 -1E20, -2*ones(1,6)]'; | ||
x_U = 2*ones(n,1); | x_U = 2*ones(n,1); | ||
Line 68: | Line 61: | ||
b_U = [1E20 -2 -2]'; | b_U = [1E20 -2 -2]'; | ||
% Must put lower and upper bounds on variables and constraints together | % Must put lower and upper bounds on variables and constraints together | ||
bl = | bl = [x_L;b_L]; | ||
bu = | bu = [x_U;b_U]; | ||
H = [ ones(1,n); 1 2 1 1 1 1 2 0 0; 1 1 3 1 1 1 -1 -1 -3; ... | H = [ ones(1,n); 1 2 1 1 1 1 2 0 0; 1 1 3 1 1 1 -1 -1 -3; ... | ||
Line 79: | Line 72: | ||
x_0 = 1./[1:n]'; | x_0 = 1./[1:n]'; | ||
% Set | % Set empty indicating default values for most variables | ||
c = []; % No linear coefficients, they are for LP/QP | c = []; % No linear coefficients, they are for LP/QP | ||
Warm = []; % No | Warm = []; % No warm start | ||
iState = []; % No | iState = []; % No warm start | ||
Upper = []; % C | Upper = []; % C is not factorized | ||
kx = []; % No | kx = []; % No warm start | ||
SpecsFile = []; % No | SpecsFile = []; % No parameter settings in a SPECS file | ||
PriLev = []; % PriLev | PriLev = []; % PriLev is not really used in LSSOL | ||
ProbName = []; % ProbName | ProbName = []; % ProbName is not really used in LSSOL | ||
optPar(1) = 50; % Set | optPar(1) = 50; % Set print level at maximum | ||
PrintFile = 'lssol.txt'; % Print result on the | PrintFile = 'lssol.txt'; % Print result on the file with name lssol.txt | ||
z0 = (y-H*x_0); | z0 = (y-H*x_0); | ||
Line 106: | Line 99: | ||
f = 0.5*z'*z; | f = 0.5*z'*z; | ||
fprintf('Optimal function value %f\n',f); | fprintf('Optimal function value %f\n',f); | ||
</ | </source> | ||
==Nonlinear Least Squares Problems== | |||
This section shows examples how to define and solve nonlinear least squares problems using the TOMLAB format. As a first illustration, the example ''ls1Demo ''in file ''lsDemo ''shows how to fit a nonlinear model of exponential type with three unknown parameters to experimental data. This problem, ''Gisela'', is also defined as problem three in '' | This section shows examples how to define and solve nonlinear least squares problems using the TOMLAB format. As a first illustration, the example ''ls1Demo ''in file ''lsDemo ''shows how to fit a nonlinear model of exponential type with three unknown parameters to experimental data. This problem, ''Gisela'', is also defined as problem three in ''ls_prob''. A weighting parameter ''K ''is sent to the residual and Jacobian routine using the ''Prob ''structure. The solver ''clsSolve ''is called directly. Note that the user only defines the routine to compute the residual vector and the Jacobian matrix of derivatives. TOMLAB has special routines ''ls f'', ''ls g ''and ''ls H ''that computes the nonlinear least squares objective function value, given the residuals, as well as the gradient and the approximative Hessian, see [[#Table: Constrained nonlinear least squares (cls) test problems.]]. The residual routine for this problem is defined in file ''ls1_r ''in the directory ''example'' with the statements | ||
< | <source lang="matlab"> | ||
function r = ls_r(x, Prob) | function r = ls_r(x, Prob) | ||
% Compute residuals to nonlinear least squares problem Gisela | % Compute residuals to nonlinear least squares problem Gisela | ||
% US_A is the | % US_A is the standard TOMLAB global parameter to be used in the | ||
% communication | % communication between the residual and the Jacobian routine | ||
global US_A | global US_A | ||
% The extra weight parameter | % The extra weight parameter K is sent as part of the structure | ||
K = Prob.user.K; | K = Prob.user.K; | ||
t = Prob.LS.t(:); % Pick up the time points | t = Prob.LS.t(:); % Pick up the time points | ||
% Exponential | % Exponential expressions to be later used when computing the Jacobian | ||
US_A.e1 = exp(-x(1) | US_A.e1 = exp(-x(1)*t); US_A.e2 = exp(-x(2)*t); | ||
r = K | r = K*x(1)*(US_A.e2 - US_A.e1) / (x(3)*(x(1)-x(2))) - Prob.LS.y; | ||
</ | </source> | ||
Note that this example also shows how to communicate information between the residual and the Jacobian routine. It is best to use any of the predefined global variables '' | Note that this example also shows how to communicate information between the residual and the Jacobian routine. It is best to use any of the predefined global variables ''US_A'' and ''US_B'', because then there will be no conflicts with respect to global variables if recursive calls are used. In this example the global variable ''US_A'' is used as structure array storing two vectors with exponential expressions. The Jacobian routine for this problem is defined in the file ''ls1_J '' in the directory ''example''. The global variable ''US_A'' is accessed to obtain the exponential expressions, see the statements below. | ||
< | <source lang="matlab"> | ||
function J = ls1_J(x, Prob) | function J = ls1_J(x, Prob) | ||
% Computes | % Computes the Jacobian to least squares problem Gisela. J(i,j) is dr_i/d_x_j | ||
% Parameter K is input in the structure | % Parameter K is input in the structure | ||
Prob a = Prob.user.K * x(1)/(x(3)*(x(1)-x(2))); | Prob a = Prob.user.K * x(1)/(x(3)*(x(1)-x(2))); | ||
Line 147: | Line 140: | ||
e1 = US_A.e1; e2 = US_A.e2; | e1 = US_A.e1; e2 = US_A.e2; | ||
% Compute | % Compute the three columns in the Jacobian, one for each of variable | ||
J = a * [ t.*e1+(e2-e1)*(1-1/b), -t.*e2+(e2-e1)/b, (e1-e2)/x(3) ]; | J = a * [ t.*e1+(e2-e1)*(1-1/b), -t.*e2+(e2-e1)/b, (e1-e2)/x(3) ]; | ||
Line 165: | Line 158: | ||
19; 17; 14; 9.5; 8.5; 7; 6; 6; 4.5; 3.6; 3; 2.2; 1.6]; | 19; 17; 14; 9.5; 8.5; 7; 6; 6; 4.5; 3.6; 3; 2.2; 1.6]; | ||
x_0 = | x_0 = [6.8729,0.0108,0.1248]'; % Initial values for unknown x | ||
% Generate the | % Generate the problem structure using the TOMLAB format (short call) | ||
% Prob = clsAssign(r, J, JacPattern, x_L, | % Prob = clsAssign(r, J, JacPattern, x_L, x_U, Name, x_0, ... | ||
% y, t, weightType, | % y, t, weightType, weightY, SepAlg, fLowBnd, ... | ||
% A, b_L, b_U, c, dc, | % A, b_L, b_U, c, dc, ConsPattern, c_L, c_U, ... | ||
% x_min, x_max, f_opt, x_opt); | % x_min, x_max, f_opt, x_opt); | ||
Prob = clsAssign('ls1_r', 'ls1_J', | Prob = clsAssign('ls1_r', 'ls1_J', [], [], [], Name, x_0, y, t); | ||
% Weighting parameter | % Weighting parameter K in model is sent to r and J computation using Prob | ||
Prob.user.K = 5; | Prob.user.K = 5; | ||
Result = tomRun('clsSolve', Prob, 2); | Result = tomRun('clsSolve', Prob, 2); | ||
</ | </source> | ||
The second example ''ls2Demo ''in file ''lsDemo ''solves the same problem as ''ls1Demo'', but using numerical differences to compute the Jacobian matrix in each iteration. To make TOMLAB avoid using the Jacobian routine, the variable ''Prob.NumDiff ''has to be set nonzero. Also in this example the flag ''Prob.optParam.IterPrint ''is set to enable one line of printing for each iteration. The changed statements are | The second example ''ls2Demo ''in file ''lsDemo ''solves the same problem as ''ls1Demo'', but using numerical differences to compute the Jacobian matrix in each iteration. To make TOMLAB avoid using the Jacobian routine, the variable ''Prob.NumDiff ''has to be set nonzero. Also in this example the flag ''Prob.optParam.IterPrint ''is set to enable one line of printing for each iteration. The changed statements are | ||
< | <source lang="matlab"> | ||
... | ... | ||
Prob.NumDiff = 1; % Use standard numerical differences | Prob.NumDiff = 1; % Use standard numerical differences | ||
Line 189: | Line 182: | ||
Result = tomRun('clsSolve',Prob,2); | Result = tomRun('clsSolve',Prob,2); | ||
</ | </source> | ||
The third example ''ls3Demo ''in file ''lsDemo ''solves the same problem as ''ls1Demo'', but six times for different values of the parameter ''K ''in the range | The third example ''ls3Demo ''in file ''lsDemo ''solves the same problem as ''ls1Demo'', but six times for different values of the parameter ''K ''in the range [3''.''8'', ''5''.''0]. It illustrates that it is not necessary to remake the problem structure ''Prob ''for each optimization, but instead just change the parameters needed. The ''Result ''structure is saved as an vector of structure arrays, to enable post analysis of the results. The changed statements are | ||
< | <source lang="matlab"> | ||
for i=1:6 | for i=1:6 | ||
Prob.user.K = 3.8 + 0.2 | Prob.user.K = 3.8 + 0.2*i; | ||
Result(i) = tomRun('clsSolve',Prob,2); | Result(i) = tomRun('clsSolve',Prob,2); | ||
end | end | ||
fprintf('\nWEIGHT PARAMETER K is %9.3f\n\n\n',Prob.user.K); | fprintf('\nWEIGHT PARAMETER K is %9.3f\n\n\n',Prob.user.K); | ||
</ | </source> | ||
Table | [[#Table: Constrained nonlinear least squares (cls) test problems.]] describes the low level routines and the initialization routines needed for the predefined constrained nonlinear least squares ('''cls''') test problems. Similar routines are needed for the nonlinear least squares ('''ls''') test problems (here no constraint routines are needed). | ||
====Table: Constrained nonlinear least squares (cls) test problems.==== | |||
{| class="wikitable" style="background-color:#ffffff;" | |||
!Function | |||
!Description | |||
|- | |- | ||
|''cls_prob''||Initialization of '''cls '''test problems. | |||
|'' | |||
|- | |- | ||
|''cls_r''||Compute the residual vector <math>r_i(x), i = 1,...,m.~ x \in \mathbb{R}^n</math> ''for '''cls ''' test problems. | |||
|'' | |||
|- | |- | ||
|''cls_J''||Compute the Jacobian matrix <math>J_{ij}(x)=\partial r_i / d x_j, | |||
|'' | i=1,...,m, j=1,...,n</math> ''for '''cls ''' test problems. | ||
|- | |- | ||
|''cls_c''||Compute the vector of constraint functions ''c''(''x'') for '''cls '''test problems. | |||
|'' | |||
|- | |- | ||
|''cls_dc''||Compute the matrix of constraint normals <math>\partial | |||
|'' | c(x)/dx</math> ''for '''cls '''test problems. | ||
|- | |- | ||
|''cls_d2c''||Compute the second part of the second derivative of the Lagrangian function for '''cls '''test problems. | |||
|'' | |||
|- | |- | ||
|''ls_f''||General routine to compute the objective function value <math>f(x) = | |||
|''ls_f''||General routine to compute the objective function value | {1\over2} r(x)^T r(x)</math> for nonlinear least squares type of problems. | ||
|- | |- | ||
|''ls_g''||General routine to compute the gradient <math>g(x) = J(x)^T r(x)</math> for nonlinear least squares type of problems. | |||
|''ls_g''||General routine to compute the gradient | |||
|- | |- | ||
|''ls_H''||General routine to compute the Hessian approximation <math>H(x) = J(x)^T * J(x)</math> for nonlinear least squares type of problems. | |||
|''ls_H''||General routine to compute the Hessian approximation | |||
|} | |} | ||
==Fitting Sums of Exponentials to Empirical Data== | |||
In TOMLAB the problem of fitting sums of positively weighted exponential functions to empirical data may be formulated either as a nonlinear least squares problem or a separable nonlinear least squares problem. Several empirical data series are predefined and artificial data series may also be generated. There are five different types of exponential models with special treatment in TOMLAB, shown in [[#Table: Exponential models treated in TOMLAB.]]. In research in cooperation with Todd Walton, Vicksburg, USA, TOMLAB has been used to estimate parameters using maximum likelihood in simulated Weibull distributions, and Gumbel and Gamma distributions with real data. TOMLAB has also been useful for parameter estimation in stochastic hydrology using real-life data. | |||
====Table: Exponential models treated in TOMLAB.==== | |||
{| class="wikitable" cellpadding="5" style="background-color:#ffffff;" | |||
|<math>f(t) = \sum\limits_i^p \alpha_i e^{-\beta_i t}</math>, || <math>\alpha_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... <\beta_p</math>. | |||
|- | |- | ||
|<math>f(t) = \sum\limits_i^p \alpha_i e^{-\beta_i t}</math>, || <math>\alpha_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... <\beta_p</math>. | |||
|<math> | |||
|- | |- | ||
|<math>f(t) = \sum\limits_i^p \alpha_i(1-e^{-\beta_i t})</math>, || <math>\alpha_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... <\beta_p</math>. | |||
|<math> | |||
|- | |- | ||
|<math>f(t) = \sum\limits_i^p t \alpha_i e^{-\beta_i t}</math>, || <math>\alpha_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... <\beta_p</math>. | |||
|<math> | |||
|- | |- | ||
|<math>f(t) = \sum\limits_i^p (t \alpha_i-\gamma_i) e^{-\beta_i t}</math>, || <math>\alpha_i,\gamma_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... >\beta_p</math>. | |||
|<math> | |||
|- | |- | ||
|<math>f(t) = \sum\limits_i^p t \alpha_i e^{-\beta_i (t - \gamma_i)}</math>, || <math>\alpha_i \geq 0</math>, || <math>0\leq\beta_1<\beta_2< ... <\beta_p</math>. | |||
|<math> | |||
|} | |} | ||
Algorithms to find starting values for different number of exponential terms are implemented. Test results show that | Algorithms to find starting values for different number of exponential terms are implemented. Test results show that these initial value algorithms are very close to the true solution for equidistant problems and fairly good for non-equidistant problems. Good initial values are extremely important when solving real life exponential fitting problems, because they are so ill-conditioned. [[#Table: Exponential fitting test problems.]] shows the relevant routines. | ||
====Table: Exponential fitting test problems.==== | |||
{| class="wikitable" style="background-color:#ffffff;" | |||
!Function | |||
!Description | |||
|- | |- | ||
|''expAssign''||Assign exponential fitting problem. | |''expAssign''||Assign exponential fitting problem. | ||
|- | |- | ||
|''exp_ArtP''||Generate artificial exponential sum problems. | |||
|'' | |||
|- | |- | ||
|''expInit''||Find starting values for the exponential parameters <math>\lambda</math>. | |||
|''expInit''||Find starting values for the exponential parameters | |||
|- | |- | ||
|''expSolve''||Solve exponential fitting problems. | |''expSolve''||Solve exponential fitting problems. | ||
|- | |- | ||
|''exp_prob''||Defines a exponential fitting type of problem, with data series (''t, y''). The file includes data from several different empirical test series. | |||
|'' | |||
|- | |- | ||
|''Helax_prob''||Defines 335 medical research problems supplied by Helax AB, Uppsala, Sweden, where an exponential model is fitted to data. The actual data series (''t, y'') are stored on one file each, i.e. 335 data files, 8MB large, and are not distributed. A sample of five similar files are part of ''exp prob''. | |||
|'' | |||
|- | |- | ||
|''exp_r''||Compute the residual vector <math>r_i(x), i = 1,...,m.~ x \in \mathbb{R}^n</math> | |||
|'' | |||
|- | |- | ||
|''exp_J''||Compute the Jacobian matrix <math>\partial r_i / d x_j, i=1,...,m, j=1,...,n</math> | |||
|'' | |||
|- | |- | ||
|''exp_d2r''||Compute the 2nd part of the second derivative for the nonlinear least squares exponential fitting problem. | |||
|'' | |||
|- | |- | ||
|''exp_c''||Compute the constraints <math>\lambda_1 < \lambda_2 < ...</math> on the exponential parameters <math>\lambda_i, i=1,...,p</math>. | |||
|'' | |||
|- | |- | ||
|''exp_dc''||Compute matrix of constraint normals for constrained exponential fitting problem. | |||
|'' | |||
|- | |- | ||
|''exp_d2c''||Compute second part of second derivative matrix of the Lagrangian function for constrained exponential fitting problem. This is a zero matrix, because the constraints are linear. | |||
|'' | |||
|- | |- | ||
|''exp_q''||Find starting values for exponential parameters <math>\lambda_i, i=1,...,p</math>. | |||
|'' | |||
|- | |- | ||
|''exp_p''||Find optimal number of exponential terms, ''p''. | |||
|'' | |||
|} | |} | ||
==Large Scale LS problems with Tlsqr== | |||
The ''Tlsqr ''MEX solver provides special parameters for advanced memory handling, enabling the user to solve extremely large linear least squares problems. | The ''Tlsqr ''MEX solver provides special parameters for advanced memory handling, enabling the user to solve extremely large linear least squares problems. | ||
Line 366: | Line 292: | ||
The normal mode of operation of ''Tlsqr ''is that memory for the ''A ''matrix is allocated and deallocated each time the solver is called. In a real-life situation with a very large ''A ''and where the solver is called repeatedly, this may become inefficient and even cause problems getting memory because of memory fragmenting. | The normal mode of operation of ''Tlsqr ''is that memory for the ''A ''matrix is allocated and deallocated each time the solver is called. In a real-life situation with a very large ''A ''and where the solver is called repeatedly, this may become inefficient and even cause problems getting memory because of memory fragmenting. | ||
The ''Tlsqr ''solver provides a parameter ''Alloc'', given as the second element of the first input parameter to control the memory handling. The possible values of ''Alloc ''and their meanings are given in Table | The ''Tlsqr ''solver provides a parameter ''Alloc'', given as the second element of the first input parameter to control the memory handling. The possible values of ''Alloc ''and their meanings are given in [[#Table: Alloc values for Tlsqr]]. | ||
====Table: ''Alloc ''values for ''Tlsqr''==== | |||
{| class="wikitable" | |||
!Alloc (m(2)) | |||
!Meaning | |||
|- | |- | ||
|0||Normal operation: allocate - solve - deallocate | |0||Normal operation: allocate - solve - deallocate | ||
|- | |- | ||
|1||Only allocate, no results returned | |1||Only allocate, no results returned | ||
|- | |- | ||
|2||Allocate and solve, no deallocate | |2||Allocate and solve, no deallocate | ||
|- | |- | ||
|3||Only solve, no allocate/deallocate | |3||Only solve, no allocate/deallocate | ||
|- | |- | ||
|4||Solve and deallocate | |4||Solve and deallocate | ||
|- | |- | ||
|5||Deallocate only, no results returned | |5||Deallocate only, no results returned | ||
|} | |} | ||
An example of the calling sequence is given below. | An example of the calling sequence is given below. | ||
Line 480: | Line 389: | ||
This function, which more often than not needs to be customized to the application in mind, should provide the following functionality: | This function, which more often than not needs to be customized to the application in mind, should provide the following functionality: | ||
< | <source lang="matlab"> | ||
function y = Tlsqrglob( mode, m, n, x, Aname, rw ) | function y = Tlsqrglob( mode, m, n, x, Aname, rw ) | ||
Line 488: | Line 397: | ||
y = A*x; | y = A*x; | ||
else | else | ||
y = A' | y = A'*x; | ||
end | end | ||
</ | </source> | ||
The purpose is to provide the possibility to define a global variable ''A ''and perform the multiplication without transferring this potentially very large matrix to the MEX function ''Tlsqr''. | The purpose is to provide the possibility to define a global variable ''A ''and perform the multiplication without transferring this potentially very large matrix to the MEX function ''Tlsqr''. | ||
If several matrices are involved, for example if ''A ''= | If several matrices are involved, for example if ''A ''= [''A''<sub>1</sub> ; ''A''<sub>2</sub> ], this approach can be used to eliminate the need to explicitly repeatedly form the composite matrix ''A ''during a run. ''Tlsqrglob.m ''should then be (copied and) modified as: | ||
< | <source lang="matlab"> | ||
function y = Tlsqrglob( mode, m, n, x, Aname, rw ) | function y = Tlsqrglob( mode, m, n, x, Aname, rw ) | ||
Line 509: | Line 418: | ||
A2' * x(M+1:end); | A2' * x(M+1:end); | ||
end | end | ||
</ | </source> | ||
To use the global approach, ''Tlsqr ''must be called with the name of the global multiplication routine, for example: | To use the global approach, ''Tlsqr ''must be called with the name of the global multiplication routine, for example: | ||
< | <source lang="matlab"> | ||
[ x, ... ] = Tlsqr(m,n,'Tlsqrglob',...); | [ x, ... ] = Tlsqr(m,n,'Tlsqrglob',...); | ||
</ | </source> |
Latest revision as of 10:18, 5 August 2014
This page is part of the TOMLAB Manual. See TOMLAB Manual. |
This section describes how to define and solve different types of linear and nonlinear least squares and parameter estimation problems. Several examples are given on how to proceed, depending on if a quick solution is wanted, or more advanced tests are needed. TOMLAB is also compatible with MathWorks Optimization TB. See TOMLAB Appendix E for more information and test examples.
All demonstration examples that are using the TOMLAB format are collected in the directory examples. The examples relevant to this section are lsDemo and llsDemo. The full path to these files are always given in the text.
Large Scale LS problems with Tlsqr contains information on solving extreme large-scale ls problems with Tlsqr.
Linear Least Squares Problems
This section shows examples how to define and solve linear least squares problems using the TOMLAB format. As a first illustration, the example lls1Demo in file llsDemo shows how to fit a linear least squares model with linear constraints to given data. This test problem is taken from the User's Guide of LSSOL.
Name='LSSOL test example';
% In TOMLAB it is best to use Inf and -Inf, not big numbers.
n = 9; % Number of unknown parameters
x_L = [-2 -2 -Inf, -2*ones(1,6)]';
x_U = 2*ones(n,1);
A = [ ones(1,8) 4; 1:4,-2,1 1 1 1; 1 -1 1 -1, ones(1,5)];
b_L = [2 -Inf -4]';
b_U = [Inf -2 -2]';
y = ones(10,1);
C = [ ones(1,n); 1 2 1 1 1 1 2 0 0; 1 1 3 1 1 1 -1 -1 -3; ...
1 1 1 4 1 1 1 1 1;1 1 1 3 1 1 1 1 1; 1 1 2 1 1 0 0 0 -1; ...
1 1 1 1 0 1 1 1 1;1 1 1 0 1 1 1 1 1; 1 1 0 1 1 1 2 2 3; ...
1 0 1 1 1 1 0 2 2];
x_0 = 1./[1:n]';
t = []; % No time set for y(t) (used for plotting)
weightY = []; % No weighting
weightType = []; % No weighting type set
x_min = []; % No lower bound for plotting
x_max = []; % No upper bound for plotting
Prob = llsAssign(C, y, x_L, x_U, Name, x_0, t, weightType, weightY, ...
A, b_L, b_U, x_min, x_max); Result = tomRun('lsei',Prob,2);
It is trivial to change the solver in the call to tomRun to a nonlinear least squares solver, e.g. clsSolve, or a general nonlinear programming solver.
Linear Least Squares Problems using the SOL Solver LSSOL
The example lls2Demo in file llsDemo shows how to fit a linear least squares model with linear constraints to given data using a direct call to the SOL solver LSSOL. The test problem is taken from the User's Guide of LSSOL.
% Note that when calling the LSSOL MEX interface directly, avoid using
% Inf and -Inf. Instead use big numbers that indicate Inf.
% The standard for the MEX interfaces is 1E20 and -1E20, respectively.
n = 9; % There are nine unknown parameters, and 10 equations
x_L = [-2 -2 -1E20, -2*ones(1,6)]';
x_U = 2*ones(n,1);
A = [ ones(1,8) 4; 1:4,-2,1 1 1 1; 1 -1 1 -1, ones(1,5)];
b_L = [2 -1E20 -4]';
b_U = [1E20 -2 -2]';
% Must put lower and upper bounds on variables and constraints together
bl = [x_L;b_L];
bu = [x_U;b_U];
H = [ ones(1,n); 1 2 1 1 1 1 2 0 0; 1 1 3 1 1 1 -1 -1 -3; ...
1 1 1 4 1 1 1 1 1;1 1 1 3 1 1 1 1 1;1 1 2 1 1 0 0 0 -1; ...
1 1 1 1 0 1 1 1 1;1 1 1 0 1 1 1 1 1;1 1 0 1 1 1 2 2 3; ...
1 0 1 1 1 1 0 2 2];
y = ones(10,1);
x_0 = 1./[1:n]';
% Set empty indicating default values for most variables
c = []; % No linear coefficients, they are for LP/QP
Warm = []; % No warm start
iState = []; % No warm start
Upper = []; % C is not factorized
kx = []; % No warm start
SpecsFile = []; % No parameter settings in a SPECS file
PriLev = []; % PriLev is not really used in LSSOL
ProbName = []; % ProbName is not really used in LSSOL
optPar(1) = 50; % Set print level at maximum
PrintFile = 'lssol.txt'; % Print result on the file with name lssol.txt
z0 = (y-H*x_0);
f0 = 0.5*z0'*z0;
fprintf('Initial function value %f\n',f0);
[x, Inform, iState, cLamda, Iter, fObj, r, kx] = ...
lssol( A, bl, bu, c, x_0, optPar, H, y, Warm, ...
iState, Upper, kx, SpecsFile, PrintFile, PriLev, ProbName );
% We could equally well call with the following shorter call:
% [x, Inform, iState, cLamda, Iter, fObj, r, kx] = ...
% lssol( A, bl, bu, c, x, optPar, H, y);
z = (y-H*x);
f = 0.5*z'*z;
fprintf('Optimal function value %f\n',f);
Nonlinear Least Squares Problems
This section shows examples how to define and solve nonlinear least squares problems using the TOMLAB format. As a first illustration, the example ls1Demo in file lsDemo shows how to fit a nonlinear model of exponential type with three unknown parameters to experimental data. This problem, Gisela, is also defined as problem three in ls_prob. A weighting parameter K is sent to the residual and Jacobian routine using the Prob structure. The solver clsSolve is called directly. Note that the user only defines the routine to compute the residual vector and the Jacobian matrix of derivatives. TOMLAB has special routines ls f, ls g and ls H that computes the nonlinear least squares objective function value, given the residuals, as well as the gradient and the approximative Hessian, see #Table: Constrained nonlinear least squares (cls) test problems.. The residual routine for this problem is defined in file ls1_r in the directory example with the statements
function r = ls_r(x, Prob)
% Compute residuals to nonlinear least squares problem Gisela
% US_A is the standard TOMLAB global parameter to be used in the
% communication between the residual and the Jacobian routine
global US_A
% The extra weight parameter K is sent as part of the structure
K = Prob.user.K;
t = Prob.LS.t(:); % Pick up the time points
% Exponential expressions to be later used when computing the Jacobian
US_A.e1 = exp(-x(1)*t); US_A.e2 = exp(-x(2)*t);
r = K*x(1)*(US_A.e2 - US_A.e1) / (x(3)*(x(1)-x(2))) - Prob.LS.y;
Note that this example also shows how to communicate information between the residual and the Jacobian routine. It is best to use any of the predefined global variables US_A and US_B, because then there will be no conflicts with respect to global variables if recursive calls are used. In this example the global variable US_A is used as structure array storing two vectors with exponential expressions. The Jacobian routine for this problem is defined in the file ls1_J in the directory example. The global variable US_A is accessed to obtain the exponential expressions, see the statements below.
function J = ls1_J(x, Prob)
% Computes the Jacobian to least squares problem Gisela. J(i,j) is dr_i/d_x_j
% Parameter K is input in the structure
Prob a = Prob.user.K * x(1)/(x(3)*(x(1)-x(2)));
b = x(1)-x(2);
t = Prob.LS.t;
% Pick up the globally saved exponential computations
global US_A
e1 = US_A.e1; e2 = US_A.e2;
% Compute the three columns in the Jacobian, one for each of variable
J = a * [ t.*e1+(e2-e1)*(1-1/b), -t.*e2+(e2-e1)/b, (e1-e2)/x(3) ];
The following statements solve the ''Gisela ''problem.
% ---------------------------------------------------------------------
function ls1Demo - Nonlinear parameter estimation with 3 unknowns
% ---------------------------------------------------------------------
Name='Gisela';
% Time values
t = [0.25; 0.5; 0.75; 1; 1.5; 2; 3; 4; 6; 8; 12; 24; 32; 48; 54; 72; 80;...
96; 121; 144; 168; 192; 216; 246; 276; 324; 348; 386];
% Observations
y = [30.5; 44; 43; 41.5; 38.6; 38.6; 39; 41; 37; 37; 24; 32; 29; 23; 21;...
19; 17; 14; 9.5; 8.5; 7; 6; 6; 4.5; 3.6; 3; 2.2; 1.6];
x_0 = [6.8729,0.0108,0.1248]'; % Initial values for unknown x
% Generate the problem structure using the TOMLAB format (short call)
% Prob = clsAssign(r, J, JacPattern, x_L, x_U, Name, x_0, ...
% y, t, weightType, weightY, SepAlg, fLowBnd, ...
% A, b_L, b_U, c, dc, ConsPattern, c_L, c_U, ...
% x_min, x_max, f_opt, x_opt);
Prob = clsAssign('ls1_r', 'ls1_J', [], [], [], Name, x_0, y, t);
% Weighting parameter K in model is sent to r and J computation using Prob
Prob.user.K = 5;
Result = tomRun('clsSolve', Prob, 2);
The second example ls2Demo in file lsDemo solves the same problem as ls1Demo, but using numerical differences to compute the Jacobian matrix in each iteration. To make TOMLAB avoid using the Jacobian routine, the variable Prob.NumDiff has to be set nonzero. Also in this example the flag Prob.optParam.IterPrint is set to enable one line of printing for each iteration. The changed statements are
...
Prob.NumDiff = 1; % Use standard numerical differences
Prob.optParam.IterPrint = 1; % Print one line each iteration
Result = tomRun('clsSolve',Prob,2);
The third example ls3Demo in file lsDemo solves the same problem as ls1Demo, but six times for different values of the parameter K in the range [3.8, 5.0]. It illustrates that it is not necessary to remake the problem structure Prob for each optimization, but instead just change the parameters needed. The Result structure is saved as an vector of structure arrays, to enable post analysis of the results. The changed statements are
for i=1:6
Prob.user.K = 3.8 + 0.2*i;
Result(i) = tomRun('clsSolve',Prob,2);
end
fprintf('\nWEIGHT PARAMETER K is %9.3f\n\n\n',Prob.user.K);
#Table: Constrained nonlinear least squares (cls) test problems. describes the low level routines and the initialization routines needed for the predefined constrained nonlinear least squares (cls) test problems. Similar routines are needed for the nonlinear least squares (ls) test problems (here no constraint routines are needed).
Table: Constrained nonlinear least squares (cls) test problems.
Function | Description |
---|---|
cls_prob | Initialization of cls test problems. |
cls_r | Compute the residual vector for cls test problems. |
cls_J | Compute the Jacobian matrix for cls test problems. |
cls_c | Compute the vector of constraint functions c(x) for cls test problems. |
cls_dc | Compute the matrix of constraint normals for cls test problems. |
cls_d2c | Compute the second part of the second derivative of the Lagrangian function for cls test problems. |
ls_f | General routine to compute the objective function value for nonlinear least squares type of problems. |
ls_g | General routine to compute the gradient for nonlinear least squares type of problems. |
ls_H | General routine to compute the Hessian approximation for nonlinear least squares type of problems. |
Fitting Sums of Exponentials to Empirical Data
In TOMLAB the problem of fitting sums of positively weighted exponential functions to empirical data may be formulated either as a nonlinear least squares problem or a separable nonlinear least squares problem. Several empirical data series are predefined and artificial data series may also be generated. There are five different types of exponential models with special treatment in TOMLAB, shown in #Table: Exponential models treated in TOMLAB.. In research in cooperation with Todd Walton, Vicksburg, USA, TOMLAB has been used to estimate parameters using maximum likelihood in simulated Weibull distributions, and Gumbel and Gamma distributions with real data. TOMLAB has also been useful for parameter estimation in stochastic hydrology using real-life data.
Table: Exponential models treated in TOMLAB.
, | , | . |
, | , | . |
, | , | . |
, | , | . |
, | , | . |
, | , | . |
Algorithms to find starting values for different number of exponential terms are implemented. Test results show that these initial value algorithms are very close to the true solution for equidistant problems and fairly good for non-equidistant problems. Good initial values are extremely important when solving real life exponential fitting problems, because they are so ill-conditioned. #Table: Exponential fitting test problems. shows the relevant routines.
Table: Exponential fitting test problems.
Function | Description |
---|---|
expAssign | Assign exponential fitting problem. |
exp_ArtP | Generate artificial exponential sum problems. |
expInit | Find starting values for the exponential parameters . |
expSolve | Solve exponential fitting problems. |
exp_prob | Defines a exponential fitting type of problem, with data series (t, y). The file includes data from several different empirical test series. |
Helax_prob | Defines 335 medical research problems supplied by Helax AB, Uppsala, Sweden, where an exponential model is fitted to data. The actual data series (t, y) are stored on one file each, i.e. 335 data files, 8MB large, and are not distributed. A sample of five similar files are part of exp prob. |
exp_r | Compute the residual vector |
exp_J | Compute the Jacobian matrix |
exp_d2r | Compute the 2nd part of the second derivative for the nonlinear least squares exponential fitting problem. |
exp_c | Compute the constraints on the exponential parameters . |
exp_dc | Compute matrix of constraint normals for constrained exponential fitting problem. |
exp_d2c | Compute second part of second derivative matrix of the Lagrangian function for constrained exponential fitting problem. This is a zero matrix, because the constraints are linear. |
exp_q | Find starting values for exponential parameters . |
exp_p | Find optimal number of exponential terms, p. |
Large Scale LS problems with Tlsqr
The Tlsqr MEX solver provides special parameters for advanced memory handling, enabling the user to solve extremely large linear least squares problems.
We'll take the problem of solving Ax = b in the least squares sense as a prototype problem for this section. Here, A ? Rm×n is a dense or sparse matrix and b ? Rm .
Controlling memory allocation in Tlsqr
The normal mode of operation of Tlsqr is that memory for the A matrix is allocated and deallocated each time the solver is called. In a real-life situation with a very large A and where the solver is called repeatedly, this may become inefficient and even cause problems getting memory because of memory fragmenting.
The Tlsqr solver provides a parameter Alloc, given as the second element of the first input parameter to control the memory handling. The possible values of Alloc and their meanings are given in #Table: Alloc values for Tlsqr.
Table: Alloc values for Tlsqr
Alloc (m(2)) | Meaning |
---|---|
0 | Normal operation: allocate - solve - deallocate |
1 | Only allocate, no results returned |
2 | Allocate and solve, no deallocate |
3 | Only solve, no allocate/deallocate |
4 | Solve and deallocate |
5 | Deallocate only, no results returned |
An example of the calling sequence is given below.
>> m = 60000; n = 1000; d = 0.01; % Size and density of A >> A = sprand(m,n,d); % Sparse random matrix >> b = ones(m,1); % Right hand side >> whos A Name Size Bytes Class A 60000x500 3584784 sparse array Grand total is 298565 elements using 3584784 bytes % ======================================================================= % Simple standard call to Tlsqr, Alloc is set to default 0 if m is scalar >> x=Tlsqr(m,n,A,[],[],b); % ======================================================================= % To solve repeatedly with e.g. the same A but different b, % the user may do: % Indicate to Tlsqr to allocate and solve the problem >> m(2) = 2 m = 60000 2 >> x = Tlsqr(m,n,A,[],[],b); % First solution % Indicate to Tlsqr that memory is already allocated, % and that no deallocation should occur on exit >> m(2) = 3 m = 60000 3 % Loop 100 times, calling Tlsqr each time - without re-allocation of memory >> for k=1:100 >> b = (...); % E.g. alter the right hand side each time >> x = Tlsqr(m,n,A,[],[],b); % Call Tlsqr, now with m(2)=3 >> end % Final call, with m(2) = 4: Solve and deallocate >> m(2) = 4 m = 60000 4 >> x=Tlsqr(m,n,A,[],[],b); % Alternatively, to just deallocate, the user could do >> m(2) = 5; >> Tlsqr(m,n,A,[],[],b); % Nothing is returned
Further Memory Control: The maxneA Parameter
If the number of non-zero elements in a sparse A matrix increases in the middle of a Tlsqr-calling loop, the initially allocated space will not be sufficient. One solution is that the user checks this prior to calling Tlsqr and reallocating if necessary. The other solution is to set m(3) to an upper limit (maxneA) of the number of nonzero elements in A in the first allocation call. For example:
>> m = [ 60000 1 1E6 ] m = 60000 1 1000000
will initiate a Tlsqr session, allocating sufficient memory to allow A matrices with up to 1.000.000 nonzeros. If the allocated memory is still insufficient, Tlsqr will try to reallocate enough space for the operation to continue.
Using Global Variables with Tlsqr and Tlsqrglob.m
For cases where it is not possible to send the A matrix to Tlsqr because it is simply too large, the user may choose to use the tomlab/mex/Tlsqrglob.m routine.
This function, which more often than not needs to be customized to the application in mind, should provide the following functionality:
function y = Tlsqrglob( mode, m, n, x, Aname, rw )
global A
if mode==1
y = A*x;
else
y = A'*x;
end
The purpose is to provide the possibility to define a global variable A and perform the multiplication without transferring this potentially very large matrix to the MEX function Tlsqr.
If several matrices are involved, for example if A = [A_{1} ; A_{2} ], this approach can be used to eliminate the need to explicitly repeatedly form the composite matrix A during a run. Tlsqrglob.m should then be (copied and) modified as:
function y = Tlsqrglob( mode, m, n, x, Aname, rw )
global A1 A2
if mode==1
y = A1*x;
y = [y ; A2*x];
else
M = size(A1,1);
y = A1' * x(1:M) + ...
A2' * x(M+1:end);
end
To use the global approach, Tlsqr must be called with the name of the global multiplication routine, for example:
[ x, ... ] = Tlsqr(m,n,'Tlsqrglob',...);