Interpreting error from computing Jordan form of 36-by-36 matrix

Question

I've been trying to compute the jordan normal form of a 36-by-36 matrix composed of only three distinct entries, 1, 1/2, and 0. The matrix is a probability transition matrix so, given these entries, the matrix is obviously sparse.

The issue I've been having is the following: whenever I try to compute

[V, J] = jordan(A),

or

[V, J] = jordan(sym(A)),

I get the following error message:

Error using mupadmex
Error in MuPAD command: Similarity matrix too large.

Error in sym/mupadmexnout (line 1546)
        out = mupadmex(fcn,args{:});

Error in sym/jordan (line 32)
        [Vsym,Jsym] = mupadmexnout('symobj::jordan',A,'All');

I've read in the MATLAB help that computation of the Jordan form is very sensitive to perturbations. I did not think my computation would be an issue, however, since all the entries of the matrix are either integers or ratios of integers.

My questions are the following:

1. How do I interpret the error output I received?
2. Are the errors I received addressable?
3. If the errors are not addressable, are there alternative methods (functions in Matlab) I could try to compute the Jordan form?

Ewald

Answer 1

1) How do I interpret the error output I received?

The point is that Matlab uses symbolic computation to evaluate the Jordan form. This is the reason that it asks you to provide rational numbers. A 36-by-36 matrix is very small when we are considering numerical programming, but (I am not sure about this) maybe this size is big for symbolic programming.

2) Why does not matlab have a toolbox to evaluate numerically the Jordan form?

The point is that this evaluation is numerically unstable. See the example in Wikipedia. Basically, any perturbation of a matrix with multiple eigenvalues (that share the same block) can cause these eigenvalues to become distinct in separated blocks of the desired Jordan Form.

3) If the errors are not addressable, are there alternative methods (functions in Matlab) I could try to compute the Jordan form?

I think Matlab does not have numerical functions to solve this task.

I dont know exactly what kind of application you are looking at... Having said that, one (very common) option is to evaluate the Schur form (both transformations transform the matrix in an upper triangular decomposition), which is numerically stable. It uses a unitary similarity transformation. Matlab's schur function implements this.

See also this Math.StackExchange question: What's the difference between Jordan and Schur decomposition?

Unify

Good answer. +1. Hope that you don't mind that I cleaned up some. – Clockwork

I think the OP updated the question so your #2 may be out of date. – Clockwork

The Jordan form is so numerically unstable that it has been proposed as the basis for uncheatable numeric benchmarks (PDF). – Clockwork

@Clockwork interesting I didn't know about that! – Unify

Answer 2

I will primarily address the third part of your question: alternative methods to compute the Jordan form.

Depending on the largest matrix you want to evaluate and possibly which Matlab version you have, yes, you can compute the Jordan form and it's similarity transformation symbolically.

I had a related question a while back so your question made me look into this again. In my question, I was not able to use the jordan function to calculate generalized eigenvectors for test matrices larger than 82-by-82 using Matlab R2013a (same in R2014b that I use now).

UPDATE: R2015a+ appears to have fixed the issue for matrices larger than 82-by-82 – at least for the 300-by-300 test matrix below. The following is a workaround for pre-R2015a.

---

The trick boils down to using a different (and probably newer) MuPAD function to calculate the Jordan form and the associated similarity transformation: linalg::jordanForm instead of symobj::jordan. Because you can't even compute a 36-by-36 Jordan form I wonder if you're using a much older version. So I can't guarantee that this will work in your version of Matlab. Here's the necessary code:

% My test matrix
A = -eye(40);
A(1,2:end-1) = -2;
A(1,end) = -1/2;
A(2,2) = 1/2;

% [V,J] = jordan(A); doesn't work for A larger than 82-by-82
P = feval(symengine,'linalg::jordanForm',sym(A),'All');
J = double(P(1));
V = double(P(2));

I was able to calculate J and V for a 300-by-300 test matrix using this, though it did take about a minute on my MacBook Pro. This is still a symbolic approach and anything over 100-by-100 can probably be considered quite large for many symbolic operations. Older versions of the Symbolic Math toolbox are even less efficient.

By the way, if you don't actually need the similarity transform, V, then you may still be able to use J = jordan(A); without getting an error.

Clockwork