Algorithme valeurs propres

[mar31]

Bonjour,
J'ai besoin de calculer toutes les valeurs propres d'une matrice. Pour cela, j'ai besoin de deux algorithmes que j'ai écrits : la méthode de la puissance qui calcule une valeur propre et la méthode de la déflation qui calcule une autre valeur propre à partir d'une valeur propre et d'un vecteur propre déjà trouvés.
Voici mon code écrit en C++:

Code: Tout sélectionner

//Fonction calculant la valeur propre et le vecteur propre d'une matrice pour un vecteur de départ et une précision dpnnée. Retourne la valeur propre.
double MatriceCarree::methodePuissances(double* vCour,double epsilon)
{
 
  double * vSuiv = new double[n];
  vSuiv = MatriceCarree::initVecteur(vSuiv,n);
  double * vPrime = new double[n];
  int nbIteration = 0;
  double lambda = 0;
  double signe = 1;
  int i = 0;
 
    if (nepsilon && nbIterationsetMat(k,l,(getMat(k,l) - getMat(i,l)*v[k]*(1/v[i])));
                }
            }
 
            //Calcul du vecteur propre et de la valeur propre de b
            mu = b->methodePuissances(vCour,EPS);
 
          //Calcul du vecteur propre de A
          for (int k=0;kafficher();
        int n=m->getDimension();
        double * v = new double[n]; MatriceCarree::initVecteur(v,n); v[0]=1; //Crée un vecteur de dimension n initialisé à 0 (sauf v[0])
        lambda = m->methodePuissances(v,10^(-6));
        for(int i=0;imethodeDeflation(v,lambda);
        system("PAUSE");
        return 0;
}

Les 2 premières valeurs propres sont bonnes mais ensuite le programme me trouve des valeurs propres identiques (ou presque) alternativement. Le résultat sera de la forme : {val1, val2, val1, val2...}.
Je n'arrive vraiment pas à régler ce problème donc si quelqu'un aurait la gentillesse de m'aiguiller ?
Merci

[fatal_error]

hello,

déjà une question évidente.
Quel algorithme utilises tu. Pe un pseudo code est-il plus avisé que du C++, qui plus est mal indenté...

Enfin attention aux tests du genre normeMax(vPrime)==0, parce que potentiellelement tu peux tomber sur du 1E-16 qui devrait etre considéré comme 0, donc du coup pe que ta condition est jamais vérifiée (sauf si normeMax tronque à 0 bien sur, mais comme on a pas sa définition...)

De même pour vCour[i]==0 et autres tests sur des valeurs calculées flottantes.

Enfin quelques remarques annexes:
- si tu fais du C++, utilises vector plutot que double[].
Ca t'éviteras de devoir gérer le delete (que tu n'as pas géré d'ailleurs..)
- prototype ton algo avec matlab/octave que directement le coder en C++, au moins pour etre sur que ton algo est correct (et c'est plus simple d'écrire du calcul avec octave qu'avec du C++)

[mar31]

Merci pour ta réponse.

Pour le C++ ce langage m'a été imposé et je ne connais pas matlab ou octave.

Les algorithmes sont les suivants :

• Méthode de la déflation :
  
  1. v0) ;
  2. v0, on dé;)nit :
     B = A - (1/v[i]) * v * A[i] où A[i] dénote le vecteur ligne formé des ´éléments de la i-ème ligne de A.
     
     => B et A ont les mêmes valeurs propres, sauf que B a 0 pour valeur propre au lieu de lambda.
     De plus, si w est un vecteur propre de B associé à mu, alors un vecteur propre associé à mu pour A est :
     v = (mu - lambda) * w + (1/v[i]) * (A[i] * w) * v

Pour les tests sur les valeurs flottantes en effet je vais essayer de corriger cela mais je ne pense pas que cela soit la source du problème ? Et pour gérer les delete en effet j'ai essayé mais j'ai eu des erreurs de pointeurs donc je vais essayer vector merci.

Aurait-tu une idée d'où se trouve mon erreur dans mon algorithme (je pense dans l'utilisation de la méthode de la déflation) ?

[fatal_error]

j'ai pas tout lu ni compris, mais

Code: Tout sélectionner

b->setMat(k,l,(getMat(k,l) - getMat(i,l)*v[k]*(1/v[i])));

ca m'a l'air faux car
j'intuite que setMat sette au coeff (k,l) la valeur getMat(k,l) - ...

