feat: Add option to eigs() to turn off eigenvector computation

  For large matrices, the eigenvector computation can be noticeably expensive
  and so it's worthwhile to have a way to turn it off if the eigenvectors
  will not be used.
  Resolves #2180.

src/expression/embeddedDocs/function/matrix/eigs.js

@@ -4,9 +4,10 @@ export const eigsDocs = {
   syntax: [
     'eigs(x)'
   ],
-  description: 'Calculate the eigenvalues and eigenvectors of a real symmetric matrix',
+  description: 'Calculate the eigenvalues and optionally eigenvectors of a square matrix',
   examples: [
-    'eigs([[5, 2.3], [2.3, 1]])'
+    'eigs([[5, 2.3], [2.3, 1]])',
+    'eigs([[1, 2, 3], [4, 5, 6], [7, 8, 9]], { precision: 1e-6, eigenvectors: false }'
   ],
   seealso: [
     'inv'

src/function/matrix/eigs.js

@@ -13,7 +13,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
   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 square matrix.
+   * Compute eigenvalues and optionally 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
@@ -31,9 +31,21 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    * in that case, however, you may still find useful information
    * in `err.values` and `err.vectors`.
    *
+   * Note that the 'precision' option does not directly specify the _accuracy_
+   * of the returned eigenvalues. Rather, it determines how small an entry
+   * of the iterative approximations to an upper triangular matrix must be
+   * in order to be considered zero. The actual accuracy of the returned
+   * eigenvalues may be greater or less than the precision, depending on the
+   * conditioning of the matrix and how far apart or close the actual
+   * eigenvalues are. Note that currently, relatively simple, "traditional"
+   * methods of eigenvalue computation are being used; this is not a modern,
+   * high-precision eigenvalue computation. That said, it should typically
+   * produce fairly reasonable results.
+   *
    * Syntax:
    *
    *     math.eigs(x, [prec])
+   *     math.eigs(x, {options})
    *
    * Examples:
    *
@@ -47,14 +59,17 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
    *     const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H if possible
    *     E[0] == UTxHxU[0][0]  // returns true always
    *
+   *     // Compute only approximate eigenvalues:
+   *     const {values} = eigs(H, {eigenvectors: false, precision: 1e-6})
+   *
    * See also:
    *
    *     inv
    *
    * @param {Array | Matrix} x  Matrix to be diagonalized
    *
-   * @param {number | BigNumber} [prec] Precision, default value: 1e-15
-   * @return {{values: Array|Matrix, eigenvectors: Array<EVobj>}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects.
+   * @param {number | BigNumber | OptsObject} [opts] Object with keys `precision`, defaulting to config.epsilon, and `eigenvectors`, defaulting to true and specifying whether to compute eigenvectors. If just a number, specifies precision.
+   * @return {{values: Array|Matrix, eigenvectors?: Array<EVobj>}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects. The eigenvectors property is undefined if eigenvectors were not requested.
    *
    */
   return typed('eigs', {
@@ -67,38 +82,46 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
     // is a roundabout way of doing type detection.
     Array: function (x) { return doEigs(matrix(x)) },
     'Array, number|BigNumber': function (x, prec) {
-      return doEigs(matrix(x), prec)
+      return doEigs(matrix(x), { precision: prec })
     },
+    'Array, Object' (x, opts) { return doEigs(matrix(x), opts) },
     Matrix: function (mat) {
-      return doEigs(mat, undefined, true)
+      return doEigs(mat, { matricize: true })
     },
     'Matrix, number|BigNumber': function (mat, prec) {
-      return doEigs(mat, prec, true)
+      return doEigs(mat, { precision: prec, matricize: true })
+    },
+    'Matrix, Object': function (mat, opts) {
+      const useOpts = { matricize: true }
+      Object.assign(useOpts, opts)
+      return doEigs(mat, useOpts)
     }
   })
 
-  function doEigs (mat, prec, matricize = false) {
-    const result = computeValuesAndVectors(mat, prec)
-    if (matricize) {
+  function doEigs (mat, opts = {}) {
+    const computeVectors = 'eigenvectors' in opts ? opts.eigenvectors : true
+    const prec = opts.precision ?? config.epsilon
+    const result = computeValuesAndVectors(mat, prec, computeVectors)
+    if (opts.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')
+      if (computeVectors) {
+        result.eigenvectors = result.eigenvectors.map(({ value, vector }) =>
+          ({ value, vector: matrix(vector) }))
       }
-    })
+    }
+    if (computeVectors) {
+      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
-    }
-
+  function computeValuesAndVectors (mat, prec, computeVectors) {
     const arr = mat.toArray() // NOTE: arr is guaranteed to be unaliased
     // and so safe to modify in place
     const asize = mat.size()
@@ -114,12 +137,12 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config,
 
       if (isSymmetric(arr, N, prec)) {
         const type = coerceTypes(mat, arr, N) // modifies arr by side effect
-        return doRealSymmetric(arr, N, prec, type)
+        return doRealSymmetric(arr, N, prec, type, computeVectors)
       }
     }
 
     const type = coerceTypes(mat, arr, N) // modifies arr by side effect
-    return doComplexEigs(arr, N, prec, type)
+    return doComplexEigs(arr, N, prec, type, computeVectors)
   }
 
   /** @return {boolean} */

src/function/matrix/eigs/realSymmetric.js

@@ -7,28 +7,31 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
    * @param {number} prec
    * @param {'number' | 'BigNumber'} type
    */
-  function main (arr, N, prec = config.epsilon, type) {
+  function main (arr, N, prec = config.epsilon, type, computeVectors) {
     if (type === 'number') {
-      return diag(arr, prec)
+      return diag(arr, prec, computeVectors)
     }
 
     if (type === 'BigNumber') {
-      return diagBig(arr, prec)
+      return diagBig(arr, prec, computeVectors)
     }
 
     throw TypeError('Unsupported data type: ' + type)
   }
 
   // diagonalization implementation for number (efficient)
-  function diag (x, precision) {
+  function diag (x, precision, computeVectors) {
     const N = x.length
     const e0 = Math.abs(precision / N)
     let psi
-    let Sij = new Array(N)
-    // Sij is Identity Matrix
-    for (let i = 0; i < N; i++) {
-      Sij[i] = Array(N).fill(0)
-      Sij[i][i] = 1.0
+    let Sij
+    if (computeVectors) {
+      Sij = new Array(N)
+      // Sij is Identity Matrix
+      for (let i = 0; i < N; i++) {
+        Sij[i] = Array(N).fill(0)
+        Sij[i][i] = 1.0
+      }
     }
     // initial error
     let Vab = getAij(x)
@@ -37,26 +40,29 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       const j = Vab[0][1]
       psi = getTheta(x[i][i], x[j][j], x[i][j])
       x = x1(x, psi, i, j)
-      Sij = Sij1(Sij, psi, i, j)
+      if (computeVectors) Sij = Sij1(Sij, psi, i, j)
       Vab = getAij(x)
     }
     const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
-    return sorting(clone(Ei), clone(Sij))
+    return sorting(clone(Ei), Sij, computeVectors)
   }
 
   // diagonalization implementation for bigNumber
-  function diagBig (x, precision) {
+  function diagBig (x, precision, computeVectors) {
     const N = x.length
     const e0 = abs(precision / N)
     let psi
-    let Sij = new Array(N)
-    // Sij is Identity Matrix
-    for (let i = 0; i < N; i++) {
-      Sij[i] = Array(N).fill(0)
-      Sij[i][i] = 1.0
+    let Sij
+    if (computeVectors) {
+      Sij = new Array(N)
+      // Sij is Identity Matrix
+      for (let i = 0; i < N; i++) {
+        Sij[i] = Array(N).fill(0)
+        Sij[i][i] = 1.0
+      }
     }
     // initial error
     let Vab = getAijBig(x)
@@ -65,15 +71,15 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       const j = Vab[0][1]
       psi = getThetaBig(x[i][i], x[j][j], x[i][j])
       x = x1Big(x, psi, i, j)
-      Sij = Sij1Big(Sij, psi, i, j)
+      if (computeVectors) Sij = Sij1Big(Sij, psi, i, j)
       Vab = getAijBig(x)
     }
     const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
     // return [clone(Ei), clone(Sij)]
-    return sorting(clone(Ei), clone(Sij))
+    return sorting(clone(Ei), Sij, computeVectors)
   }
 
   // get angle
@@ -234,13 +240,15 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
   }
 
   // sort results
-  function sorting (E, S) {
+  function sorting (E, S, computeVectors) {
     const N = E.length
     const values = Array(N)
-    const vecs = Array(N)
-
-    for (let k = 0; k < N; k++) {
-      vecs[k] = Array(N)
+    let vecs
+    if (computeVectors) {
+      vecs = Array(N)
+      for (let k = 0; k < N; k++) {
+        vecs[k] = Array(N)
+      }
     }
     for (let i = 0; i < N; i++) {
       let minID = 0
@@ -252,11 +260,14 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
         }
       }
       values[i] = E.splice(minID, 1)[0]
-      for (let k = 0; k < N; k++) {
-        vecs[i][k] = S[k][minID]
-        S[k].splice(minID, 1)
+      if (computeVectors) {
+        for (let k = 0; k < N; k++) {
+          vecs[i][k] = S[k][minID]
+          S[k].splice(minID, 1)
+        }
       }
     }
+    if (!computeVectors) return { values }
     const eigenvectors = vecs.map((vector, i) => ({ value: values[i], vector }))
     return { values, eigenvectors }
   }

test/unit-tests/function/matrix/eigs.test.js

@@ -94,12 +94,18 @@ describe('eigs', function () {
         [1.0, 1.0, 1.0]]).values,
     [0, 0, 3]
     )
-    approx.deepEqual(eigs(
+    const sym4 =
       [[0.6163396801190624, -3.8571699139231796, 2.852995822026198, 4.1957619745869845],
         [-3.8571699139231796, 0.7047577966772156, 0.9122549659760404, 0.9232933211541949],
         [2.852995822026198, 0.9122549659760404, 1.6598316026960402, -1.2931270747054358],
-        [4.1957619745869845, 0.9232933211541949, -1.2931270747054358, -4.665994662426116]]).values,
-    [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432]
+        [4.1957619745869845, 0.9232933211541949, -1.2931270747054358, -4.665994662426116]]
+    const fullValues = eigs(sym4).values
+    approx.deepEqual(fullValues,
+      [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432]
+    )
+    assert.deepStrictEqual(
+      fullValues,
+      eigs(sym4, { eigenvectors: false }).values
     )
   })
 
@@ -147,6 +153,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 justvalues = eigs(H, { eigenvectors: false })
+    assert.deepStrictEqual(E, justvalues.values)
     testEigenvectors(ans,
       (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v))
     )
@@ -251,6 +259,8 @@ describe('eigs', function () {
       [4.24, -4.68, -3.33, 1.67, 2.80, 2.73],
       [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]])
     const ans = eigs(H)
+    const justvalues = eigs(H, { eigenvectors: false })
+    assert.deepStrictEqual(ans.values, justvalues.values)
     const E = ans.values
     const Vcols = ans.eigenvectors.map(obj => obj.vector)
     const V = matrixFromColumns(...Vcols)