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);


// On voit que la somme des éléments sur les colonnes, diagonales et lignes
// vaut 0 : cela s'explique par le fait que Dpred = GMp + GM0
// On voit bien qu'ici M0 vaut 0, et donc ne contribue en rien dans le modèle.

// (7)
Rm = Vp*Vp';
sumcolsRm = sum(Rm, 'c'); // =1
sumrowsRm = sum(Rm, 'r'); // =1
// D_true = GM_true
// M_est = inv(G_t G) G_t D_obs
// Quand on suppose que D_obs = D_true on a donc
// M_est = inv(G_t G) G_t G M_true
// M_est = Rm M_true
// Si M_est = M_true, alors la matrice Rm est diagonale.
// Dans notre exemple, ce n'est pas le cas et on a une dispersion.
// Mais la somme des colonnes est égale à 1, c'est logique.

Rmdiag = matrix(diag(Rm), 3, 3);
// On voit que, pour cette matrice, seule le modèle m33 est parfaitement défini
// cela s'explique par le fait que le dernier rayon passe UNIQUEMENT à travers
// M33, ce qui permet de le définir de manière parfaite.

// (8)
Rd = Up*Up';
sumcolsRd = sum(Rd, 'c'); // = 1
sumrowsRd = sum(Rd, 'r'); // = 1
// La somme des colonnes et des lignes vaut 1 :
// Dans le genre : La dispersion sur les données, au même titre que Rm
// fait qu'on a pas une diagonale parfaite :
// d_obs != d_true et on a donc une dispersion autour de la diagonale

//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 ;