|
|
Ublas : |
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2007-06-01 11:14:51
On Friday 01 June 2007 15:10, Marchesi Mauro wrote:
> Hello.
>
> I tried to solve a linear system Ax=b by using the following code:
>
> ublas::lu_factorize( M, pm );
> ublas::lu_substitute( M, pm, x );
>
> This works fine, unless A is singular.
> However it turns out that, because of the way the matrix is built, the
> rank of the coefficient matrix A equals the rank of the augmented matrix
> Ab.
The LU-decomposition (or Gaussian elimination) only works for regular
matrices.
> So, even in this case, the linear system admits solutions, but I can't
> figure out how to get them by using the uBLAS library.
>
> I'd like to get the set of solutions in any form: for example as a
> vector x0 along with a set of vectors {x1,...,xN} such that every
> solution of the system can be expressed in the form x0 + k1*x1 + ... +
> kN*xN (clearly the set {x1,...,xN} is empty iff A is non-singular).
>
> Can anyone help me?
If you are looking for solutions to over/underdetermined systems you have to
solve an auxiliary problem that has a unique solution.
||Ax-b|| -> min (if Ax=b has no solution)
||x||->min, s.t. Ax=b (if Ax=b has many solutions)
both can be done by using a QR decomposition and solving the so called "normal
equation": A^T A x = A^T b
(QR is not included in ublas, yet)
If you have A = QR with Q^T Q = I and R upper triangular then the normal
equation is equivalent to
R^T R x = R^T Q^T b
because (Q^T R)^T QR = R^T R.
If R is like
[ R' ]
[ 0 ]
with a regular square matrix R' you finally have to solve
R'^T R x = R^T b' <==> R x = b'
(b' are the first k rows of Q^T b)
If R is like
[ R' 0 ]
with a regular square matrix R' you finally have to solve
R x' = Q^T b
(where x' are the first k entries of x)
Summary: The solution x of R x = Q^T b (with appropriate zero padding) is
always the best possible solution of Ax=b in the above sense.
If you need more information you should do a singular value decomposition of
A.
mfg
Gunter
PS: read Golub van Loan (or similar book on numerical linear algebra)