core/decompositions/LU.js

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

/**
 * Calculates the LUP decomposition of the Matrix,
 * where L is lower triangular matrix which diagonal entries are always 1,
 * U is upper triangular matrix, and P is permutation matrix.<br><br>
 * 
 * It is implemented using Gaussian Elimination with Partial Pivoting in order to
 * reduce the error caused by floating-point arithmetic.<br><br>
 * 
 * Note that if optimized is true, P is a Permutation Array and both L and U are merged
 * into one matrix in order to improve performance.
 * @memberof Matrix
 * @static
 * @param {Matrix} A - Any matrix
 * @param {boolean} [optimized=false] - Returns [P, LU] if it is true, [P, L, U] if it is false
 * @returns {Matrix[]} The LUP decomposition of Matrix
 */
function LU(A, optimized = false) {
  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 permutation = initPermutation(row);
  const copy = this.clone(A)._matrix;

  for (let i = 0; i < row - 1; i++) {
    const currentCol = Math.min(i, col);

    // apply Partial Pivoting
    PartialPivoting(copy, permutation, currentCol, row, col);

    const ith = permutation[i];
    const pivot = copy[ith][currentCol];

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

    for (let j = i + 1; j < row; j++) {
      const jth = permutation[j];
      const entry = copy[jth][currentCol];

      if (Math.abs(entry) >= EPSILON) {
        const factor = entry / pivot;
        for (let k = currentCol; k < col; k++) {
          copy[jth][k] -= factor * copy[ith][k];
        }
        copy[jth][currentCol] = factor;
      }
    }
  }

  const result = new Array(row);
  for (let i = 0; i < row; i++) {
    result[i] = copy[permutation[i]];
  }

  if (optimized) {
    return [permutation, new this(result)];
  }

  const P = this.generate(row, row, (i, j) => {
    const idx = permutation[i];
    if (j === idx) {
      return 1;
    }
    return 0;
  });
  const L = this.generate(row, size, (i, j) => {
    if (i === j) {
      return 1;
    }
    if (i < j) {
      return 0;
    }
    return result[i][j];
  });
  const U = this.generate(size, col, (i, j) => {
    if (i > j) {
      return 0;
    }
    return result[i][j];
  });

  return [P, L, U];
};

function initPermutation(size) {
  const permutation = new Array(size);
  for (let i = 0; i < size; i++) {
    permutation[i] = i;
  }
  return permutation;
}

function PartialPivoting(matrix, permutation, pos, row, col) {
  const currentCol = Math.min(pos, col);
  let maxIdx = pos;
  let max = Math.abs(matrix[permutation[pos]][currentCol]);
  for (let i = pos + 1; i < row; i++) {
    const value = Math.abs(matrix[permutation[i]][currentCol]);
    if (value > max) {
      maxIdx = i;
      max = value;
    }
  }
  const t = permutation[pos];
  permutation[pos] = permutation[maxIdx];
  permutation[maxIdx] = t;
}

module.exports = LU;