// Minimisation de fonctionnelles quadratique

exec('optim_macros.sci'); 

N=10;

A=defpos_sym(N,2,3);
b=rand(N,1);

// le pas de gradient maximal 

rho_max = rmax(A)
// méthode de gradient 
x0= rand(N,1);
rho= rho_max/10;

// iterations de gradient 

[c,xn,fxn]=gradient(x0,f,df,rho,1.e-3,500,A,b);

// Calcul par résolution d'un système linéaire 

xs= A\(-b) 
fxs=f(xs,A,b);

// test des resultats 

plot(c-fxs) 

// Rendre la fonction gradient plus générique 

[c,xn,fxn]=gradient_1(x0,f,df,rho,1.e-3,500,A,b);

plot(c-fxs) 

// Rajout d'une contrainte
// On rajoute une contrainte. de la forme T*x = 0 
// T est une matrice NxN aléatoire dont le noyau est de dimension N/2 
// ---------------------------

T=contrainte(N,N/2) ; 

// verifions 
rank(T)

// la contrainte est de la forme T*x = 0 
// on utilise une penalisation 

eps = 0.1 ;
rho_max = rmax(A+ (1/eps)*T'*T); 

// méthode de gradient 

x0= rand(N,1);
rho= rho_max/10;

[c,xn,fxn]=gradient_1(x0,fc,dfc,rho,1.e-4,500,A,b,T,eps);

//  Un peu d'algèbre linéaire
//  pour résoudre le problème avec contraintes 
// ---------------------------- 

[W,rk]=colcomp(T);
// les rk premieres colonnes de W donne une base du noyau de T 
// x = W(:,1:rk)* y 

W1= W(:,1:rk); 

// On minimise sans contraintes avec un nouveau A 

A1= W1'*A*W1; 
b1 = W1'*b; 

// contraintes sur le pas de gradient 

rho_max  = rmax(A1); 

// méthode de gradient 
N1= size(A1,1) 
y0= rand(N1,1); yn=y0;
rho= rho_max/10;

[c,yn,fyn]=gradient_1(x0,f,df,rho,1.e-8,500,A,b);

ys= A1\(-b1) 
f(ys,A1,b1)
// on revient à xs 
xs= W1*ys ;
// norm(T*xs) 

// On change les contraintes pour rendres les lignes linéairement inépendantes 

T1= contrainte_plein_rang(T);
X= [ (A+A')/2, T1'; T1, zeros(N/2,N/2)] \ [ -b ; zeros(N/2,1)];

norm(X(1:N) -xs) 

// avec T 

X2 = [ (A+A')/2, T'; T, zeros(N,N)] \ [ -b ; zeros(N,1)];

// le systeme est singulier, la partie en lambda est 
// non unique. 

norm(X(1:N) -xs)