add shifted eigenvalue procedure and unit test

doc/software_manual.org

@@ -1035,22 +1035,22 @@ double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance,
   return lambda;
 }
 #+END_SRC
-*** ~least_dominant_eigenvalue~
+*** ~shift_inverse_power_eigenvalue~
 + Author: Elizabeth Hunt
 + Name: ~least_dominant_eigenvalue~
 + Location: ~src/eigen.c~
 + Input: a pointer to an invertible matrix ~m~, an initial eigenvector guess ~v~ (that is non
-  zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~tolerance~ and
+  zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~shift~ to act as the
- ~max_iterations~ that act as stop conditions
+  shifted \delta, and ~tolerance~ and ~max_iterations~ that act as stop conditions.
-+ Output: the least dominant eigenvalue with the lowest magnitude, approximated with the Inverse
++ Output: the eigenvalue closest to ~shift~ with the lowest magnitude closest to 0, approximated
-  Power Iteration Method
+  with the Inverse Power Iteration Method
 #+BEGIN_SRC c
-double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
-                                 double tolerance, size_t max_iterations) {
+                                      double shift, double tolerance,
+                                      size_t max_iterations) {
   assert(m->rows == m->cols);
   assert(m->rows == v->size);
 
-  double shift = 0.0;
   Matrix_double *m_c = copy_matrix(m);
   for (size_t y = 0; y < m_c->rows; ++y)
     m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift;
@@ -1065,22 +1065,38 @@ double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
     Array_double *normalized_eigenvector_2 =
         scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2));
     free_vector(eigenvector_2);
-    eigenvector_2 = normalized_eigenvector_2;
 
-    Array_double *mx = m_dot_v(m, eigenvector_2);
+    Array_double *mx = m_dot_v(m, normalized_eigenvector_2);
     double new_lambda =
-        v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2);
+        v_dot_v(mx, normalized_eigenvector_2) /
+        v_dot_v(normalized_eigenvector_2, normalized_eigenvector_2);
 
     error = fabs(new_lambda - lambda);
     lambda = new_lambda;
     free_vector(eigenvector_1);
-    eigenvector_1 = eigenvector_2;
+    eigenvector_1 = normalized_eigenvector_2;
   }
 
   return lambda;
 }
 #+END_SRC
 
+*** ~least_dominant_eigenvalue~
++ Author: Elizabeth Hunt
++ Name: ~least_dominant_eigenvalue~
++ Location: ~src/eigen.c~
++ Input: a pointer to an invertible matrix ~m~, an initial eigenvector guess ~v~ (that is non
+  zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~tolerance~ and
+ ~max_iterations~ that act as stop conditions.
++ Output: the least dominant eigenvalue with the lowest magnitude closest to 0, approximated
+  with the Inverse Power Iteration Method.
+#+BEGIN_SRC c
+double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+                                 double tolerance, size_t max_iterations) {
+  return shift_inverse_power_eigenvalue(m, v, 0.0, tolerance, max_iterations);
+}
+#+END_SRC
+
 *** ~leslie_matrix~
 + Author: Elizabeth Hunt
 + Name: ~leslie_matrix~

homeworks/hw-7.org

@@ -21,7 +21,23 @@ finds the least dominant eigenvalue on the matrix:
 0 & 2 & 6 
 \end{bmatrix}
 
-which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should produce $\sqrt{17}$.
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should thus produce $5 - \sqrt{17}$.
 
 See also the entry ~Eigen-Adjacent -> least_dominant_eigenvalue~ in the LIZFCM API
 documentation.
+* Question Four
+See ~UTEST(eigen, shifted_eigenvalue)~ in ~test/eigen.t.c~ which
+finds the least dominant eigenvalue on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should thus produce $2.0$.
+
+With the initial guess: $[0.5, 1.0, 0.75]$.
+
+See also the entry ~Eigen-Adjacent -> shift_inverse_power_eigenvalue~ in the LIZFCM API
+documentation.

inc/lizfcm.h

@@ -76,6 +76,9 @@ extern double fixed_point_secant_bisection_method(double (*f)(double),
 
 extern double dominant_eigenvalue(Matrix_double *m, Array_double *v,
                                   double tolerance, size_t max_iterations);
+extern double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
+                                             double shift, double tolerance,
+                                             size_t max_iterations);
 extern double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
                                         double tolerance,
                                         size_t max_iterations);

src/eigen.c

@@ -49,12 +49,12 @@ double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance,
   return lambda;
 }
 
-double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
-                                 double tolerance, size_t max_iterations) {
+                                      double shift, double tolerance,
+                                      size_t max_iterations) {
   assert(m->rows == m->cols);
   assert(m->rows == v->size);
 
-  double shift = 0.0;
   Matrix_double *m_c = copy_matrix(m);
   for (size_t y = 0; y < m_c->rows; ++y)
     m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift;
@@ -69,17 +69,22 @@ double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
     Array_double *normalized_eigenvector_2 =
         scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2));
     free_vector(eigenvector_2);
-    eigenvector_2 = normalized_eigenvector_2;
 
-    Array_double *mx = m_dot_v(m, eigenvector_2);
+    Array_double *mx = m_dot_v(m, normalized_eigenvector_2);
     double new_lambda =
-        v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2);
+        v_dot_v(mx, normalized_eigenvector_2) /
+        v_dot_v(normalized_eigenvector_2, normalized_eigenvector_2);
 
     error = fabs(new_lambda - lambda);
     lambda = new_lambda;
     free_vector(eigenvector_1);
-    eigenvector_1 = eigenvector_2;
+    eigenvector_1 = normalized_eigenvector_2;
   }
 
   return lambda;
 }
