« Previous « Start » Next »
8 Solving Least Squares and Parameter Estimation Problems
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 Appendix
E for
more information and test examples.
All demonstration examples that are using the TOMLAB (TQ) format
are collected in the directory
examples . Running the menu
program
tomMenu , it is possible to run all demonstration
examples. It is also possible to run each example separately. The
examples relevant to this section are
lsDemo and
llsDemo . All files that show how to use the Init File
format are collected in the directory
usersguide . The full
path to these files are always given in the text.
Section
8.5
contains information on
solving extreme large-scale
ls problems with
Tlsqr .
8.1 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 Users Guide of
LSSOL [
32].
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.
8.2 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 Users Guide of
LSSOL
[
32].
% 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);
8.3 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
8.3. 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 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);
fprintf('\nWEIGHT PARAMETER K is %9.3f\n\n\n',Prob.user.K);
end
Table
8.3 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).
Constrained nonlinear least squares (cls) test problems. |
|
| Function |
Description |
|
cls_prob |
Initialization of cls test problems. |
| cls_r |
Compute the residual vector ri(x), i = 1,...,m. x  Rn
for cls test problems. |
| cls_J |
Compute the Jacobian matrix Jij(x)=∂ ri / d xj,
i=1,...,m, j=1,...,n 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 ∂ c(x)/dx for
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 f(x) =
1/2 r(x)T r(x) for nonlinear least squares type of problems. |
| ls_g |
General routine to compute the gradient g(x) = J(x)T r(x) for
nonlinear least squares type of problems. |
| ls_H |
General routine to compute the Hessian approximation H(x) = J(x)T * J(x)
for nonlinear least squares type of problems. |
|
8.4 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 [
69].
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
11.
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 11: Exponential models treated in TOMLAB.
|
| f(t) = Σip αi e−βi t, |
αi ≥ 0, |
0≤β1<β2< ... <βp. |
| f(t) = Σip αi(1−e−βi t), |
αi ≥ 0, |
0≤β1<β2< ... <βp. |
| f(t) = Σip t αi e−βi t, |
αi ≥ 0, |
0≤β1<β2< ... <βp. |
| f(t) = Σip (t αi−γi) e−βi t, |
αi,γi ≥ 0, |
0≤β1<β2< ... <βp. |
| f(t) = Σip t αi e−βi (t − γi), |
αi ≥ 0, |
0≤β1<β2< ... <βp. |
|
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,
see the thesis by Petersson [
64]. Good initial values are
extremely important when solving real life exponential fitting
problems, because they are so ill-conditioned. Table
8.4 shows the relevant routines. The best way to define
new problems of the predefined exponential type is to edit the
exp_prob.m Init File as described in
Appendix
D.9.
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 ri(x), i = 1,...,m. x Rn |
| exp_J |
Compute the Jacobian matrix ∂ ri / d xj, i=1,...,m, j=1,...,n. |
| exp_d2r |
Compute the 2nd part of the second derivative for the nonlinear least
squares exponential fitting problem. |
| exp_c |
Compute the constraints λ1 < λ2 < ... on the exponential
parameters λi, i=1,...,p. |
| 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 λi, i=1,...,p. |
| exp_p |
Find optimal number of exponential terms, p. |
|
The algorithmic development implemented in TOMLAB is further discussed
in [
52]. An overview of the field is also given in this
reference.
8.5 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
12.
Table 12:
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=[
A1 ;
A2], 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',...);
« Previous « Start » Next »
TOMLAB
