core/properties/eigenvalues.js

/* eslint-disable no-param-reassign */
// reference: https://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

const Complex = require('@rayyamhk/complex');
const Matrix = require('../..');
const { INVALID_SQUARE_MATRIX } = require('../../Error');

/**
 * Calculates the eigenvalues of any square Matrix using QR Algorithm.<br><br>
 * 
 * The eigenvalues can be either real number or complex number.
 * Note that all eigenvalues are instance of Complex,
 * for more details please visit [Complex.js]{@link https://rayyamhk.github.io/Complex.js}.<br><br>
 * 
 * The eigenvalues are cached.
 * @memberof Matrix
 * @instance
 * @returns {Complex[]} Array of eigenvalues
 */

function eigenvalues() {
  if (!this.isSquare()) {
    throw new Error(INVALID_SQUARE_MATRIX);
  }

  if (this._eigenvalues !== undefined) {
    return this._eigenvalues;
  }

  const size = this.size()[0];
  const values = [];
  const digit = this._digit;
  const EPSILON = 1 / ((10 ** digit) * 2);

  const clone = Matrix.clone(this)._matrix;
  let isConvergent = true; // flag
  let skip = false;

  // Transform matrix to Hessenberg matrix
  HouseholderTransform(clone, digit);

  for (let i = size - 1; i > 0; i--) {
    let divergenceCount = 0;
    let prev; // used to determine convergence

    // if obtains complex eigenvalues pair in previous iteration, skip current round
    if (skip) {
      skip = false;
      continue;
    }

    const shift = clone[size - 1][size - 1];

    // eslint-disable-next-line no-constant-condition
    while (true) {
      if (!isConvergent) { // if the current eigenvalue is not real
        prev = size2Eigenvalues(
          clone[i - 1][i - 1],
          clone[i - 1][i],
          clone[i][i - 1],
          clone[i][i],
        ).metric;
      } else { // if the current eigenvalue is real
        prev = Math.abs(clone[i][i - 1]);
      }

      // apply single shift
      for (let j = 0; j < size; j++) {
        clone[j][j] -= shift;
      }

      // Apply QR Algorithm
      HessenbergQR(clone, digit);

      for (let j = 0; j < size; j++) {
        clone[j][j] += shift;
      }

      if (isConvergent && prev < Math.abs(clone[i][i - 1])) {
        divergenceCount++;
      }

      // if the current eigenvalue is real and the entry is almost ZERO => break;
      if (isConvergent && Math.abs(clone[i][i - 1]) < EPSILON) {
        values[i] = new Complex(clone[i][i]);
        break;
      }

      // if the current eigenvalues pair is complex, if the difference of the previous eiganvalues and the
      // eigenvalues of submatrix is almost ZERO => break
      const {
        metric,
        eigen1,
        eigen2,
      } = size2Eigenvalues(clone[i - 1][i - 1], clone[i - 1][i], clone[i][i - 1], clone[i][i]);
      if (!isConvergent && Math.abs(prev - metric) < EPSILON) {
        isConvergent = true; // re-initialize
        skip = true;
        const { re: re1, im: im1 } = eigen1;
        const { re: re2, im: im2 } = eigen2;
        values[i] = new Complex(re1, im1);
        values[i - 1] = new Complex(re2, im2);
        break;
      }

      // if the entry doesn't converge => complex eigenvalues pair
      if (divergenceCount > 3) {
        isConvergent = false;
      }
    }
  }

  if (!skip) {
    values[0] = new Complex(clone[0][0]);
  }

  this._eigenvalues = values;
  return values;
};

