fix: Find eigenvectors of defective matrices #3037
Conversation
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.
|
Looks like I need to add tests in which the precision argument is used to pass the codecov/patch target. I will try that now. |
|
Oh, and I forgot one more excuse/justification for changing nearlyEqual from |
There was a problem hiding this comment.
Thanks a lot Glen, this looks good! I had expected that it would require more changes in the logic to fix this.
I'm glad with your extra refactoring steps. I think the change in behavior around DBL_EPSILON is no problem.
Just good to mention: this PR contains a breaking change in the API of eigs and will be scheduled for v12.
…ergence Although we might want to use better shifts in the future, we might just use a library instead. But for now I think this: Resolves josdejong#2178. Also responds to the review feedback provided in PR josdejong#3037.
|
OK, I think I've covered the review feedback, and I have incorporated an improvement to #2178 that I think is sufficient for the moment. |
|
Thanks for the updates Glen, this looks good to go. |
|
I'll merge the PR into a new |
|
This fix has been published in |
Previously, attempting to take the
eigsof any defective matrixwas 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 #3014, this
PR also changes the return value of
eigsso that the eigenvectorsare passed back in a property
eigenvectorswhich is an array ofplain 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 thesechanges 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 diffdoes properly track thefile 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 #3036.
Resolves #2879.
Resolves #2927.
Resolves #3014.