Je pense que ta double boucle a pour but de setter chacun des coeffs de B.
et que tu utilises B = A - (1/v[i]) * v * A[i]
mais si tu regardes getMat(i,l)*v[k]*(1/v[i]))
ca ne correspond pas, car dans l'énoncé tu as A[i] un vecteur, et v un vecteur, alors que dans ton implem, tu as un coeff (v[k] ) et un autre coeff A(i,l)==getMat(i,l)

De manière générale, tu devrais t'écrire une petite méthode pour faire le produit de deux vecteurs.
Et également tester ta méthode de puissance pour etre sur que c'est déflation qui marche pas.

Je pense pas qu'il n'y ai que ca, mais déjà ca ca sent mauvais!

[mar31]

Je ne pense pas que cela soit une erreur : j'ai redéfini l'équation à la main pour qu'il soit plus facile à implanter algorithmiquement et je pense que ça correspond.
Et j'ai testé ma méthode des puissances qui a l'air de très bien fonctionner.

[fatal_error]

si t'es sûr de toi je peux pas t'aider plus.

Il faut voir si tu penses te tromper dans l'algorithme ou dans l'implémentation de l'algo.
Comme tu as l'air serein sur l'algorithme, ca veut dire que t'as glissé une coquille dans l'algo.

Généralement quand t'as ce genre de choses, ben tu regardes à coup de std::cout ou de gdb les différentes valeurs au fur et à mesure.

si j'ai rien à faire je te proposerai une deflation dans la soirée, mais c'est pas super car on saura pas à priori pourquoi ca marchait pas pour toi..

edit: cela dit pe que quelqu'un d'autre prendra la peine de tester ton code et de le debugger directement..

[fatal_error]

En partant depuis ce pdf

L'algo octave est assez immédiat:

Code: Tout sélectionner

clear;
%http://lmlm6-62.univ-lille1.fr/lml/perso/calzavarini/Enrico_Calzavarini_Home/MNI_files/C4-Valeur-propres-Enrico.pdf
function [lambda,v]=iteratedPower(A)
  z=zeros(size(A,1),1);
  zPrime=z;
  z(1)=1;
  count=0;
  while count++1e-15
    zPrime=z;
    z=A*z;
    z=z/norm(z,2);
  end
  lambda=z'*A*z/(z'*z);
  v=z;
%  isEigenValue=A*v-lambda*v
endfunction

% A needs to be defined positive and symmetric
function [P,D]=deflate(A)
  D=A*0;
  for i=1:size(A,1)
    [lambda,v]=iteratedPower(A);
    P(:,i)=v;
    D(i,i)=lambda;
    A=A-lambda*v*v';
  end
endfunction
A=[
3 2 3
2 6 7
3 7 11
];
[P,D]=deflate(A);
P
D
expectZeros = P*D*inv(P) - A


output:

Code: Tout sélectionner

P =

   0.251526   0.963929   0.087040
   0.556314  -0.217580   0.801981
   0.791991  -0.153298  -0.590974

D =

   16.86974    0.00000    0.00000
    0.00000    2.07145    0.00000
    0.00000    0.00000    1.05881

expectZeros =

   0.0000e+00   5.3291e-15  -3.5527e-15
   0.0000e+00   8.8818e-15  -6.2172e-15
  -1.7764e-15   1.1546e-14  -1.0658e-14




S'ensuit le code C++

Code: Tout sélectionner

