Lesson of 26 February 2010¶
Lesson file in a zip file
The header of Greville algorithm for pseudoinverse
File greville.hh
1#ifndef GREVILLE_HH
2#define GREVILLE_HH
3
4namespace Greville {
5 typedef double valueType ;
6 typedef unsigned indexType ;
7
8 void
9 pseudoinverse( indexType nRow, // number of rows of A
10 indexType nCol, // number of columns of A
11 valueType const A[], // the marix A
12 indexType ldA, // row dimension of A
13 valueType iA[], // computed pseudoinverse
14 indexType ldiA, // row dimension of A
15 valueType work[] ) ;
16}
17
18#endif
The code of Greville algorithm for pseudoinverse
File greville.cc
1#include "greville.hh"
2#include "alglin.hh"
3
4namespace Greville {
5
6 void
7 pseudoinverse( indexType nRow, // number of rows of A
8 indexType nCol, // number of columns of A
9 valueType const A[], // the marix A
10 indexType ldA, // row dimension of A
11 valueType Aplus[], // computed pseudoinverse
12 indexType ldAplus, // row dimension of A
13 valueType work[] ) {
14
15 valueType epsi = 1E-10 ;
16 valueType *dk = work ;
17 valueType *ck = dk + nCol ;
18
19 // init recursion, evalute pseudo inverse first column
20 alglin::copy( nRow, A, 1, Aplus, ldAplus ) ; // copy first column of A to the first row of iA
21 valueType len = alglin::nrm2( nRow, A, 1 ) ; // norm of first column of A
22 alglin::scal( nRow, 1/len/len, Aplus, ldAplus ) ;
23
24 // Greville recursion
25 for ( indexType k = 1 ; k < nCol ; ++k ) {
26 valueType const *ak = A + k*ldA ; // k-th column of A
27 valueType *bk = Aplus + k ; // k-th row of A^+
28 // evaluate dk = A^+_{k-1} * ak
29 alglin::gemv( "No Transpose", k, nRow,
30 1.0, Aplus, ldAplus,
31 ak, 1,
32 0.0, dk, 1 ) ;
33 // evaluate ck = ak - A_{k-1} * dk
34 alglin::copy( nRow, ak, 1, ck, 1 );
35 alglin::gemv( "No Transpose", nRow, k,
36 -1.0, A, ldA,
37 dk, 1,
38 1.0, ck, 1 ) ;
39 valueType lenck = alglin::nrm2( nRow, ck, 1 ) ;
40 if ( lenck < epsi ) {
41 // bk = (A^+)^T dk / (1+dk^2)
42 valueType bf = alglin::nrm2( k, dk, 1 ) ;
43 alglin::gemv( "Transpose", k, nRow,
44 1.0/(1+bf*bf), Aplus, ldAplus, dk, 1,
45 0.0, bk, ldAplus ) ;
46 } else {
47 // bk = ck^+
48 alglin::copy( nRow, ck, 1, bk, ldAplus ) ;
49 alglin::scal( nRow, 1/lenck/lenck, bk, ldAplus ) ;
50 }
51 // A^+_{k-1} -= dk^T*bk
52 alglin::ger( k, nRow, -1.0, dk, 1, bk, ldAplus, Aplus, ldAplus ) ;
53 }
54 }
55}
A driver code
File testGreville.cc
1#include <cmath>
2#include <iostream>
3#include <iomanip>
4#include <vector>
5#include <algorithm>
6#include "greville.hh"
7#include "alglin.hh"
8
9using namespace std ;
10using namespace Greville ;
11
12#define NROW 3
13#define NCOL 2
14
15int
16main() {
17
18 valueType iA[NROW*NCOL], work[NROW+NCOL], B[NROW*NROW+NCOL*NCOL] ;
19 valueType A[NROW*NCOL] = { 2, 1, 0, 1, 1, 1 } ;
20
21 pseudoinverse( NROW, NCOL, A, NROW, iA, NCOL, work ) ;
22
23 alglin::gemm( "No Transpose", "No Transpose", NROW, NROW, NCOL, 1.0, A, NROW, iA, NCOL, 0.0, B, NROW ) ;
24 for ( indexType i = 0 ; i < NROW ; ++i ) {
25 for ( indexType j = 0 ; j < NROW ; ++j )
26 cout << "(A*iA)(" << i << "," << j << ") = " << B[i+NROW*j] << '\n' ;
27 cout << '\n' ;
28 }
29
30 alglin::gemm( "No Transpose", "No Transpose", NCOL, NCOL, NROW, 1.0, iA, NCOL, A, NROW, 0.0, B, NCOL ) ;
31 for ( indexType i = 0 ; i < NCOL ; ++i ) {
32 for ( indexType j = 0 ; j < NCOL ; ++j )
33 cout << "(iA*A)(" << i << "," << j << ") = " << B[i+NCOL*j] << '\n' ;
34 cout << '\n' ;
35 }
36
37}