Pseudo inverse of complex matrix

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

Pseudo inverse of complex matrix

Postby nawras » Thu Sep 22, 2011 11:30 am

Dear LAPACK users,

I need some help please. I need a subroutine (FORTRAN 77) to find the Pseudo-inverse of a complex rectangular matrix. I am wondering if anyone can guide me.

Your time and help are very appreciated.

Regards

Nawras
nawras
 
Posts: 4
Joined: Thu Sep 22, 2011 10:39 am

Re: Pseudo inverse of complex matrix

Postby admin » Fri Sep 23, 2011 11:46 am

You can have a look at our forum. There are a few discussions on how to compute the pseudoinverse using LAPACK.
See for example viewtopic.php?f=2&t=160
Julien.
admin
Site Admin
 
Posts: 488
Joined: Wed Dec 08, 2004 7:07 pm

Re: Pseudo inverse of complex matrix

Postby nawras » Fri Sep 23, 2011 6:15 pm

But all discussions were made for real matrices and I am interested in the Pseudoinverse of complex matrix.
Thanks very much.
nawras
 
Posts: 4
Joined: Thu Sep 22, 2011 10:39 am

Re: Pseudo inverse of complex matrix

Postby Julien Langou » Mon Sep 26, 2011 11:33 am

Hi Nawras, Henc gave it a shot and wrote a quick example of computation of the pseudoinverse. First he computes the SVD of A. ( A = U S V^H .) Then he computes the pseudoinverse of A with A = V * S^(-1) * U^H. The codes are in Fortran and in C using double complex precision. See below. Henc is using DGESVD for the SVD, you can also use DGESDD. Henc is assuming that the matrix is full rank. (So all singular values of A are different from zero so the 1/S(i) is safe.) If some singular values are zero (or very small), (then A is not numerically nonsingular), and you should set 1/S(i) to zero in Henc's code. In the case of A full rank, another way (faster) is to go with a QR factorization instead of an SVD. Cheers, Julien.

{ Thanks Henc! }

Code: Select all
   Program PseudoInverse

   Implicit none
   
   external ZLANGE
   double precision ZLANGE
   
   integer i, j, M, N, K, L, LWORK, INFO
   parameter (M=15)
   parameter (N=10)

   parameter (K = MIN(M,N))
   parameter (L = MAX(M,N))

   parameter (LWORK = MAX(1,2*K+L))

   double complex, dimension(M,N) :: A1, A2, SIGMA
   double complex, dimension(N,M) :: PINV
   double complex, dimension(M,K) :: U
   double complex, dimension(K,N) :: VT
   double complex, dimension(N,N) :: BUFF
   double complex, dimension(LWORK) :: WORK
   double precision, dimension(5*K) :: RWORK
   double precision, dimension(K) :: S
   integer, dimension(4) :: ISEED

   double precision :: normA, normAPA, normPAP

   data ISEED/0,0,0,1/

c  Fill A1 with random values and copy into A2
   call ZLARNV( 1, ISEED, M*N, A1 )
   do i=1,M
      do j=1,N
         A2(i,j) = A1(i,j)
      end do
   end do

c  Compute the SVD of A1
   call ZGESVD( 'S', 'S', M, N, A1, M, S, U, M, VT, K, WORK, LWORK,
     $   RWORK, INFO)

c  Compute PINV = VT**T * SIGMA * U**T in two steps
   do j = 1, K
      call ZSCAL( M, dcmplx(1 / S( j )), U( 1, j ), 1 )
   end do
   call ZGEMM( 'C', 'C', N, M, K, dcmplx(1.0), VT, K, U, M,
     $   dcmplx(0.0), PINV, N)

c   check the result
   normA = ZLANGE( 'F', M, N, A2, M, NULL() )
   call ZGEMM( 'N', 'N', N, N, M, dcmplx(1.0), PINV, N, A2, M,
     $ dcmplx(0.0), BUFF, N )

   call ZGEMM( 'N', 'N', M, N, N, dcmplx(-1.0), A2, M, BUFF, N,
     $ dcmplx(1.0), A2, M );
   normAPA = ZLANGE( 'F', M, N, A2, M, NULL() )

   call ZGEMM( 'N', 'N', N, M, N, dcmplx(-1.0), BUFF, N, PINV, N,
     $ dcmplx(1.0), PINV, N );
   normPAP = ZLANGE( 'F', N, M, PINV, N, NULL() )

   write(*,"(A, e10.4)") '|| A - A*P*A || = ', normAPA/normA
   write(*,"(A, e10.4)") '|| P - P*A*P || = ', normPAP/normA
   end


