Power iteration method implementation

kennethassogba · Jan 1, 2025 · 583e25e · 583e25e
1 parent ed7f6c2
commit 583e25e
Show file tree

Hide file tree

Showing 3 changed files with 34 additions and 26 deletions.
diff --git a/examples/main.cpp b/examples/main.cpp
@@ -31,10 +31,14 @@ int main() {
   // Test power iteration solver
   // 1st problem, vp : (5, 2)
   Eigen::MatrixXd A1(2,2);
+  Eigen::VectorXd  x1(2);
+  Eigen::VectorXd  y1(2);
+
   A1 << 4, 1,
         1, 3;
-  // PowerIteration PowIte(A);
-  // PowIte.solve(); 
+
+  PowerIteration PowIteA(A1);
+  PowIteA.solve(); 
 
   // 2nd problem, vp : (11.25, 7.34, 5.11, 4.30)
   Eigen::MatrixXd B1(4,4);
@@ -43,6 +47,9 @@ int main() {
          0, 1, 6, 1,
          0, 0, 1, 4;
 
+  PowerIteration PowIteB(B1);
+  PowIteB.solve(); 
+
   // 3rd problem
   // vp : (11.34, 11.28, 8.56, 6.68, 4.84, 2.20)
   Eigen::MatrixXd C1(6,6);
@@ -53,6 +60,9 @@ int main() {
          0,  0, 0, 1, 5, 1,
          0,  0, 0, 0, 1, 3;
 
+  PowerIteration PowIteC(C1);
+  PowIteC.solve(); 
+
   // Test linear solvers
   Eigen::Matrix2f A, b;
   A << 2, -1, -1, 3; 

diff --git a/src/demeter/solve/solver.cpp b/src/demeter/solve/solver.cpp
@@ -22,35 +22,33 @@ void Solver::solve() {
 void PowerIteration::solve() {
 
   size_t size;
+  int nb_ite;
   double lambda_k1;
   double lambda_k2;
   double error_eigvalue;
 
   // Initialization
-  lambda_k1 = 0;
+  nb_ite = 0;
+  lambda_k1 = 1;
   lambda_k2 = 0;
+  error_eigvalue = 10;
   size = A_.rows();
-  Eigen::ArrayXd x1 = Eigen::ArrayXd::Ones(size);
-  Eigen::ArrayXd x2 = Eigen::ArrayXd::Ones(size);
-
-  // Iterations
-  while (error_eigvalue >= tol_eigvalue_ ) {
-
-     // Vector update 
-     // x2 = A_ * x1;
-     // x2 = x2 / |x2|
-
-     // Value update 
-     //lambda_k1 = 1;     
-     //lambda_k2 = <x2,x1> / (<x1,x1>;
-     //x2 = x2 / |x2|
-
-     // Convergence test
-     error_eigvalue = lambda_k2-lambda_k1;
-     std::cout << "Lambda: " << lambda_k2 << std::endl;
-     } 
-
+  Eigen::VectorXd x1 = Eigen::ArrayXd::Random(size);
+  Eigen::VectorXd x2 = Eigen::ArrayXd::Ones(size);
+
+  while (error_eigvalue >= tol_eigvalue_) {
+    // Update
+    nb_ite++;
+    x2 = A_ * x1;
+    lambda_k2 = x2.dot(x1) / x1.dot(x1);
+    x1 = x2;
+
+    // Convergence test
+    error_eigvalue = std::abs((lambda_k2 - lambda_k1) / lambda_k2);
+    lambda_k1 = lambda_k2;
+  }
+
   std::cout << "Eigenvalue: " << lambda_k2 << std::endl;
-
+  std::cout << "Iterations: " << nb_ite   << std::endl;
 }
 }  // namespace Demeter
diff --git a/src/demeter/solve/solver.hpp b/src/demeter/solve/solver.hpp
@@ -65,9 +65,9 @@ class PowerIteration {
   void solve();
 
  private:
-  
+
   double lambda_;
-  double tol_eigvector_ = 1e-5;
+  double tol_eigvector_ = 1e-10;
   double tol_eigvalue_  = 1e-5;
 
   Eigen::MatrixXd A_;