The topic of eigenvalues and eigenvectors is extremely complex. So I won’t try to explain other than to say that if you have an n x n square matrix, then there will be n eigenvalues (single scalar values) and n eigenvectors, each with n elements. For example, the Wikipedia article on the topic gives a 3×3 example where:

A = [[2, 0, 0],
     [0, 3, 4],
     [0, 4, 9]]

and shows how to compute the three eigenvalues: (2.0, 1.0, 11.0) and the three eigenvectors: [1, 0, 0], [0, -2, 1], [0, 1, 2]. The order in which the eigenvalues are computed is arbitrary, but it’s common practice to order them from largest to smallest (which I didn’t do in the demo, but probably should have). The eigenvectors can be multiplied by any constant (including -1) and still be the eigenvectors.

In most scenarios, eigenvectors are normalized by dividing by their norm (vector length — square root of sum of elements squared) so that they have length = 1. For the example above, the normalized eigenvectors are [1, 0, 0], [0, 0.8944, -0.4472], [0, 0.4472, 0.8944].

The eigenvectors are calculated from the eigenvalues and the source matrix A. A few days ago I coded up a C# implementation that shows how to compute eigenvalues. See the post on Monday, Oct. 9 at https://jamesmccaffreyblog.com/2023/10/09/matrix-eigenvalues-from-scratch-using-csharp/. This new blog post shows how to use the eigenvalues to compute the associated eigenvectors.

It is possible to compute the eigenvalues using one algorithm, and then compute the eigenvectors from the eigenvalues using a second algorithm. But in my opinion, in most scenarios, it’s best to compute the eigenvalues and eigenvectors at the same time:

public static void Eigen(double[][] M,
  out double[] eigenVals, out double[][] eigenVecs)
{
  // compute eigenvalues and eigenvectors at the same time
 
  int n = M.Length;
  double[][] X = MatCopy(M);  // mat must be square
  double[][] Q; double[][] R;
  double[][] pq = MatIdentity(n);
  int maxCt = 10000;

  int ct = 0;
  while (ct "lt" maxCt)
  {
    MatDecomposeQR(X, out Q, out R, false);
    pq = MatProduct(pq, Q);
    X = MatProduct(R, Q);  // note order
    ++ct;

    if (MatIsUpperTri(X, 1.0e-8) == true)
      break;
  }

  // eigenvalues are diag elements of X
  double[] evals = new double[n];
  for (int i = 0; i "lt" n; ++i)
    evals[i] = X[i][i];

  // eigenvectors are columns of pq
  double[][] evecs = MatCopy(pq);

  eigenVals = evals;
  eigenVecs = evecs;
}

The Eigen() function is deceptive because there are hundreds of lines of helper code that aren’t visible. In particular, MatDecomposeQR() does most of the work but it’s complex and calls several non-trivial helper functions.

The Main() method of the demo program is:

using System;
using System.IO;

namespace Eigen
{
  internal class EigenProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nEigenvectors from scratch ");

      // en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
      double[][] A = MatCreate(3, 3);
      A[0] = new double[] { 2, 0, 0 };
      A[1] = new double[] { 0, 3, 4 };
      A[2] = new double[] { 0, 4, 9 };

      Console.WriteLine("\nSource matrix: ");
      MatShow(A, 1, 6);

      double[] eigVals;
      double[][] eigVecs;
      Eigen(A, out eigVals, out eigVecs);

      Console.WriteLine("\nEigenvalues: ");
      VecShow(eigVals, 4, 9);

      Console.WriteLine("\nEigenvectors (columns): ");
      MatShow(eigVecs, 4, 9);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main

    // all the functions and helper functions go in here
  } // Program class
} // ns

I enjoy implementing algorithms from scratch. The process gives me insights into how library versions work, which in turn allows me to use library functions more accurately/efficiently. In addition, implementing from scratch keeps my coding skills sharp — coding, like any skill, must be practiced.

Eigenvalues and eigenvectors are used in many applications, including the Google page rank algorithm. But you still have to get the search phrase correct. Left: A search for “JabaScript” instead of “JavaScript”. Right: A search for “machine leaning” instead of “machine learning”.

Demo code. Long! Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;
using System.IO;

namespace Eigen
{
  internal class EigenProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nEigenvectors from scratch ");

      // en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
      double[][] A = MatCreate(3, 3);
      A[0] = new double[] { 2, 0, 0 };
      A[1] = new double[] { 0, 3, 4 };
      A[2] = new double[] { 0, 4, 9 };

      Console.WriteLine("\nSource matrix: ");
      MatShow(A, 1, 6);

      double[] eigVals;
      double[][] eigVecs;
      Eigen(A, out eigVals, out eigVecs);

      Console.WriteLine("\nEigenvalues: ");
      VecShow(eigVals, 4, 9);

      Console.WriteLine("\nEigenvectors (columns): ");
      MatShow(eigVecs, 4, 9);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main

    // ------------------------------------------------------

