fix: use relative cutoffs for singular values and residuals

matrix.d.ts

@@ -1393,7 +1393,7 @@ export interface ILinearDependenciesOptions {
   thresholdValue?: number;
 
   /**
-   * If the error is inferior to that threshold, the linear combination found is accepted and the row is dependent from other rows.
+   * If the error, relative to the magnitude of the row being explained, is inferior to that threshold, the linear combination found is accepted and the row is dependent from other rows.
    * @default `10e-10`
    */
   thresholdError?: number;
@@ -1415,7 +1415,7 @@ export function linearDependencies(
 /**
  * Returns inverse of a matrix if it exists or the pseudoinverse.
  * @param matrix
- * @param threshold - Threshold for taking inverse of singular values. Default: `Number.EPSILON`.
+ * @param threshold - Relative threshold for taking inverse of singular values. Singular values smaller than `threshold * max(rows, columns) * largestSingularValue` are treated as zero. Default: `Number.EPSILON`.
  * @returns - The (pseudo)inverted matrix.
  */
 export function pseudoInverse(matrix: MaybeMatrix, threshold?: number): Matrix;

src/__tests__/matrix/linearDependencies.test.js

@@ -26,4 +26,40 @@ describe('Linear Dependencies', () => {
       3,
     );
   });
+
+  it('should not depend on the magnitude of the matrix', () => {
+    // row 3 is 2 * row 1 whatever the scale, so the result must not change
+    const rows = [
+      [1, 2, 3],
+      [4, 5, 6],
+      [2, 4, 6],
+    ];
+    for (const k of [1e-9, 1e-3, 1, 1e3, 1e6, 1e9, 1e12]) {
+      const dependencies = linearDependencies(new Matrix(rows).mul(k));
+      expect(dependencies.to2DArray()).toBeDeepCloseTo(
+        [
+          [0, 0, 0.5],
+          [0, 0, 0],
+          [2, 0, 0],
+        ],
+        6,
+      );
+    }
+  });
+
+  it('should not report dependencies for a full rank matrix', () => {
+    const rows = [
+      [1, 2, 3],
+      [4, 5, 6],
+      [7, 8, 10],
+    ];
+    for (const k of [1e-9, 1e-3, 1, 1e3, 1e6, 1e9, 1e12]) {
+      const dependencies = linearDependencies(new Matrix(rows).mul(k));
+      expect(dependencies.to2DArray()).toStrictEqual([
+        [0, 0, 0],
+        [0, 0, 0],
+        [0, 0, 0],
+      ]);
+    }
+  });
 });

src/__tests__/matrix/utility.test.js

@@ -578,6 +578,122 @@ describe('utility methods', () => {
     expect(result).toStrictEqual([[], []]);
   });
 
