Schur decomposition: What is the Mathematica equivalent to MATLAB's rsf2csf function?

Question

MATLAB has a function rsf2csf that can convert real Schur form to complex Schur form.

The description of the function can be found here.

How can I do this in Mathematica, and is there a Mathematica function that can do the job?

Answer 1

I think I may have found a solution:

Mathematica's SchurDecomposition routine comes with an option RealBlockDiagonalForm which we can set to False to obtain the complex Schur form of the matrix.

Consider the matrix $$ X = \begin{bmatrix} 1 & 1 & 1 & 3 \\ 1 & 2 & 1 & 1 \\ -2 & 1 & 1 & 4 \end{bmatrix} $$

In Mathematica, running

X = {{1, 1, 1, 3}, {1, 2, 1, 1}, {1, 1, 3, 1}, {-2, 1, 1, 4}};
{U, T} = SchurDecomposition[N[X]];

gives the following two matrices

U = {{-0.491624, -0.489955, -0.633055, -0.342771}, {-0.497979, 
  0.2403, -0.232494, 0.800137}, {-0.675062, 0.428792, 
  0.423036, -0.425992}, {-0.233698, -0.719956, 0.605168, 0.246617}}

T {{4.81213, 1.19723, -2.22729, -1.00674}, {0., 
  1.92022, -3.04846, -1.83806}, {0., 0.712929, 1.92022, 
  0.256561}, {0., 0., 0., 1.34743}}

and on the other hand running

X = {{1, 1, 1, 3}, {1, 2, 1, 1}, {1, 1, 3, 1}, {-2, 1, 1, 4}};
{Unew, Tnew} = 
  SchurDecomposition[N[X], RealBlockDiagonalForm -> False];

gives us

Unew = {{-0.172637 + 0.460315 I, -0.443908 - 0.271038 I, 
  0.258268 + 0.550996 I, -0.342771 + 0. I}, {-0.174869 + 0.466266 I, 
  0.00633486 + 0.238756 I, -0.0875173 + 0.217011 I, 
  0.800137 + 0. I}, {-0.237052 + 0.632071 I, 
  0.337498 + 0.262733 I, -0.216582 - 0.364666 I, -0.425992 + 
   0. I}, {-0.0820646 + 0.218815 I, -0.054562 - 0.697517 I, 
  0.268798 - 0.568161 I, 0.246617 + 0. I}}

Tnew = {{4.81213 + 0. I, -1.4439 + 0.130783 I, -2.03655 + 
   0.380418 I, -0.353522 - 0.942624 I}, {0. + 0. I, 
  1.92022 + 1.47422 I, 
  2.16546 + 0.874913 I, -0.640886 + 1.52965 I}, {0. + 0. I, 0. + 0. I,
   1.92022 - 1.47422 I, 0.779147 + 0.294321 I}, {0. + 0. I, 0. + 0. I,
   0. + 0. I, 1.34743 + 0. I}}

Now if we compare it with Matlab, we can can see that the matrices $U_{new}$ and $T_{new}$ are different.

But it is important to note that in both cases, the diagonal of the matrix $T_{new}$ is the same, and this is crucial because it represents the eigenvalues of the original matrix $X$.

We can verify this easily in Mathematica, viz.,

In[]:= Diagonal[Tnew] // Chop

Out[]= {4.81213, 1.92022 + 1.47422 I, 1.92022 - 1.47422 I, 1.34743}

In[]:= Eigenvalues[X] // N

Out[]= {4.81213, 1.92022 + 1.47422 I, 1.92022 - 1.47422 I, 1.34743}