    public static void Eigen(double[][] M,
      out double[] eigenVals, out double[][] eigenVecs)
    {
      // compute eigenvalues, eigenvectors at same time
 
      int n = M.Length;
      double[][] X = MatCopy(M);  // mat must be square
      double[][] Q; double[][] R;
      double[][] pq = MatIdentity(n);
      int maxCt = 10000;

      int ct = 0;
      while (ct "lt" maxCt)
      {
        MatDecomposeQR(X, out Q, out R, false);
        pq = MatProduct(pq, Q);
        X = MatProduct(R, Q);  // note order
        ++ct;

        if (MatIsUpperTri(X, 1.0e-8) == true)
          break;
      }

      // eigenvalues are diag elements of X
      double[] evals = new double[n];
      for (int i = 0; i "lt" n; ++i)
        evals[i] = X[i][i];

      // eigenvectors are columns of pq
      double[][] evecs = MatCopy(pq);

      eigenVals = evals;
      eigenVecs = evecs;
    }

    // ------------------------------------------------------

    public static double[] Eigenvalues(double[][] A)
    {
      // just eigenvalues only
      int n = A.Length;
      double[][] X = MatCopy(A);  // mat must be square
      double[][] Q; double[][] R;
      int maxCt = 10000;
      
      int ct = 0;
      while (ct "lt" maxCt)
      {
        MatDecomposeQR(X, out Q, out R, false);
        X = MatProduct(R, Q);  // note order
        ++ct;

        if (MatIsUpperTri(X, 1.0e-8) == true)
          break;
      }

      if (ct "gte" maxCt)
        Console.WriteLine("Warn: Eigenvalues convergence ");

      double[] result = new double[n];
      for (int j = 0; j "lt" n; ++j)
        result[j] = X[j][j];
      return result;
    }

    // ------------------------------------------------------

    public static double[] MatEigenVecFromEigenVal(
      double[][] A, double eigenVal)
    {
      // I suspect this funcion has bugs -- do not use!!
      // A is the source matrix that produced eigenvalue
      double tol = 1.0e-6; // stopping criteria
      int maxIter = 1000;
      double noise = 1.0e-8; // avoid inverse failure
      int n = A.Length;
      double[] curr = new double[n];
      for (int i = 0; i "lt" n; ++i)
        curr[i] = 1.0;
      double[] prev;

      double[][] AminusEigenvec = MatCopy(A);
      for (int i = 0; i "lt" n; ++i)
        AminusEigenvec[i][i] -= eigenVal + noise;

      int iter = 0;
      while (iter "lt" maxIter)
      {
        prev = curr;
        curr = MatLinSolveQR(AminusEigenvec, prev);

        // normalize curr sign and length
        if (curr[0] "lt" 0.0)  // OK but not necessary
          curr = VecNegation(curr);
        
        double norm = VecNorm(curr);
        for (int j = 0; j "lt" n; ++j)
          curr[j] /= norm;

        ++iter;
        if (VecEqual(prev, curr, tol) == true)
          break;

        double[] negatedCurr = VecNegation(curr);
        if (VecEqual(prev, negatedCurr, tol) == true)
          break;
      }

      if (iter "gte" maxIter)
        Console.WriteLine("Warn: Eigenvec max iterations ");

      return curr;
    }

    // ------------------------------------------------------

    static double[] MatLinSolveQR(double[][] A, double[] b)
    {
      // solve general system of equations using QR
      // A * x = b
      // Q*R * x = b
      // inv(Q*R) * Q*R * x = inv(Q*R) * b
      // x = inv(R) * inv(Q) * b
      // x = inv(R) * tr(Q) * b where R is upper tri
      double[][] Q; double[][] R;
      MatDecomposeQR(A, out Q, out R, false);
      double[][] Ri = MatInverseUpperTri(R);
      double[][] Qt = MatTranspose(Q);
      double[][] B = VecToMat(b, b.Length, 1);
      double[][] RiQt = MatProduct(Ri, Qt);
      double[][] RiQtB = MatProduct(RiQt, B);
      return MatToVec(RiQtB);
    }

    // ------------------------------------------------------

