Méthode de la puissance

Loly

Salut à tous
Voici le sujet d'un tp que je dois réaliser, sous Matlab.

Ce que j'ai fait : (bien que les tests effectués ne retournent pas ce que je souhaite obtenir)
Des idées... ? HELP

function [lambda, x, niter]=puissance1(S,x0,kmax,tol)
     %init
     x=x0;
     lambda=0;
     niter=1;
     err=tol+1;
     %Boucle
     while err>tol && niter<kmax
               z=S*x;
               lambda=(x')*z;
               x=z/norm(z,1);
               niter=niter+1;
      end
end

tatof

Pour commencer tu fais la même erreur que dans tes autres sujets, la ligne $\text{err}=\text{tol}+1$ ne sert à rien car elle n'est pas dans la boucle while. De plus, l'expression pour l'erreur ne peut pas être celle là (l'erreur doit normalement être calculée sur la valeur propre). Autre chose, dans ton code tu écris la normalisation après le calcul de la valeur propre et du vecteur propre. Or cela doit être fait tout de suite après le début de la boucle.

Voici les codes avec la fonction et un fichier test:

function [lambda,x,niter]=puissance1(S,x0,kmax,tol)

%% Entrées
% S: matrice stochastique
% x0: vecteur d'initialisation
% kmax: nombre d'itération maximale
% tol: seuil de l'erreur

%% Sorties
% lambda: approximation de la plus grande valeur propre en module
% x: vecteur propre associé à lambda
% niter: nombre d'itération réalisée si la condition sur l'erreur est atteinte

%%
b=x0;
lambda=0;
for i=1:kmax
    lambda_ac=lambda;
    x=b/norm(b,1); %normalisation
    b=S*x;
    lambda=x'*b;
    
    err=abs(lambda-lambda_ac); %expression de l'erreur
    if (err<tol)
        break;
    end
end
niter=i;
end %fin fonction

Fichier test:

%test algo puissance
clear all;close all;clc;

%% Paramètres
S=[0.5 0.3 0.2;0.2 0.8 0;0.3 0.3 0.4];
x0=ones(size(S,1),1);
kmax=10;
tol=1e-3;

%% Algo puissance
[lambda,x,niter]=puissance1(S,x0,kmax,tol);

Loly

En effet, je n'ai pas fait attention à cette boucle... Je vais essayer de comprendre l'idée.
Par ailleurs, saurais-tu construire une fonction $L=matriceL(n)$ définie par la photo ci-jointe. I.e la matrice $L$ sera creuse, et de grande taille.

tatof

Je suppose que tu étudies l'algo PageRank. Pour ta question :

function [L]=matriceL(n,seuil)

%% Entrée
% n: taille de la matrice L
% seuil: variable dans [0,1] permettant d'avoir une matrice creuse. 
% Plus seuil proche de 1, plus L aura des 0

%% Sortie
% L: matrice creuse

%%
L=(rand(n)>seuil);

end %fin fonction

Loly

Oui c'est bien ça.
Ne devrais-je pas utiliser la fonction sprand ?

tatof

Je ne connais pas la commande sprand. Cependant, ma fonction te donne une matrice creuse qui dépend du seuil fixé. Lorsque ce seuil est proche de $1$, tu peux voir que la matrice a beaucoup de $0$.

tatof

J'ai compris ton problème. En fait, c'est une question de mémoire. Dans ma fonction je garde la matrice entière (zéro compris) en mémoire, ce qui n'est pas efficace si $n$ est très grand. Lorsque c'est le cas, il vaut mieux utiliser le code ci-dessous qui permet d'exploiter la spécificité de la matrice creuse en ne gardant que les indices des éléments non nuls. De plus, tu peux également remarquer avec ce nouveau code que si $n$ est petit par exemple $n=4$. Nous avons des répétitions dans les liens que tu peux voir par les valeurs de la matrice $A$ du code. Cela justifie le fait de prendre $n$ assez grand, par exemple $n=20$.
Code.

function [L]=connec(n,nb_con,i,j)

