EliminateConstraints

Eliminates variables from a problem with linear equality constraints to give an unconstrained proble

Description

function [C,d]=eliminateConstraints(A,b) eliminates variables from a problem with linear equality constraints to give an unconstrained problem. This is useful e.g. when solving a problem with linear constraints and a nonlinear objective or further nonlinear constraints; eliminating the linear constraints makes the problem easier.

Where the original constrained problem has:
 - variable vector x (length n).
 - linear constraint A*x=b, where the matrix A and vector b are
     eliminateConstraints's arguments.
 
And the final unconstrained problem has:
 - variable vector y (length p).
 - an x satisfying the original constraint can be generated by x=C*y+d, where the matrix C and vector d are returned by eliminateConstraints.
 
Inputs:
 - A (an m*n matrix of real numbers)
 - b (a m*1 vector of real numbers)

 - C (an n*p matrix of real numbers)
 - d (a n*1 vector of real numbers)

Hopefully, p<n, so the unconstrained problem is smaller than the constrained problem.

The functionality provided by eliminateConstraints is trivial for small matricies. For larger matrices, it can be done using MATLAB's various matrix decompositions. However, as far as the author is aware, there is no trivial MATLAB solution for the case where A is very large and sparse, so we require a C which is also sparse without using any very large full matrices.

The case for which this was developed automatically creates a well-defined problem, so only minimal validation and error checking has been implemented at this stage.

• MATLAB release: MATLAB 7.6 (R2008a)

[Dimitri Shvorob]

I'll be happy to apologize if I'm being dumb, but: 'solution' of what? In what sense does this finding of C and d 'eliminate linear constraints'? What is it for?

[John D'Errico]

Dimitri is not being dumb here. The author has failed to provide a tool that has any useful documentation. Your job as a programmer does not stop when you write the last line of code. If you think so, then you should be fired. You should document your code. Provide help. Otherwise, that code is just a bunch of random bits, useful to nobody else in the world.

This code has a couple of uses that I might think of. (The fact that I need to document the author's code is a major failure on this author's part. Does my irritation show through?)

You MIGHT choose to use this tool to enable an unconstrained nonlinear least squares solver to solve a linearly equality constrained problem, with a nonlinear objective function. For example, suppose you wish to use fminsearch to solve a problem where the sum of the parameters must be 1. You could use this function to help. You would still need to know what to do, but the author failed to provide any hints. Too bad.

Another possible use might be to solve a linear least squares problem with linear equality constraints. Suppose you wished to solve the problem

  G*x = h

subject to the equality constraint(s)

  A*x = b

This code provides a transformation x = C*y + d, so that the linear least squares problem reduces to

  A*(C*y+d) = b

or

  (A*C)*y = b - A*d

Solve it as

  y = (A*C)\(b - A*d)

and then recover the value of x from the supplied relation

  x = C*y + d

So, in theory, this code has some utility. Without the comments I've just provided, it is useless though to anybody who does not know enough about the problem that what it does is trivial. My point is that while I personally might in theory have found a tool like this of some use, I also know how to do what it does in only a few lines, so it still offers no value to me.

How about the help? Poor. It does not tell you what the inputs should be, in terms of their order, shape, class, or size. The outputs are similarly undescribed.

There is error checking, but since there is no real description of the inputs, how can you check for errors? There is no error check on the number of columns of b either (should be 1).

I feel there is also a flaw in the error checks. If you pass in an inconsistent system, the code will error out, but the error generated will be a very confusing one.

There are some internal comments, and the author has chosen descriptive variable names, so you can read the code.

How about the code itself? This uses loops to do what can be done in one line of code, at least if your matrix is full. I've not checked the algorithm itself to decide if it is as stable, as say, a call to null, but on a few tests I found that a numerically nonsingular (although poorly conditioned) matrix caused this code to fail.

An interesting point is, if A has many columns (i.e., x has many elements) but few rows, then the resulting C matrix will be very large, but if A is sparse, then C will be stored in a sparse format, and it will be sparse. The null solution would yield a full matrix here.

Had this code full documentation, coupled with complete examples on how to use it for at least the two classes of problems I've described, and if it had good help, then it would have earned a better rating from me. As it is, I can't justify a better rating than to say it needs improvement.

[John D'Errico]

If I could recommend something...

I'd add a published m-file to this that has several examples of how to use it, and what features it has. When you publish the file, it will be seen as a list of examples that clearly show how to use your code. In fact, I'll even give you a start on it. Publish the code I've included below, then zip it all up with your m-file and update the zipped file here on the FEX.

%% Minimize an equality constrained optimization problem
A = [1 1 1];
b = 1;
[C,d] = eliminateConstraints(A,b)

%%
% set up a quadratic form to minimize.
%
% With no constraints, the minimizer is at [0 0 0].
H = [2 -1 0;-1 2 -1;0 -1 2];

fun = @(Y) (C*Y+d)'*H*(C*Y+d);
Y = fminsearch(fun,[1 1]',optimset('tolx',1.e-7))

%%
% The solution should be [.3 .4 .3]
X = C*Y+d

%% Solve a constrained linear least squares
% G*X = h, aubject to the constraint A*x = b.
G = rand(10,5);
h = rand(10,1);

A = rand(2,5);
b = rand(2,1);

%%
% Solve it using lsqlin (in the optimization toolbox)
X0 = lsqlin(G,h,[],[],A,b)

%%
% Use eliminateConstraints instead
[C,d] = eliminateConstraints(A,b)

%%
% and solve it using backslash
Y1 = (G*C)\(h - G*d)
X1 = C*Y1 + d

%%
% Show that X0 and X1 were the same
[X0,X1]

%% Show that for a sparse matrix A that C is also sparse
A = sprand(10,100,.05);
spy(A)
title 'Spy of A matrix'
b = rand(10,1);
[C,d] = eliminateConstraints(A,b);
figure
spy(C)
title 'Spy for C matrix'