optimtips_21_36.m

New users and old of optimization in MATLAB will find useful tips and tricks in this document, as we

Consider the very simple (linear) regression model:

  y = a0 + a1*x + error.

The sum of squares surface as a function of the parameters
(a0,a1) will have ellipsoidal contours. Theory tells us this,
but its always nice to back up theory with an example. Then
I'll look at what happens in a nonlinear regression.

%}

%% 

% A simple linear model
n = 20;
x = randn(n,1);
y = 1 + x + randn(n,1);
close
figure
plot(x,y,'o')
title 'Linear data with noise'
xlabel 'x'
ylabel 'y'

%%

% linear regression estimates
M = [ones(n,1),x];
a0a1 = M\y

% look at various possible values for a0 & a1
v = -1:.1:3;
nv = length(v);
[a0,a1] = meshgrid(v);
a0=a0(:)';
a1=a1(:)';
m = length(a0);

SS = sum(((ones(n,1)*a0 + x*a1) - repmat(y,1,m)).^2,1);
SS = reshape(SS,nv,nv);

surf(v,v,SS)
title 'Sum of squares error surface'
xlabel 'a0'
ylabel 'a1'
figure
contour(v,v,SS)
hold on
plot(a0a1(1),a0a1(2),'rx')
plot([[-1;3],[1;1]],[[1;1],[-1;3]],'g-')
hold off
title 'Linear model: Sum of squares (elliptical) error contours'
xlabel 'a0'
ylabel 'a1'
axis equal
axis square
grid on

% Note: the min SSE will occur roughly at (1,1), as indicated
% by the green crossed lines. The linear regression estimates
% are marked by the 'x'.

% As predicted, the contours are elliptical.

%%

% next, we will do the same analysis for a nonlinear model
% with two parameters:
%
%   y = a0*exp(a1*x.^2) + error
%
% First, we will look at the error surface

% A simple nonlinear model
n = 20;
x = linspace(-1.5,1.5,n)';
y = 1*exp(1*x.^2) + randn(n,1);
figure
plot(x,y,'o')
title 'Nonlinear data with noise'
xlabel 'x'
ylabel 'y'

%%

% Nonlinear regression estimates. Use pleas for this one.
a1_start = 2;
[a1,a0] = pleas({@(a1,xdata) exp(a1*xdata.^2)},a1_start,x,y)
a0a1 = [a0,a1];

%%

% look at various possible values for a0 & a1
v = .5:.01:1.5;
nv = length(v);
[a0,a1] = meshgrid(v);
a0=a0(:)';
a1=a1(:)';
m = length(a0);

SS = sum((((ones(n,1)*a0).*exp(x.^2*a1)) - repmat(y,1,m)).^2,1);
SS = reshape(SS,nv,nv);
minSS = min(SS(:));

surf(v,v,SS)
title 'Nonlinear model: Sum of squares error surface'
xlabel 'a0'
ylabel 'a1'

figure
contour(v,v,SS,minSS + (0:.25:2))
hold on
plot(a0a1(1),a0a1(2),'rx')
plot([[-1;3],[1;1]],[[1;1],[-1;3]],'g-')
hold off
title 'Nonlinear model: Sum of squares error contours'
xlabel 'a0'
ylabel 'a1'
axis equal
axis square
grid on

% Note: This time, the sums of squares surface is not as
% neatly elliptical. 