+
+double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+                                 double tolerance, size_t max_iterations) {
+  return shift_inverse_power_eigenvalue(m, v, 0.0, tolerance, max_iterations);
+}

test/eigen.t.c

@@ -1,50 +1,7 @@
 #include "lizfcm.test.h"
 
-UTEST(eigen, leslie_matrix) {
+Matrix_double *eigen_test_matrix() {
-  Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8});
+  // produces a matrix that has eigenvalues [5 + sqrt{17}, 2, 5 - sqrt{17}]
-  Array_double *survivor_ratios = InitArray(double, {0.8, 0.55});
-
-  Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0);
-  m->data[0]->data[0] = 0.0;
-  m->data[0]->data[1] = 1.5;
-  m->data[0]->data[2] = 0.8;
-  m->data[1]->data[0] = 0.8;
-  m->data[2]->data[1] = 0.55;
-
-  Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity);
-
-  EXPECT_TRUE(matrix_equal(leslie, m));
-
-  free_matrix(leslie);
-  free_matrix(m);
-  free_vector(felicity);
-  free_vector(survivor_ratios);
-}
-
-UTEST(eigen, leslie_matrix_dominant_eigenvalue) {
-  Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8});
-  Array_double *survivor_ratios = InitArray(double, {0.8, 0.55});
-  Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity);
-  Array_double *v_guess = InitArrayWithSize(double, 3, 2.0);
-  double tolerance = 0.0001;
-  uint64_t max_iterations = 64;
-
-  double expect_dominant_eigenvalue = 1.22005;
-
-  double approx_dominant_eigenvalue =
-      dominant_eigenvalue(leslie, v_guess, tolerance, max_iterations);
-
-  EXPECT_NEAR(expect_dominant_eigenvalue, approx_dominant_eigenvalue,
-              tolerance);
-
-  free_vector(v_guess);
-  free_vector(survivor_ratios);
-  free_vector(felicity);
-  free_matrix(leslie);
-}
-
-UTEST(eigen, least_dominant_eigenvalue) {
-
   Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0);
   m->data[0]->data[0] = 2.0;
   m->data[0]->data[1] = 2.0;
@@ -54,12 +11,17 @@ UTEST(eigen, least_dominant_eigenvalue) {
   m->data[1]->data[2] = 7.0;
   m->data[2]->data[1] = 2.0;
   m->data[2]->data[2] = 6.0;
+  return m;
+}
+
+UTEST(eigen, least_dominant_eigenvalue) {
+  Matrix_double *m = eigen_test_matrix();
 
   double expected_least_dominant_eigenvalue = 0.87689; // 5 - sqrt(17)
   double tolerance = 0.0001;
   uint64_t max_iterations = 64;
 
-  Array_double *v_guess = InitArrayWithSize(double, 3, 2.0);
+  Array_double *v_guess = InitArrayWithSize(double, 3, 1.0);
   double approx_least_dominant_eigenvalue =
       least_dominant_eigenvalue(m, v_guess, tolerance, max_iterations);
 
@@ -88,3 +50,63 @@ UTEST(eigen, dominant_eigenvalue) {
   free_matrix(m);
   free_vector(v_guess);
 }
+
+UTEST(eigen, shifted_eigenvalue) {
+  Matrix_double *m = eigen_test_matrix();
+
+  double least_dominant_eigenvalue = 0.87689; // 5 - sqrt{17}
+  double dominant_eigenvalue = 9.12311;       // 5 + sqrt{17}
+  double expected_middle_eigenvalue = 2.0;
+  double shift = (dominant_eigenvalue + least_dominant_eigenvalue) / 2.0;
+
+  double tolerance = 0.0001;
+  uint64_t max_iterations = 64;
+  Array_double *v_guess = InitArray(double, {0.5, 1.0, 0.75});
+
+  double approx_middle_eigenvalue = shift_inverse_power_eigenvalue(
+      m, v_guess, shift, tolerance, max_iterations);
+
+  EXPECT_NEAR(approx_middle_eigenvalue, expected_middle_eigenvalue, tolerance);
+}
+
+UTEST(eigen, leslie_matrix_dominant_eigenvalue) {
+  Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8});
+  Array_double *survivor_ratios = InitArray(double, {0.8, 0.55});
+  Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity);
+  Array_double *v_guess = InitArrayWithSize(double, 3, 2.0);
+  double tolerance = 0.0001;
+  uint64_t max_iterations = 64;
+
+  double expect_dominant_eigenvalue = 1.22005;
+
+  double approx_dominant_eigenvalue =
+      dominant_eigenvalue(leslie, v_guess, tolerance, max_iterations);
+
+  EXPECT_NEAR(expect_dominant_eigenvalue, approx_dominant_eigenvalue,
+              tolerance);
+
+  free_vector(v_guess);
+  free_vector(survivor_ratios);
+  free_vector(felicity);
+  free_matrix(leslie);
+}
+UTEST(eigen, leslie_matrix) {
+  Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8});
+  Array_double *survivor_ratios = InitArray(double, {0.8, 0.55});
+
+  Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0);
+  m->data[0]->data[0] = 0.0;
+  m->data[0]->data[1] = 1.5;
+  m->data[0]->data[2] = 0.8;
+  m->data[1]->data[0] = 0.8;
+  m->data[2]->data[1] = 0.55;
+
+  Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity);
+
+  EXPECT_TRUE(matrix_equal(leslie, m));
+
+  free_matrix(leslie);
+  free_matrix(m);
+  free_vector(felicity);
+  free_vector(survivor_ratios);
+}