diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js
index cdf4172fe0..472d544719 100644
--- a/src/function/matrix/eigs.js
+++ b/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} */
diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js
index 2e337fc207..5bfa6944ef 100644
--- a/src/function/matrix/eigs/complexEigs.js
+++ b/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,15 +439,11 @@ 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)
       }
 
@@ -460,14 +451,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul
       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))
     }
diff --git a/src/function/matrix/eigs/realSymetric.js b/src/function/matrix/eigs/realSymmetric.js
similarity index 86%
rename from src/function/matrix/eigs/realSymetric.js
rename to src/function/matrix/eigs/realSymmetric.js
index 57c215f2c6..7e24d406da 100644
--- a/src/function/matrix/eigs/realSymetric.js
+++ b/src/function/matrix/eigs/realSymmetric.js
@@ -27,7 +27,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     let Sij = new Array(N)
     // Sij is Identity Matrix
     for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
+      Sij[i] = Array(N).fill(0)
       Sij[i][i] = 1.0
     }
     // initial error
@@ -40,7 +40,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       Sij = Sij1(Sij, psi, i, j)
       Vab = getAij(x)
     }
-    const Ei = createArray(N, 0) // eigenvalues
+    const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
@@ -55,7 +55,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     let Sij = new Array(N)
     // Sij is Identity Matrix
     for (let i = 0; i < N; i++) {
-      Sij[i] = createArray(N, 0)
+      Sij[i] = Array(N).fill(0)
       Sij[i][i] = 1.0
     }
     // initial error
@@ -68,7 +68,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       Sij = Sij1Big(Sij, psi, i, j)
       Vab = getAijBig(x)
     }
-    const Ei = createArray(N, 0) // eigenvalues
+    const Ei = Array(N).fill(0) // eigenvalues
     for (let i = 0; i < N; i++) {
       Ei[i] = x[i][i]
     }
@@ -101,8 +101,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const N = Sij.length
     const c = Math.cos(theta)
     const s = Math.sin(theta)
