choleskyDecomposition method

@override List<RealMatrix> choleskyDecomposition()

Uses the the Cholesky decomposition algorithm to factor the matrix into the product of a lower triangular matrix and its conjugate transpose. In particular, this method returns the L and L^T matrices of the

• A = L x L^T

relation. The algorithm might fail in case the square root of a negative number were encountered.

The returned list contains L at index 0 and L^T at index 1.

A MatrixException object is thrown if the matrix isn't square.

Implementation

@override
List<RealMatrix> choleskyDecomposition() {
  // Making sure that the matrix is squared
  if (!isSquareMatrix) {
    throw const MatrixException(
      'LU decomposition only works with square matrices!',
    );
  }

  // Exit immediately because if [0,0] is a negative number, the algorithm
  // cannot even start since the square root of a negative number in R is not
  // allowed.
  if (this(0, 0) <= 0) {
    throw const SystemSolverException('The matrix is not positive-definite.');
  }

  // Creating L and Lt matrices
  final L = List.generate(
    rowCount * columnCount,
    (_) => 0.0,
    growable: false,
  );
  final transpL = List.generate(
    rowCount * columnCount,
    (_) => 0.0,
    growable: false,
  );

  // Computing the L matrix so that A = L * Lt (where 'Lt' is L transposed)
  for (var i = 0; i < rowCount; i++) {
    for (var j = 0; j <= i; j++) {
      var sum = 0.0;
      if (j == i) {
        for (var k = 0; k < j; k++) {
          sum += _getDataAt(L, j, k) * _getDataAt(L, j, k);
        }
        _setDataAt(L, j, j, sqrt(this(i, j) - sum));
      } else {
        for (var k = 0; k < j; k++) {
          sum += _getDataAt(L, i, k) * _getDataAt(L, j, k);
        }
        _setDataAt(L, i, j, (this(i, j) - sum) / _getDataAt(L, j, j));
      }
    }
  }

  // Computing Lt, the transposed version of L
  for (var i = 0; i < rowCount; i++) {
    for (var j = 0; j < rowCount; j++) {
      _setDataAt(transpL, i, j, _getDataAt(L, j, i));
    }
  }

  return [
    RealMatrix.fromFlattenedData(
      rows: rowCount,
      columns: columnCount,
      data: L,
    ),
    RealMatrix.fromFlattenedData(
      rows: rowCount,
      columns: columnCount,
      data: transpL,
    ),
  ];
}