diff --git a/matrix.d.ts b/matrix.d.ts index a418fe8..b58b1e1 100644 --- a/matrix.d.ts +++ b/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; diff --git a/src/__tests__/matrix/linearDependencies.test.js b/src/__tests__/matrix/linearDependencies.test.js index 019cbd0..2255379 100644 --- a/src/__tests__/matrix/linearDependencies.test.js +++ b/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], + ]); + } + }); }); diff --git a/src/__tests__/matrix/utility.test.js b/src/__tests__/matrix/utility.test.js index c3aef3c..b1893a3 100644 --- a/src/__tests__/matrix/utility.test.js +++ b/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, + ); + } + } +} diff --git a/src/linearDependencies.js b/src/linearDependencies.js index abb3494..c2024c8 100644 --- a/src/linearDependencies.js +++ b/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), diff --git a/src/pseudoInverse.js b/src/pseudoInverse.js index faed62f..d603437 100644 --- a/src/pseudoInverse.js +++ b/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;