RE: Computing a Projection Matrix

["Mikael Adlers" <mikael dot adlers at mathcore dot com>]

Hi!
It depends on what you want do with the matrix. If you only
want to apply the matrix on vector then you shouldn't form the matrix since
it costs very much and you can loos a lot of accuracy (the projection matrix
may become very large and costly to apply)

I would compute the result in the following way (if the matrix is not ill-conditioned)

Assume that you want to compute 

                 T  -1 T
    y = (I - A (A A)  A ) b 

Compute:

1) v = A^Tb

2) solve system (A^TA)x = v (you store the factorization instead of the inverse)

3) y = x - Ax

If you have computed the factorization of A^TA and A is a m-by-n matrix.
Then to form y costs 2nm + 2n^2 + 2nm + 2m mult/additions. 

Compare that to what it costs you if you form the projection matrix : 2m^2

m=100
n=10
the vector approach: 2000 + 200 + 2000 + 200 = 4400
the matrix approach: 20000

If is also very! expensive to form the projection matrix.


Note that A^TA is symmetric and it is more efficient to compute the Cholesky 
factorization of the matrix (B = L L^T) (see http://sources.redhat.com/gsl/ref/gsl-ref_13.html#SEC216)


If A is ill-conditioned the matrix A^TA has the condition number of A squared
and you can get a big loss of accuracy. However, there is a way around.

Look at the step 1) and 2) above: (A^TA)x = A^Tb
This is the normal equation that solves the least squares problem
min_x ||Ax - b||_2. So if you would like to solve this part
with high accuracy you can use the QR decomposition of A to compute x.
(see http://sources.redhat.com/gsl/ref/gsl-ref_13.html#SEC213)

Hope this helps.

Regards,
Mikael Adlers


------------------------------------------------------------------ 
 Mikael Adlers, Ph.D.          email: mikael@mathcore.com 
 MathCore AB                   phone: +4613 32 85 07 
 Wallenbergsgata 4             fax:         21 27 01
 SE-583 35 Linköping, Sweden   http://www.mathcore.com



> -----Original Message-----
> From: banerjee@wernicke.ccn.upenn.edu 
> [mailto:banerjee@wernicke.ccn.upenn.edu] 
> Sent: den 13 mars 2002 00:50
> To: gsl-discuss@sources.redhat.com
> Subject: Computing a Projection Matrix
> 
> 
> Folks,
> 
> I need to compute the projector:
> 
>                 T  -1 T
>        (I - A (A A)  A )
> 
> Here's how I am doing it:
> 
> Call gsl_blas_dgemm() to compute the matrix product:
> 
>                  T
>                (A A)
> 
> Call gsl_permutation_calloc(), gsl_linalg_LU_decomp(), and
> gsl_linalg_LU_invert() to compute the matrix inverse:
> 
>                 T  -1    
>               (A A)  
> 
> Call gsl_blas_dgemm() to compute the matrix product:
> 
>                 T  -1 T
>               (A A)  A 
> 
> Finally call gsl_blas_dgemm() to compute the matrix product:
> 
>                 T  -1 T
>             A (A A)  A  
> 
> (I won't actually go to the trouble of subtracting the
> above matrix from the identity matrix I.)
> 
> Is this the best approach? Can I avoid computing the matrix
> inverse?
> 
> Thanks.
> 
> K. Banerjee
> 
>

Attachment: GNU Scientific Library -- Reference Manual - Linear Algebra.url
Description: GNU Scientific Library -- Reference Manual - Linear Algebra.url