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)
Verify that the constraints are satisfied
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.
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.
How large are the matrices it can handle?
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.
why d needs to be a vector (mx1), shouldnt it allow (mxq)?
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!
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.
Great piece of code! I do have a question? Is there a way to bound the results to be positive (i.e. x>0)?
Great code! Do you have any references on the approach you have implemented? (maybe you could put them in the documentation for release 3?)
Thank you. It's a very good work!
Bug fix - Single constraint problems on R2009 releases failed due to a qr issue - this fix repairs that bug.
Version 2.0 - Used a better solution for when pinv was specified. Also allowed the code to work when no constraints at all are supplied.