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 josdejong#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 josdejong#3036.

  Resolves josdejong#2879.
  Resolves josdejong#2927.
  Resolves josdejong#3014.

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,63 +54,69 @@ 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<EVobj>}} 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 computeValuesAndVectors(matrix(x)) },
     'Array, number|BigNumber': function (x, prec) {
-      const mat = matrix(x)
-      return computeValuesAndVectors(mat, prec)
+      return computeValuesAndVectors(matrix(x), prec)
     },
-
     Matrix: function (mat) {
-      const { values, vectors } = computeValuesAndVectors(mat)
-      return {
-        values: matrix(values),
-        vectors: matrix(vectors)
-      }
+      return computeValuesAndVectors(mat, undefined, true)
     },
-
     'Matrix, number|BigNumber': function (mat, prec) {
-      const { values, vectors } = computeValuesAndVectors(mat, prec)
-      return {
-        values: matrix(values),
-        vectors: matrix(vectors)
-      }
+      return computeValuesAndVectors(mat, prec, true)
     }
   })
 
-  function computeValuesAndVectors (mat, prec) {
+  function computeValuesAndVectors (mat, prec, matricize = false) {
     if (prec === undefined) {
       prec = config.epsilon
     }
 
-    const size = mat.size()
+    let result
+    const arr = mat.toArray()
+    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)
 
       if (isSymmetric(arr, N, prec)) {
         const type = coerceTypes(mat, arr, N)
-        return doRealSymetric(arr, N, prec, type)
+        result = doRealSymmetric(arr, N, prec, type)
       }
     }
-
-    const type = coerceTypes(mat, arr, N)
-    return doComplexEigs(arr, N, prec, type)
+    if (!result) {
+      const type = coerceTypes(mat, arr, N)
+      result = doComplexEigs(arr, N, prec, type)
+    }
+    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
   }
 
   /** @return {boolean} */

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
@@ -46,14 +46,13 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     // (So U = C⁻¹ 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 }
   }
 
   /**
@@ -350,7 +349,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))
           }
@@ -433,9 +431,6 @@ 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
@@ -444,30 +439,19 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
       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)
       }
 
       // Transform back into original array coordinates
       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: λ, vector: flatten(v) })))
     }
 
     return vectors
@@ -514,19 +498,17 @@ 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 - λI
     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]]
+    // N⃗₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
+    // N⃗₁ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
+
+    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]]
     }
   }
 
@@ -614,12 +596,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,7 +653,9 @@ 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) {
+    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∥² w
       v = subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w))
     }