File Exchange

LSE

version 1.1 (4.56 KB) by

A linear least squares solver, subject to linear equality constraints

Updated

This submission was written by request - as a tool to handle linear least squares problems, subject to linear equality constraints that may potentially be rank deficient. (It handles problems with full rank constraints of course too.) In the event of a rank deficient constraint system, it tests for consistency of the constraints.

I added a few other features to LSE:

- It allows multiple right hand sides to the least squares problem, fully vectorized of course.
- Weights may be supplied.
- You are offered a choice of least squares solvers, either backslash or pinv.

LSE solves the problem (for an unknown vector x)

argmin norm(A*x - b)

subject to the constraints

C*x = d

As an example, consider the random system
A = rand(10,3);
b = rand(10,1);

With a rank deficient constraint set
C = [1 1 1;1 1 1];
d = [1;1];

X = lse(A,b,C,d)
X =
0.5107
0.57451
-0.085212

Verify that the constraints are satisfied

C*X
ans =
1
1

Column pivoting is used to eliminate variables from the constraint system when \ is specified, and when pinv is specified, an svd is used for the final solution.

John D'Errico

John D'Errico (view profile)

Royi - I'm sorry, but that question is impossible to answer, as it depends on too many factors. There are virtually no hard, fast limits, but this is true of any of the linear algebra calls in MATLAB.

How much memory do you have? At some point, matrix operations are no longer done strictly in the available RAM, but are forced to start swapping using virtual memory. So your available RAM is a limiting factor. The speed of your computer, your system bus, hard drive, etc., all are factors.

As well, your patience is a factor. At some point, you would decide that it was taking too long, and I can never know what you would consider too much time for a given size problem.

Finally, it is not only the number of unknowns that are a factor. Other factors over which I have no control are the number of data points in the least squares problem, and the number of equality constraints.

I would only point out that the QR based (default) option will probably be considerably faster than the SVD/PINV based alternative solver in this code.

Royi Avital

John D'Errico

John D'Errico (view profile)

Sorry, but making d, the right hand side of the CONSTRAINT system seems to provide little value at the cost of making the code a bit more complicated for little gain.

If you absolutely need separate right hand sides for the constraint, you can simply put a loop around the call. Surely that is the simple solution.

Erd

Erd (view profile)

why d needs to be a vector (mx1), shouldnt it allow (mxq)?

Michael Bevis

Michael Bevis (view profile)

I recently encountered a least squares problem with equality constraints that caused lsqlin to crash, but lse solved it nicely producing an excellent fit to the data vector, and obeying the constraints. Thanks, John!

John D'Errico

John D'Errico (view profile)

Greg - If you want to apply inequality or bound constraints, the simplest way is to use lsqlin from the optimization toolbox. That tool implements all forms of linear constraints.

Greg

Greg (view profile)

Great piece of code! I do have a question? Is there a way to bound the results to be positive (i.e. x>0)?

Matt Ferrara

Great code! Do you have any references on the approach you have implemented? (maybe you could put them in the documentation for release 3?)

Lukas Malec

Thank you. It's a very good work!