Sorted Schur form

[atisharma]

As I understand it, es! gives the Schur decomposition of a matrix by calling MKL gees. Unfortunately there is no way to get sorted real Schur form, which is sometimes important (e.g. subspace calculations). MKL has a 'sort' option which achieves this. Is it possible for es! to allow the call to specify sorted form?

This is in relation to a Riccati equation solver I'm writing for a maths library project

[blueberry]

Hi Ati,

Currently, es! does not offer the sort option, but that can be changed in one of the future releases of Neanderthal. What I'd need is a clear and explicit test case that demonstrates how it would be used (with input numbers and the expected results).

[atisharma]

That would be great, thank you.

It's a bit tricky because Schur form is not unique. I think it can be done as follows, though.

(def M (dge 3 3 [8 3 4 1 5 9 6 7 2]))

#RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →       8.00    1.00    6.00
   →       3.00    5.00    7.00
   →       4.00    9.00    2.00
   ┗                               ┛

(def T (copy M))

gives real Schur form (unordered) of

(uncomplicate.neanderthal.linalg/es! T w U)

#RealGEMatrix[double, mxn:3x2, layout:column, offset:0]
   ▥       ↓       ↓       ┓
   →      15.00    0.00
   →       4.90    0.00
   →      -4.90    0.00
   ┗                       ┛

matlib.core=> T
#RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →      15.00    0.00   -0.00
   →       0.00    4.90   -3.46
   →       0.00    0.00   -4.90
   ┗                               ┛

matlib.core=> U
#RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →      -0.58   -0.81   -0.07
   →      -0.58    0.47   -0.67
   →      -0.58    0.34    0.74
   ┗                               ┛

My (loose) understanding of MKL based on the documentation, is that you have to specify which eigenvalue or eigenvalues are to appear in the top left. Let's say we specify only the one that is close to -4.8990, in which case we have the ordered form (according to Matlab):

U-ordered:
   -0.3416   -0.5774   -0.7416
   -0.4714   -0.5774    0.6667
    0.8131   -0.5774    0.0749

T-ordered:
   -4.8990   -0.0000    3.4641
         0   15.0000    0.0000
         0         0    4.8990

Since the other two eigenvalues could appear in any order, I think the test should verify (1) that (mm U T (trans U)) reproduces M and that (2) the first element of T-ordered is the eigenvalue we specify. Then (nrm1 (axpy -1 M (mm U T (trans U)))) should be less that some tolerance (I'm not yet sure what you use). We also have to check that (mm U (trans U)) is the identity matrix.

Is that approach suitable for you? If so I can have a go at drafting a test.

[blueberry]

Q:

is that you have to specify which eigenvalue or eigenvalues are to appear in the top left. Let's say we specify only the one that is close to -4.8990,

How do we specify that in code?

Is that approach suitable for you? If so I can have a go at drafting a test.

Yes, please. However, please mind that I have a pile of work already in the queue, so realistically, it will be months rather than days until I can work on this, since it also requires that I extend the API in a non-intrusive way.

[atisharma]

How do we specify that in code?

Perhaps specify the last element of the Schur form (the last eigenvalue) to be the first element. This does require doing the unordered Schur form first.

I'll work on it, it should clarify things.

[atisharma]

I drafted a test in this gist.
It's a simple modification of your existing es! test. Please let me know if it's suitable.

PS - I fixed the link in the original post.

[blueberry]

Thank you. I hope so.

…

On Fri, Jul 3, 2020 at 13:39 Ati Sharma ***@***.***> wrote: I drafted a test in this gist <https://gist.github.com/atisharma/5b36df662d9897557748c59d9b0f6c3f>. It's a simple modification of your existing es! test. Please let me know if it's suitable. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#93 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAZ2JEORTGUGXSGVT66ECTRZW7PPANCNFSM4OPA2QCQ> .