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;
}
Posted: Thu Sep 22, 2011 11:30 am