From b54aefb3f3709ae3e267aad677a81a5fa32beb07 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Thu, 5 Oct 2023 02:24:51 -0700 Subject: [PATCH] fix: Find eigenvectors of defective matrices (#3037) * fix: Find eigenvectors of defective matrices Previously, attempting to take the `eigs` of any defective matrix was doomed to fail in an attempt to solve a singular linear system. This PR detects the situation (as best as it can given the inherent numerical instability of the current methods used) and handles it. Note that in such cases, it's not possible to return a square matrix whose columns are the eigenvectors corresponding to the returned eigenvalues. In light of that fact and issue #3014, this PR also changes the return value of `eigs` so that the eigenvectors are passed back in a property `eigenvectors` which is an array of plain objects `{value: e, vector: v}`. Note that this PR makes the ancillary changes of correcting the spelling of the filename which was "realSymetric.js," and replacing the now-unnecessary auxiliary function "createArray" therein with `Array(size).fill(element)`. The rationale for performing these changes not strictly related to the issues at hand is that this file is rarely touched and with the level of maintenance hours we have at hand, it's more efficient to do these small refactorings in parallel with the actual bugfixes, which are orthogonal and so will not be obfuscated by this refactor. Note `git diff` does properly track the file name change. However, it also makes a potentially more pervasive change: in order for the numerically-sensitive algorithm to work, it changes the condition on when two very close (double) numbers are "nearlyEqual" from differing by less than DBL_EPSILON to differing by less than or equal to DBL_EPSILON. Although this may change other behaviors than the ones primarily being addressed, I believe it is an acceptable change because (a) It preserves all tests. (b) DBL_EPSILON is well below the standard config.epsilon anyway (c) I believe there are extant issues noting the odd/inconsistent behavior of nearlyEqual near 0 anyway, so I believe this will be overhauled in the future in any case. If so, the eigenvector computation will make a good test that a future nearlyEqual algorithm is working well. To be clear, the direct motivation for the change is that there are multiple cases in the eigenvector computation in which a coefficient that is "supposed" to be zero comes out to precisely DBL_EPSILON, which is fairly unsurprising given that these coefficients are produced by subtracting an eigenvalue from a diagonal entry of a matrix, which is likely to be essentially equal to that eigenvalue. As many tests of defective matrices as I could readily find by web searching have been added as unit tests (and one more in the typescript type testing). An additional case I found still fails, but in the _eigenvalue_ computation rather than the _eigenvector_ search, so that was deemed beyond the scope of this PR and has been filed as issue #3036. Resolves #2879. Resolves #2927. Resolves #3014. * refactor: remove comma that lint now doesn't like * test: add a test for eigs with a precision argument * feat: Use simple shifts in QR eigenvalue iterations that improve convergence Although we might want to use better shifts in the future, we might just use a library instead. But for now I think this: Resolves #2178. Also responds to the review feedback provided in PR #3037. --- AUTHORS | 1 - src/function/matrix/eigs.js | 110 +++++++++++------- src/function/matrix/eigs/complexEigs.js | 107 +++++++++-------- .../{realSymetric.js => realSymmetric.js} | 51 +++----- src/utils/number.js | 2 +- test/typescript-tests/testTypes.ts | 11 ++ test/unit-tests/function/matrix/eigs.test.js | 103 +++++++++++++--- types/index.d.ts | 13 ++- 8 files changed, 242 insertions(+), 156 deletions(-) rename src/function/matrix/eigs/{realSymetric.js => realSymmetric.js} (86%) diff --git a/AUTHORS b/AUTHORS index 5ee87d1522..6c5e9274d7 100644 --- a/AUTHORS +++ b/AUTHORS @@ -233,6 +233,5 @@ BuildTools Anik Patel <74193405+Bobingstern@users.noreply.github.com> Vrushaket Chaudhari <82214275+vrushaket@users.noreply.github.com> Praise Nnamonu <110940850+praisennamonu1@users.noreply.github.com> -vrushaket # Generated by tools/update-authors.js diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js index cdf4172fe0..50edd5cb70 100644 --- a/src/function/matrix/eigs.js +++ b/src/function/matrix/eigs.js @@ -1,22 +1,34 @@ import { factory } from '../../utils/factory.js' import { format } from '../../utils/string.js' import { createComplexEigs } from './eigs/complexEigs.js' -import { createRealSymmetric } from './eigs/realSymetric.js' +import { createRealSymmetric } from './eigs/realSymmetric.js' import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is.js' const name = 'eigs' // The absolute state of math.js's dependency system: -const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot'] -export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolve, usolveAll, im, re, smaller, matrixFromColumns, dot }) => { - const doRealSymetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add }) - const doComplexEigs = createComplexEigs({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) +const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'size', 'reshape', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot'] +export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, size, reshape, qr, usolve, usolveAll, im, re, smaller, matrixFromColumns, dot }) => { + const doRealSymmetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add }) + const doComplexEigs = createComplexEigs({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, size, reshape, qr, inv, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) /** - * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending. - * An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix – - * the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`). - * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information + * Compute eigenvalues and eigenvectors of a square matrix. + * The eigenvalues are sorted by their absolute value, ascending, and + * returned as a vector in the `values` property of the returned project. + * An eigenvalue with algebraic multiplicity k will be listed k times, so + * that the returned `values` vector always has length equal to the size + * of the input matrix. + * + * The `eigenvectors` property of the return value provides the eigenvectors. + * It is an array of plain objects: the `value` property of each gives the + * associated eigenvalue, and the `vector` property gives the eigenvector + * itself. Note that the same `value` property will occur as many times in + * the list provided by `eigenvectors` as the geometric multiplicity of + * that value. + * + * If the algorithm fails to converge, it will throw an error – + * in that case, however, you may still find useful information * in `err.values` and `err.vectors`. * * Syntax: @@ -25,14 +37,15 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, * * Examples: * - * const { eigs, multiply, column, transpose } = math + * const { eigs, multiply, column, transpose, matrixFromColumns } = math * const H = [[5, 2.3], [2.3, 1]] - * const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]} + * const ans = eigs(H) // returns {values: [E1,E2...sorted], eigenvectors: [{value: E1, vector: v2}, {value: e, vector: v2}, ...] * const E = ans.values - * const U = ans.vectors - * multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0)) - * const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H - * E[0] == UTxHxU[0][0] // returns true + * const V = ans.eigenvectors + * multiply(H, V[0].vector)) // returns multiply(E[0], V[0].vector)) + * const U = matrixFromColumns(...V.map(obj => obj.vector)) + * const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H if possible + * E[0] == UTxHxU[0][0] // returns true always * * See also: * @@ -41,62 +54,71 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, * @param {Array | Matrix} x Matrix to be diagonalized * * @param {number | BigNumber} [prec] Precision, default value: 1e-15 - * @return {{values: Array|Matrix, vectors: Array|Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns. + * @return {{values: Array|Matrix, eigenvectors: Array}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects. * */ return typed('eigs', { - Array: function (x) { - const mat = matrix(x) - return computeValuesAndVectors(mat) - }, - + // The conversion to matrix in the first two implementations, + // just to convert back to an array right away in + // computeValuesAndVectors, is unfortunate, and should perhaps be + // streamlined. It is done because the Matrix object carries some + // type information about its entries, and so constructing the matrix + // is a roundabout way of doing type detection. + Array: function (x) { return doEigs(matrix(x)) }, 'Array, number|BigNumber': function (x, prec) { - const mat = matrix(x) - return computeValuesAndVectors(mat, prec) + return doEigs(matrix(x), prec) }, - Matrix: function (mat) { - const { values, vectors } = computeValuesAndVectors(mat) - return { - values: matrix(values), - vectors: matrix(vectors) - } + return doEigs(mat, undefined, true) }, - 'Matrix, number|BigNumber': function (mat, prec) { - const { values, vectors } = computeValuesAndVectors(mat, prec) - return { - values: matrix(values), - vectors: matrix(vectors) - } + return doEigs(mat, prec, true) } }) + function doEigs (mat, prec, matricize = false) { + const result = computeValuesAndVectors(mat, prec) + if (matricize) { + result.values = matrix(result.values) + result.eigenvectors = result.eigenvectors.map(({ value, vector }) => + ({ value, vector: matrix(vector) })) + } + Object.defineProperty(result, 'vectors', { + enumerable: false, // to make sure that the eigenvectors can still be + // converted to string. + get: () => { + throw new Error('eigs(M).vectors replaced with eigs(M).eigenvectors') + } + }) + return result + } + function computeValuesAndVectors (mat, prec) { if (prec === undefined) { prec = config.epsilon } - const size = mat.size() + const arr = mat.toArray() // NOTE: arr is guaranteed to be unaliased + // and so safe to modify in place + const asize = mat.size() - if (size.length !== 2 || size[0] !== size[1]) { - throw new RangeError('Matrix must be square (size: ' + format(size) + ')') + if (asize.length !== 2 || asize[0] !== asize[1]) { + throw new RangeError(`Matrix must be square (size: ${format(asize)})`) } - const arr = mat.toArray() - const N = size[0] + const N = asize[0] if (isReal(arr, N, prec)) { - coerceReal(arr, N) + coerceReal(arr, N) // modifies arr by side effect if (isSymmetric(arr, N, prec)) { - const type = coerceTypes(mat, arr, N) - return doRealSymetric(arr, N, prec, type) + const type = coerceTypes(mat, arr, N) // modifies arr by side effect + return doRealSymmetric(arr, N, prec, type) } } - const type = coerceTypes(mat, arr, N) + const type = coerceTypes(mat, arr, N) // modifies arr by side effect return doComplexEigs(arr, N, prec, type) } diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index 2e337fc207..fc376577f1 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -1,6 +1,6 @@ import { clone } from '../../../utils/object.js' -export function createComplexEigs ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) { +export function createComplexEigs ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, size, reshape, inv, qr, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) { /** * @param {number[][]} arr the matrix to find eigenvalues of * @param {number} N size of the matrix @@ -23,9 +23,9 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const R = balance(arr, N, prec, type, findVectors) // R is the row transformation matrix - // arr = A' = R A R⁻¹, A is the original matrix + // arr = A' = R A R^-1, A is the original matrix // (if findVectors is false, R is undefined) - // (And so to return to original matrix: A = R⁻¹ arr R) + // (And so to return to original matrix: A = R^-1 arr R) // TODO if magnitudes of elements vary over many orders, // move greatest elements to the top left corner @@ -35,7 +35,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // updates the transformation matrix R with new row operationsq // MODIFIES arr by side effect! reduceToHessenberg(arr, N, prec, type, findVectors, R) - // still true that original A = R⁻¹ arr R) + // still true that original A = R^-1 arr R) // find eigenvalues const { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors) @@ -43,17 +43,16 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // values is the list of eigenvalues, C is the column // transformation matrix that transforms arr, the hessenberg // matrix, to upper triangular - // (So U = C⁻¹ arr C and the relationship between current arr + // (So U = C^-1 arr C and the relationship between current arr // and original A is unchanged.) - let vectors + let eigenvectors if (findVectors) { - vectors = findEigenvectors(arr, N, C, R, values, prec, type) - vectors = matrixFromColumns(...vectors) + eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type) } - return { values, vectors } + return { values, eigenvectors } } /** @@ -255,7 +254,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // The Francis Algorithm // The core idea of this algorithm is that doing successive - // A' = Q⁺AQ transformations will eventually converge to block- + // A' = QtAQ transformations will eventually converge to block- // upper-triangular with diagonal blocks either 1x1 or 2x2. // The Q here is the one from the QR decomposition, A = QR. // Since the eigenvalues of a block-upper-triangular matrix are @@ -277,7 +276,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // N×N matrix describing the overall transformation done during the QR algorithm let Qtotal = findVectors ? diag(Array(N).fill(one)) : undefined - // n×n matrix describing the QR transformations done since last convergence + // nxn matrix describing the QR transformations done since last convergence let Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined // last eigenvalue converged before this many steps @@ -290,7 +289,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // Perform the factorization - const k = 0 // TODO set close to an eigenvalue + const k = arr[n - 1][n - 1] // TODO this is apparently a somewhat + // old-fashioned choice; ideally set close to an eigenvalue, or + // perhaps better yet switch to the implicit QR version that is sometimes + // specifically called the "Francis algorithm" that is alluded to + // in the following TODO. (Or perhaps we switch to an independently + // optimized third-party package for the linear algebra operations...) for (let i = 0; i < n; i++) { arr[i][i] = subtract(arr[i][i], k) @@ -350,7 +354,6 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul )) inflateMatrix(Qpartial, N) Qtotal = multiply(Qtotal, Qpartial) - if (n > 2) { Qpartial = diag(Array(n - 2).fill(one)) } @@ -413,18 +416,18 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const uniqueValues = [] const multiplicities = [] - for (const λ of values) { - const i = indexOf(uniqueValues, λ, equal) + for (const lambda of values) { + const i = indexOf(uniqueValues, lambda, equal) if (i === -1) { - uniqueValues.push(λ) + uniqueValues.push(lambda) multiplicities.push(1) } else { multiplicities[i] += 1 } } - // find eigenvectors by solving U − λE = 0 + // find eigenvectors by solving U − lambdaE = 0 // TODO replace with an iterative eigenvector algorithm // (this one might fail for imprecise eigenvalues) @@ -433,26 +436,19 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const b = Array(N).fill(zero) const E = diag(Array(N).fill(one)) - // eigenvalues for which usolve failed (due to numerical error) - const failedLambdas = [] - for (let i = 0; i < len; i++) { - const λ = uniqueValues[i] - const S = subtract(U, multiply(λ, E)) // the characteristic matrix + const lambda = uniqueValues[i] + const S = subtract(U, multiply(lambda, E)) // the characteristic matrix let solutions = usolveAll(S, b) solutions.shift() // ignore the null vector // looks like we missed something, try inverse iteration + // But if that fails, just presume that the original matrix truly + // was defective. while (solutions.length < multiplicities[i]) { const approxVec = inverseIterate(S, N, solutions, prec, type) - - if (approxVec == null) { - // no more vectors were found - failedLambdas.push(λ) - break - } - + if (approxVec === null) { break } // no more vectors were found solutions.push(approxVec) } @@ -460,14 +456,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const correction = multiply(inv(R), C) solutions = solutions.map(v => multiply(correction, v)) - vectors.push(...solutions.map(v => flatten(v))) - } - - if (failedLambdas.length !== 0) { - const err = new Error('Failed to find eigenvectors for the following eigenvalues: ' + failedLambdas.join(', ')) - err.values = values - err.vectors = vectors - throw err + vectors.push( + ...solutions.map(v => ({ value: lambda, vector: flatten(v) }))) } return vectors @@ -478,7 +468,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul * @return {[number,number]} */ function eigenvalues2x2 (a, b, c, d) { - // λ± = ½ trA ± ½ √( tr²A - 4 detA ) + // lambda_+- = 1/2 trA +- 1/2 sqrt( tr^2 A - 4 detA ) const trA = addScalar(a, d) const detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c)) const x = multiplyScalar(trA, 0.5) @@ -489,7 +479,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul /** * For an 2x2 matrix compute the transformation matrix S, - * so that SAS⁻¹ is an upper triangular matrix + * so that SAS^-1 is an upper triangular matrix * @return {[[number,number],[number,number]]} * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html @@ -514,24 +504,22 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul } // matrix is not diagonalizable - // compute off-diagonal elements of N = A - λI - // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ ) - // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ ) - + // compute diagonal elements of N = A - lambdaI const na = subtract(a, l1) - const nb = subtract(b, l1) - const nc = subtract(c, l1) const nd = subtract(d, l1) - if (smaller(abs(nb), prec)) { - return [[na, one], [nc, zero]] + // col(N,2) = 0 implies S = ( col(N,1), e_1 ) + // col(N,2) != 0 implies S = ( col(N,2), e_2 ) + + if (smaller(abs(b), prec) && smaller(abs(nd), prec)) { + return [[na, one], [c, zero]] } else { - return [[nb, zero], [nd, one]] + return [[b, zero], [nd, one]] } } /** - * Enlarge the matrix from n×n to N×N, setting the new + * Enlarge the matrix from nxn to NxN, setting the new * elements to 1 on diagonal and 0 elsewhere */ function inflateMatrix (arr, N) { @@ -614,12 +602,19 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // you better choose a random vector before I count to five let i = 0 - while (true) { + for (; i < 5; ++i) { b = randomOrthogonalVector(N, orthog, type) - b = usolve(A, b) - + try { + b = usolve(A, b) + } catch { + // That direction didn't work, likely because the original matrix + // was defective. But still make the full number of tries... + continue + } if (larger(norm(b), largeNum)) { break } - if (++i >= 5) { return null } + } + if (i >= 5) { + return null // couldn't find any orthogonal vector in the image } // you better converge before I count to ten @@ -664,8 +659,10 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul * Project vector v to the orthogonal complement of an array of vectors */ function orthogonalComplement (v, orthog) { - for (const w of orthog) { - // v := v − (w, v)/∥w∥² w + const vectorShape = size(v) + for (let w of orthog) { + w = reshape(w, vectorShape) // make sure this is just a vector computation + // v := v − (w, v)/|w|^2 w v = subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w)) } diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymmetric.js similarity index 86% rename from src/function/matrix/eigs/realSymetric.js rename to src/function/matrix/eigs/realSymmetric.js index 57c215f2c6..7e24d406da 100644 --- a/src/function/matrix/eigs/realSymetric.js +++ b/src/function/matrix/eigs/realSymmetric.js @@ -27,7 +27,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c let Sij = new Array(N) // Sij is Identity Matrix for (let i = 0; i < N; i++) { - Sij[i] = createArray(N, 0) + Sij[i] = Array(N).fill(0) Sij[i][i] = 1.0 } // initial error @@ -40,7 +40,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c Sij = Sij1(Sij, psi, i, j) Vab = getAij(x) } - const Ei = createArray(N, 0) // eigenvalues + const Ei = Array(N).fill(0) // eigenvalues for (let i = 0; i < N; i++) { Ei[i] = x[i][i] } @@ -55,7 +55,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c let Sij = new Array(N) // Sij is Identity Matrix for (let i = 0; i < N; i++) { - Sij[i] = createArray(N, 0) + Sij[i] = Array(N).fill(0) Sij[i][i] = 1.0 } // initial error @@ -68,7 +68,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c Sij = Sij1Big(Sij, psi, i, j) Vab = getAijBig(x) } - const Ei = createArray(N, 0) // eigenvalues + const Ei = Array(N).fill(0) // eigenvalues for (let i = 0; i < N; i++) { Ei[i] = x[i][i] } @@ -101,8 +101,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const N = Sij.length const c = Math.cos(theta) const s = Math.sin(theta) - const Ski = createArray(N, 0) - const Skj = createArray(N, 0) + const Ski = Array(N).fill(0) + const Skj = Array(N).fill(0) for (let k = 0; k < N; k++) { Ski[k] = c * Sij[k][i] - s * Sij[k][j] Skj[k] = s * Sij[k][i] + c * Sij[k][j] @@ -118,8 +118,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const N = Sij.length const c = cos(theta) const s = sin(theta) - const Ski = createArray(N, bignumber(0)) - const Skj = createArray(N, bignumber(0)) + const Ski = Array(N).fill(bignumber(0)) + const Skj = Array(N).fill(bignumber(0)) for (let k = 0; k < N; k++) { Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j])) Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j])) @@ -138,8 +138,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const s = bignumber(sin(theta)) const c2 = multiplyScalar(c, c) const s2 = multiplyScalar(s, s) - const Aki = createArray(N, bignumber(0)) - const Akj = createArray(N, bignumber(0)) + const Aki = Array(N).fill(bignumber(0)) + const Akj = Array(N).fill(bignumber(0)) // 2cs Hij const csHij = multiply(bignumber(2), c, s, Hij[i][j]) // Aii @@ -174,8 +174,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const s = Math.sin(theta) const c2 = c * c const s2 = s * s - const Aki = createArray(N, 0) - const Akj = createArray(N, 0) + const Aki = Array(N).fill(0) + const Akj = Array(N).fill(0) // Aii const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j] const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j] @@ -237,10 +237,10 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c function sorting (E, S) { const N = E.length const values = Array(N) - const vectors = Array(N) + const vecs = Array(N) for (let k = 0; k < N; k++) { - vectors[k] = Array(N) + vecs[k] = Array(N) } for (let i = 0; i < N; i++) { let minID = 0 @@ -253,29 +253,12 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c } values[i] = E.splice(minID, 1)[0] for (let k = 0; k < N; k++) { - vectors[k][i] = S[k][minID] + vecs[i][k] = S[k][minID] S[k].splice(minID, 1) } } - - return { values, vectors } - } - - /** - * Create an array of a certain size and fill all items with an initial value - * @param {number} size - * @param {number} value - * @return {number[]} - */ - function createArray (size, value) { - // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it) - const array = new Array(size) - - for (let i = 0; i < size; i++) { - array[i] = value - } - - return array + const eigenvectors = vecs.map((vector, i) => ({ value: values[i], vector })) + return { values, eigenvectors } } return main diff --git a/src/utils/number.js b/src/utils/number.js index 730f2cf7ac..fd72b57c04 100644 --- a/src/utils/number.js +++ b/src/utils/number.js @@ -623,7 +623,7 @@ export function nearlyEqual (x, y, epsilon) { if (isFinite(x) && isFinite(y)) { // check numbers are very close, needed when comparing numbers near zero const diff = Math.abs(x - y) - if (diff < DBL_EPSILON) { + if (diff <= DBL_EPSILON) { return true } else { // use relative error diff --git a/test/typescript-tests/testTypes.ts b/test/typescript-tests/testTypes.ts index 1a777c71c4..23a1a34fd9 100644 --- a/test/typescript-tests/testTypes.ts +++ b/test/typescript-tests/testTypes.ts @@ -1282,6 +1282,17 @@ Matrices examples assert.strictEqual(math.matrix([1, 2, 3]) instanceof math.Matrix, true) } + // Eigenvalues and eigenvectors + { + const D = [ + [1, 1], + [0, 1] + ] + const eig = math.eigs(D) + assert.ok(math.deepEqual(eig.values, [1, 1])) + assert.deepStrictEqual(eig.eigenvectors, [{ value: 1, vector: [1, 0] }]) + } + // Fourier transform and inverse { assert.ok( diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index f907b1bd61..f77862d651 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -1,9 +1,14 @@ import assert from 'assert' import math from '../../../../src/defaultInstance.js' import approx from '../../../../tools/approx.js' -const { eigs, add, complex, divide, exp, fraction, matrix, multiply, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math +const { eigs, add, complex, divide, exp, fraction, matrix, matrixFromColumns, multiply, abs, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math describe('eigs', function () { + // helper to examine eigenvectors + function testEigenvectors (soln, predicate) { + soln.eigenvectors.forEach((ev, i) => predicate(ev.vector, i)) + } + it('only accepts a square matrix', function () { assert.throws(function () { eigs(matrix([[1, 2, 3], [4, 5, 6]])) }, /Matrix must be square/) assert.throws(function () { eigs([[1, 2, 3], [4, 5, 6]]) }, /Matrix must be square/) @@ -16,23 +21,30 @@ describe('eigs', function () { it('follows aiao-mimo', function () { const realSymArray = eigs([[1, 0], [0, 1]]) assert(Array.isArray(realSymArray.values) && typeof realSymArray.values[0] === 'number') - assert(Array.isArray(realSymArray.vectors) && typeof realSymArray.vectors[0][0] === 'number') + testEigenvectors(realSymArray, vector => assert(Array.isArray(vector))) + assert(typeof realSymArray.eigenvectors[0].vector[0] === 'number') const genericArray = eigs([[0, 1], [-1, 0]]) assert(Array.isArray(genericArray.values) && genericArray.values[0] instanceof Complex) - assert(Array.isArray(genericArray.vectors) && genericArray.vectors[0][0] instanceof Complex) + testEigenvectors(genericArray, + vector => assert(Array.isArray(vector) && vector[0] instanceof Complex) + ) const realSymMatrix = eigs(matrix([[1, 0], [0, 1]])) assert(realSymMatrix.values instanceof Matrix) assert.deepStrictEqual(size(realSymMatrix.values), matrix([2])) - assert(realSymMatrix.vectors instanceof Matrix) - assert.deepStrictEqual(size(realSymMatrix.vectors), matrix([2, 2])) + testEigenvectors(realSymMatrix, vector => { + assert(vector instanceof Matrix) + assert.deepStrictEqual(size(vector), matrix([2])) + }) const genericMatrix = eigs(matrix([[0, 1], [-1, 0]])) assert(genericMatrix.values instanceof Matrix) assert.deepStrictEqual(size(genericMatrix.values), matrix([2])) - assert(genericMatrix.vectors instanceof Matrix) - assert.deepStrictEqual(size(genericMatrix.vectors), matrix([2, 2])) + testEigenvectors(genericMatrix, vector => { + assert(vector instanceof Matrix) + assert.deepStrictEqual(size(vector), matrix([2])) + }) }) it('only accepts a matrix with valid element type', function () { @@ -117,8 +129,7 @@ describe('eigs', function () { // inverse iteration is stochastic, check it multiple times for (let i = 0; i < 5; i++) { - const { vectors } = eigs(m) - const eigenRows = transpose(vectors) + const eigenRows = eigs(m).eigenvectors.map(obj => obj.vector) // if we scale each row to the expected scale, they should match for (let j = 0; j < 5; j++) { approx.deepEqual(divide(eigenRows[i], eigenRows[i][oneIndex[i]]), @@ -136,7 +147,11 @@ describe('eigs', function () { [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]] const ans = eigs(H) const E = ans.values - const V = ans.vectors + testEigenvectors(ans, + (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v)) + ) + const Vcols = ans.eigenvectors.map(obj => obj.vector) + const V = matrixFromColumns(...Vcols) const VtHV = multiply(transpose(V), H, V) const Ei = Array(H.length) for (let i = 0; i < H.length; i++) { @@ -151,12 +166,10 @@ describe('eigs', function () { const cnt = 0.1 const Ath = multiply(exp(multiply(complex(0, 1), -cnt)), A) const Hth = divide(add(Ath, transpose(Ath)), 2) - const { values, vectors } = eigs(Hth) - const R = transpose(vectors) // rows are eigenvectors - for (const i of [0, 1, 2]) { - const v = R[i] - approx.deepEqual(multiply(Hth, v), multiply(values[i], v)) - } + const example = eigs(Hth) + testEigenvectors(example, (v, i) => + approx.deepEqual(multiply(Hth, v), multiply(example.values[i], v)) + ) }) it('supports fractions', function () { @@ -169,6 +182,61 @@ describe('eigs', function () { ) }) + it('handles some 2x2 defective matrices', function () { + const check = eigs([[2.0, 1.0], [0.0, 2.0]]) // Test case from #2879 + assert.deepStrictEqual(check, { + values: [2, 2], + eigenvectors: [{ value: 2, vector: [1, 0] }] + }) + const fromWeb = eigs([[-2, 1], [-1, 0]]) // https://ocw.mit.edu/courses/18-03sc-differential-equations-fall-2011/051316d5fa93f560934d3e410f8d153d_MIT18_03SCF11_s33_8text.pdf + assert.strictEqual(fromWeb.eigenvectors.length, 1) + const vec = fromWeb.eigenvectors[0].vector + approx.equal(vec[0], vec[1]) + }) + + it('handles a 3x3 defective matrix', function () { + const fromWeb = eigs([[2, -5, 0], [0, 2, 0], [-1, 4, 1]]) // https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_for_Engineers_(Lebl)/3%3A_Systems_of_ODEs/3.7%3A_Multiple_Eigenvalues + assert.strictEqual(fromWeb.eigenvectors.length, 2) + const ev = fromWeb.eigenvectors + approx.equal(ev[0].value, 1) + approx.equal(ev[1].value, 2) + approx.equal(ev[0].vector[0], 0) + approx.equal(ev[0].vector[1], 0) + assert.ok(abs(ev[0].vector[2]) > math.config.epsilon) + approx.equal(ev[1].vector[0], -ev[1].vector[2]) + approx.equal(ev[1].vector[1], 0) + const web2 = eigs([[1, 1, 0], [0, 1, 2], [0, 0, 3]]) // https://www2.math.upenn.edu/~moose/240S2013/slides7-31.pdf + assert.strictEqual(web2.eigenvectors.length, 2) + const ev2 = web2.eigenvectors + assert.strictEqual(ev2[0].value, 1) + assert.strictEqual(ev2[1].value, 3) + assert.strictEqual(ev2[0].vector[1], 0) + assert.strictEqual(ev2[0].vector[2], 0) + assert.ok(abs(ev2[0].vector[0]) > math.config.epsilon) + assert.strictEqual(ev2[1].vector[1], ev2[1].vector[2]) + approx.equal(ev2[1].vector[1], 2 * ev2[1].vector[0]) + }) + + it('accepts a precision argument', function () { + // The following is a matrix with an algebraically triple eigenvalue + // equal to 2 which has a unique eigenvector (up to scale, of course). + // It is from https://web.uvic.ca/~tbazett/diffyqs/sec_multeigen.html + // The iterative eigenvalue calculation currently being used has a + // great deal of difficulty converging. We can use a fine precision, + // but it still doesn't produce good eigenvalues. Hopefully someday + // we'll be able to get closer. + const difficult = [[2, 0, 0], [-1, -1, 9], [0, -1, 5]] + const poor = eigs(difficult, 1e-14) + assert.strictEqual(poor.values.length, 3) + approx.deepEqual(poor.values, [2, 2, 2], 6e-6) + // Note the eigenvectors are junk, so we don't test them. The function + // eigs thinks there are three of them, for example. Hopefully some + // future iteration of mathjs will be able to discover there is really + // only one. + const poorm = eigs(matrix(difficult), 1e-14) + assert.deepStrictEqual(poorm.values.size(), [3]) + }) + it('diagonalizes matrix with bigNumber', function () { const x = [[bignum(1), bignum(0)], [bignum(0), bignum(1)]] approx.deepEqual(eigs(x).values, [bignum(1), bignum(1)]) @@ -184,7 +252,8 @@ describe('eigs', function () { [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]]) const ans = eigs(H) const E = ans.values - const V = ans.vectors + const Vcols = ans.eigenvectors.map(obj => obj.vector) + const V = matrixFromColumns(...Vcols) const VtHV = multiply(transpose(V), H, V) const Ei = Array(H.length) for (let i = 0; i < H.length; i++) { diff --git a/types/index.d.ts b/types/index.d.ts index c14a05d3f9..a013d1c2d5 100644 --- a/types/index.d.ts +++ b/types/index.d.ts @@ -1753,9 +1753,8 @@ declare namespace math { * Compute eigenvalues and eigenvectors of a matrix. * The eigenvalues are sorted by their absolute value, ascending. * An eigenvalue with multiplicity k will be listed k times. - * The eigenvectors are returned as columns of a matrix – the eigenvector - * that belongs to the j-th eigenvalue in the list (eg. values[j]) is the - * j-th column (eg. column(vectors, j)). If the algorithm fails to converge, + * The eigenvectors are returned as an array of objects, each with a + * `value` and a `vector`. If the algorithm fails to converge, * it will throw an error – in that case, however, you may still find useful * information in err.values and err.vectors * @param x Matrix to be diagonalized @@ -1765,7 +1764,13 @@ declare namespace math { eigs( x: MathCollection, prec?: number | BigNumber - ): { values: MathCollection; vectors: MathCollection } + ): { + values: MathCollection + eigenvectors: { + value: number | BigNumber + vector: MathCollection + }[] + } /** * Compute the matrix exponential, expm(A) = e^A. The matrix must be