Jordan form of a matrix

Specific Domains
1

Is there an equivalent of Matlab’s jordan function, to compute the Jordan canonical form of a matrix? I didn’t find it by tab-completion and in the julia observer didn’t find anything relevant. A version over exact rings would be fine.

2

jordan_form is available through SymPy.jl, eg. in this example:

    M = Sym[ 6  5 -2 -3;
            -3 -1  3  3;
             2  1 -2 -3;
            -1  1  5  5]
    P, J = jordan_form(M)
    @test J == [2 1 0 0;
                0 2 0 0;
                0 0 2 1;
                0 0 0 2]
@test J == inv(P) * M * P
3

As noted in the Matlab manual, a floating-point jordan function is nearly useless, because the tiniest roundoff errors will generally break up Jordan blocks (by breaking eigenvalue degeneracies). The built-in linear algebra functions in Julia are mostly focused on floating-point computations. Computing the Jordan form over finite fields is quite a different problem, and something like SymPy may be more appropriate for that.

(In numerical computation, it is rare to want the Jordan form precisely because it is so unstable. Sometimes you can use the Schur factorization instead. In some applications, you have a matrix that is close to a defective matrix with a known Jordan form, e.g. close to a matrix with a 2x2 Jordan block, and in this case there are efficient algorithms to find the generalized eigenvectors of the “nearby” defective matrix.)

7 Likes
4

Thanks for the insights and the reference! I’ll have a look.

While reading your answer I remembered a passage on Moler and Van Loan’s Nineteen dubious ways to compute the matrix exponential. In fact the kind of problem that i’m interested in is very much the one depicted in the “Method 18. Block Diagonal” (page 25): the aim is to do some sort of “interpolation” between the Jordan form and the Schur form. We thought – naively – computing with Julia both the Jordan form and a Schur form, then taking a linear combination of these, which minimizes some objective norm. Anyways, in that section some strategies are outlined, but i haven’t thought too much (yet) if this is something that can be readily done with existing functions in Julia, or it would require a lot of new development.

5

You could try using the NORMFORM package with Reduce.jl, which is a julia pkg I made to interface the Reduce computer algebra system within Julia, the NORMFORM package can be called from it.

using Reduce
load_package(:normform)

If I have some time available, I might also create a built-in feature to interface the Jordan form into a julia function directly, there are still many opportunities for extended that pacakge.

6

Okay, so this feature does not quite work out of the box in Julia yet, since I have not implemented parsing of matrices through the Reduce pipe. However, I will consider this as my next task for moving forward.

julia> using Reduce
Reduce (Free CSL version, revision 4015), 05-May-17 ...

julia> load_package(:normform)

reduce> A:=mat((1,y),(y^2,3));

     [1   y]
     [     ]
a := [ 2   ]
     [y   3]


reduce> first jordansymbolic(A);

As of now, the matrix output is not properly detected and parsed yet, but this is a feature I will add in next.

Then you will be able to do Jordan form using this package as well (which is compatible with 0.7 unlike SymPy).

1 Like
7

If you want to study the Jordan structure of dense matrices numerically, you might find a generalized nullspace analysis helpful. To that end, I packaged up an implementation here.

Note that I said “study the Jordan structure” (i.e. characterize the invariant subspaces of nearby matrices) and not “perform a Jordan decomposition”, since (as @stevengj indicated above), the latter is usually ill-advised in numerics.

1 Like
8

Alright, so the symbolic matrix parsing is now implemented on the master branch,

julia> :([1 y; y^2 1]) |> RExpr

[1   y]
[     ]
[ 2   ]
[y   1]


julia> ans |> parse
:([1 y; y ^ 2 1])

So the parser now handles matrices in both directions (from Reduce into Julia and from Julia into Reduce)

Now it is only a matter of making the NORMFORM package play nicely with the Pipe

9

Thanks for the updates and your efforts to use Reduce’s functionality in Julia :slight_smile:
What is the Pipe?

10

“Generalized Null Space Decomposition” – Okay, looking at your description in the README I’m excited about trying this out in our research problem; basically we have some control systems (LTI) and perform a set-based analysis in the state-space through block decompositions.

If you need more tests, I can send you some matrices or write test scripts (PS: we’ve been working with the SLICOT benchmarks models).