[FEATURE REQUEST] `linalg.eig` for asymmetric matrices (and SVD?)

[JamesTev]

I need to be able to find eigenvalues on an asymmetric square matrix (as offered by numpy), although it seems that ulab.linalg.eig only works on symmetric square matrices at the moment. Is there a reason for this? I assume it's for computational reasons.

Also, is implementation of SVD and/or QR factorisation in the pipeline at all? I see you use Givens rotations for linalg.eig so QR shouldn't be too bad.

Thanks for all your amazing work! This is an awesome project.

[v923z]

@JamesTev

I need to be able to find eigenvalues on an asymmetric square matrix (as offered by numpy), although it seems that ulab.linalg.eig only works on symmetric square matrices at the moment. Is there a reason for this? I assume it's for computational reasons.

That was the simplest case to implement. To be honest, you are the first, who has asked a question about eig, so I wasn't even sure, whether people use it at all, and there seemed to be no reason for investing too much time into something that is useless.

The snag is, each type of matrix has its own numerical difficulties, and if you look at serious numerical packages, different methods are used for the different types. In other words, there is no silver bullet.

Also, is implementation of SVD and/or QR factorisation in the pipeline at all?

No, but this doesn't mean that it can't be.

I see you use Givens rotations for linalg.eig so QR shouldn't be too bad.

Do you have a particular algorithm in mind, or would you like to take a stab at the implementation?

[JamesTev]

I see - that makes sense.

In other words, there is no silver bullet.

I'm not too familiar with your particular eig implementation but do you think handling asymmetric matrices would be difficult? Perhaps Schur decomp (see here) could be an alternative? It works for arbitrary square matrices. I'm not too familiar with it but I believe it can be solved using QR factorisation which itself can be computed using Givens rotations.

Do you have a particular algorithm in mind, or would you like to take a stab at the implementation?

At the moment, I'm trying to compute CCA in real time on an ESP32 which requires SVD or being able to find eigen values. Unfortunately, my C skills aren't great so taking a stab at the implementation would take quite some time!

As an aside: I think having fast numpy-like linear algebra functions is really useful and convenient and far easier to use than some of the standard C/C++ implementations. Please keep up the good work.

[v923z]

I see - that makes sense.

In other words, there is no silver bullet.

I'm not too familiar with your particular eig implementation but do you think handling asymmetric matrices would be difficult? Perhaps Schur decomp (see here) could be an alternative? It works for arbitrary square matrices. I'm not too familiar with it but I believe it can be solved using QR factorisation which itself can be computed using Givens rotations.

I actually use Jacobi rotations, but that is off the point. I think QR would be relatively straightforward. Give me a couple of days. But that might not take you to the finish line.

I forgot to mention yesterday that one of the reasons for sticking with symmetric matrices was that their eigenvalues are real, and I didn't want to deal with complex numbers. (I have nothing against them, but they add too much to the firmware.)

At the moment, I'm trying to compute CCA in real time on an ESP32 which requires SVD or being able to find eigen values.

In that case, would an implementation of #188 help you? The only thing missing there is the addition of the keepdims keyword to some of the numerical functions. I believe that is not terribly difficult. I have meant to fix that, but I am a bit tied down with another feature.

Unfortunately, my C skills aren't great so taking a stab at the implementation would take quite some time!

That's fine, but I might ask you to help with the testing.

As an aside: I think having fast numpy-like linear algebra functions is really useful and convenient and far easier to use than some of the standard C/C++ implementations. Please keep up the good work.

This project lives by the interest of the users.

[v923z]

@JamesTev Here is a new issue, you could comment on the specific questions there: #364

[v923z]

@JamesTev I have opened a draft under #366. With that, we should be able to insert a generic eigenvalue solver in linalg. Could you, please, spell out once more, what exactly you need? Matrix manipulations in linear algebra are really a huge topic, and we shouldn't squander resources on an implementation that is of no interest to anybody.

[JamesTev]

@v923z Awesome, thank you! I agree that it's important to focus resources on the most useful functionality. To this end, a generic eigenvalue solver would be useful for countless algorithms and applications. My need (to perform CCA for EEG analysis) is just one. Although dimensionality reduction type problems are often fine because covariance matrices are symmetric, I've encountered several problems that produce asymmetric ones.