%% Entrée
% n: nombre de site web
% nb_con: nombre de connexion entre les sites
% (i,j): liens entre les sites web

%% Sortie
% L: matrice de connectivité

%% Vérification
if (nb_con>n^2)
    disp('Erreur dans nb_con'); %le nombre de connexion ne peut pas être supérieur au nombre d'élement de la matrice
end

%% Matrice de connectivité
L=sparse(i,j,1,n,n);

end %fin fonction

Fichier test.

% matrice creuse
clear all;close all;clc;

%% Paramètres
n=20; 
nb_con=9; 
i=round(1+rand(1,nb_con)*(n-1));
j=round(1+rand(1,nb_con)*(n-1)); 
A=[i;j]; %matrice permettant de voir les répétitions de liaison

%% 
[L]=connec(n,nb_con,i,j);
L2=full(L); %affichage de la matrice complète

Loly

Ton code me parait un peu trop complexe pour le coup. Dans le sujet, la seule variable d'entrée est n. Je pense qu'il faut passer par plusieurs boucles for, peut-être .. Ne t'embête pas plus. Ce n'est pas grave.
Voici la question sur la régression linéaire en lien avec la factorisation QR, en pdf, si tu as le temps d'y jeter un œil.
L'utilisation de la fonction polyfit, de Matlab, est autorisée.

