Added functions: schur, sylvester, and lyap

[egidioln]

Hi there,

I added these 3 functions: schur, sylvester, and lyap. Respectively, they perform a real Schur decomposition, solve a Sylvester equation and a continuous-time Lyapunov equation.

I notice that there are no javascript implementations of these online and they are pretty useful linear-algebra operations for signal processing and control systems.

In the future, these could be extended to deal with complex matrices and the discrete-time version of lyap (lyapd) could be implemented.

I followed the syntax from Matlab and the values from the tests were compared to the outputs from there.

Tests are passing locally. I hope all is good :)

[josdejong]

Thanks Lucas, that looks good at first sight. I will review your PR more in-depth soon.

Some first notes:

1. Can you write some unit tests to test whether the functions work correctly with Array and Matrix, with number and BigNumber? The current tests only test an Array with numbers.
2. I'm glad that you wrote the TypeScript definitions. In these defintions, there is also a chained version of every function inside interface MathJsChain<TValue> { ... }, can you add the method calls there too?
3. Maybe the category "Algebra" would be a better fit than "Matrix"?

[egidioln]

You're welcome and thanks for the review! :)

1. Can you write some unit tests to test whether the functions work correctly with Array and Matrix, with number and BigNumber? The current tests only test an Array with numbers.

Yes, I'll try to do it in the next few days (I am a bit overwhelmed this week)...

2. I'm glad that you wrote the TypeScript definitions. In these defintions, there is also a chained version of every function inside `interface MathJsChain<TValue> { ... }`, can you add the method calls there too?

Sure thing!

3. Maybe the category "Algebra" would be a better fit than "Matrix"?

Yes, I was hesitating between both but I agree that makes sense.

By the way, CI is failing but I cannot check if the error was due to my changes...

[josdejong]

Thanks, take it easy 👍

Don't worry about the failing build-and-test, I can look into it if needed. It could indeed be unrelated to your changes.

As for the types: so there are two definitions: a static function like math.add(a, b), and a chained method like math.chain(a).add(b).done().

[egidioln]

Hey @josdejong, I added the tests for Matrix and Array. I didn't do it for BigNumbers because I am using qr() and apparently it does not support it (am I wrong?).

I also added the docs for the chained versions (only lyap and schur because sylvester requires 3 inputs and would be a bit confusing to use it chained, but I can do it if you think it makes sense)

As you suggested I moved the functions to algebra and schur to the subfolder /decomposition

Finally, CI was failing because of the doc tests. I added schur and sylvester to knownProblems in doc.test.js.

Let me know if there's something else to be done :)

[josdejong]

Thanks for the updates Lucas 👍 , I'll review your updates somewhere this week

[egidioln]

Hello,

I just solved the conflict in the author's file.

[josdejong]

Thanks for the ping @egidioln , for some reason I totally forgot about this PR, I'm sorry 😳 .

You're right about the BigNumbers, that is not supported by qr unfortunately, so it's not possible to make this work for these three new functions right now.

All in all this PR looks very good. Thanks for writing the documentation and unit tests. Two more questions:

1. When I evaluate the matrices in the unit tests using Octave, I get differing results for schur and sylvester (lyap gives the same results though with inverted sign). How can I know that these outcomes are indeed valid? Are there multiple valid results possible?
2. Thanks for adding the type definitions, including the chained calls. I think you forgot to define the chained function call for sylvester? Can you also add a (small) usage example in types/index.ts to verify that the definitions are valid?

[egidioln]

No Worries @josdejong, thanks for doing this check :)

1. When I evaluate the matrices in the unit tests using Octave, I get differing results for `schur` and `sylvester` (`lyap` gives the same results though with inverted sign). How can I know that these outcomes are indeed valid? Are there multiple valid results possible?

Unfortunately the schur decomposition is not unique. One can verify the correctness of the decomposition by comparing $A$ (input) and $UTU^\top$ (outputs).

For sylvester and lyap, there are some mismatching definitions for these equations in the literature. For instance, the Sylvester equation is defined as $AX+XB=C$ in octave/matlab, as $AX-XB=C$ in Bartels-Stewart's paper (which I used to code the algorithm) and as $AX+sXB^\top=C$ in LAPACK, where s is an input 1 or -1.

For the Lyapunov equation in octave/matlab it is $AX+XA'=-Q$ but I actually coded it as $AX+XA'=Q$, which is a format I also see often and matches better the Sylverster equation.

However, thinking better now, maybe it is a good idea to have a notation for these that matches octave/matlab, do you agree?

2. Thanks for adding the type definitions, including the chained calls. I think you forgot to define the chained function call for `sylvester`? Can you also add a (small) usage example in types/index.ts to verify that the definitions are valid?

Actually, maybe I wasn't clear but I had mentioned above that I didn't implement it because sylvester requires 3 inputs $(A,B,C)$ and would be a bit confusing to use it in a chained call as I am not sure which of the arguments should be declared/taken from the chained operation (maybe $A$?).

[josdejong]

ah, that makes sense. Your explanation of the different definitions in differing math applications helps me a lot understanding what is going on. That complicates things.

So first thing is to clearly document this, which you did. It may be good to mention the existence of an alternative definition in the documentation so the user is aware that he may be assuming one or the other definition.

In general I try to use the "most commonly" used definition: what do Matlab, Mathematica, and Pythons NumPi for example use? Do you have an idea what is the most common definition for these three functions?

As for the missing chained definition for sylvester: chained usage like chain(A).sylvester(B, C) indeed looks a bit odd. It is fine with me not to define this.

[egidioln]

Hey @josdejong, Sorry, it took me so long to go through this...

So I changed sylvester's definition to match those of Matlab, Octave, Scipy and Julia. I also changed lyap's definition to match those of Matlab, Octave, and Julia (Scipy uses a different one, which I had implemented before, but I'd go with following the majority here, especially Matlab because I think they may have the largest community interested in using these).

I changed the docs and tests accordingly. Let me know if there's something else to be done here 😉

[josdejong]

Thanks @egidioln , this looks very good and well documented, thanks for all your work 👍

Looking at the code, I think the implementations also support bignumbers out of the box, which is great. Can you write some (small) unit tests to validate that?

[egidioln]

Thanks @egidioln , this looks very good and well documented, thanks for all your work +1

Looking at the code, I think the implementations also support bignumbers out of the box, which is great. Can you write some (small) unit tests to validate that?

Thanks @josdejong 😄! So is qr() working now for bignumbers? That was the reason why I didn't do the tests for them before

https://github.com/josdejong/mathjs/pulls#issuecomment-1274438181

[josdejong]

o wow, of course. Sorry I forgot about the dependency on qr. No there are no changes in that regard.

OK all is ready then, thanks a ton for all your effort Lucas!

[egidioln]

o wow, of course. Sorry I forgot about the dependency on qr. No there are no changes in that regard.

OK all is ready then, thanks a ton for all your effort Lucas!

Cool!! Thank you for your guidance Jos! I'll try to code the discrete-time version of lyap soon :)

[josdejong]

ow, that sounds promising! Thanks.

[josdejong]

Published now in v11.4.0