Code: Select all
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <complex.h>

int cprintmatrix( char *matname, int m, int n, double complex *A);
int dprintmatrix( char *matname, int m, int n, double *A);

void zlarnv_( int* idist, int* iseed, int* n, double complex* x);
void zgesvd_( char* jobu, char* jobvt, int* m, int* n,
                    double complex* a, int* lda, double* s,
                    double complex* u, int* ldu,
                    double complex* vt, int* ldvt,
                    double complex* work, int* lwork,
                    double* rwork, int *info );
void zgemm_( char* transa, char* transb, int* m, int* n, int* k,
                    double complex* alpha, double complex* a, int* lda,
                    double complex* b, int* ldb, double complex* beta,
                    double complex* c, int* ldc );
double zlange_( char* norm, int* m, int* n, double complex* a,
                    int* lda, double* work );
void zscal_( int* n, double complex* alpha, double complex* x, int* incx );

void print_help(){
   printf( "Usage: ./example_pinv   -h             : help (this text)\n" );
   printf( "                        -n <width>   : nb of columns\n" );
   printf( "                        -m <height>  : nb of rows\n" );
   printf( "                        -print       : output in matlab format\n" );
}

int main(int argc, char **argv){

   int i, j;
   int IONE = 1;
   int M, N, K, L, LDA, LDAxN;
   int LDU, LDVT, LWORK, INFO;
   int out2matlab;
   char JOBU, JOBVT;
   char NOTRANS, TRANS, CONJTRANS;
   char NORM;
   double complex zpone, znone, zzero, tempS;
   double normA, normAPA, normPAP;

   M     = 15;
   N     = 10;
   out2matlab = 0;

   for ( i = 1; i < argc; i++ ) {
      if( strcmp( argv[i], "-h" ) == 0 ) {print_help(); return EXIT_SUCCESS; };
      if( strcmp( argv[i], "-m" ) == 0 ) { M = atoi(argv[i+1]); i++; }
      if( strcmp( argv[i], "-n" ) == 0 ) { N = atoi(argv[i+1]); i++; }
      if( strcmp( argv[i], "-print" ) == 0 ) { out2matlab = 1; }
   }

   K     = M < N ? M : N;
   L     = M > N ? M : N;
   LDA   = M;
   LDAxN = LDA * N;
   LDU   = M;
   LDVT  = K;
   JOBU  = 'S';
   JOBVT = 'S';
   LWORK = 1 > 2*K+L ? 1 : 2*K+L;

   NOTRANS   = 'N';
   TRANS     = 'T';
   CONJTRANS = 'C';

   NORM = 'F';

   zpone = 1.0; znone = -1.0; zzero = 0.0;

   int            *ISEED  = ( int            * )malloc(        4 * sizeof( int            ));
   double complex *A1     = ( double complex * )malloc( LDA  * N * sizeof( double complex ));
   double complex *A2     = ( double complex * )malloc( LDA  * N * sizeof( double complex ));
   double complex *PINV   = ( double complex * )malloc( LDA  * N * sizeof( double complex ));
   double complex *BUFF   = ( double complex * )malloc( N    * N * sizeof( double complex ));
   double complex *U      = ( double complex * )malloc( LDU  * K * sizeof( double complex ));
   double complex *VT     = ( double complex * )malloc( LDVT * N * sizeof( double complex ));
   double complex *WORK   = ( double complex * )malloc(    LWORK * sizeof( double complex ));
   double         *S      = ( double         * )malloc(        K * sizeof( double         ));
   double         *RWORK  = ( double         * )malloc( 5    * K * sizeof( double         ));

   ISEED[0] = 0; ISEED[1] = 0; ISEED[2] = 0; ISEED[3] = 1;
   /* Check if unable to allocate memory */
   if ((!A1)||(!A2)||(!PINV)){
      printf("Out of Memory \n ");
      exit(0);
   }
   if ((!U)||(!VT)||(!WORK)||(!S)||(!RWORK)){
      printf("Out of Memory \n ");
      exit(0);
   }

   /* Initialize A1 and A2 */
   zlarnv_( &IONE, ISEED, &LDAxN, A1);
   for (i = 0; i < M; i++)
      for (j = 0; j < N; j++)
         A2[LDA*j+i] = A1[LDA*j+i] ;

   /* Compute the SVD of A1 */
   zgesvd_( &JOBU, &JOBVT, &M, &N, A1, &LDA, S, U, &LDU, VT, &LDVT, WORK, &LWORK, RWORK, &INFO );

   /* Compute the pseudo inverse */
   for (i = 0; i < K; i++){
      tempS = (double complex)(1 / S[i]);
      zscal_( &M, &tempS, &(U[i*LDU]), &IONE );
   }
   zgemm_( &CONJTRANS, &CONJTRANS, &N, &M, &K, &zpone, VT, &LDVT, U, &M, &zzero, PINV, &N );

   if (out2matlab){
      printf("clear;\n");
      cprintmatrix("A1",LDA,N,A1);
      cprintmatrix("A2",LDA,N,A2);
      cprintmatrix("US",LDU,K,U);
      cprintmatrix("VT",LDVT,N,VT);
      dprintmatrix("S",K,1,S);
      cprintmatrix("P",N,M,PINV);
      printf("PINV = VT'*US';\n");
      printf("fprintf('|| A - A*pinv(A)*A || = %%1.4e\\n', norm(A2 - A2*PINV*A2,'fro'))\n");
      printf("fprintf('|| pinv(A) - pinv(A)*A*pinv(A) || = %%1.4e\\n',norm(PINV-PINV*A2*PINV,'fro'))\n");
      printf("fprintf('|| A - A*P*A || = %%1.4e\\n', norm(A2 - A2*P*A2,'fro'))\n");
      printf("fprintf('|| P - P*A*P || = %%1.4e\\n',norm(P-P*A2*P,'fro'))\n");
   }

   /* check the result */
   normA = zlange_( &NORM, &M, &N, A2, &M, NULL );
   zgemm_( &NOTRANS, &NOTRANS, &N, &N, &M, &zpone, PINV, &N, A2, &LDA, &zzero, BUFF, &N );

   /* || A - A * pinv(A) * A || */
   zgemm_( &NOTRANS, &NOTRANS, &M, &N, &N, &znone, A2, &LDA, BUFF, &N, &zpone, A2, &LDA );
   normAPA = zlange_( &NORM, &M, &N, A2, &M, NULL );

   /* || pinv(A) - pinv(A) * A * pinv(A) || */
   zgemm_( &NOTRANS, &NOTRANS, &N, &M, &N, &znone, BUFF, &N, PINV, &N, &zpone, PINV, &N );
   normPAP = zlange_( &NORM, &N, &M, PINV, &N, NULL );

   printf( "%% || A - A*pinv(A)*A || / || A ||             = %1.4e\n", normAPA / normA );
   printf( "%% || pinv(A) - pinv(A)*A*pinv(A) || / || A || = %1.4e\n", normPAP / normA );

   free( A1    );
   free( A2    );
   free( PINV  );
   free( BUFF   );
   free( U     );
   free( VT    );
   free( WORK  );
   free( S     );
   free( RWORK );

   return 0;

}