    public static double[][]
      MatInverseUpperTri(double[][] upper)
    {
      int n = upper.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[j][k] -= result[j][i] * upper[i][k];
          }
          result[j][k] /= upper[k][k];
        }
      }
      return result;
    }

    // ------------------------------------------------------

    static double[][] MatTranspose(double[][] m)
    {
      int nr = m.Length;
      int nc = m[0].Length;
      double[][] result = MatCreate(nc, nr);  // note
      for (int i = 0; i "lt" nr; ++i)
        for (int j = 0; j "lt" nc; ++j)
          result[j][i] = m[i][j];
      return result;
    }

    // ------------------------------------------------------

    static double[] MatToVec(double[][] m)
    {
      int rows = m.Length;
      int cols = m[0].Length;
      double[] result = new double[rows * cols];
      int k = 0;
      for (int i = 0; i "lt" rows; ++i)
        for (int j = 0; j "lt" cols; ++j)
        {
          result[k++] = m[i][j];
        }
      return result;
    }

    // ------------------------------------------------------

    public static void MatDecomposeQR(double[][] mat,
      out double[][] q, out double[][] r,
      bool standardize)
    {
      // QR decomposition, Householder algorithm.
      // assumes square matrix

      int n = mat.Length;  // assumes mat is nxn
      int nCols = mat[0].Length;
      if (n != nCols) Console.WriteLine("M not square ");

      double[][] Q = MatIdentity(n);
      double[][] R = MatCopy(mat);
      for (int i = 0; i "lt" n - 1; ++i)
      {
        double[][] H = MatIdentity(n);
        double[] a = new double[n - i];
        int k = 0;
        for (int ii = i; ii "lt" n; ++ii)  // last part col [i]
          a[k++] = R[ii][i];

        double normA = VecNorm(a);
        if (a[0] "lt" 0.0) { normA = -normA; }
        double[] v = new double[a.Length];
        for (int j = 0; j "lt" v.Length; ++j)
          v[j] = a[j] / (a[0] + normA);
        v[0] = 1.0;

        double[][] h = MatIdentity(a.Length);
        double vvDot = VecDot(v, v);
        double[][] alpha = VecToMat(v, v.Length, 1);
        double[][] beta = VecToMat(v, 1, v.Length);
        double[][] aMultB = MatProduct(alpha, beta);

        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            h[ii][jj] -= (2.0 / vvDot) * aMultB[ii][jj];

        // copy h into lower right of H
        int d = n - h.Length;
        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            H[ii + d][jj + d] = h[ii][jj];

        Q = MatProduct(Q, H);
        R = MatProduct(H, R);
      } // i

      if (standardize == true)
      {
        // standardize so R diagonal is all positive
        double[][] D = MatCreate(n, n);
        for (int i = 0; i "lt" n; ++i)
        {
          if (R[i][i] "lt" 0.0) D[i][i] = -1.0;
          else D[i][i] = 1.0;
        }
        Q = MatProduct(Q, D);
        R = MatProduct(D, R);
      }

      q = Q;
      r = R;

    } // QR decomposition

    // ------------------------------------------------------
    
    public static double[][] MatCreate(int rows, int cols)
    {
      double[][] result = new double[rows][];
      for (int i = 0; i "lt" rows; ++i)
        result[i] = new double[cols];
      return result;
    }

    // ------------------------------------------------------
    
    public static double[][] MatIdentity(int n)
    {
      double[][] result = MatCreate(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      return result;
    }

    // ------------------------------------------------------
    
    public static double[][] MatCopy(double[][] m)
    {
      int nRows = m.Length; int nCols = m[0].Length;
      double[][] result = MatCreate(nRows, nCols);
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = m[i][j];
      return result;
    }

    // ------------------------------------------------------
    
    public static double[][] MatProduct(double[][] matA,
      double[][] matB)
    {
      int aRows = matA.Length;
      int aCols = matA[0].Length;
      int bRows = matB.Length;
      int bCols = matB[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

      double[][] result = MatCreate(aRows, bCols);

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += matA[i][k] * matB[k][j];

      return result;
    }

    // ------------------------------------------------------
    
    public static bool MatIsUpperTri(double[][] mat,
      double tol)
    {
      int n = mat.Length;
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lt" i; ++j)
        {  // check lower vals
          if (Math.Abs(mat[i][j]) "gt" tol)
          {
            return false;
          }
        }
      }
      return true;

    }

    // ------------------------------------------------------
    
    public static void MatShow(double[][] m,
      int dec, int wid)
    {
      for (int i = 0; i "lt" m.Length; ++i)
      {
        for (int j = 0; j "lt" m[0].Length; ++j)
        {
          double v = m[i][j];
          if (Math.Abs(v) "lt" 1.0e-8) v = 0.0;  // hack
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

    // ------------------------------------------------------
    
    public static double VecDot(double[] v1, double[] v2)
    {
      double result = 0.0;
      int n = v1.Length;
      for (int i = 0; i "lt" n; ++i)
        result += v1[i] * v2[i];
      return result;
    }

    // ------------------------------------------------------
    
    public static double VecNorm(double[] vec)
    {
      int n = vec.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += vec[i] * vec[i];
      return Math.Sqrt(sum);
    }

    // ------------------------------------------------------

    public static bool VecEqual(double[] v1, double[] v2,
      double tol)
    {
      for (int i = 0; i "lt" v1.Length; ++i)
        if (Math.Abs(v1[i] - v2[i]) "gt" tol)
          return false;
      return true;
    }

    // ------------------------------------------------------
    
    public static double[] VecNegation(double[] vec)
    {
      int n = vec.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = -vec[i];
      return result;
    }

    // ------------------------------------------------------

    public static double[][] VecToMat(double[] vec,
      int nRows, int nCols)
    {
      double[][] result = MatCreate(nRows, nCols);
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = vec[k++];
      return result;
    }

    // ------------------------------------------------------
    
    public static void VecShow(double[] vec,
      int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
      {
        double x = vec[i];
        if (Math.Abs(x) "lt" 1.0e-8) x = 0.0;
        Console.Write(x.ToString("F" + dec).
          PadLeft(wid));
      }
      Console.WriteLine("");
    }

    // ------------------------------------------------------

  } // Program
} // ns