diff --git a/README.md b/README.md index 0af64f0..574a92d 100644 --- a/README.md +++ b/README.md @@ -168,6 +168,7 @@ const { LuDecomposition, CholeskyDecomposition, EigenvalueDecomposition, + SingularValueDecomposition, } = require('ml-matrix'); ``` #### Inverse and Pseudo-inverse @@ -285,6 +286,22 @@ var real = e.realEigenvalues; var imaginary = e.imaginaryEigenvalues; var vectors = e.eigenvectorMatrix; ``` +##### Singular Value Decomposition +```js +var A = new Matrix([ + [2, 3, 5], + [4, 1, 6], + [1, 3, 0], +]); + +var svd = new SingularValueDecomposition(A); +var U = svd.leftSingularVectors; +var s = svd.diagonal; +var V = svd.rightSingularVectors; +// U * diag(s) * V.transpose() gives A back +``` +The decomposition settles one singular value at a time through repeated sweeps. `maxIterations`, 100 by default, caps how many sweeps any one value gets. If another sweep would exceed the cap, the decomposition reports that it did not converge instead of continuing without a bound. + #### Linear dependencies ```js var A = new Matrix([ diff --git a/matrix.d.ts b/matrix.d.ts index c6f6153..b1bb03e 100644 --- a/matrix.d.ts +++ b/matrix.d.ts @@ -1478,6 +1478,13 @@ export interface ISVDOptions { * @default `false` */ autoTranspose?: boolean; + + /** + * Maximum number of sweeps spent settling any single singular value. If another sweep + * would exceed the limit, the decomposition throws instead of continuing without a bound. + * @default `100` + */ + maxIterations?: number; } /** diff --git a/src/__tests__/decompositions/svdMaxIterations.test.js b/src/__tests__/decompositions/svdMaxIterations.test.js new file mode 100644 index 0000000..68b02e1 --- /dev/null +++ b/src/__tests__/decompositions/svdMaxIterations.test.js @@ -0,0 +1,68 @@ +import { XSadd } from 'ml-xsadd'; +import { describe, it, expect } from 'vitest'; + +import { Matrix, SingularValueDecomposition } from '../..'; + +const wellBehaved = new Matrix([ + [1, 2, 3], + [4, 5, 6], + [7, 8, 10], + [2, 9, 4], +]); + +function randomMatrix(rows, columns, seed) { + return Matrix.rand(rows, columns, { random: new XSadd(seed).random }); +} + +describe('SVD sweep ceiling', () => { + it('a decomposition that converges is untouched by the default', () => { + const svd = new SingularValueDecomposition(wellBehaved); + const product = svd.leftSingularVectors + .mmul(Matrix.diag(svd.diagonal)) + .mmul(svd.rightSingularVectors.transpose()); + for (let i = 0; i < wellBehaved.rows; i++) { + for (let j = 0; j < wellBehaved.columns; j++) { + expect(product.get(i, j)).toBeCloseTo(wellBehaved.get(i, j), 10); + } + } + }); + + it('a ceiling of one sweep is reported rather than passed over', () => { + expect( + () => + new SingularValueDecomposition(randomMatrix(40, 25, 42), { + maxIterations: 1, + }), + ).toThrow( + /^SVD did not converge after 1 iteration on a single singular value$/, + ); + }); + + it('the ceiling is counted per singular value, not over the whole run', () => { + expect( + () => + new SingularValueDecomposition(randomMatrix(200, 150, 43), { + maxIterations: 20, + }), + ).not.toThrow(); + }); + + it('a generous ceiling behaves like the default', () => { + const bounded = new SingularValueDecomposition(wellBehaved, { + maxIterations: 1000, + }); + const plain = new SingularValueDecomposition(wellBehaved); + expect(bounded.diagonal).toStrictEqual(plain.diagonal); + }); + + it('rejects a ceiling that is not a positive integer', () => { + for (const value of [0, -1, 2.5, '10', null]) { + expect( + () => + new SingularValueDecomposition(wellBehaved, { + maxIterations: value, + }), + ).toThrow(/^maxIterations must be a positive integer$/); + } + }); +}); diff --git a/src/dc/svd.js b/src/dc/svd.js index 6b78606..3d44f5e 100644 --- a/src/dc/svd.js +++ b/src/dc/svd.js @@ -18,8 +18,13 @@ export default class SingularValueDecomposition { computeLeftSingularVectors = true, computeRightSingularVectors = true, autoTranspose = false, + maxIterations = 100, } = options; + if (!Number.isInteger(maxIterations) || maxIterations < 1) { + throw new RangeError('maxIterations must be a positive integer'); + } + let wantu = Boolean(computeLeftSingularVectors); let wantv = Boolean(computeRightSingularVectors); @@ -302,6 +307,15 @@ export default class SingularValueDecomposition { break; } case 3: { + // iter resets when the current singular value settles. Refuse the + // next sweep once this value has consumed its full allowance. + if (iter >= maxIterations) { + const iterationWord = + maxIterations === 1 ? 'iteration' : 'iterations'; + throw new Error( + `SVD did not converge after ${maxIterations} ${iterationWord} on a single singular value`, + ); + } const scale = Math.max( Math.abs(s[p - 1]), Math.abs(s[p - 2]),