[Merci d'écrire les mots en entier. :-) AD]

tatof

Voici un autre code pour la matrice creuse avec $1$ argument en entrée. J'ai mis les autres variables dans le code. J'ai aussi ajouté une ligne permettant de ne pas avoir de doublon entre les sites.

function [L,nb_con,A]=connec_v2(n)

%% Entrée
% n: nombre de site web

%% Sorties
% L: matrice de connectivité
% nb_con: nombre de connexion entre les sites
% A=[i;j]: matrice représentant les liaisons entre les sites web

%% Création du réseau 
nb_con=round(rand(1)*(n^2-1))+1;
i=round(1+rand(1,nb_con)*(n-1));
j=round(1+rand(1,nb_con)*(n-1)); %(i,j): liens entre les sites web
A=[i;j];
[~,ind]=unique(A','rows','first'); %on supprime les colonnes identiques de A
A=A(:,sort(ind));

%% Matrice de connectivité
L=sparse(A(1,:),A(2,:),1,n,n);

end %fin fonction

Loly

Pas de problème, je vais voir ce que je peux faire avec cette fonction.
As-tu regardé la question du pdf ? Je ne sais pas trop par où commencer...
Je tenais à préciser que dans ce pdf, t1 correspond à la deuxième colonne du tableau et t2 à la troisième colonne. Si cela peut aider.

tatof

Pour ce type de question, il vaut mieux demander dans la rubrique statistique. Personnellement, je vois pas trop ce qu'il faut faire.

tatof

Pour ton problème, je pense que tu cherches une fonction polynomiale qui passe à peu près par tous les points de données. Dans le code ci-dessous qui correspond à la Table 2 méthode Jacobi, tu peux remarquer que plus la variable deg est grande, plus la courbe correspondante sera lisse.

Code:

%test interpolation
clear all;close all;clc;

%% Paramètres
x=10:2:26;
y=[0.159 0.296 0.524 0.865 1.351 2.016 2.892 4.030 5.495];
deg=3; %degré du polynôme d'interpolation

%% Interpolation
p=polyfit(x,y,deg);
u=polyval(p,x);

plot(x,u);hold on;scatter(x,y);xlabel('Taille');ylabel('Temps');
legend('interpolation','données')

Loly

Oui c'est bien ça, je dois tester du degrés 1 au degrés 3.
Par contre, je suis censé utiliser la méthode qr...
Ce n'est pas utile ?

tatof

Lorsque je regarde la documentation de la fonction polyfit d'Octave, il est mentionné que le polynôme est trouvé au sens des moindres carrés. Du coup, je ne vois pas l'intérêt de la méthode QR.

Loly

Salut,

Je me permets de revenir vers toi concernant la matrice $L$ de la PageRank, l'utilisation des fonctions sparse, sprand, min et max de Matlab est autorisée. Il ne serait pas plus judicieux de passer par sparse et sprand ?
C'est peut être pour ça que ta fonction est plus complexe ? (Du moins elle me paraît complexe)

tatof

C'est exactement ce que je fais dans la fonction connec_v2. J'utilise sparse.

Loly

Salut
Je dois créer une fonction qui génère une matrice $A$ (voir fichier joint) en fonction des paramètres a,b,m et p.
Je bloque sur la fin, aurais-tu une idée ?

function A = mat(a,b,m,p)
  
  deltax = a/(m+1);
  deltay= b/(p+1);
  alpha = (2/(deltax)^2) + (2/(deltay)^2) ;
  beta = -(1/(deltax)^2);
  gamma = -(1/(deltay)^2);
  v = ones(m,1).*alpha;
  u = ones(m-1,1).*beta;
  G = diag(v) + diag(u,1) + diag(u,-1); 
endfunction

tatof

Quelle est la taille de la matrice A?

Loly

Elle est de taille n

tatof

Voici ce que je te propose avec ma fonction matri_bloc. J'ai tout d'abord remplacé tes multiplications d'un scalaire par un vecteur colonne avec des $1$ par la fonction repmat. Puis, j'ai construit la matrice $A$ en partant de la matrice $G$ et en rajoutant les diagonales avec une somme. J'ai mis un fichier test pour vérifier. Tu peux voir que la construction correspond.

Fonction:

function [A]=matri_bloc(a,b,m,p)

%% Entrées
% a,b,m,p: nombre réels

%% Sortie
% A: matrice triadiagonale par bloc

%%
deltax=a/(m+1);
deltay=b/(p+1);
alpha=(2/(deltax)^2)+(2/(deltay)^2) ;
beta=-(1/(deltax)^2);
gamma=-(1/(deltay)^2);
v=repmat(alpha,1,m);
u=repmat(beta,1,m-1);

G=diag(v)+diag(u,1)+diag(u,-1);
G2=repmat({G},1,2);
A1=blkdiag(G2{:}); %bloc de la diagonale
A=diag(repmat(gamma,1,m),m)+A1+diag(repmat(gamma,1,m),-m);
end %fin fonction

Fichier test:

%test matrice bloc
clear all;close all;clc;

%% Données
a=1;
b=1;
m=3;
p=3;

%%
[A]=matri_bloc(a,b,m,p);

Loly

J'ai pu constater que ça marchait bien ! Merci.

Par ailleurs, on me demande ensuite de calculer la plus petite valeur propre $\lambda_1$ de cette matrice et son vecteur propre associé. Je propose d'utiliser la méthode de la puissance inverse mais je ne suis pas certain de mon code...

function[lambda,x,iter,L]=puissinv(A,mu,x0,maxiter,tol,kmax)
  %init
  x=x0;
  lambda=0;
  iter=1;
  L=zeros(maxiter,1);
  n=length(x0);
  I=eye(n,n);
  for i = 1:kmax
      z=(A-mu*I)./x
      lambda = (x').*z;
      x=z/norm(z,2);
      L(iter)=lambda;
      if (err<tol)
        break;
      endif
      if (iter>1)
        err=abs(L(iter-1)-L(iter));
      endif  
  endfor
  lambda=(1/lambda) + mu;
endfunction

Ensuite on me demande de représenter graphiquement ce vecteur associé, à l'aide de ce code. Les paramètres ne sont pas imposés.

xt = linspace(0,a,m+2); yt = linspace(0,b,p+2); [Xt,Yt]=meshgrid(xt,yt);
figure
Uplot = [zeros(1,m+2); zeros(p,1) U' zeros(p,1); zeros(1,m+2)];
surf(Xt,Yt,Uplot);
shading interp

tatof

Tu as déjà codé la fonction pour la puissance itérée. En regardant sur le net, j'ai vu que la puissance inverse revient à appliquer la puissance itérée avec la matrice $A^{-1}$.

Loly

Merci pour ton aide.