core/decompositions/QR.js

const { INVALID_MATRIX } = require('../../Error');

/**
 * Calculates the QR decomposition of the Matrix
 * where Q is orthogonal matrix, R is upper triangular matrix.<br><br>
 * 
 * The algorithm is implemented using Householder Transform instead of Gram–Schmidt process
 * because the Householder Transform is more numerically stable.
 * @memberof Matrix
 * @static
 * @param {Matrix} A - Any matrix
 * @returns {Matrix[]} The QR decomposition of matrix in the form of [Q, R]
 */
function QR(A) {
  if (!(A instanceof this)) {
    throw new Error(INVALID_MATRIX);
  }

  const [row, col] = A.size();
  const size = Math.min(row, col);
  const EPSILON = 1 / ((10 ** A._digit) * 2);

  const matrixR = this.clone(A)._matrix;
  const matrixQ = this.identity(row)._matrix;

  for (let j = 0; j < size; j++) {
    // if all entries below main diagonal are considered as zero, skip this round
    let skip = true;
    for (let i = j + 1; i < row; i++) {
      if (Math.abs(matrixR[i][j]) >= EPSILON) {
        skip = false;
        break;
      }
    }

    if (!skip) {
      // Apply Householder transform
      let norm = 0;
      for (let i = j; i < row; i++) {
        norm += matrixR[i][j] ** 2;
      }
      norm = Math.sqrt(norm);

      // reduce floating point arithmatic error
      let s = -1;
      if (matrixR[j][j] < 0) {
        s = 1;
      }
      const u1 = matrixR[j][j] - s * norm;

      let w = new Array(row - j);
      for (let i = 0; i < row - j; i++) {
        w[i] = matrixR[i + j][j] / u1;
      }
      w[0] = 1;

      const tau = (-1 * s * u1) / norm;

      const subR = new Array(row - j);
      for (let i = 0; i < row - j; i++) {
        const newRow = new Array(col);
        for (let k = 0; k < col; k++) {
          newRow[k] = matrixR[j + i][k];
        }
        subR[i] = newRow;
      }

      for (let i = j; i < row; i++) {
        for (let k = 0; k < col; k++) {
          let summation = 0;
          for (let m = 0; m < row - j; m++) {
            summation += subR[m][k] * w[m];
          }
          matrixR[i][k] = subR[i - j][k] - tau * w[i - j] * summation;
        }
      }

      const subQ = new Array(row);
      for (let i = 0; i < row; i++) {
        const newRow = new Array(row - j);
        for (let k = 0; k < row - j; k++) {
          newRow[k] = matrixQ[i][j + k];
        }
        subQ[i] = newRow;
      }

      for (let i = 0; i < row; i++) {
        for (let k = j; k < row; k++) {
          let summation = 0;
          for (let m = 0; m < row - j; m++) {
            summation += subQ[i][m] * w[m];
          }
          matrixQ[i][k] = subQ[i][k - j] - tau * w[k - j] * summation;
        }
      }
    }
  }

  for (let i = 0; i < row; i++) {
    for (let j = 0; j < col; j++) {
      if (i > j) {
        matrixR[i][j] = 0;
      }
    }
  }

  return [new this(matrixQ), new this(matrixR)];
};

module.exports = QR;