load(lapack)$ fpprintprec : 4;

G: matrix( [0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 1, 0], [0, 0, 1, 1, 0, 2 ] );

 dgesvd (G, true, true);

rank(G);

load (eigen)$

gramschmidt (G);

GR: -matrix( [0, 1, 0, 1], [1, 0, 1, 0], [0, 0, 0, 2 ] );

dgesvd (GR, true, true);

GR: -matrix( [1, 1], [ 0,  0], [0, 0 ] );

dgesvd (GR, true, true);

ejemplo

GR: matrix(  [ 0,  0, -1], [-1, -1, 0], [0, 0, 0 ] );

[d, U, VT] : dgesvd (GR, true, true);

d;

Sigma: genmatrix(lambda ([i, j], if  i=j  then d[i] else 0), 3, 3);

iSigma: genmatrix(lambda ([i, j], if  i=j  then ( if is(equal(d[i],0)) then 0 else 1/d[i] ) else 0), 3, 3);

U . Sigma . VT;

pINVc: transpose(VT) . iSigma . transpose( U );

pINVc . GR;

test U, Sigma, VT ok

V: transpose(VT);

V . VT;

VT . V;

test V, VT ok

Ur: submatrix( U, 3);

transpose( Ur ) . Ur;

test Ur ok

Sigmar: submatrix( 3, Sigma, 3);

iSigmar: invert( Sigmar);

iSigmar . Sigmar;

test iSigmar, Sigmar ok

tVr: submatrix( 3, VT);

Vr: transpose( tVr);

tVr . transpose( tVr);

Vr : submatrix(  V, 3);

 transpose( Vr) . Vr;

Ur . Sigmar . transpose( Vr );

test Ur.Sigmar.tVr ok

tUr: transpose( Ur);

tUr . Ur;

 transpose( Vr) . Vr;

pINV: Vr . iSigmar. tUr;

tVr . transpose (iSigmar . tUr . GR );

pINV . GR;

GR . pINV;

axiomas de la pseudoinversa

GR . pINV . GR - GR;

I: ok

pINV . GR . pINV   - pINV;

II: ok

transpose( GR . pINV ) - GR . pINV;

III: ok

transpose(  pINV . GR ) - pINV . GR;

IV: ok

pINVc - pINV;

pINV;

submatrix(V, 1, 2);

 R: F * transpose( matrix( [1, 1, 0 ] ) );

-F * pINV . transpose( matrix( [1, 1, 0 ] ) ) - m*submatrix(V, 1, 2);

-F * pINV . transpose( matrix( [1, 1, 0 ] ) ) + (diagmatrix(3,1)-pINV . GR) . transpose( [a, b, c] );