function HouseholderTransform(A, digit) {
  const size = A.length;
  const EPSILON = 1 / ((10 ** digit) * 2);

  for (let j = 0; j < size - 2; j++) {
    let xNorm = 0;
    const u = new Array(size - j - 1);
    for (let i = j + 1; i < size; i++) {
      const entry = A[i][j];
      xNorm += entry ** 2;
      u[i - j - 1] = entry;
    }
    xNorm = Math.sqrt(xNorm);

    if (Math.abs(xNorm) < EPSILON) {
      continue;
    }

    if (u[0] >= 0) {
      u[0] += xNorm;
    } else {
      u[0] -= xNorm;
    }

    // Make 'u' unit vector
    let uNorm = 0;
    for (let i = 0; i < u.length; i++) {
      uNorm += u[i] ** 2;
    }
    uNorm = Math.sqrt(uNorm);

    for (let i = 0; i < u.length; i++) {
      u[i] /= uNorm;
    }

    // update the matrix, multiply P from left
    for (let n = j; n < size; n++) { // column
      const v = new Array(size - j - 1);
      for (let m = j + 1; m < size; m++) {
        v[m - j - 1] = A[m][n];
      }

      let scaler = 0;
      for (let m = 0; m < v.length; m++) {
        scaler += v[m] * u[m];
      }
      scaler *= 2;

      for (let m = j + 1; m < size; m++) { // row
        if (n === j && m !== j + 1) {
          A[m][n] = 0;
        } else {
          A[m][n] = v[m - j - 1] - scaler * u[m - j - 1];
        }
      }
    }

    // update the matrix, multiply P from right
    for (let m = 0; m < size; m++) { // row
      const v = new Array(size - j - 1);
      for (let n = j + 1; n < size; n++) {
        v[n - j - 1] = A[m][n];
      }

      let scaler = 0;
      for (let n = 0; n < v.length; n++) {
        scaler += v[n] * u[n];
      }
      scaler *= 2;

      for (let n = j + 1; n < size; n++) { // column
        A[m][n] = v[n - j - 1] - scaler * u[n - j - 1];
      }
    }
  }
}

function HessenbergQR(H, digit) {
  const size = H.length;
  const EPSILON = 1 / ((10 ** digit) * 2);
  const sincos = new Array(size - 1);

  for (let i = 0; i < size - 1; i++) {
    const a = H[i][i];
    const c = H[i + 1][i];

    const norm = Math.sqrt(a ** 2 + c ** 2);
    if (norm < EPSILON) {
      continue;
    }
    const cos = a / norm;
    const sin = (c * -1) / norm;
    sincos[i] = [sin, cos];

    const row1 = new Array(size - i);
    const row2 = new Array(size - i);

    for (let j = i; j < size; j++) {
      row1[j - i] = H[i][j];
      row2[j - i] = H[i + 1][j];
    }

    for (let j = i; j < size; j++) {
      H[i][j] = cos * row1[j - i] + sin * -1 * row2[j - i];
      if (i === j) {
        H[i + 1][j] = 0;
      } else {
        H[i + 1][j] = sin * row1[j - i] + cos * row2[j - i];
      }
    }
  }

  for (let j = 0; j < size - 1; j++) {
    if (!sincos[j]) {
      continue;
    }
    const [sin, cos] = sincos[j];

    const col1 = new Array(j + 2);
    const col2 = new Array(j + 2);
    for (let i = 0; i <= j + 1; i++) {
      col1[i] = H[i][j];
      col2[i] = H[i][j + 1];
    }

    for (let i = 0; i <= j + 1; i++) {
      H[i][j] = col1[i] * cos - col2[i] * sin;
      H[i][j + 1] = col1[i] * sin + col2[i] * cos;
    }
  }
}

// find the eigenvalues of 2x2 matrix
function size2Eigenvalues(e11, e12, e21, e22) {
  const b = (e11 + e22) * -1;
  const c = (e11 * e22 - e21 * e12);
  const delta = b ** 2 - 4 * c;
  let re1;
  let im1;
  let re2;
  let im2;
  if (delta >= 0) {
    im1 = 0;
    im2 = 0;
    if (b >= 0) {
      re1 = (b * -1 - Math.sqrt(delta)) / 2;
    } else {
      re1 = (b * -1 + Math.sqrt(delta)) / 2;
    }
    re2 = c / re1;
  } else {
    re1 = -b / 2;
    re2 = re1;
    im1 = Math.sqrt(delta * -1) / 2;
    im2 = im1 * -1;
  }
  return {
    metric: Math.sqrt(re1 ** 2 + im1 ** 2),
    eigen1: { re: re1, im: im1 },
    eigen2: { re: re2, im: im2 },
  };
}

module.exports = eigenvalues;