feat: bound the number of sweeps the SVD spends on one singular value

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([

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;
 }
 
 /**

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$/);
+    }
+  });
+});

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]),