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

* 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.

* fix: Add test for precision in options arg of eigs

  And also a fix for a small bug that the new test uncovered.

* test: check eigs with matrix and options

* refactor: remove dead code from complexEigs.js

* fix: add new signatures of eigs to typescript

* test: ensure eigenvectors property not present with eigenvectors: false option

* fix: correct balancing code in complexEigs

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/complexEigs.js

@@ -10,11 +10,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
    *
    * @returns {{ values: number[], vectors: number[][] }}
    */
-  function complexEigs (arr, N, prec, type, findVectors) {
-    if (findVectors === undefined) {
-      findVectors = true
-    }
-
+  function complexEigs (arr, N, prec, type, findVectors = true) {
     // TODO check if any row/col are zero except the diagonal
 
     // make sure corresponding rows and columns have similar magnitude
@@ -46,13 +42,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
     // (So U = C^-1 arr C and the relationship between current arr
     // and original A is unchanged.)
 
-    let eigenvectors
-
     if (findVectors) {
-      eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type)
+      const eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type)
+      return { values, eigenvectors }
     }
 
-    return { values, eigenvectors }
+    return { values }
   }
 
   /**
@@ -95,9 +90,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
 
         for (let j = 0; j < N; j++) {
           if (i === j) continue
-          const c = abs(arr[i][j]) // should be real
-          colNorm = addScalar(colNorm, c)
-          rowNorm = addScalar(rowNorm, c)
+          colNorm = addScalar(colNorm, abs(arr[j][i]))
+          rowNorm = addScalar(rowNorm, abs(arr[i][j]))
         }
 
         if (!equal(colNorm, 0) && !equal(rowNorm, 0)) {
@@ -136,21 +130,21 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
               if (i === j) {
                 continue
               }
-              arr[i][j] = multiplyScalar(arr[i][j], f)
-              arr[j][i] = multiplyScalar(arr[j][i], g)
+              arr[i][j] = multiplyScalar(arr[i][j], g)
+              arr[j][i] = multiplyScalar(arr[j][i], f)
             }
 
             // keep track of transformations
             if (findVectors) {
-              Rdiag[i] = multiplyScalar(Rdiag[i], f)
+              Rdiag[i] = multiplyScalar(Rdiag[i], g)
             }
           }
         }
       }
     }
 
     // return the diagonal row transformation matrix
-    return diag(Rdiag)
+    return findVectors ? diag(Rdiag) : null
   }
 
   /**

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/typescript-tests/testTypes.ts

@@ -1291,6 +1291,13 @@ Matrices examples
     const eig = math.eigs(D)
     assert.ok(math.deepEqual(eig.values, [1, 1]))
     assert.deepStrictEqual(eig.eigenvectors, [{ value: 1, vector: [1, 0] }])
+    const eigvv = math.eigs(D, { precision: 1e-6 })
+    assert.ok(math.deepEqual(eigvv.values, [1, 1]))
+    assert.deepStrictEqual(eigvv.eigenvectors, [{ value: 1, vector: [1, 0] }])
+    const eigv = math.eigs(D, { eigenvectors: false })
+    assert.ok(math.deepEqual(eigv.values, [1, 1]))
+    //@ts-expect-error  ...verify that eigenvectors not expected to be there
+    eigv.eigenvectors
   }
 
   // Fourier transform and inverse