Adding eigs function

[arkajitmandal]

This function computes the eigenvalues and eigenvectors of real symmetric matrix.

[arkajitmandal]

Hi !
Some initial points:

• for numbers : I think most of the people want a fast diagonalization function for 2d array with numbers. I think (I could be wrong) the mathjs add or substract, sin, cos function is more expensive than the native Math function due to the type checking. Since diagonalization itself is a very computationally expensive task, for the purpose of practicality having a faster diagonalization for number would be very important.
• for fractions: I am slightly confused on how to deal with fractions. The algorithm used here tries to numerically converge the result and heavily relies on the sin and cos functions. As far as I understand, the sin and cos function did not work with fraction type and therefore I had no choice but to convert it to a number.
• for bigNumber: The only way to perform diagonalization is off-course by using the mathjs add, sin etc functions.
• for mixed : It seems bigNumber and fractions are not compatible, while number and bigNumber and fraction and number are. So diagonalization having mixed elements are possible as long as its either Bignumber and number or number and fraction. But again since fractions don't work for cos and sin, fraction + number implementation would be by converting fraction to numbers. Yet tobe implemented.

Given this, let me know how you want me to go ahead. If I only have bignumber implementation, and remove number implementation then the fraction type might give some trouble. Mixed type is definitely going to be an issue.

[arkajitmandal]

How about these choices:

1. I have two implementation, one for number and bignumber. I take an array and convert to matrix. If the
   datatype is number I use the fast method, if its "fraction" I convert to number and use fast diagonalization, if its 'bignumber' I use the slow method, if its mixed I use slow method by converting whole matrix to bignumber. (Pro: fast for numbers, Con: having more code)
2. Delete number implementation, convert mixed type to BigNumber. (Pro: having less code, Con: slow for numbers) . (Not sure if its possible though)
   Let me know if you want me to do either of them.

[josdejong]

I think you have to add one more pro/con to the choices 1. and 2. : Option 1 will not work for new and/or custom data types out of the box, option 2 will.

I'm quite sure option 2 with mixed content be much slower out of the box, and you're right that regular numbers is the main use case so it's important to optimize for that. However, option 2 does not exclude optimizing for numbers: if eigs is constructed via a factory function:

const createEigs = factory(name, dependencies, ({ add, subtract, ... }) => { ... }

and you pass number-only implementations like add: (x, y) => x + y, the created function should be just as fast as the current implementation for numbers. Except that you can also pass and implementation for add etc which does support BigNumbers. Or other numeric types. So, what I try to say is: this option 2 is much more powerful and flexible than option 1. It will require some work to set this up so that it indeed loads a version which is optimized for numbers.

We're talking about optimization without doing measurements though, so it's all based on assumptions. I think we can go in the following directions:

1. leave it as it is now (2 implementations), it works. In a next iteration we can have another look.
2. drop the number optimized version for now, only keep the mixed implementation (which we assume is slower), and some day in a next iteration have another look at optimization.
3. write out two versions and a (small) benchmark to get a feeling of what the performance of the different approaches is, and decide based on the outcome is. We have three approaches:
   1. a generic version supporting mixed content (basically the current "Big" implementation)
   2. the number optimized version that you have implemented
   3. the idea I propose of instantiating the generic version with number-only implementations for it's dependencies, which in theory should be as fast as the number optimized version (ii).

Would you find it interesting to do a bit of benchmarking? Or shall we leave it as it is now? (Both ok with me)

[arkajitmandal]

How about lets keep the current version. So that people get to use it now. Meanwhile by next week I clean up the code and do some benchmarking as well as code up the option iii. Then, with the results of the benchmarking we can talk again to decide which direction we should move ahead.

What do you think?

[josdejong]

Yes agree, let's be pragmatic here and leave the benchmarking and optimization for a next stage.

I've done some more testing with the function. It works like a charm 👍.

Some more thoughts about naming and API.

1. Most languages I'm familiar with (Matlab, Python, Mathematica) call this function eig. Python even has a separate function eigh for real, symmetric matrices (which is what we have right now, right?). Is there a special reason to call the function eigs? Or is it maybe handier to call it eig since most people will expect that? Or even take the same approach as Python?
2. There are not many functions in mathjs which return multiple results (in this case eigen values and eigen vectors). The function math.qr does, and there we choose to utilize objects instead of arrays, so you get named results back: { Q: ..., R: ...}. I quite like this API since it's very explicit in the output what is what (compared to having to pick the right index from an array). It may be good to choose a similar API for eigs, so it could return an object { values: ..., vectors: ... }. What do you think?
3. So far, all functions in mathjs follow the behavior Array in -> Array out, Matrix in -> Matrix out. It would be neat if eigs also returns Matrix as results when the input is a Matrix. This basically means that you wrap the resulting eigenvalues and eigenvectors in math.matrix(...) before returning.
4. I'm still curious if there is a specific reason to implement a custom sorting function inside eigs. Would it be possible to reuse the existing sort function?

[arkajitmandal]

I have tried to address most of your comments.

1. name : eigh stands for eigenvalue eigenvector for hermitian matrix, the current implementation only applies to real symmetric matrix only (and not complex) and thats why I named it eigs. That said, I plan to add the functionality for complex types of matrix next and at that point we could rename this as eigh. For now lets keep it in that way.

2. The results are now returned as an object type : { values: ..., vectors: ... }

3. Now the function handles Array in -> {values: Array, vectors: Array} out and Matrix in -> {values: Array, vectors: Matrix} out.

4. For sorting I would require the sorted ids since I need to sort the corresponding eigenvalues in the same way. I dont know if I could do that with mathjs sort. If I can do that let me know I would do that in the next iteration.

5. Used addScalar, multiplyScalar etc. wherever I could.

[josdejong]

1. Ah, now I get it. Makes sense 👍

Thanks for all the updates (2, 3, 5), I love it. Agree about sorting, let's see if there are options for improvement there in a later stage, I'm not sure.

Will merge your work now, thanks again Arkajit!

[josdejong]

I've done some minor tweaks in your code: it was failing on IE which doesn't support Array.fill. I've also replaced all occurrences of var with let and const, and I've changed the returned type of .values to be a Matrix too when the input is a Matrix (so both .values and .vectors are a Matrix).

[arkajitmandal]

Great ! I am nearly done with the complex part will add that soon.

[josdejong]

The function eigs is now available in v6.6.0 🎉 . Thanks again Arkajit.

[arkajitmandal]

Thats great... got busy with work... will soon come back to finish the other parts...

[josdejong]

👍 take it easy