From 327aa782a913e9953a2827b6e1c38df354f9e26c Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Sat, 22 Jun 2024 14:56:08 +0300 Subject: [PATCH] Implement matrix inversion Implement matmul and eye to verify implementation. Implement det. Implement transpose. --- src/__tests__/basis.spec.ts | 236 ++++++++++++++++++++ src/basis.ts | 423 ++++++++++++++++++++++++++++++++++++ 2 files changed, 659 insertions(+) diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts index f99892f..195d62f 100644 --- a/src/__tests__/basis.spec.ts +++ b/src/__tests__/basis.spec.ts @@ -1,11 +1,22 @@ import {describe, expect, it} from 'vitest'; import { + det, + eye, + fractionalDet, + fractionalEye, fractionalGram, + fractionalInv, fractionalLenstraLenstraLovasz, + fractionalMatmul, + fractionalTranspose, gram, + inv, lenstraLenstraLovasz, + matmul, + transpose, } from '../basis'; import { + FractionalMonzo, applyWeights, fractionalDot, fractionalMonzosEqual, @@ -322,3 +333,228 @@ describe('Precise LLL basis reduction', () => { } }); }); + +describe('Matrix multiplication', () => { + it('multiplies two matrices', () => { + const A = [ + [1, 0, 1], + [2, 1, 1], + [0, 1, 1], + [1, 1, 2], + ]; + const B = [ + [1, 2, 1], + [2, 3, 1], + [4, 2, 2], + ]; + const C = matmul(A, B); + expect(C).toEqual([ + [5, 4, 3], + [8, 9, 5], + [6, 5, 3], + [11, 9, 6], + ]); + }); + + it('multiplies two matrices (fractions)', () => { + const A = [ + [1, 0, 1], + [2, 1, 1], + [0, 0.5, 1], + [1, 1, 2], + ]; + const B = [ + [1, 2, 1], + [2, 3, 1], + ['4/3', 2, 2], + ]; + const C = fractionalMatmul(A, B); + expect(C.map(row => row.map(f => f.toFraction()))).toEqual([ + ['7/3', '4', '3'], + ['16/3', '9', '5'], + ['7/3', '7/2', '5/2'], + ['17/3', '9', '6'], + ]); + }); +}); + +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('computes another 3x3 inverse', () => { + const mat = [ + [-2, -1, 2], + [2, 1, 4], + [-3, 3, -1], + ]; + const inverse = inv(mat).map(row => row.map(x => Math.round(100000 * x))); + expect(inverse).toEqual([ + [-24074, 9259, -11111], + [-18519, 14815, 22222], + [16667, 16667, 0], + ]); + }); + + 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('computes another 3x3 inverse with fractional result', () => { + const mat = [ + [-2, -1, 2], + [2, 1, 4], + [-3, 3, -1], + ]; + const inverse = fractionalInv(mat).map(row => row.map(x => x.toFraction())); + expect(inverse).toEqual([ + ['-13/54', '5/54', '-1/9'], + ['-5/27', '4/27', '2/9'], + ['1/6', '1/6', '0'], + ]); + }); + + 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'); + }); + + it.runIf(FUZZ)('fuzzes for random inverses', () => { + for (let k = 0; k < 10000; ++k) { + const mat: number[][] = []; + const N = Math.ceil(1 + Math.random() * 10); + for (let i = 0; i < N; ++i) { + const row: number[] = []; + for (let j = 0; j < N; ++j) { + row.push(Math.random() * 10 - 5); + } + mat.push(row); + } + const inverse = inv(mat); + const I = matmul(mat, inverse).map(row => + row.map(x => Math.round(x * 1024) / 1024 || 0) + ); + expect(I).toEqual(eye(N)); + } + }); + + it.runIf(FUZZ)('fuzzes for random inverses (fractional)', () => { + for (let k = 0; k < 1000; ++k) { + const mat: number[][] = []; + const N = Math.ceil(1 + Math.random() * 7); + for (let i = 0; i < N; ++i) { + const row: number[] = []; + for (let j = 0; j < N; ++j) { + row.push(Math.round(Math.random() * 10 - 5)); + } + mat.push(row); + } + let inverse: FractionalMonzo[] | undefined; + try { + inverse = fractionalInv(mat); + } catch { + /** empty */ + } + if (inverse) { + const I = fractionalMatmul(mat, inverse); + expect(I).toEqual(fractionalEye(N)); + } + } + }); +}); + +describe('Determinant', () => { + it('computes the determinant of a 3x3 matrix', () => { + const mat = [ + [-2, -1, 2], + [2, 1, 4], + [-3, 3, -1], + ]; + const d = det(mat); + expect(d).toBe(54); + }); + + it('computes the determinant of a 3x3 matrix with fractional entries', () => { + const mat = [ + [-2, -1, 2], + [2, 0.5, 4], + [-3, 3, '-1/3'], + ]; + const d = fractionalDet(mat); + expect(d.toFraction()).toBe('152/3'); + }); +}); + +describe('Transpose', () => { + it('transposes a 3x2 matrix', () => { + const mat = [[1, 2], [3], [4, 5]]; + expect(transpose(mat)).toEqual([ + [1, 3, 4], + [2, 0, 5], + ]); + }); + + it('transposes a 3x2 matrix with rational entries', () => { + const mat = [[1, 0.5], [3], ['2/7', 5]]; + expect( + fractionalTranspose(mat).map(row => row.map(f => f.toFraction())) + ).toEqual([ + ['1', '3', '2/7'], + ['1/2', '0', '5'], + ]); + }); +}); diff --git a/src/basis.ts b/src/basis.ts index a8ed0f4..88b16e3 100644 --- a/src/basis.ts +++ b/src/basis.ts @@ -218,3 +218,426 @@ export function fractionalLenstraLenstraLovasz( }, }; } + +/** + * Return a 2-D array with ones on the diagonal and zeros elsewhere. + * @param N Number of rows in the output. + * @param M Number of columns in the output. + * @param k Index of the diagonal. + * @returns An array where all elements are equal to zero, except for the `k`-th diagonal, whose values are equal to one. + */ +export function eye(N: number, M?: number, k = 0) { + M ??= N; + const result: number[][] = []; + for (let i = 0; i < N; ++i) { + result.push(Array(M).fill(0)); + if (i + k < M) { + result[i][i + k] = 1; + } + } + return result; +} + +/** + * Return a 2-D array with ones on the diagonal and zeros elsewhere. + * @param N Number of rows in the output. + * @param M Number of columns in the output. + * @param k Index of the diagonal. + * @returns An array where all elements are equal to zero, except for the `k`-th diagonal, whose values are equal to one. + */ +export function fractionalEye(N: number, M?: number, k = 0): FractionalMonzo[] { + M ??= N; + const result: FractionalMonzo[] = []; + for (let i = 0; i < N; ++i) { + const row: FractionalMonzo = []; + for (let j = 0; j < M; ++j) { + if (j === i + k) { + row.push(new Fraction(1)); + } else { + row.push(new Fraction(0)); + } + } + result.push(row); + } + return result; +} + +// 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; + } + } + // 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'); + } + } + 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); + } + } + // 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) { + // Row echelon form + for (let y = x + 1; y < height; ++y) { + if (matrix_[y][x].n) { + 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.n) { + throw new Error('Matrix is singular'); + } + } + if (!s.isUnity()) { + matrix_[x] = matrix_[x].map(a => a.div(s)); + result[x] = result[x].map(a => a.div(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; +} + +/** + * Compute the matrix product of two arrays of arrays of numbers. + * @param A The left operand. + * @param B The right operand. + * @returns The matrix product of the operands. + */ +export function matmul(A: number[][], B: number[][]) { + const height = A.length; + let width = 0; + for (const row of B) { + width = Math.max(width, row.length); + } + let n = 0; + for (const row of A) { + n = Math.max(n, row.length); + } + B = [...B]; + while (B.length < n) { + B.push([]); + } + const result: number[][] = []; + for (let i = 0; i < height; ++i) { + const row = Array(width).fill(0); + const rowA = A[i]; + for (let j = 0; j < width; ++j) { + for (let k = 0; k < rowA.length; ++k) { + row[j] += rowA[k] * (B[k][j] ?? 0); + } + } + result.push(row); + } + return result; +} + +/** + * Compute the matrix product of two arrays of arrays of fractions. + * @param A The left operand. + * @param B The right operand. + * @returns The matrix product of the operands. + */ +export function fractionalMatmul( + A: ProtoFractionalMonzo[], + B: ProtoFractionalMonzo[] +): FractionalMonzo[] { + const height = A.length; + let width = 0; + for (const row of B) { + width = Math.max(width, row.length); + } + let n = 0; + for (const row of A) { + n = Math.max(n, row.length); + } + const B_ = B.map(row => row.map(f => new Fraction(f))); + while (B_.length < n) { + B_.push([]); + } + for (const row of B_) { + while (row.length < width) { + row.push(new Fraction(0)); + } + } + const result: FractionalMonzo[] = []; + for (let i = 0; i < height; ++i) { + const row: FractionalMonzo = []; + while (row.length < width) { + row.push(new Fraction(0)); + } + const rowA = A[i].map(f => new Fraction(f)); + for (let j = 0; j < width; ++j) { + for (let k = 0; k < rowA.length; ++k) { + row[j] = row[j].add(rowA[k].mul(B_[k][j])); + } + } + result.push(row); + } + return result; +} + +/** + * Compute the determinant of a matrix. + * @param matrix Array of arrays of numbers to calculate determinant for. + * @returns The determinant. + */ +export function det(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'); + } + matrix = matrix.map(row => [...row]); + let result = 1; + // Put zeros in the lower triangle + 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; + } + } + result *= d; + d = 1 / d; + for (let y = x + 1; y < height; ++y) { + const s = matrix[y][x] * d; + if (s) { + matrix[y] = matrix[y].map((a, i) => a - s * (matrix[x][i] ?? 0)); + } + } + } + return result; +} + +/** + * Compute the determinant of a matrix with rational entries. + * @param matrix Array of arrays of fractions to calculate determinant for. + * @returns The determinant. + */ +export function fractionalDet(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 matrix_ = matrix.map(row => row.map(f => new Fraction(f))); + let result = new Fraction(1); + // Put zeros in the lower triangle + for (let x = 0; x < width; ++x) { + let d = matrix_[x][x]; + if (d === undefined || !d.n) { + // 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.neg(); + break; + } + } + d = matrix_[x][x]; + if (d === undefined || !d.n) { + return new Fraction(0); + } + } + result = result.mul(d); + for (let y = x + 1; y < height; ++y) { + const s = matrix_[y][x].div(d); + if (s.n) { + matrix_[y] = matrix_[y].map((a, i) => a.sub(s.mul(matrix_[x][i] ?? 0))); + } + } + } + return result; +} + +/** + * Transpose a 2-D matrix (swap rows and columns). + * @param matrix Matrix to transpose. + * @returns The transpose. + */ +export function transpose(matrix: number[][]) { + let width = 0; + for (const row of matrix) { + width = Math.max(row.length, width); + } + const result: number[][] = []; + for (let i = 0; i < width; ++i) { + const row: number[] = []; + for (let j = 0; j < matrix.length; ++j) { + row.push(matrix[j][i] ?? 0); + } + result.push(row); + } + return result; +} + +/** + * Transpose a 2-D matrix with rational entries (swap rows and columns). + * @param matrix Matrix to transpose. + * @returns The transpose. + */ +export function fractionalTranspose(matrix: ProtoFractionalMonzo[]) { + let width = 0; + for (const row of matrix) { + width = Math.max(row.length, width); + } + const result: FractionalMonzo[] = []; + for (let i = 0; i < width; ++i) { + const row: FractionalMonzo = []; + for (let j = 0; j < matrix.length; ++j) { + row.push(new Fraction(matrix[j][i] ?? 0)); + } + result.push(row); + } + return result; +}