Implement matrix inversion

src/__tests__/basis.spec.ts

@@ -1,8 +1,10 @@
 import {describe, expect, it} from 'vitest';
 import {
   fractionalGram,
+  fractionalInv,
   fractionalLenstraLenstraLovasz,
   gram,
+  inv,
   lenstraLenstraLovasz,
 } from '../basis';
 import {
@@ -322,3 +324,71 @@ describe('Precise LLL basis reduction', () => {
     }
   });
 });
+
+describe('Matrix inverse', () => {
+  it('computes a 3x3 inverse', () => {
+    const mat = [
+      [2, -1], // Missing entry interpreted as 0
+      [-1, 2, -1],
+      [0, -1, 2],
+    ];
+    const inverse = inv(mat).map(row => row.map(x => Math.round(4 * x) / 4));
+    expect(inverse).toEqual([
+      [0.75, 0.5, 0.25],
+      [0.5, 1, 0.5],
+      [0.25, 0.5, 0.75],
+    ]);
+  });
+
+  it('throws for non-square matrix', () => {
+    expect(() =>
+      inv([
+        [1, 2],
+        [3, 4],
+        [5, 6],
+      ])
+    ).toThrow('Non-square matrix');
+  });
+
+  it('throws for singular matrix', () => {
+    expect(() =>
+      inv([
+        [1, 0],
+        [0, 0],
+      ])
+    ).toThrow('Matrix is singular');
+  });
+
+  it('computes a 3x3 with fractional entries', () => {
+    const mat = [
+      [2, -1], // Missing entry interpreted as 0
+      [-1, '2/1', -1],
+      [0, -1, '2'],
+    ];
+    const inverse = fractionalInv(mat).map(row => row.map(x => x.toFraction()));
+    expect(inverse).toEqual([
+      ['3/4', '1/2', '1/4'],
+      ['1/2', '1', '1/2'],
+      ['1/4', '1/2', '3/4'],
+    ]);
+  });
+
+  it('throws for non-square matrix with fractional entries', () => {
+    expect(() =>
+      fractionalInv([
+        [1, 2],
+        [3, 4],
+        [5, 6],
+      ])
+    ).toThrow('Non-square matrix');
+  });
+
+  it('throws for singular matrix with fractional entries', () => {
+    expect(() =>
+      fractionalInv([
+        [1, 0],
+        [0, 0],
+      ])
+    ).toThrow('Matrix is singular');
+  });
+});

src/basis.ts

@@ -218,3 +218,172 @@ export function fractionalLenstraLenstraLovasz(
     },
   };
 }
