fix: Find eigenvectors of defective matrices

AUTHORS

@@ -233,6 +233,5 @@ BuildTools <anikpatel1322@gmail.com>
 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 <vrushu00@gmail.com>
 
 # Generated by tools/update-authors.js

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))
     }