Merge pull request #68 from jbrown1618/cholesky-decomposition

Implemented Cholesky Decomposition

Author: jbrown1618

docs/VECTOR.api.md

@@ -378,7 +378,7 @@ export class NumberOperations extends ScalarOperations<number> {
     getAdditiveInverse(x: number): number;
     getMultiplicativeIdentity(): number;
     getMultiplicativeInverse(x: number): number | undefined;
-    getPrincipalSquareRoot(x: number): number;
+    getPrincipalSquareRoot(x: number): number | undefined;
     multiply(first: number, second: number): number;
     norm(x: number): number;
     prettyPrint(x: number): string;
@@ -453,7 +453,7 @@ export abstract class ScalarOperations<S> {
     abstract getAdditiveInverse(x: S): S;
     abstract getMultiplicativeIdentity(): S;
     abstract getMultiplicativeInverse(x: S): S | undefined;
-    abstract getPrincipalSquareRoot(x: S): S;
+    abstract getPrincipalSquareRoot(x: S): S | undefined;
     abstract multiply(first: S, second: S): S;
     negativeOne(): S;
     abstract norm(x: S): number;

docs/api/vector.numberoperations.getprincipalsquareroot.md

@@ -9,7 +9,7 @@ Returns the principal square root of a scalar.
 <b>Signature:</b>
 
 ```typescript
-getPrincipalSquareRoot(x: number): number;
+getPrincipalSquareRoot(x: number): number | undefined;
 ```
 
 ## Parameters
@@ -20,7 +20,7 @@ getPrincipalSquareRoot(x: number): number;
 
 <b>Returns:</b>
 
-`number`
+`number | undefined`
 
 The square root
 

docs/api/vector.scalaroperations.getprincipalsquareroot.md

@@ -9,7 +9,7 @@ Returns the principal square root of a scalar.
 <b>Signature:</b>
 
 ```typescript
-abstract getPrincipalSquareRoot(x: S): S;
+abstract getPrincipalSquareRoot(x: S): S | undefined;
 ```
 
 ## Parameters
@@ -20,7 +20,7 @@ abstract getPrincipalSquareRoot(x: S): S;
 
 <b>Returns:</b>
 
-`S`
+`S | undefined`
 
 The square root
 

package.json

@@ -1,7 +1,7 @@
 {
   "name": "@josh-brown/vector",
   "author": "Joshua Brown <[email protected]>",
-  "version": "0.0.3",
+  "version": "0.1.0",
   "description": "A linear algebra library written in TypeScript that focuses on generality, extensibility, and ease of use.",
   "license": "ISC",
   "repository": {

src/algorithms/CholeskyDecomposition.spec.ts

@@ -0,0 +1,39 @@
+import { expect } from 'chai';
+import { describe, it } from 'mocha';
+import { NumberMatrix } from '../types/matrix/NumberMatrix';
+import { calculateCholeskyDecomposition } from './CholeskyDecomposition';
+
+describe('CholeskyDecomposition', () => {
+  const matrixBuilder = NumberMatrix.builder();
+
+  describe('calculateCholeskyDecomposition', () => {
+    it('calculates the Cholesky Decomposition of a matrix A', () => {
+      const A = matrixBuilder.fromArray([[4, 12, -16], [12, 37, -43], [-16, -43, 98]]);
+      const expectedL = matrixBuilder.fromArray([[2, 0, 0], [6, 1, 0], [-8, 5, 3]]);
+
+      const decomposition = calculateCholeskyDecomposition(A);
+      expect(decomposition).not.to.be.undefined;
+      const { L } = decomposition as any;
+      expect(L).to.deep.equal(expectedL);
+
+      // Check that LL* equals A
+      const LLstar = L.multiply(L.adjoint());
+      expect(LLstar.equals(A)).to.be.true;
+    });
+
+    it('returns undefined for a non-square matrix', () => {
+      const A = matrixBuilder.fromArray([[1, 2]]);
+      expect(calculateCholeskyDecomposition(A)).to.be.undefined;
+    });
+
+    it('returns undefined for a non-symmetric matrix', () => {
+      const A = matrixBuilder.fromArray([[4, 12, -16], [12, 37, -43], [-16, -46, 98]]);
+      expect(calculateCholeskyDecomposition(A)).to.be.undefined;
+    });
+
+    it('returns undefined for a non-positive-definite matrix', () => {
+      const A = matrixBuilder.fromArray([[-1]]);
+      expect(calculateCholeskyDecomposition(A)).to.be.undefined;
+    });
+  });
+});

src/algorithms/CholeskyDecomposition.ts

@@ -0,0 +1,55 @@
+import { Matrix } from '../types/matrix/Matrix';
+import { isHermitian } from '../utilities/MatrixProperties';
+
+export interface CholeskyDecomposition<S> {
+  L: Matrix<S>;
+}
+
+/**
+ * Uses the serial version of the Cholesky algorith to calculate the Cholesky
+ * decomposition of a matrix `A`.
+ *
+ * @remarks
+ * A Cholesky decomposition of a matrix `A` consists of a lower-triangular
+ * matrix `L` such that _LL* = A_.
+ *
+ * A Cholesky decomposition only exists if `A` is symmetrix and positive-definite.
+ * @param A - The matrix to decompose
+ */
+export function calculateCholeskyDecomposition<S>(
+  A: Matrix<S>
+): CholeskyDecomposition<S> | undefined {
+  if (!isHermitian(A)) return undefined;
+
+  const ops = A.ops();
+  const builder = A.builder();
+  const dim = A.getNumberOfColumns();
+
+  let L = builder.zeros(dim);
+
+  for (let j = 0; j < dim; j++) {
+    const Ajj = A.getEntry(j, j);
+    let entrySquared = Ajj;
+    for (let p = 0; p < j; p++) {
+      const Ljp = L.getEntry(j, p);
+      entrySquared = ops.subtract(entrySquared, ops.multiply(Ljp, ops.conjugate(Ljp)));
+    }
+    const Ljj = ops.getPrincipalSquareRoot(entrySquared);
+    if (Ljj === undefined) return undefined;
+    L = L.set(j, j, Ljj);
+
+    for (let i = j + 1; i < dim; i++) {
+      let difference = A.getEntry(i, j);
+      for (let p = 0; p < j; p++) {
+        const Lip = L.getEntry(i, p);
+        const Ljp = L.getEntry(j, p);
+        difference = ops.subtract(difference, ops.multiply(Lip, ops.conjugate(Ljp)));
+      }
+      const Lij = ops.divide(difference, Ljj);
+      if (Lij === undefined) return undefined;
+
+      L = L.set(i, j, Lij);
+    }
+  }
+  return { L };
+}

src/types/scalar/NumberOperations.spec.ts

@@ -100,7 +100,7 @@ describe('NumberOperations', () => {
     it('returns the positive square root', () => {
       testNumbers(n => {
         if (n < 0) {
-          expect(ops.getPrincipalSquareRoot(n)).to.be.NaN;
+          expect(ops.getPrincipalSquareRoot(n)).to.be.undefined;
         } else {
           expect(ops.getPrincipalSquareRoot(n)).to.equal(Math.sqrt(n));
         }

src/types/scalar/NumberOperations.ts

@@ -62,7 +62,8 @@ export class NumberOperations extends ScalarOperations<number> {
   /**
    * {@inheritdoc ScalarOperations.getPrincipalSquareRoot}
    */
-  public getPrincipalSquareRoot(x: number): number {
+  public getPrincipalSquareRoot(x: number): number | undefined {
+    if (x < 0) return undefined;
     return Math.sqrt(x);
   }
 

src/types/scalar/ScalarOperations.ts

@@ -170,7 +170,7 @@ export abstract class ScalarOperations<S> {
    * @returns The square root
    * @public
    */
-  public abstract getPrincipalSquareRoot(x: S): S;
+  public abstract getPrincipalSquareRoot(x: S): S | undefined;
 
   /**
    * Returns the norm (absolute value or magnitude) of a scalar