ublas and matrix expression
I have an expression a = (X^T * X)^-1 * X^T * y where x is an n * k matrix and y is an n row vector How do I implement this expression in ublas? I can create and populate the matrices and can do the product of matrices and the transpose use prod and trans, but I'm unsure on the ^-1. Can I do this with ublas? (I don't know much about matrices, but have been given the above expression as a way of finding a polynomial fit for a set of data). Thanks in advance Russell
Russell Hind wrote:
I have an expression
a = (X^T * X)^-1 * X^T * y
where x is an n * k matrix and y is an n row vector
How do I implement this expression in ublas? I can create and populate the matrices and can do the product of matrices and the transpose use prod and trans, but I'm unsure on the ^-1. Can I do this with ublas?
(I don't know much about matrices, but have been given the above expression as a way of finding a polynomial fit for a set of data).
Hi, Russell. You are trying to solve a matrix expression. A = B^-1 * C*D is equivalent to B*A = B*B^-1 * C*D which collapses to B*A = C*D So, you can reprase your problem as solve the system B*A = E where B,E are known, E=C*D, to give the unknown A. Have a look at the (undocumented :-() lu_factorize and lu_substitute routines. They should enable you to complete the task. HTH, -- Angus
Angus Leeming wrote:
Hi, Russell.
You are trying to solve a matrix expression.
A = B^-1 * C*D is equivalent to B*A = B*B^-1 * C*D which collapses to B*A = C*D
So, you can reprase your problem as solve the system B*A = E where B,E are known, E=C*D, to give the unknown A.
Have a look at the (undocumented :-() lu_factorize and lu_substitute routines. They should enable you to complete the task.
Thanks, where are these methods defined? I'm using 1.30.2 and have searched the ublas directory for files with lu_factorize in and can't find it. Thanks Russell
Russell Hind wrote:
Thanks, where are these methods defined? I'm using 1.30.2 and have searched the ublas directory for files with lu_factorize in and can't find it.
I do not know when these routines arrived, but they are certainly
present in boost 1.31. Example usage:
#include
Angus Leeming wrote:
Russell Hind wrote:
Thanks, where are these methods defined? I'm using 1.30.2 and have searched the ublas directory for files with lu_factorize in and can't find it.
I do not know when these routines arrived, but they are certainly present in boost 1.31. Example usage:
Yep, found lu.hpp in CVS. This is what I think I have
matrix<double> X;
vector<double> y;
These are known and I have created them.
I'm trying to solve
Xa = y
So should I be doing this? (I'm getting compile errors with bcc so
can't check anything at the moment).
permutation_matrix
Russell Hind wrote:
Yep, found lu.hpp in CVS. This is what I think I have
matrix<double> X; vector<double> y;
These are known and I have created them.
I'm trying to solve
Xa = y
So should I be doing this? (I'm getting compile errors with bcc so can't check anything at the moment).
permutation_matrix
pm(X.size1()); lu_factorize(X, pm) lu_substitute(X, pm, y); so on the way out of lu_substitute, y will contain the values for a?
Correct. Note that X is modified by the call to lu_factorize and y is modified by the call to lu_substitute.
Or am I completely on the wrong track?
If not, is there another way of solving this without using lu_*?
Sure. Use the matrix solution method of your choice. Does the MTL not interoperate with uBlas? Matrix solvers aren't anything to do with BLAS and the fact that they're in uBlas is probably an 'embarassment'. (By which I mean that they extend the remit of the library beyond its original intent.)
the initial expression I gave came from re-arranging Xa = y to give it in terms of a. I'm currently stuck with 1.30.2 but would like to be able to solve this using boost, rather than having to find another library.
Well, if matrix X isn't too big, then the lu_ methods are just what you need. -- Angus
Angus Leeming wrote:
Correct. Note that X is modified by the call to lu_factorize and y is modified by the call to lu_substitute.
Thanks. Unfortunately I'm stuck with boost-1.30.2 and it doesn't have lu_* in it. I have lu.hpp but I get compile errors with bcc-5.6.4 in lu.hpp hence why I wondered if there was a way to do it using ublas, but not lu_*.
If not, is there another way of solving this without using lu_*?
Sure. Use the matrix solution method of your choice. Does the MTL not interoperate with uBlas?
Don't know what MTL is. I have a look. Thanks
Well, if matrix X isn't too big, then the lu_ methods are just what you need.
X is typically ~500 by 5 (don't know if that counts as big or not, last time I did matrix work, we only used 3x3 and 4x4 so to me it seems pretty huge, but we were doing stuff by hand). I'll have to have a look around for another library if there is no way to do it with ublas without using lu_*. Thanks for your help Russell
Russell Hind wrote:
Don't know what MTL is. I have a look. Thanks
Matrix Template Library. Jeremy Siek's master's thesis.
X is typically ~500 by 5 (don't know if that counts as big or not, last time I did matrix work, we only used 3x3 and 4x4 so to me it seems pretty huge, but we were doing stuff by hand).
It isn't too terrible, but to invert it it must be square... -- Angus
Angus Leeming wrote:
Russell Hind wrote:
Don't know what MTL is. I have a look. Thanks
Matrix Template Library. Jeremy Siek's master's thesis.
I've found the link. I'll see if I have any luck with that
X is typically ~500 by 5 (don't know if that counts as big or not, last time I did matrix work, we only used 3x3 and 4x4 so to me it seems pretty huge, but we were doing stuff by hand).
It isn't too terrible, but to invert it it must be square...
Well, the bit I need to invert is prod(trans(X), X) so X is 500 x 5. Is this possible? Thanks Russell
Russell Hind wrote:
It isn't too terrible, but to invert it it must be square...
Well, the bit I need to invert is prod(trans(X), X) so X is 500 x 5. Is this possible?
X is 500x5 trans(X) is 5x500 prod(trans(X),X) is 5x5, which is a square matrix. So, yes, it can be inverted (all other requirements assumed satisfied...). And if it's 5x5 the lu_ functions are perfect. In fact, you could probably invert it by hand ;-) -- Angus
Angus Leeming wrote:
X is 500x5 trans(X) is 5x500 prod(trans(X),X) is 5x5, which is a square matrix. So, yes, it can be inverted (all other requirements assumed satisfied...). And if it's 5x5 the lu_ functions are perfect. In fact, you could probably invert it by hand ;-)
Is there no method in ublas that can invert it? I haven't done matrices for years, and don't have any docs on doing it. The initial problem I'm trying to solve is a = (X^T X)^-1 X^T y where y is a 500 row vector and X is the 500 x 5 matrix. But I can't use lu_* with the compiler/boost version I have. ublas::solve says it does A^-1 * b but what is the tag parameter passed to it? It can be lower_tag or upper_tag but I can't find what these mean in the docs, or the difference between them or which I should be using. Thanks Russell
seeing the dimension of your matrix, the easiest way to do that is: 1. compute z = X^T y: z is 5 x 1 2. compute A = X^T X with prod (trans (X), X): A is 5 x 5 3. use LAPACK to solve: A a = z (5 x 5 linear system) ublas::solve only works for particular types of matrices (for example, triangular matrices - upper or lower... - like the tag...). regards, lorenzo. On Thu, Apr 08, 2004 at 04:36:17PM +0100, Russell Hind wrote:
Angus Leeming wrote:
X is 500x5 trans(X) is 5x500 prod(trans(X),X) is 5x5, which is a square matrix. So, yes, it can be inverted (all other requirements assumed satisfied...). And if it's 5x5 the lu_ functions are perfect. In fact, you could probably invert it by hand ;-)
Is there no method in ublas that can invert it? I haven't done matrices for years, and don't have any docs on doing it.
The initial problem I'm trying to solve is
a = (X^T X)^-1 X^T y
where y is a 500 row vector and X is the 500 x 5 matrix. But I can't use lu_* with the compiler/boost version I have.
ublas::solve says it does A^-1 * b but what is the tag parameter passed to it? It can be lower_tag or upper_tag but I can't find what these mean in the docs, or the difference between them or which I should be using.
Thanks
Russell
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
-- CH3 | N / \ N----C C==O || || | || || | CH C N--CH3 \ / \ / N C | || CH3 O
Lorenzo Bolla wrote:
seeing the dimension of your matrix, the easiest way to do that is: 1. compute z = X^T y: z is 5 x 1 2. compute A = X^T X with prod (trans (X), X): A is 5 x 5 3. use LAPACK to solve: A a = z (5 x 5 linear system)
Thats where I've go to, I just need the routines to solve the Aa=z. As mentioned, I don't have boost-1.31.0 so I can't use the lu_* methods from there. I had a quick look at the matrix template library, but didn't have much luck getting it to compile with Borland bcc32. Are there any small, light-weight libraries out there that can solve this? I'd rather something that was source only (as uBlas is), so I don't have to build and link with another library / dll.
ublas::solve only works for particular types of matrices (for example, triangular matrices - upper or lower... - like the tag...).
Thanks, I've just found out what triangular matrices are :) Cheers Russell
Russell Hind wrote: [...]
I'll have to have a look around for another library if there is no way to do it with ublas without using lu_*.
You can try to use bindings library which `binds' ublas with ATLAS and LAPACK (so you must install either ATLAS or LAPACK). You can download bindings library either from files section of ublas-dev mailing list: http://groups.yahoo.com/group/ublas-dev/ (and I think you must subscribe to access it), or from boost-sandbox CVS (headers are in boost/numeric/bindings, and docs and examples are in libs/numeric/bindings). Regards, fres PS. Home page of the MTL is http://www.osl.iu.edu/research/mtl/ And AFAIK MTL doesn't work with ublas (or vice versa).
depending on the shape and sparsity of the matrix you have, there are a lot of ways you can solve your problem. simplifying, you can use: - direct solvers: for example LU decomposition. i think the best library to do that is UMFPACK 4.3 if the matrix is sparse and any known LAPACK implementation if it's dense: uBlas has bindings for both of them (not included in the library - but look for "uBlas bindings" in google...). use direct solvers only if the matrix is not too big, because they need a lot of memory. - iterative solvers: here the choice is much more difficult because the "best" algorithm strongly depends on the matrix you have. None of them are implemented in uBlas (even if it would be nice to...). to try, you can have a look at a very basic implementation of the most common iterative algorithms in the IML++ (Iterative Methods Library). regards, Lorenzo. On Thu, Apr 08, 2004 at 02:33:35PM +0100, Russell Hind wrote:
Angus Leeming wrote:
Russell Hind wrote:
Thanks, where are these methods defined? I'm using 1.30.2 and have searched the ublas directory for files with lu_factorize in and can't find it.
I do not know when these routines arrived, but they are certainly present in boost 1.31. Example usage:
Yep, found lu.hpp in CVS. This is what I think I have
matrix<double> X; vector<double> y;
These are known and I have created them.
I'm trying to solve
Xa = y
So should I be doing this? (I'm getting compile errors with bcc so can't check anything at the moment).
permutation_matrix
pm(X.size1()); lu_factorize(X, pm) lu_substitute(X, pm, y); so on the way out of lu_substitute, y will contain the values for a?
Or am I completely on the wrong track?
If not, is there another way of solving this without using lu_*? the initial expression I gave came from re-arranging Xa = y to give it in terms of a. I'm currently stuck with 1.30.2 but would like to be able to solve this using boost, rather than having to find another library.
Thanks
Russell
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
-- CH3 | N / \ N----C C==O || || | || || | CH C N--CH3 \ / \ / N C | || CH3 O
participants (4)
-
Angus Leeming -
Kresimir Fresl -
Lorenzo Bolla -
Russell Hind