Fix: Allow QR decomposition (and therefore eigenvalue finding) of singular matrices. Fixes #129 (#132)

Author: jbrown1618

src/decompositions/QRDecomposition.ts

@@ -50,10 +50,7 @@ export function calculateQRDecomposition<S>(A: Matrix<S>): QRDecomposition<S> {
   // If any columns are the zero vector, then A was not full-rank to begin with.
   const qColumns: Vector<S>[] = [];
   for (let i = 0; i < dim; i++) {
-    const qi = normalize(uColumns[i]);
-    if (qi === undefined) {
-      throw Error('A is singular; no QR decomposition exists');
-    }
+    const qi = normalize(uColumns[i]) ?? uColumns[i]; // Default to the zero-vector
     qColumns.push(qi);
   }
   const Q = matrixBuilder.fromColumnVectors(qColumns);

src/decompositions/__tests__/QRDecomposition.spec.ts

@@ -22,25 +22,29 @@ describe('QRDecomposition', () => {
         [0, 0, 35],
       ]);
 
-      const result = calculateQRDecomposition(A);
-      expect(result.Q.equals(expectedQ)).toBe(true);
-      expect(result.R.equals(expectedR)).toBe(true);
+      const { Q, R } = calculateQRDecomposition(A);
+      expect(Q.equals(expectedQ)).toBe(true);
+      expect(R.equals(expectedR)).toBe(true);
 
       // Check that QR equals A
-      expect(result.Q.multiply(result.R).equals(A)).toBe(true);
+      expect(Q.multiply(R).equals(A)).toBe(true);
     });
 
     test('rejects a non-square matrix', () => {
       const A = mat([[1, 2]]);
       expect(() => calculateQRDecomposition(A)).toThrow();
     });
 
-    test('rejects a matrix with linearly dependent columns', () => {
+    test('handles a singular matrix', () => {
       const A = mat([
         [1, 1],
         [1, 1],
       ]);
-      expect(() => calculateQRDecomposition(A)).toThrow();
+
+      const { Q, R } = calculateQRDecomposition(A);
+
+      // The decomposition is not unique in this case, so just verify that it reconstitutes A
+      expect(Q.multiply(R).equals(A)).toBe(true);
     });
   });
 });

src/eigenvalues/__tests__/Eigenvalues.spec.ts

@@ -38,6 +38,22 @@ describe('Eigenvalues', () => {
         expect(pair.eigenvector.equals(expectedVectors[i])).toBe(true);
       });
     });
+
+    test('calculates eigenvalue-eigenvector pairs of a singular matrix', () => {
+      const A = mat([
+        [2, -1, -1],
+        [-1, 2, -1],
+        [-1, -1, 2],
+      ]);
+      const expectedValues = [3, 3, 0];
+      const expectedVectors = [vec([-2, 1, 1]), vec([-2, 1, 1]), vec([1, 1, 1])];
+
+      const pairs = eig(A);
+      pairs.forEach((pair, i) => {
+        expect(pair.eigenvalue).toBeCloseTo(expectedValues[i], 5);
+        expect(pair.eigenvector.equals(expectedVectors[i])).toBe(true);
+      });
+    });
   });
 
   describe('calculateEigenvalues', () => {

src/operations/RowOperations.ts

@@ -12,6 +12,8 @@ export interface RowOperationResult<S> {
   operator: Matrix<S>;
 }
 
+type VectorComparator<S> = (v1: Vector<S>, v2: Vector<S>) => number;
+
 /**
  * A wrapper for static methods representing the elementary row operations
  * @public
@@ -118,7 +120,7 @@ export class RowOperations {
     const ops = matrix.ops();
     // We will sort the rows of an identity matrix according to the number
     // of leading zeros in the corresponding row of `matrix`
-    const comparator: (r1: Vector<S>, r2: Vector<S>) => number = (iRow1, iRow2) => {
+    const comparator: VectorComparator<S> = (iRow1, iRow2) => {
       const rowIndex1 = getNumberOfLeadingZeros(iRow1);
       const rowIndex2 = getNumberOfLeadingZeros(iRow2);
       const mRow1 = matrix.getRow(rowIndex1);
@@ -127,6 +129,11 @@ export class RowOperations {
       const leadingZeros2 = getNumberOfLeadingZeros(mRow2);
 
       if (leadingZeros1 === leadingZeros2) {
+        const rowsAreAllZero = leadingZeros1 === mRow1.getDimension();
+        if (rowsAreAllZero) {
+          // The rows are equivalent - the sort order does not matter
+          return 0;
+        }
         // If they have the same number of leading zeros, put
         // the entry with the greatest magnitude on top
         const firstEntry1 = mRow1.getEntry(leadingZeros1);

src/operations/__tests__/RowOperations.spec.ts

@@ -1,4 +1,4 @@
-import { mat } from '../../utilities/aliases';
+import { mat, zeros } from '../../utilities/aliases';
 import { RowOperations } from '../RowOperations';
 
 describe('RowOperations', () => {
@@ -74,5 +74,33 @@ describe('RowOperations', () => {
       expect(result.result).toStrictEqual(sorted);
       expect(result.operator).toStrictEqual(permutation);
     });
+
+    test('sorts by the first entry when the number of leading zeros is the same', () => {
+      const unsorted = mat([
+        [0, 5, 5],
+        [5, 5, 5],
+        [0, 4, 5],
+      ]);
+      const sorted = mat([
+        [5, 5, 5],
+        [0, 5, 5],
+        [0, 4, 5],
+      ]);
+      const permutation = mat([
+        [0, 1, 0],
+        [1, 0, 0],
+        [0, 0, 1],
+      ]);
+
+      const result = RowOperations.pivot(unsorted);
+      expect(result.result).toStrictEqual(sorted);
+      expect(result.operator).toStrictEqual(permutation);
+    });
+
+    test('handles all-zero rows', () => {
+      const unsorted = zeros([3, 3]);
+      const { result } = RowOperations.pivot(unsorted);
+      expect(result).toStrictEqual(unsorted);
+    });
   });
 });