Computing a Projection Matrix

[banerjee at wernicke dot ccn dot upenn dot edu]

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