+  it('pseudoinverse of rank-deficient matrices', () => {
+    // Actual values calculated by the Numpy library
+    let result = pseudoInverse(
+      new Matrix([
+        [10, 20, 30],
+        [40, 50, 60],
+        [70, 80, 90],
+      ]),
+    ).to2DArray();
+
+    expect(result[0][0]).toBeCloseTo(-6.38888889e-2, 8);
+    expect(result[0][1]).toBeCloseTo(-1.66666667e-2, 8);
+    expect(result[0][2]).toBeCloseTo(3.05555556e-2, 8);
+
+    expect(result[1][0]).toBeCloseTo(-5.55555556e-3, 8);
+    expect(result[1][1]).toBeCloseTo(0, 8);
+    expect(result[1][2]).toBeCloseTo(5.55555556e-3, 8);
+
+    expect(result[2][0]).toBeCloseTo(5.27777778e-2, 8);
+    expect(result[2][1]).toBeCloseTo(1.66666667e-2, 8);
+    expect(result[2][2]).toBeCloseTo(-1.94444444e-2, 8);
+
+    result = pseudoInverse(
+      new Matrix([
+        [1, 2],
+        [2, 4],
+        [3, 6],
+        [4, 8],
+      ]),
+    ).to2DArray();
+
+    expect(result[0][0]).toBeCloseTo(6.66666667e-3, 8);
+    expect(result[0][1]).toBeCloseTo(1.33333333e-2, 8);
+    expect(result[0][2]).toBeCloseTo(2.0e-2, 8);
+    expect(result[0][3]).toBeCloseTo(2.66666667e-2, 8);
+
+    expect(result[1][0]).toBeCloseTo(1.33333333e-2, 8);
+    expect(result[1][1]).toBeCloseTo(2.66666667e-2, 8);
+    expect(result[1][2]).toBeCloseTo(4.0e-2, 8);
+    expect(result[1][3]).toBeCloseTo(5.33333333e-2, 8);
+  });
+
+  it('pseudoinverse is scale invariant', () => {
+    // pinv(k*A) = pinv(A)/k exactly, at any k
+    const matrices = [
+      [
+        [1, 2, 3],
+        [4, 5, 6],
+        [7, 8, 9],
+      ],
+      [
+        [2, 4],
+        [7, 1],
+      ],
+      [
+        [1, 2],
+        [2, 4],
+        [3, 6],
+        [4, 8],
+      ],
+      [
+        [1, 2, 3, 4],
+        [2, 4, 6, 8],
+      ],
+    ];
+
+    for (const rows of matrices) {
+      const reference = pseudoInverse(new Matrix(rows)).to2DArray();
+      for (const k of [1e-17, 1e-12, 1e-6, 1e-2, 1e2, 1e6, 1e12, 1e17]) {
+        const scaled = pseudoInverse(new Matrix(rows).mul(k)).to2DArray();
+        expectCloseRelative(
+          scaled.map((row) => row.map((value) => value * k)),
+          reference,
+        );
+      }
+    }
+  });
+
+  it('pseudoinverse satisfies the Moore-Penrose conditions', () => {
+    const matrices = [
+      [
+        [1, 2, 3],
+        [4, 5, 6],
+        [7, 8, 9],
+      ],
+      [
+        [4, 7],
+        [2, 6],
+      ],
+      [
+        [1, 2],
+        [2, 4],
+        [3, 6],
+        [4, 8],
+      ],
+      [
+        [1, 2, 3, 4],
+        [2, 4, 6, 8],
+      ],
+    ];
+
+    for (const rows of matrices) {
+      for (const k of [1e-8, 1, 1e8]) {
+        const A = new Matrix(rows).mul(k);
+        const P = pseudoInverse(A);
+        const AP = A.mmul(P);
+        const PA = P.mmul(A);
+
+        expectCloseRelative(AP.mmul(A).to2DArray(), A.to2DArray());
+        expectCloseRelative(PA.mmul(P).to2DArray(), P.to2DArray());
+        expectCloseRelative(AP.transpose().to2DArray(), AP.to2DArray());
+        expectCloseRelative(PA.transpose().to2DArray(), PA.to2DArray());
+      }
+    }
+  });
+
   it('isEchelonForm', () => {
     const matrix = new Matrix([
       [1, 0],
@@ -808,3 +924,21 @@ describe('utility methods', () => {
     ]);
   });
 });
+
+// Compares entries relative to the magnitude of the expected matrix, so the
+// same assertion holds for results spanning many orders of magnitude.
+function expectCloseRelative(actual, expected) {
+  let scale = 0;
+  for (const row of expected) {
+    for (const value of row) {
+      scale = Math.max(scale, Math.abs(value));
+    }
+  }
+  for (let i = 0; i < expected.length; i++) {
+    for (let j = 0; j < expected[i].length; j++) {
+      expect(Math.abs(actual[i][j] - expected[i][j])).toBeLessThan(
+        1e-12 * scale,
+      );
+    }
+  }
+}

src/linearDependencies.js

@@ -43,7 +43,9 @@ export function linearDependencies(matrix, options = {}) {
     let Abis = matrix.subMatrixRow(xrange(n, i)).transpose();
     let svd = new SingularValueDecomposition(Abis);
     let x = svd.solve(b);
-    let error = Matrix.sub(b, Abis.mmul(x)).abs().max();
+    // The residual scales with the row it explains; the coefficients don't.
+    let scale = Matrix.abs(b).max() || 1;
+    let error = Matrix.sub(b, Abis.mmul(x)).abs().max() / scale;
     results.setRow(
       i,
       dependenciesOneRow(error, x, i, thresholdValue, thresholdError),

src/pseudoInverse.js

@@ -15,8 +15,12 @@ export function pseudoInverse(matrix, threshold = Number.EPSILON) {
   let V = svdSolution.rightSingularVectors;
   let s = svdSolution.diagonal;
 
+  // Singular values scale with the matrix, so the cutoff must too. Same
+  // tolerance as SVD.rank, and the `rcond * max(s)` rule used by LAPACK.
+  const cutoff = threshold * Math.max(matrix.rows, matrix.columns) * s[0];
+
   for (let i = 0; i < s.length; i++) {
-    if (Math.abs(s[i]) > threshold) {
+    if (Math.abs(s[i]) > cutoff) {
       s[i] = 1.0 / s[i];
     } else {
       s[i] = 0.0;