#include
#include
#include
#include
class Matrix;
// Matrix and Vector are highly inefficient, just provided for no dependencies
class Vector{
  public:
  Vector(const std::vector &v):_v(v){}
  Vector(size_t i):_v(i){}
  double operator[](size_t i)const{return _v[i];}
  double& operator[](size_t i){return _v[i];}
  Matrix matProd(const Vector &other)const;
  Vector operator*(double scal)const;
  Vector operator/(double scal)const;
  Vector operator-(const Vector& other)const;
  double operator*(const Vector& other)const;
  size_t size()const{return _v.size();}
  std::string toString()const{
    std::stringstream oss;
    for(size_t i=0;i _v;
};

 
 
class Matrix{
  public:
  Matrix(const std::vector &v):_v(v), _width(sqrt(v.size())){}
  Matrix(size_t width):_v(width*width),_width(width){}
  double& at(size_t i, size_t j){return _v[i*_width+j];}
  double at(size_t i, size_t j)const{return _v[i*_width+j];}
  Vector getVector(size_t i)const{
    std::vector v(_width);
    for(size_t j=0;j result(v.size());
    for(size_t i=0;i v(_v.size());
    for(size_t i=0;i v(_v.size());
    for(size_t i=0;i _v;
  size_t _width;
};

Vector Vector::operator*(double scal)const{
  Vector v(_v);
  for(size_t i=0;ioperator*(1./scal);
}
double Vector::operator*(const Vector& other)const{
  double result=0;
  for(size_t i=0;i v(other.size());
  for(size_t i=0;i

struct EigenData{
  Matrix P;
  Matrix D;
  EigenData(size_t i):P(i),D(i){}
};
struct EigenCouple{
  double lambda;
  Vector v;
  EigenCouple(size_t i):v(i){}
};
double norm(const Vector& v,size_t dummy){
  if(dummy || true){
    return sqrt(v*v);
  }
  return 0;
}
void iteratedPower(const Matrix& A, struct EigenCouple& out){
  Vector z(A.width());
  z=z*0;
  Vector zPrime=z;
  z[0]=1;
  int count=0;
  while( count++1e-15){
    zPrime=z;
    z=A*z;
    z=z/norm(z,2);
  }
  out.lambda=(A*z)*z/(z*z);
  out.v=z;
}
void deflate(const Matrix& m, struct EigenData &out){
  Matrix A=m;
  out.D=out.D*0;
  struct EigenCouple eigCouple(A.width());
  for(size_t i=0;i v(values, values+sizeof(values)/sizeof(double));
  Matrix A(v);
  EigenData eig(A.size());
  deflate(A, eig);
  std::cout<<eig.D.toString()<<std::endl;
/*
16.8697 0 0
0 2.07145 0
0 0 1.05881
*/

  std::cout<<eig.P.toString()<<std::endl;
/*
0.251526 0.963929 0.0870397
0.556314 -0.21758 0.801981
0.791991 -0.153298 -0.590974
*/
  return 0;
}



[mar31]

:doh: Mille mercis je regarde ça !

[mar31]

J'ai une erreur dans mon "Vector.cpp": Pour la fonction

Code: Tout sélectionner

Matrix Vector::matProd(const Vector &other)const{
  Matrix result(_v.size());
  for(size_t i=0;i<_v.size();++i){
    Vector v_i = other*_v[i];
    result.setVectorAt(i, v_i);
  }
  return result;
}

A la compilation j'ai l'erreur suivante : 27|error: return type 'class Matrix' is incomplete|

[fatal_error]

ben inclus Matrix.hpp dans Vector.cpp
Et si t'as pas Matrix.hpp, ben je peux pas t'aider plus si tu montres pas l'arborescence de tes fichiers.

Mais c'est juste un problème de forward declaration, rien de bien palpitant :lol3:

[mar31]

Ah oui en effet je pensais que la déclaration anticipée de la classe "Matrix" dans le header de la classe "Vector" suffisait. Du coup il n'est plus utile de faire cette déclaration anticipée ? Ne risque-t-il pas d'y avoir d'imbrications infinies ?

En tous cas merci tu m'as beaucoup aidée je ne m'en sortais plus avec mes classes très "brouillones" je peux repartir avec des structures de classes beaucoup plus propres pour essayer de faire l'algorithme de départ. Car le but est justement de comparer l'algorithme de déflation que j'ai donné à celui d'un logiciel de type R.

[fatal_error]

Du coup il n'est plus utile de faire cette déclaration anticipée ? Ne risque-t-il pas d'y avoir d'imbrications infinies ?

Dans le header, soit tu inclues Matrix, soit tu fais une forward declaration.
Comme Vector dépend de Matrix et Matrix dépend de Vector, si tu inclues Matrix dans Vector.hpp, il faut que tu forward declare Vector dans le header de Matrix.hpp et que tu inclues Vector dans Matrix.cpp (ou vice versa)

Quant aux imbrications infinies, tu les gères avec les classiques ifndef ...

mes classes très "brouillones"
que j'ai donné à celui d'un logiciel de type R.

Attention, ce que je t'ai proposé, c'est intuitif, mais c'est pas optimisé du tout. Il y a des librairies qui font ca bien mieux que cette implem de 100 lignes...
De fait, tirer des conclusions par rapport à R, c'est pas super si tu te bases sur ton implem C++ (sauf si tu obtiens déjà mieux bien sûr... mais j'en doute )