-    const Ski = createArray(N, 0)
-    const Skj = createArray(N, 0)
+    const Ski = Array(N).fill(0)
+    const Skj = Array(N).fill(0)
     for (let k = 0; k < N; k++) {
       Ski[k] = c * Sij[k][i] - s * Sij[k][j]
       Skj[k] = s * Sij[k][i] + c * Sij[k][j]
@@ -118,8 +118,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const N = Sij.length
     const c = cos(theta)
     const s = sin(theta)
-    const Ski = createArray(N, bignumber(0))
-    const Skj = createArray(N, bignumber(0))
+    const Ski = Array(N).fill(bignumber(0))
+    const Skj = Array(N).fill(bignumber(0))
     for (let k = 0; k < N; k++) {
       Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
       Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
@@ -138,8 +138,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const s = bignumber(sin(theta))
     const c2 = multiplyScalar(c, c)
     const s2 = multiplyScalar(s, s)
-    const Aki = createArray(N, bignumber(0))
-    const Akj = createArray(N, bignumber(0))
+    const Aki = Array(N).fill(bignumber(0))
+    const Akj = Array(N).fill(bignumber(0))
     // 2cs Hij
     const csHij = multiply(bignumber(2), c, s, Hij[i][j])
     //  Aii
@@ -174,8 +174,8 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
     const s = Math.sin(theta)
     const c2 = c * c
     const s2 = s * s
-    const Aki = createArray(N, 0)
-    const Akj = createArray(N, 0)
+    const Aki = Array(N).fill(0)
+    const Akj = Array(N).fill(0)
     //  Aii
     const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
     const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
@@ -237,10 +237,10 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
   function sorting (E, S) {
     const N = E.length
     const values = Array(N)
-    const vectors = Array(N)
+    const vecs = Array(N)
 
     for (let k = 0; k < N; k++) {
-      vectors[k] = Array(N)
+      vecs[k] = Array(N)
     }
     for (let i = 0; i < N; i++) {
       let minID = 0
@@ -253,29 +253,12 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c
       }
       values[i] = E.splice(minID, 1)[0]
       for (let k = 0; k < N; k++) {
-        vectors[k][i] = S[k][minID]
+        vecs[i][k] = S[k][minID]
         S[k].splice(minID, 1)
       }
     }
-
-    return { values, vectors }
-  }
-
-  /**
-   * Create an array of a certain size and fill all items with an initial value
-   * @param {number} size
-   * @param {number} value
-   * @return {number[]}
-   */
-  function createArray (size, value) {
-    // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
-    const array = new Array(size)
-
-    for (let i = 0; i < size; i++) {
-      array[i] = value
-    }
-
-    return array
+    const eigenvectors = vecs.map((vector, i) => ({ value: values[i], vector }))
+    return { values, eigenvectors }
   }
 
   return main
diff --git a/src/utils/number.js b/src/utils/number.js
index 730f2cf7ac..fd72b57c04 100644
--- a/src/utils/number.js
+++ b/src/utils/number.js
@@ -623,7 +623,7 @@ export function nearlyEqual (x, y, epsilon) {
   if (isFinite(x) && isFinite(y)) {
     // check numbers are very close, needed when comparing numbers near zero
     const diff = Math.abs(x - y)
-    if (diff < DBL_EPSILON) {
+    if (diff <= DBL_EPSILON) {
       return true
     } else {
       // use relative error
diff --git a/test/typescript-tests/testTypes.ts b/test/typescript-tests/testTypes.ts
index 1a777c71c4..e55b5d2c75 100644
--- a/test/typescript-tests/testTypes.ts
+++ b/test/typescript-tests/testTypes.ts
@@ -1282,6 +1282,17 @@ Matrices examples
     assert.strictEqual(math.matrix([1, 2, 3]) instanceof math.Matrix, true)
   }
 
+  // Eigenvalues and eigenvectors (should this be extended?)
+  {
+    const D = [
+      [1, 1],
+      [0, 1],
+    ]
+    const eig = math.eigs(D)
+    assert.ok(math.deepEqual(eig.values, [1, 1]))
+    assert.deepStrictEqual(eig.eigenvectors, [{ value: 1, vector: [1, 0] }])
+  }
+
   // Fourier transform and inverse
   {
     assert.ok(
diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js
index f907b1bd61..d9815d3a27 100644
--- a/test/unit-tests/function/matrix/eigs.test.js
+++ b/test/unit-tests/function/matrix/eigs.test.js
@@ -1,7 +1,7 @@
 import assert from 'assert'
 import math from '../../../../src/defaultInstance.js'
 import approx from '../../../../tools/approx.js'
-const { eigs, add, complex, divide, exp, fraction, matrix, multiply, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math
+const { eigs, add, complex, divide, exp, fraction, matrix, matrixFromColumns, multiply, abs, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math
 
 describe('eigs', function () {
   it('only accepts a square matrix', function () {
@@ -16,23 +16,36 @@ describe('eigs', function () {
   it('follows aiao-mimo', function () {
     const realSymArray = eigs([[1, 0], [0, 1]])
     assert(Array.isArray(realSymArray.values) && typeof realSymArray.values[0] === 'number')
-    assert(Array.isArray(realSymArray.vectors) && typeof realSymArray.vectors[0][0] === 'number')
+    for (let ix = 0; ix < 2; ++ix) {
+      assert(Array.isArray(realSymArray.eigenvectors[ix].vector))
+    }
+    assert(typeof realSymArray.eigenvectors[0].vector[0] === 'number')
 
     const genericArray = eigs([[0, 1], [-1, 0]])
     assert(Array.isArray(genericArray.values) && genericArray.values[0] instanceof Complex)
-    assert(Array.isArray(genericArray.vectors) && genericArray.vectors[0][0] instanceof Complex)
+    for (let ix = 0; ix < 2; ++ix) {
+      assert(Array.isArray(genericArray.eigenvectors[ix].vector))
+    }
+    assert(genericArray.eigenvectors[0].vector[0] instanceof Complex)
 
     const realSymMatrix = eigs(matrix([[1, 0], [0, 1]]))
     assert(realSymMatrix.values instanceof Matrix)
     assert.deepStrictEqual(size(realSymMatrix.values), matrix([2]))
-    assert(realSymMatrix.vectors instanceof Matrix)
-    assert.deepStrictEqual(size(realSymMatrix.vectors), matrix([2, 2]))
+    for (let ix = 0; ix < 2; ++ix) {
+      assert(realSymMatrix.eigenvectors[ix].vector instanceof Matrix)
+      assert.deepStrictEqual(
+        size(realSymMatrix.eigenvectors[ix].vector),
+        matrix([2]))
+    }
 
     const genericMatrix = eigs(matrix([[0, 1], [-1, 0]]))
     assert(genericMatrix.values instanceof Matrix)
     assert.deepStrictEqual(size(genericMatrix.values), matrix([2]))
-    assert(genericMatrix.vectors instanceof Matrix)
-    assert.deepStrictEqual(size(genericMatrix.vectors), matrix([2, 2]))
+    for (let ix = 0; ix < 2; ++ix) {
+      assert(genericMatrix.eigenvectors[ix].vector instanceof Matrix)
+      assert.deepStrictEqual(
+        size(genericMatrix.eigenvectors[ix].vector), matrix([2]))
+    }
   })
 
   it('only accepts a matrix with valid element type', function () {
@@ -117,8 +130,7 @@ describe('eigs', function () {
 
     // inverse iteration is stochastic, check it multiple times
     for (let i = 0; i < 5; i++) {
-      const { vectors } = eigs(m)
-      const eigenRows = transpose(vectors)
+      const eigenRows = eigs(m).eigenvectors.map(obj => obj.vector)
       // if we scale each row to the expected scale, they should match
       for (let j = 0; j < 5; j++) {
         approx.deepEqual(divide(eigenRows[i], eigenRows[i][oneIndex[i]]),
@@ -136,7 +148,12 @@ describe('eigs', function () {
       [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]]
     const ans = eigs(H)
     const E = ans.values
-    const V = ans.vectors
+    for (let j = 0; j < 6; j++) {
+      const v = ans.eigenvectors[j].vector
+      approx.deepEqual(multiply(E[j], v), multiply(H, v))
+    }
+    const Vcols = ans.eigenvectors.map(obj => obj.vector)
+    const V = matrixFromColumns(...Vcols)
     const VtHV = multiply(transpose(V), H, V)
     const Ei = Array(H.length)
     for (let i = 0; i < H.length; i++) {
@@ -151,10 +168,9 @@ describe('eigs', function () {
     const cnt = 0.1
     const Ath = multiply(exp(multiply(complex(0, 1), -cnt)), A)
     const Hth = divide(add(Ath, transpose(Ath)), 2)
-    const { values, vectors } = eigs(Hth)
-    const R = transpose(vectors) // rows are eigenvectors
+    const { values, eigenvectors } = eigs(Hth)
     for (const i of [0, 1, 2]) {
-      const v = R[i]
+      const v = eigenvectors[i].vector
       approx.deepEqual(multiply(Hth, v), multiply(values[i], v))
     }
   })
@@ -169,6 +185,41 @@ describe('eigs', function () {
     )
   })
 
+  it('handles some 2x2 defective matrices', function () {
+    const check = eigs([[2.0, 1.0], [0.0, 2.0]]) // Test case from #2879
+    assert.deepStrictEqual(check, {
+      values: [2, 2],
+      eigenvectors: [{ value: 2, vector: [1, 0] }]
+    })
+    const fromWeb = eigs([[-2, 1], [-1, 0]]) // https://ocw.mit.edu/courses/18-03sc-differential-equations-fall-2011/051316d5fa93f560934d3e410f8d153d_MIT18_03SCF11_s33_8text.pdf
+    assert.strictEqual(fromWeb.eigenvectors.length, 1)
+    const vec = fromWeb.eigenvectors[0].vector
+    approx.equal(vec[0], vec[1])
+  })
+
+  it('handles a 3x3 defective matrix', function () {
+    const fromWeb = eigs([[2, -5, 0], [0, 2, 0], [-1, 4, 1]]) // https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_for_Engineers_(Lebl)/3%3A_Systems_of_ODEs/3.7%3A_Multiple_Eigenvalues
+    assert.strictEqual(fromWeb.eigenvectors.length, 2)
+    const ev = fromWeb.eigenvectors
+    approx.equal(ev[0].value, 1)
+    approx.equal(ev[1].value, 2)
+    approx.equal(ev[0].vector[0], 0)
+    approx.equal(ev[0].vector[1], 0)
+    assert.ok(abs(ev[0].vector[2]) > math.config.epsilon)
+    approx.equal(ev[1].vector[0], -ev[1].vector[2])
+    approx.equal(ev[1].vector[1], 0)
+    const web2 = eigs([[1, 1, 0], [0, 1, 2], [0, 0, 3]]) // https://www2.math.upenn.edu/~moose/240S2013/slides7-31.pdf
+    assert.strictEqual(web2.eigenvectors.length, 2)
+    const ev2 = web2.eigenvectors
+    assert.strictEqual(ev2[0].value, 1)
+    assert.strictEqual(ev2[1].value, 3)
+    assert.strictEqual(ev2[0].vector[1], 0)
+    assert.strictEqual(ev2[0].vector[2], 0)
+    assert.ok(abs(ev2[0].vector[0]) > math.config.epsilon)
+    assert.strictEqual(ev2[1].vector[1], ev2[1].vector[2])
+    approx.equal(ev2[1].vector[1], 2 * ev2[1].vector[0])
+  })
+
   it('diagonalizes matrix with bigNumber', function () {
     const x = [[bignum(1), bignum(0)], [bignum(0), bignum(1)]]
     approx.deepEqual(eigs(x).values, [bignum(1), bignum(1)])
@@ -184,7 +235,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 V = ans.vectors
+    const Vcols = ans.eigenvectors.map(obj => obj.vector)
+    const V = matrixFromColumns(...Vcols)
     const VtHV = multiply(transpose(V), H, V)
     const Ei = Array(H.length)
     for (let i = 0; i < H.length; i++) {
diff --git a/types/index.d.ts b/types/index.d.ts
index c14a05d3f9..a013d1c2d5 100644
--- a/types/index.d.ts
+++ b/types/index.d.ts
@@ -1753,9 +1753,8 @@ declare namespace math {
      * 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,
+     * The eigenvectors are returned as an array of objects, each with a
+     * `value` and a `vector`. 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
      * @param x Matrix to be diagonalized
@@ -1765,7 +1764,13 @@ declare namespace math {
     eigs(
       x: MathCollection,
       prec?: number | BigNumber
-    ): { values: MathCollection; vectors: MathCollection }
+    ): {
+      values: MathCollection
+      eigenvectors: {
+        value: number | BigNumber
+        vector: MathCollection
+      }[]
+    }
 
     /**
      * Compute the matrix exponential, expm(A) = e^A. The matrix must be