+
+// XXX: I'm sorry. This matrix inversion algorithm is not particularly sophisticated. Existing solutions just come with too much bloat.
+
+/**
+ * Compute the (multiplicative) inverse of a matrix.
+ * @param matrix Matrix to be inverted.
+ * @returns The multiplicative inverse.
+ * @throws An error if the matrix is not square or not invertible.
+ */
+export function inv(matrix: number[][]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  const result: number[][] = [];
+  for (let i = 0; i < height; ++i) {
+    result.push(Array(width).fill(0));
+    result[i][i] = 1;
+  }
+  // Don't modify input
+  matrix = matrix.map(row => [...row]);
+  // Coerce missing entries to zeros
+  for (let y = 0; y < height; ++y) {
+    for (let x = matrix[y].length; x < width; ++x) {
+      matrix[y][x] = 0;
+    }
+  }
+  // Row echelon form
+  for (let x = 0; x < width; ++x) {
+    for (let y = x; y < height; ++y) {
+      if (matrix[y][x]) {
+        if (x === y) {
+          break;
+        }
+        let temp = matrix[y];
+        matrix[y] = matrix[x];
+        matrix[x] = temp;
+
+        temp = result[y];
+        result[y] = result[x];
+        result[x] = temp;
+        break;
+      }
+    }
+  }
+  // Put ones along the diagonal, zeros in the lower triangle
+  for (let x = 0; x < width; ++x) {
+    let s = matrix[x][x];
+    if (!s) {
+      throw new Error('Matrix is singular');
+    }
+    if (s !== 1) {
+      s = 1 / s;
+      matrix[x] = matrix[x].map(a => a * s);
+      result[x] = result[x].map(a => a * s);
+    }
+    for (let y = x + 1; y < height; ++y) {
+      s = matrix[y][x];
+      if (s) {
+        matrix[y] = matrix[y].map((a, i) => a - s * matrix[x][i]);
+        result[y] = result[y].map((a, i) => a - s * result[x][i]);
+      }
+    }
+  }
+  // Eliminate remaining entries in the upper triangle
+  for (let x = width - 1; x > 0; --x) {
+    for (let y = x - 1; y >= 0; --y) {
+      const s = matrix[y][x];
+      if (s) {
+        // No need to keep track of these entries anymore.
+        // matrix[y] = matrix[y].map((a, i) => a - s * matrix[x][i]);
+        result[y] = result[y].map((a, i) => a - s * result[x][i]);
+      }
+    }
+  }
+  return result;
+}
+
+/**
+ * Compute the (multiplicative) inverse of a matrix.
+ * @param matrix Matrix to be inverted.
+ * @returns The multiplicative inverse.
+ * @throws An error if the matrix is not square or not invertible.
+ */
+export function fractionalInv(matrix: ProtoFractionalMonzo[]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  const result: FractionalMonzo[] = [];
+  for (let i = 0; i < height; ++i) {
+    const row: FractionalMonzo = [];
+    for (let j = 0; j < width; ++j) {
+      if (i === j) {
+        row.push(new Fraction(1));
+      } else {
+        row.push(new Fraction(0));
+      }
+    }
+    result.push(row);
+  }
+  // Don't modify input
+  const matrix_: FractionalMonzo[] = matrix.map(row =>
+    row.map(f => new Fraction(f))
+  );
+  // Coerce missing entries to zeros
+  for (let y = 0; y < height; ++y) {
+    for (let x = matrix_[y].length; x < width; ++x) {
+      matrix_[y][x] = new Fraction(0);
+    }
+  }
+  // Row echelon form
+  for (let x = 0; x < width; ++x) {
+    for (let y = x; y < height; ++y) {
+      if (matrix_[y][x].n) {
+        if (x === y) {
+          break;
+        }
+        let temp = matrix_[y];
+        matrix_[y] = matrix_[x];
+        matrix_[x] = temp;
+
+        temp = result[y];
+        result[y] = result[x];
+        result[x] = temp;
+        break;
+      }
+    }
+  }
+  // Put ones along the diagonal, zeros in the lower triangle
+  for (let x = 0; x < width; ++x) {
+    let s = matrix_[x][x];
+    if (!s.n) {
+      throw new Error('Matrix is singular');
+    }
+    if (!s.isUnity()) {
+      s = s.inverse();
+      matrix_[x] = matrix_[x].map(a => a.mul(s));
+      result[x] = result[x].map(a => a.mul(s));
+    }
+    for (let y = x + 1; y < height; ++y) {
+      s = matrix_[y][x];
+      if (s.n) {
+        matrix_[y] = matrix_[y].map((a, i) => a.sub(s.mul(matrix_[x][i])));
+        result[y] = result[y].map((a, i) => a.sub(s.mul(result[x][i])));
+      }
+    }
+  }
+  // Eliminate remaining entries in the upper triangle
+  for (let x = width - 1; x > 0; --x) {
+    for (let y = x - 1; y >= 0; --y) {
+      const s = matrix_[y][x];
+      if (s.n) {
+        // No need to keep track of these entries anymore.
+        // matrix_[y] = matrix_[y].map(...);
+        result[y] = result[y].map((a, i) => a.sub(s.mul(result[x][i])));
+      }
+    }
+  }
+  return result;
+}