a = sqrt(2);
G = [
1 0 0 1 0 0 1 0 0;
0 1 0 0 1 0 0 1 0;
0 0 1 0 0 1 0 0 1;
1 1 1 0 0 0 0 0 0;
0 0 0 1 1 1 0 0 0;
0 0 0 0 0 0 1 1 1;
a 0 0 0 a 0 0 0 a;
0 0 0 0 0 0 0 0 a;]

// (3)
[U, S, V] = svd(G);

// (4)
r=rank(G); // rang = 7

// (5)
// Les deux dernieres colonnes de V sont les deux vecteurs
// du null space

Up = U(:,1:r);
Sp = S(1:r, 1:r);
Vp = V(:,1:r);
V0 = V(:,r+1:r+2);
Ginvg = Vp * inv(Sp) * Up'

// (6)
M1 = V(:, r+1); M2 = V(:, r+2);
M1 = matrix(M1, 3, 3); M2 = matrix(M2, 3, 3);

// (7)
Rm = Vp*Vp';
Rmdiag = matrix(diag(Rm), 3, 3);

// (8)
Rd = Up*Up';
sumcols = sum(Rd, 'c');
sumrows = sum(Rd, 'r');

//calcul de l'inverse généralisé
Ginvg = Vp*inv(Sp)*Up';

//test d'un modèle ( 0 0 0 0 1 0 0 0 0 )
M = [0 0 0 0 1 0 0 0 0];
D = G*M';
Mapprox = Ginvg * D ;

// (10) ajout de bruit gaussien
Dnoisy = D+ rand(8, 1, 'uniform')*0.1;

Mapprox2 = Ginvg*Dnoisy ;