solutions for linear systems

I found a page with an implementation of an algorithm that sounds usefull, called "Bounded Variable Least Squares": http://lib.stat.cmu.edu/general/bvls

As it turns out, the projection formula can actually be used to find the solution to a system of equations that minimizes the least-squared error. Put all the 'a' vectors into the columns of a matrix, and call it A. The solution is then (A^TA)^-1Ab (assuming the dimensions work out). I'm not sure if there is a simpler approach when the values of A are constrained, but this equation works in the general case so it should also work in the constrained case.

Lord Crc - Thanks for the link, I'll check it out after work.

Zipster - That's a very nice solution. At the moment I don't think the dimensions will work out in the majority of the cases (I haven't quite figured out the dimensions of the vectors or the amount of terms yet), but as far as I understand it, I could pad the vectors with zeros until the dimensions do work out. Would that still give me reasonable results?

edit: Actually I don't think that would work. Padding with zeros would give me a determinant of 0, which means I couldn't invert the matrix. I think...

[Edited by - Thergothon on September 1, 2008 12:47:13 AM]

Quote:

Original post by Thergothon Zipster - That's a very nice solution. At the moment I don't think the dimensions will work out in the majority of the cases (I haven't quite figured out the dimensions of the vectors or the amount of terms yet), but as far as I understand it, I could pad the vectors with zeros until the dimensions do work out. Would that still give me reasonable results?

The formula Zipster gave is for the Moore-Penrose pseudoinverse. It is not sensitive to the dimensions of the matrix. Any sized matrix will work, inasmuch as all operations will be well-defined.

To expand on what Sneftel said, the only requirement for using this pseudo-inverse is that the quantity A^TA be invertible. Note that this quantity is always a square matrix, regardless of the dimensions of A. But in order for it to be invertible, A has to have linearly independent rows and/or columns -- in other words, full rank in at least one dimension. You can't really think in terms of the determinant, because that's only valid for square matrices and A isn't necessarily square. Adding "pad" zeroes isn't necessarily a bad thing as long as you don't suddenly add any linear dependencies.

Thanks for all the help guys, I'd never even heard of a pseudo-inverse before, that will be a handy thing to remember.

As for the link posted by lord crc, it looks perfect, but seriously, they need to republish it in something that isn't Fortran. Once (if) I get it working, I'll definitely post the source here to save others the pain.

On Fortran, does it use 1 or 0 based indexing? i.e. 'do 10 i=1, n' equivalent to 'for(size_t i = 1; i <= n; ++i)'. Also does 'sign(sqrt(sq), a(j,j))' mean 'the value of sqrt(sq) with the sign of a(j, j)'? And what does this: 'Relies on FORTRAN77 do-loop conventions!' comment refer to?

I can tell this is going to be a nightmare to debug...

Once again, your help is much appreciated.
s.

edit: And could a mod please change the thread title to something like 'finding solutions to linear systems with bounded variables'? Thanks again.

Glad my googles helped you :) I'd be interested in the source as well if you manage to port it.

As for your questions:

"The sign transfer function SIGN(X,Y) takes the sign of the second argument and puts it on the first argument, ABS(X) if Y >= 0 and -ABS(X) if Y < 0." (from this nice page).

As for the loop questions, I found this page to be nice.

If you use the conjugate gradient method http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf. It will give you a solution even if the matrix is ill conditioned or does not have an inverse. In those cases the answer you will get will be the soultion with the "highest energy". It also uses a set number of operations which can be useful in a situation where the timing is critical.