Idea: linear algebra methods for decomposition, eigenvalues and etc..

[JasonShin]

Hi again =)

Throwing another idea of implementing some useful methods of linalg, mostly inspired by Numpy. Again, anyone who is doing Machine Learning in JS will benefit from Math.js if we can support APIs like:
https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.svd.html#numpy.linalg.svd
https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig

[josdejong]

yeah, this would be awesome of course 😎

[josdejong]

Related: #672

[honeybar]

Hi, just letting you know I am currently working on eigenvalue/vector function.

[josdejong]

awesome @honeybar !

[backspaces]

Any progress?

[josdejong]

I haven't had any update about it.

[m4recek]

Any progress now? We have just "found out" that numeric.js has a bug in calculating eig values. And I was not able to find any trustworthy javascript library.

+ can this be takend out of Math.js docs? Its recommending numeric.js for eig calculations
https://mathjs.org/docs/core/extension.html

Issues with numeric.js
https://github.com/sloisel/numeric/issues?utf8=%E2%9C%93&q=eig

[honeybar]

@josdejong Hi, I am so sorry. It is been a while and I am busy with work. I dont think I am able to complete it. Feel free to pass it over to anyone.

[josdejong]

Thanks for getting back @honeybar, no problem at all!

Anyone interested in implementing eigenvalues and eigenvectors?

[arkajitmandal]

I am interested. I have already written a small diagonalization library : https://github.com/arkajitmandal/DiagonalizeJS
check the demo: https://arkajitmandal.github.io/DiagonalizeJS/

[josdejong]

That is exactly what we're looking for Arkajit! Thanks for your offer. Would be great if you can implement this in mathjs. Note that right now mathjs supports numbers and BigNumbers. An implementation for both would be great but it makes sense to me to start simple with just support for numbers. What do you think is a good approach to integrated this functionality in mathjs?

[arkajitmandal]

I agree, I would start with simple. I will begin with the eigenvalue eigenvectors for real symmetric matrix.
Lets call it "eigs.js" in the matrix folder. What do you think? Once this is done I will do a function for generalized eigenvalue problem (we will call it "eigg.js"?) then one for hermitian complex ("eigh"?). Or do you think I should make a folder under matrix called LA or LinearAlgebra ?

[josdejong]

I think /src/function/matrix/eig.js is indeed the best place to put the new function. Sounds good to do implement the functions one at a time 👍 . If you have questions about the architecture of mathjs please don't hesitate to ask.

[arkajitmandal]

I have started coding. I will ask you if I get some problem. Thanks

[arkajitmandal]

Hi !
I am trying to add the eigs function. And I have some questions regarding this. So far, this is what I have done:

1. created a file src/function/matrix/eigs.js , which is exporing createEigs with name = 'eigs' and some dependencies. Right now the this file only checks whether the input is square or not.
2. Added test/unit-tests/function/matrix/eigs.test.js which now has some preliminary unit tests that passes on npm tests
3. Added "export { createEigs } from './function/matrix/eigs' " to factoriesAny.js

After running build-and-test my unit-tests pass, how ever I get this error:

35 passing (3s)
1 failing

1. mainAny
   should export all functions and constants:

   AssertionError [ERR_ASSERTION]: Missing entry in bundle. Path: ["docs","eigs"], expected type: Object, actual type: undefined

   • expected - actual

   -undefined
   +Object

Could you tell me what am I missing thanks in advance.

[arkajitmandal]

Update:

eigs is now working but I still have "the Missing entry in bundle" error. Thanks

[josdejong]

That looks promising Arkajit 👍

The issue Missing entry in bundle. Path: ["docs","eigs"] is probably because mathjs expects embedded documentation for every function. To solve this add an entry eigs in src/expression/embeddedDocs/embeddedDocs.js. You can look at docs of existing functions to get the idea (and if still unclear just ask me).

[cshaa]

@arkajitmandal This looks awesome, any updates? :)

[arkajitmandal]

The present version of MathJS now has a math.eigs function that can compute eigenvalue/eigenvector for real symmetrix matrix.

