Improve inv() and det() numerical stability

ref #42
xenharmonic-devs · Jun 24, 2024 · c9cbdcd · c9cbdcd
1 parent a0b30ca
commit c9cbdcd
Show file tree

Hide file tree

Showing 3 changed files with 92 additions and 36 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,7 @@
 # Change log
 
+## 0.10.2
+ * Tweak: Make matrix inversion and determinant more stable numerically. [#42](https://github.com/xenharmonic-devs/xen-dev-utils/issues/42)
 ## 0.10.1
  * Feature: Accept vectors in `matmul` and `fractionalMatmul` and adjust return types accordingly.
  * Feature: `matscale`, `matadd`, `matsub` and fractional variants. [#41](https://github.com/xenharmonic-devs/xen-dev-utils/issues/41)

diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts
@@ -33,7 +33,8 @@ import {
   unapplyWeights,
 } from '../monzo';
 import {dot} from '../number-array';
-import {LOG_PRIMES} from '../primes';
+import {LOG_PRIMES, PRIMES} from '../primes';
+import {Fraction} from '../fraction';
 
 const FUZZ = 'FUZZ' in process.env;
 
@@ -479,6 +480,22 @@ describe('Matrix inverse', () => {
     ]);
   });
 
+  it.each([2, 3, 4, 5, 6, 7, 8, 9, 10, 11])(
+    'computes inverse of a Vandermonde matrix %s',
+    (N: number) => {
+      const mat: number[][] = [];
+      for (const p of PRIMES.slice(0, N)) {
+        const row = [...Array(N).keys()].map(i => p ** -i);
+        mat.push(row);
+      }
+      expect(
+        matmul(mat, inv(mat)).map(row =>
+          row.map(x => Math.round(10000 * x) / 10000 || 0)
+        )
+      ).toEqual(eye(N));
+    }
+  );
+
   it('throws for non-square matrix', () => {
     expect(() =>
       inv([
@@ -526,6 +543,22 @@ describe('Matrix inverse', () => {
     ]);
   });
 
+  it.each([2, 3, 4, 5, 6])(
+    'computes exact inverse of a Vandermonde matrix %s',
+    (N: number) => {
+      const mat: Fraction[][] = [];
+      for (const p of PRIMES.slice(0, N)) {
+        const row = [...Array(N).keys()].map(i => new Fraction(p).pow(-i)!);
+        mat.push(row);
+      }
+      expect(
+        fractionalMatmul(mat, fractionalInv(mat)).map(row =>
+          row.map(x => x.valueOf())
+        )
+      ).toEqual(eye(N));
+    }
+  );
+
   it('throws for non-square matrix with fractional entries', () => {
     expect(() =>
       fractionalInv([
@@ -605,6 +638,24 @@ describe('Determinant', () => {
     expect(determinant).toBe(54);
   });
 
+  it.each([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])(
+    'computes the determinant of a Vandermonde matrix %s',
+    (N: number) => {
+      const mat: number[][] = [];
+      for (const p of PRIMES.slice(0, N)) {
+        const row = [...Array(N).keys()].map(i => p ** -i);
+        mat.push(row);
+      }
+      let analytic = 1;
+      for (let i = 0; i < N; ++i) {
+        for (let j = i + 1; j < N; ++j) {
+          analytic *= 1 / PRIMES[j] - 1 / PRIMES[i];
+        }
+      }
+      expect(det(mat) / analytic).toBeCloseTo(1, 1);
+    }
+  );
+
   it('computes 0 for the origin', () => {
     const mat = [[0, 0], []];
     const determinant = det(mat);

diff --git a/src/basis.ts b/src/basis.ts
@@ -296,26 +296,27 @@ export function inv(matrix: number[][]) {
   }
   // Put ones along the diagonal, zeros in the lower triangle
   for (let x = 0; x < width; ++x) {
-    let s = matrix[x][x];
-    if (!s) {
-      // Row echelon form
-      for (let y = x + 1; y < height; ++y) {
-        if (matrix[y][x]) {
-          let temp = matrix[y];
-          matrix[y] = matrix[x];
-          matrix[x] = temp;
-
-          temp = result[y];
-          result[y] = result[x];
-          result[x] = temp;
-          break;
-        }
-      }
-      s = matrix[x][x];
-      if (!s) {
-        throw new Error('Matrix is singular');
+    // Maintain row echelon form by pivoting on the most dominant row.
+    let pivot: number | undefined;
+    let s = 0;
+    for (let y = x; y < height; ++y) {
+      if (Math.abs(matrix[y][x]) > Math.abs(s)) {
+        pivot = y;
+        s = matrix[y][x];
       }
     }
+    if (pivot === undefined) {
+      throw new Error('Matrix is singular');
+    }
+    if (x !== pivot) {
+      let temp = matrix[pivot];
+      matrix[pivot] = matrix[x];
+      matrix[x] = temp;
+
+      temp = result[pivot];
+      result[pivot] = result[x];
+      result[x] = temp;
+    }
     if (s !== 1) {
       s = 1 / s;
       matrix[x] = matrix[x].map(a => a * s);
@@ -389,7 +390,8 @@ export function fractionalInv(matrix: ProtoFractionalMonzo[]) {
   for (let x = 0; x < width; ++x) {
     let s = matrix_[x][x];
     if (!s.n) {
-      // Row echelon form
+      // Row echelon form (pivoting makes no difference over rationals)
+      // TODO: Figure out if there's a strategy to avoid blowing safe limits during manipulation.
       for (let y = x + 1; y < height; ++y) {
         if (matrix_[y][x].n) {
           let temp = matrix_[y];
@@ -598,24 +600,25 @@ export function det(matrix: number[][]) {
   matrix = matrix.map(row => [...row]);
   let result = 1;
   for (let x = 0; x < width; ++x) {
-    let d = matrix[x][x];
-    if (!d) {
-      // Row echelon form
-      for (let y = x + 1; y < height; ++y) {
-        if (matrix[y][x]) {
-          const temp = matrix[y];
-          matrix[y] = matrix[x];
-          matrix[x] = temp;
-
-          result = -result;
-          break;
-        }
-      }
-      d = matrix[x][x];
-      if (!d) {
-        return 0;
+    // Maintain row echelon form by pivoting on the most dominant row.
+    let pivot: number | undefined;
+    let d = 0;
+    for (let y = x; y < height; ++y) {
+      if (Math.abs(matrix[y][x]) > Math.abs(d)) {
+        pivot = y;
+        d = matrix[y][x];
       }
     }
+    if (pivot === undefined) {
+      return 0;
+    }
+    if (x !== pivot) {
+      const temp = matrix[pivot];
+      matrix[pivot] = matrix[x];
+      matrix[x] = temp;
+
+      result = -result;
+    }
     result *= d;
     d = 1 / d;
     for (let y = x + 1; y < height; ++y) {