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 josdejong#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)