[cshaa]

Cool! In my opinion the best way to implement a general-case math.eigs would be to implement symbols (#1732) and an equation solver (#38) for polynomials, and then use the classic algorithm:

function eigs_general(M: Matrix)
{
    let E = math.identity( math.size(M) );
    let l = math.symbols('l');

    let charEqn = math.evaluate( 'det( M - l*E )', {M, E, l} );
    let eigVals = math.solve( charEqn, l );

    let eigVecs = [];
    for (ev of eigVals) {
        mat = math.evaluate('M - ev*E', {M, E});
        eigVecs.push( ...math.kernel(mat) );
    }
    
    return [eigVals, eigVecs];
}

I know it sounds daunting at first, but imho it would benefit math.js as a whole. If you don't have time for that, I should be able to do that some time later this month.

But if you think another approach would be better, go ahead!

[arkajitmandal]

I think I would focus on other parts, .... JosDeJong would be able to let you know if this is something he would like to have for mathjs. :)

[josdejong]

@m93a do you mean like first use symbolic computation to turn det( M - l*E ) into a polynomial equation, and then use a solver to solve the polynomial equation, resulting in the eigenvalues?

I think implementing a solver for polynomial equations will be great to have on it's own, but as far as I know you can't "just" solve every polynomial symbolically, or is it? And how will this scale if you have a large real world matrix with ugly values? Will it feasible computationally to write it out as a polynomial and evaluate it? I have the feeling that it will explode.

Probably going the numerical way is probably more robust and performant (and there are proven algorithms and libraries we can use for inspiration). I don't know for sure though. What do you guys think?

[cshaa]

@josdejong I think that the problem of finding eigenvalues of a matrix is equivalent to that of solving a polynomial. So yes, the polynomial would grow with O(n²), but so would the problem itself. But I haven't studied this problem for long enough, so it's possible I'm missing something important. However, there's a great article about eigenvalue algorithms on Wikipedia. It says, for example, that:

Any monic polynomial is the characteristic polynomial of its companion matrix. Therefore, a general algorithm for finding eigenvalues could also be used to find the roots of polynomials. [It was shown that] for dimensions greater than 4 [it's impossible to do this in finite steps]. For this reason algorithms that exactly calculate eigenvalues in a finite number of steps only exist for a few special classes of matrices. For general matrices, algorithms are iterative, producing better approximate solutions with each iteration.

So it seems that we can take the naïve approach I suggested and move the hard numerical algorithms to the polynomial solver.

I'll take more time to look at it after my QM exam next week 😅️

[josdejong]

That is an excellent Wikipedia page 👍 . What sounds like the easiest solution to me is pick one of the iterative, numeric methods listed on that page, one that works for generic matrices.

It could be interesting to experiment with your idea of writing out the polynomial and solve that (numerically with a root-finding algorithm?). My gut feeling though is that that will be more complicated and very slow. But maybe you can prove me wrong 😄 .

[cshaa]

I started reading Numerical recipes in Fortran 77 and there they say:

Root-searching in the characteristic equation is usually a very poor computational method for finding eigenvalues. We will learn much better ways in this chapter, as well as efficient ways for finding corresponding eigenvectors.

So you were right, there are much better ways to tackle the issue. According to the book, a good library will provide methods that return:

• all eigenvalues and no eigenvectors
• all eigenvalues and some corresponding eigenvectors
• all eigenvalues and all corresponding eigenvectors

Computing eigenvectors takes a lot of time, even when you know the eigenvalues. They also say that the library should have algorithms for these specific cases:

• real, symmetric, tridiagonal
• real, symmetric, banded (only a small number of sub- and superdiagonals are nonzero)
• real, symmetric
• real, nonsymmetric
• complex, Hermitian
• complex, non-Hermitian

I'll try to implement an algorithm that computes all eigenvalues and all eigenvectors for a complex, non-Hermitian matrix, I'm starting issue #1741 for that. But even after I'm done with the algorithm, there will still be a lot to improve.

[josdejong]

Thanks for your research and the update @m93a , sounds like a good start 👍