From c9cbdcd503d9a70e0988d7fdc62f71b700cf4525 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Mon, 24 Jun 2024 18:43:47 +0300 Subject: [PATCH] Improve inv() and det() numerical stability ref #42 --- CHANGELOG.md | 2 + src/__tests__/basis.spec.ts | 53 ++++++++++++++++++++++++++- src/basis.ts | 73 +++++++++++++++++++------------------ 3 files changed, 92 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e5bdfb1..68506dd 100644 --- 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 index cc93e94..aeec2d7 100644 --- 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 index db5ce09..fd045a5 100644 --- 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) {