Extend eigs to a general matrix

[cshaa]

As discussed in #1175, right now mathjs can only find eigenvalues of real symmetric matrices. It should also be able to find the eigenvalues and eigenvectors for general complex matrices. I'm consulting Numerical recipes in Fortran 77 when working on this issue.

The steps that need to be implemented in order to solve the eigenproblem are the following:

• Balancing (for numeric stability; perform a similarity transform to make the norms of corresponding row-column pairs close to each other)
• Reduce to Hessenberg form (upper triangular plus one subdiagonal; using Gaussian elimination)
• Perform the QR algorithm for Hessenberg matrices

[josdejong]

👍

[atiyabzafar]

I have been searching for a solution to a technical problem I have encountered. I want to calculate eigenvalues for a non-symmetric matrix. Any help would be appreciated as I am working on a web application that requires me to fund out the eigenvalue of a matrix whose size would go to 50x50. I was using power method to find the largest eigenvalue. But that fails if I have degenerate largest eigenvalue.

Do you have any updates?

[cshaa]

Hey Atiya,
sadly I'm on a tight schedule until Friday the 18th because of my uni. After that I will return to the PR and hopefully finish it in no time, there doesn't seem to be much work left to be done. If you need the solution earlier, you can try fiddling with the code yourself, there's a good chance you'll get it working.

EDIT: To get it working, this line needs to be updated to work with usolveAll. If there aren't any bugs, that's it.

EDIT 2: Sorry, I misread that you only need eigenvalues. Than the code in my PR should be good to go, just modify it to skip the finding of eigenvectors.

[atiyabzafar]

Cool. Thanks I will try fiddling around with it. For now I want the eigenvalues. I would require the eigenvectors as well but not as of yet. I will update here if I am able to get it working. Thanks

[atiyabzafar]

@m93a Hey. So I was fiddling around the code and updated the comlex.js function to only return values and not find vectors. I compiled the build. And tried using it. I got an error in console that larger is not defined. Looking at the source, I can see a comment added
TODO proper comparison of bignum and frac
I am from a physics background and haven't really worked with bignumbers before and I am a beginner in javascript as well. My guess is if I remove the "larger" function that you have used and remove the bignumber facility completely and instead use greater than (>) with numbers. But I am sure this would require too much work.

Just wondering if It is worth jumping into.

[cshaa]

Hey Atiya, try replacing the whole function isSymmetric with let isSymmetric = () => false, maybe it will help 🤷‍♂️️ 😁️
If you find the code too broken, just wait 2 weeks and I'll fix it! Or, if you don't have time for that, try using a different library: ml-matrix should have an eigenproblem algorithm.

[josdejong]

I got an error in console that larger is not defined.

@atiyabzafar this can very well be caused by a merge conflict that still needs to be resolved. If that is so, it could be a matter of adding the function larger to the list with dependencies of the function eigs:

const dependencies = [ ... ]  // add 'larger' to the array

export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ ... }) => { // add larger to the destructuring object
  // ...
}

For reference see this comment and the following few: #1743 (comment).

[stur86]

Just wondering, how is this going? I'm developing a project that needs stuff like this (at the moment I think I can make do with eigs, but I'm not sure about future requirements) and had started by relying on Numeric.js. I am considering switching to Math.js as that project is dead and I just found out it has a fatal bug when dealing with matrices that have degenerate eigenvalues. I can't figure out what algorithm exactly they're using as I don't know them all and the code is poorly commented (it uses an upper Hessenberg transformation followed by QR decomposition, but that's all I can gleam without looking through it in detail). If I knew in the future support for all sorts of diagonalisations was guaranteed I'd be a lot happier to switch over. Also I could potentially give a hand - though I haven't ever implemented personally a diagonalisation algorithm beyond basic power iterations with shifts.

[cshaa]

@stur86 I am sorry for my late response. The algorithm is very near completion, see the PR for details. I also believe that my code is well-written and well documented, so you'll have a fine time reading it, if your comment is still current.

[atiyabzafar]

Hey, replying to this as I haven't been able to solve eigenvalue for a simple matrix for some reason. As you can see above, I have been trying to get eigenvalues and vectors for general matrices. I saw that this was closed and merged in #1743. So I used the latest build, i.e. v9.4.1, and tried calculating the values for a simple matrix. I ran into a complication right away. Earlier I had written a code from scratch that was doing the same. I thought of getting rid of the bulk code from my app by using math.js.

https://js.do/atizaf/eigen1

The above link would lead you to the example. It seems to me that eigs() is not able to calculate eigenvectors for complex eigenvalues. Is that the case?

EDIT: I just checked #1743 , that @josdejong has mentioned the algorithm would still fail 50% of the time for the vector. My own implementation for QR failed miserably when i was working on it last year. And the project just halted. I would like to know what exactly is the issue here.

[cshaa]

Hey! 😁️
Sorry for the inconvenience, this is an artifact of the current eigenvector finder.

Since we don't have an iterative eigenvector finder yet, we search for eigenvectors using the equation (A − λE) v = 0. But sometimes the numeric error builds up, so that v is no longer a precise solution to that equation. This problem will be solved once somebody implements an iterative eigenvalue finder (see issue #2179).

If you only needed the eigenvalues, the issue could be side-stepped by not computing the eigenvectors... if we had an API for that. Sadly, we still don't (see issue #2180).

(On a completely unrelated note, you might find the math.zeros method interesting – maybe it will save you some time setting up matrices. Here's how I'd refactor your example code: https://runkit.com/embed/2ggld1mjham3 .)

If you really want to know the eigenvalues and can't wait for one of the two issues to be resolved, you can either use a different library, or use code like this: https://runkit.com/embed/qwm8i4qz5okm . This takes advantage of the fact that even if the algorithm fails, all the eigenvalues and eigenvectors that have been found so far will be passed as properties on the Error.

EDIT: I've just merged #2237. Expect your code to work in the next minor version.

[atiyabzafar]

Thanks for the reply @m93a I look forward to checking out the new merger (#2237 ) for sure. Perhaps I'll throw it matrices to test the new reverse iterative algorithms. I will update here once I have it in the next few days.
I have actually been using math.zeros in the main code. Since I am simply adapting from a python script that I wrote for my thesis where I use NumPy's zeros function.