int dprintmatrix( char *matname, int m, int n, double *A){
   int i,j;

   printf("%s = [\n", matname);
   for( i = 0; i < m; i++){
      for( j = 0; j < n; j++ )
         printf("%1.16e ",A[i+j*m]);
      printf("\n");
   }
   printf("]; \n");

   return 0;
}

int cprintmatrix( char *matname, int m, int n, double complex *A){
   int i,j;

   printf("%s = [\n", matname);
   for( i = 0; i < m; i++){
      for( j = 0; j < n; j++ )
         printf("%1.16e + %1.16ei ",creal(A[i+j*m]),cimag(A[i+j*m]));
      printf("\n");
   }
   printf("]; \n");

   return 0;
}
Julien Langou
 
Posts: 734
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Pseudo inverse of complex matrix

Postby nawras » Mon Sep 26, 2011 12:21 pm

Thanks very much for Henc and also for Julien. I will try it soon!!
nawras
 
Posts: 4
Joined: Thu Sep 22, 2011 10:39 am

Re: Pseudo inverse of complex matrix

Postby MM2011 » Wed Dec 28, 2011 3:35 am

What are the size limits for matrices that LAPACK can handle? Can the C code provided in this thread calculate the pseudoinverse of a 300k X 300k matrix?
MM2011
 
Posts: 2
Joined: Tue Dec 27, 2011 8:52 pm


Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 1 guest