Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

orthogonalComplement in complexEigs.js should return a vector but returns a matrix #2927

Closed
vjancik opened this issue Mar 30, 2023 · 3 comments · Fixed by #3037
Closed

orthogonalComplement in complexEigs.js should return a vector but returns a matrix #2927

vjancik opened this issue Mar 30, 2023 · 3 comments · Fixed by #3037
Labels

Comments

@vjancik
Copy link

vjancik commented Mar 30, 2023

Describe the bug
I'm using math.eigs by passing a 2D row-major matrix array.

function arrMatToRows(arrMat, res = [new Array(4), new Array(4), new Array(4), new Array(4)]) {
  const dim = Math.sqrt(arrMat.length);
  for (let i = 0; i < dim; ++i) {
    for (let j = 0; j < dim; ++j) {
      res[i][j] = arrMat[j*dim + i];
    }
  }
  return res;
}

eigs = math.eigs(arrMatToRows(mZoomInvert, rowBuff));

Line numbers are from https://cdnjs.cloudflare.com/ajax/libs/mathjs/11.7.0/math.js
It throws with

Uncaught RangeError: Expected a column vector, instead got a matrix of size (4, 4)
    o dot.js:46
    Array | DenseMatrix, Array | DenseMatrix dot.js:64
    me typed-function.js:1794
    C complexEigs.js:681
    M complexEigs.js:696
    A complexEigs.js:660
    E complexEigs.js:618
    I complexEigs.js:448
    M complexEigs.js:439
    z complexEigs.js:52
    q eigs.js:100
    Array eigs.js:51
    me typed-function.js:1794
    <anonymous> pen.js:258
    EventListener.handleEvent* pen.js:235
dot.js:46:12

I've traced the problem to orthogonalComplement at complexEigs.js:672.
Here is argument v at the beginning of the function call:
bugBeforeMathjs
and here it is at the end (which ends up throwing in the dot function from the norm call), including the ortho argument :
bugAfterMathjs

To Reproduce
The error is thrown in this demo app but it requires a specific combination of zoom + pan + change rotation origin + rotate while un-zooming to occur. It appears the function only gets invoked when round-off errors accumulate and distort the transformation matrix.
EDIT: I've since changed the code, possibly removing use of math.eigs, you can copy paste the code (from the bug occurrence) from here. You also have to add the previously quoted math.js script in Settings > JS

Sidenote
I've also tried to call math.eigs with a matrix returned from

const temp = math.matrix(mZoomInvert);
temp.resize([4, 4]);
eigs = math.eigs(temp);

but that threw every time. Is that also worth exploring or is that normal?

@josdejong
Copy link
Owner

Thanks for reporting Victor. Can you create a minimal example demonstrating the issue? Something that can be tried out straight in the developer console of https://mathjs.org/ for example, like:

const matrix = [...] // array or nested array containing numbers
const eigs = math.eigs(matrix)

In the examples you post, there are arrays like mZoomInvert and rowBuff from which I do not know what they contain and what type they are.

@gwhitney
Copy link
Collaborator

gwhitney commented Oct 2, 2023

This will be fixed by the fix to #2879.

gwhitney added a commit to gwhitney/mathjs that referenced this issue Oct 2, 2023
  Previously, attempting to take the `eigs` of any defective matrix
  was doomed to fail in an attempt to solve a singular linear system.
  This PR detects the situation (as best as it can given the
  inherent numerical instability of the current methods used) and
  handles it. Note that in such cases, it's not possible to return
  a square matrix whose columns are the eigenvectors corresponding to
  the returned eigenvalues. In light of that fact and issue josdejong#3014, this
  PR also changes the return value of `eigs` so that the eigenvectors
  are passed back in a property `eigenvectors` which is an array of
  plain objects `{value: e, vector: v}`.

  Note that this PR makes the ancillary changes of correcting the
  spelling of the filename which was "realSymetric.js," and replacing
  the now-unnecessary auxiliary function "createArray" therein with
  `Array(size).fill(element)`. The rationale for performing these
  changes not strictly related to the issues at hand is that this
  file is rarely touched and with the level of maintenance hours we have
  at hand, it's more efficient to do these small refactorings in parallel
  with the actual bugfixes, which are orthogonal and so will not be
  obfuscated by this refactor. Note `git diff` does properly track the
  file name change.

  However, it also makes a potentially more pervasive change: in order for
  the numerically-sensitive algorithm to work, it changes the condition
  on when two very close (double) numbers are "nearlyEqual" from differing by
  less than DBL_EPSILON to differing by less than or equal to DBL_EPSILON.
  Although this may change other behaviors than the ones primarily being
  addressed, I believe it is an acceptable change because

  (a) It preserves all tests.
  (b) DBL_EPSILON is well below the standard config.epsilon anyway
  (c) I believe there are extant issues noting the odd/inconsistent
      behavior of nearlyEqual near 0 anyway, so I believe this will
      be overhauled in the future in any case. If so, the eigenvector
      computation will make a good test that a future nearlyEqual
      algorithm is working well.

  To be clear, the direct motivation for the change is that there are
  multiple cases in the eigenvector computation in which a coefficient
  that is "supposed" to be zero comes out to precisely DBL_EPSILON, which
  is fairly unsurprising given that these coefficients are produced by
  subtracting an eigenvalue from a diagonal entry of a matrix, which is
  likely to be essentially equal to that eigenvalue.

  As many tests of defective matrices as I could readily find by web
  searching have been added as unit tests (and one more in the typescript
  type testing). An additional case I found still fails, but in the
  _eigenvalue_ computation rather than the _eigenvector_ search, so that
  was deemed beyond the scope of this PR and has been filed as issue josdejong#3036.

  Resolves josdejong#2879.
  Resolves josdejong#2927.
  Resolves josdejong#3014.
@gwhitney gwhitney added the bug label Oct 3, 2023
@josdejong
Copy link
Owner

Fixed in v12.0.0 now via #3037

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants