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