[EmDash00]

I'm going to chime in here saying that an SVD implementation would let me implement a generic pseudoinverse solver for elliptical fitting. I'm a PhD student at UW and currently my "SVD implementation" uses a rudimentary implementation based on power iteration, but a more robust and efficient approach such as one using the QR algorithm would be helpful.

Seeing that the QR decomposition was implemented by #427, I don't believe it would be too bad to implement it now. Unfortunately numerical algorithms just aren't my area and I'm mainly interested in the applications. If it would be possible to implement generic SVD, that would be great.

I can certainly try to help if need be with the SVD algorithm if time allows.

[v923z]

There are a couple of features in the pipeline, but I might be able to work on this in the next one or two weeks. When this issue came up, I looked into possible implementations of the SVD algorithm, and it's definitely doable with moderate effort.
Once I have a reasonable prototype, I'd like to ask you to test it.

[kwagyeman]

@v923z - I'd like to support you to get SVD support added for ulab along with everything necessary to implement a homography.

E.g. https://medium.com/analytics-vidhya/using-homography-for-pose-estimation-in-opencv-a7215f260fdd

While I can easily implement a method to do this in OpenMV thanks to AprilTag's matrix library in our code base I think it should be in ulab instead of image library directly.

Please let me know how I can help.

...

All the matrix math you need: https://github.com/openmv/apriltag/blob/master/common/matd.c#L980

And how to do a homography: https://github.com/openmv/apriltag/blob/master/common/homography.c

...

We use the homography_to_pose method to produce what we need given the camera matrix.

[v923z]

@kwagyeman SVD in ulab would be a great addition, I completely agree. (You have to look no further than this feature request.) It seems to me that most of the work has already been done in matd.c, thanks for the pointer! I can add the required interface, that would be absolutely no problem.

As for the homography routine, do you want to add that ulab? If so, can we place it in the utils module, or something similar, given that it seems to have no corresponding numpy/scipy function?

Also, would you be interested in supporting blob searching and similar? There is some code for that e.g., here: https://github.com/teuler/make-thermal-cam/blob/main/blob_ulab2_c_code/user.c, and I've been flirting with the idea for a long time. Having some image processing functionality could be useful beyond openmv, and I would be happy to add those functions, but I just don't know if there is interest.

[kwagyeman]

I'm wouldn't add anything to ulab that needs performance. Pretty much anything for image processing needs to go into SIMD land if you really want to push the speed.

I'd keep it more high level and focused on doing mathematics for arrays.

Note, The AprilTag folks released an AprilTag 3 library that may have an updated version of that library with faster/better code.

As for the homography... I could add something to our image library to handle this. However, it may be more efficient to have in ulab since it needs the SVD code base. So, there would just be duplicates of the same code all over. Given that, adding the methods in their homography c file to a utils module that's generic for everyone to use would save code space I think in the long run.

If you could make it possible to call these methods from C still that would be great. Then I could add a method to our image lib if we needed to for some reason without having two copies of the same base library.

[v923z]

I'm wouldn't add anything to ulab that needs performance. Pretty much anything for image processing needs to go into SIMD land if you really want to push the speed.

I'd keep it more high level and focused on doing mathematics for arrays.

OK.

Note, The AprilTag folks released an AprilTag 3 library that may have an updated version of that library with faster/better code.

As far as I see, the last commit to the SVD routine was two years ago: https://github.com/AprilRobotics/apriltag/blob/master/common/svd22.c

As for the homography... I could add something to our image library to handle this. However, it may be more efficient to have in ulab since it needs the SVD code base. So, there would just be duplicates of the same code all over. Given that, adding the methods in their homography c file to a utils module that's generic for everyone to use would save code space I think in the long run.
If you could make it possible to call these methods from C still that would be great. Then I could add a method to our image lib if we needed to for some reason without having two copies of the same base library.

One the one hand, I guess it would be enough, if we made the functions non-static, and then you can call them from everywhere. But on the other hand, if you want to implement another interface, that would definitely lead to code duplication, because the argument parsing etc. would be at two places. I might be missing the point, though.