%% 28. Confidence limits on a regression model
%{
There are various things that people think of when the phrase
"confidence limits" arises.

- We can ask for confidence limits on the regression parameters
themselves. This is useful to decide if a model term is statistically
significant, if not, we may choose to drop it from a model.

- We can ask for confidence limits on the model predictions (yhat)

These goals are related of course. They are obtained from the
parameter covariance matrix from the regression. (The reference to
look at here is again Draper and Smith.)

If our regression problem is to solve for x, such that A*x = y,
then we can compute the covariance matrix of the parameter vector
x by the simple

  V_x = inv(A'*A)*s2

where s2 is the error variance. Typically the error variance is
unknown, so we would use a measure of it from the residuals.

  s2 = sum((y-yhat).^2)/(n-p);

Here n is the number of data points, and p the number of parameters
to be estimated. This presumes little or no lack of fit in the model.
Of course, significant lack of fit would invalidate any confidence
limits of this form.

Often only the diagonal of the covariance matrix is used. This
provides simple variance estimates for each parameter, assuming
independence between the parameters. Large (in absolute value)
off-diagonal terms in the covariance matrix will indicate highly
correlated parameters.

In the event that only the diagonal of the covariance matrix is
required, we can compute it without an explicit inverse of A'*A,
and in way that is both computationally efficient and as well
conditioned as possible. Thus if one solves for the solution to
A*x=y using a qr factorization as 

  x = R\(Q*y)

then recognize that 

  inv(A'*A) = inv(R'*R) = inv(R)*inv(R') = inv(R)*inv(R)'

If we have already computed R, this will be more stable
numerically. If A is sparse, the savings will be more dramatic.
There is one more step to take however. Since we really want
only the diagonal of this matrix, we can get it as:

  diag(inv(A'*A)) = sum(inv(R).^2,2)
%}

%%

% Compare the two approaches on a random matrix.
A=rand(10,3);

diag(inv(A'*A))

[Q,R]=qr(A,0);
sum(inv(R).^2,2)

% Note that both gave the same result.

%%

% As a test, run 1000 regressions on random data, compare how
% many sets of 95% confidence intervals actually contained the
% true slope coefficient.

m = 1000;  % # of runs to verify confidence intervals
n = 100;   % # of data points
% ci will contain the upper and lower limits on the slope
% parameter in columns 1 and 2 respectively.
ci = zeros(m,2);
p = 2;
for i=1:m
  x = randn(n,1);
  y = 1+2*x + randn(n,1);

  % solve using \
  M = [ones(n,1),x];
  coef = M\y;
  
  % the residual vector
  res = M*coef - y;
  
  % yes, I'll even be lazy and use inv. M is well conditioned
  % so there is no fear of numerical problems, and there are
  % only 2 coefficients to estimate so there is no time issue.
  s2 = sum(res.^2)/(n-p);
  
  sigma = s2*diag(inv(M'*M))';
  
  % a quick excursion into the statistics toolbox to get
  % the critical value from tinv for 95% limits:
  % tinv(.975,100)
  % ans =
  %    1.9840
  ci(i,:) = coef(2) + 1.984*sqrt(sigma(2))*[-1 1];
end

% finally, how many of these slope confidence intervals
% actually contained the true value (2). In a simulation
% with 1000 events, we expect to see roughly 950 cases where
% the bounds contained 2.
sum((2>=ci(:,1)) & (2<=ci(:,2)))

% 950 may not have been too bad a prediction after all.

%%

% Given a complete covariance matrix estimate, we can also compute
% uncertainties around any linear combination of the parameters.
% Since a linear regression model is linear in the parameters,
% confidence limits on the predictions of the model at some point
% or set of points are easily enough obtained, since the prediction
% at any point is merely a linear combination of the parameters.

% An example of confidence limits around the curve on a linear
% regression is simple enough to do.
n = 10;
x = sort(rand(n,1))*2-1;
y = 1 + 2*x + 3*x.^2 + randn(n,1);

M = [ones(n,1),x,x.^2];
coef = M\y;

% Covariance matrix of the parameters
yhat = M*coef;
s2 = sum((y - yhat).^2)/(n-3);
sigma = s2*inv(M'*M)

% Predictions at a set of equally spaced points in [0,1].
x0 = linspace(-1,1,101)';
M0 = [ones(101,1),x0,x0.^2];
y0 = M0*coef;
V = diag(M0*sigma*M0');

% Use a (2-sided) t-statistic for the intervals at a 95% level
% tinv(.975,100)
% ans =
%    1.9840
tcrit = 1.9840;
close all
plot(x,y,'ro')
hold on
plot(x0,y0,'b-')
plot(x0,[y0-tcrit*V,y0+tcrit*V],'g-')
hold off

%% 

% Equality constraints are easily dealt with, since they can be
% removed from the problem (see section 25.) However, active
% inequality constraints are a problem! For a confidence limit,
% one cannot assume that an active inequality constraint was truly
% an equality constraint, since it is permissible to fall any tiny
% amount inside the constraint.
% 
% If you have this problem, you may wish to consider the jackknife
% or bootstrap procedures.

%% 29. Confidence limits on the parameters in a nonlinear regression
%{
A nonlinear regression is not too different from a linear one. Yes,
it is nonlinear. The (approximate) covariance matrix of the parameters 
is normally derived from a first order Taylor series. Thus, of J is
the (nxp) Jacobian matrix at the optimum, where n is the number of
data points, and p the number of parameters to estimate, the
covariance matrix is just

  S = s2*inv(J'*J)

where s2 is the error variance. See section 27 for a description
of how one may avoid some of the problem with inv and forming J'*J,
especially if you only wish to compute the diagonals of this matrix.

What you do not want to use here is the approximate Hessian matrix
returned by an optimizer. These matrices tend to be formed from
updates. As such, they are often only approximations to the true
Hessian at the optimum, at the very least lagging behind a few
iterations. Since they are also built from rank 1 updates to the
Hessian, they may lack the proper behavior in some dimensions.
(For example, a lucky starting value may allow an optimization
to converge in only one iteration. This would also update the
returned Hessian approximation in only 1 dimension, so probably a
terrible estimate of the true Hessian at that point.)

%}

%% 30. Quadprog example, unrounding a curve
%{
A nice, I'll even claim elegant, example of the use of quadprog
is the problem of unrounding a vector. Code for this in the form of
a function can be found on the file exchange as unround.m

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8719&objectType=FILE

Consider a vector composed of some smooth function evaluated
at consecutive (equally spaced) values of its independent variable.
Then the function values are rounded to the nearest integer. Can
we recover the original smooth function values? Clearly we cannot
do so exactly for an arbitrary function, but we can attempt to find
the smoothest set of function values that are consistent with the
given rounded set.

We will treat each element of the set as an unknown, so a vector
of length n will have n unknowns to solve for.
%}

%% 

% A simple function of x
n=101;
x = linspace(0,1,n);
v0 = (10*x.^2);
% rounded vector
vr = round(v0);
close
plot(x,v0,'r-',x,vr,'bo')
title 'Base functional form, with rounding'
xlabel 'x'
ylabel 'y'

%%

% bound constraints. since vr is composed of all integers, then
% the origin, unrounded values must lie within simple bounds
% from our rounded set.
LB = vr - 0.5;
UB = vr + 0.5;

% sparse system representing smoothness. Consider any row of S.
% When multiplied by the vector of unknowns, it computes a
% second order finite difference of the vector around some point.
% This, by itself cannot be used in quadprog, because S is not
% an appropriate quadratic form.
S = spdiags(repmat([1 -2 1],n-2,1),[0 1 2],n-2,n);

% H is positive semi-definite. Quadprog will be happy.
H = S'*S;
f = zeros(n,1);

% Solve the quadratic programming problem. Make sure we use the
% largescale solver. Note that we have carefully constructed
% H to be sparse. The largescale solver will allow simple bound
% constraints. As important as the constraints, when LargeScale
% is 'on', quadprog can handle quite large problems, numbering
% into many thousands of variables. Note that I have supplied H
% as a sparse matrix. It is tightly banded.
options = optimset('quadprog');
options.LargeScale = 'on';
options.Display = 'final';

vur = quadprog(H,f,[],[],[],[],LB,UB,[],options);

plot(x,v0,'r-',x,vr,'bo',x,vur,'k--')
legend('Base curve','Rounded data','Un-rounded curve','Location','NorthWest')
xlabel 'x'
ylabel 'y'

% Clearly the two curves overlay quite nicely, only at the ends
% do they deviate from each other significantly. That deviation
% is related to the behavior of a natural spline, despite the
% fact that we never explicitly defined a spline model. In effect,
% we minimized the approximate integral of the square of the second
% derivative of our function. It is this integral from which cubic
% splines are derived.

%% 31. R^2
%{
R^2 is a measure of the goodness of fit of a regression model. What
does it mean? We can write R^2 in terms of the vectors y (our data)
and yhat (our model predictions) as

  R^2 = 1 - ( norm(yhat-y) / norm(y-mean(y)) )^2

Intuitively, we might describe the term "norm(y-mean(y))" as the 
total information content of our data. Likewise, the term "norm(yhat-y)"
is the amount of signal that remains after removing the predictions
of our model. So when the model describes the data perfectly, this
ratio will be zero. Squared and subtracted from 1, a perfect fit will
yield an R^2 == 1.

At the other end of the spectrum, a fit that describes the data no
better than does the mean will yield a zero value for R^2. It also
suggests that a model with no mean term can produce a negative
value for R^2. The implication in the event of a negative R^2 is the
model is a poorer predictor of the data than would be a simple mean.

Finally, note that the formula above is not specific to a linear
regression model. R^2 may apply as well to a nonlinear regression
model.

When do we use R^2? In the words of Draper & Smith, R^2 is the first
thing they might look at when perusing the output of a regression
analysis. It is also clearly not the only thing they look at. My point
is that R^2 is only one measure of a model. There is no critical value
of R^2 such that a model is acceptable. The following pair of examples
may help to illustrate this point.

%}

%% 

% This curve fit will yield an R^2 of roughly 0.975 from my tests.
% Is it a poor model? It depends upon how much you need to fit the
% subtle curvature in that data.

n=25;
x = sort(rand(n,1));
y = exp(x/2)+randn(size(x))/30;
M = [ones(n,1),x];

coef = M\y;

yhat = M*coef;
Rsqr = 1 - (norm(yhat-y)/norm(y-mean(y))).^2
close
plot(x,y,'bo',x,yhat,'r-')
title 'R^2 is approximately 0.975'
xlabel x
ylabel y

% Was the linear model chosen adequate? It depends on your data and
% your needs. My personal opinion of the data in this example is
% that if I really needed to know the amount of curvature in this
% functional relationship, then I needed better data or much more
% data.

%%

% How about this data? It yields roughly the same value of R^2 as did
% the data with a linear model above.

n = 25;
x = linspace(0,1,n)';
y = x + ((x-.5).^2)*0.6;
M = [ones(n,1),x];

coef = M\y;

yhat = M*coef;
Rsqr = 1 - (norm(yhat-y)/norm(y-mean(y))).^2

plot(x,y,'bo',x,yhat,'r-')
title 'R^2 is approximately 0.975'
xlabel x
ylabel y

% Is there any doubt that there is lack of fit in this model? You
% may decide that a linear model is not adequate here. That is a
% decision that you must make based on your goals for a model.

%%

%{
To sum up, the intrepid modeler should look at various statistics
about their model, not just R^2. They should look at plots of the
model, plots of residuals in various forms. They should consider
if lack of fit is present. They should filter all of this information
based on their needs and goals for the final model and their knowledge
of the system being modeled.
%}

%% 32. Estimation of the parameters of an implicit function
%{
How can we estimate the parameters for implicit functions, one
simple example might be:

  y = a*(y - b*x)^2

The simplistic approach is to formulate this as a direct least
squares problem for lsqnonlin. 

  y - a*(y - b*x)^2 = 0