# Small linear matrix(5×5) equation taking too long that can only abort

The story revolves around something similar to "Jordan normal form":

$$P^{-1}AP=J.$$

First, to generate two 4×4 lower triangular matrices — $A,\ J$ (5×5 in the end):

n = 4;

rules = {{i_, i_} -> i - 1, {i_, j_} /; i >= j -> j};
lowT[n_, range__: All] := SparseArray[rules[[range]], {n, n}]

{A, J} = lowT[n, #] & /@ {All, 1} // Normal;
MatrixForm /@ {A, J}


$\{ {\left( {\matrix{0 & 0 & 0 & 0 \cr 1 & 1 & 0 & 0 \cr 1 & 2 & 2 & 0 \cr 1 & 2 & 3 & 3 \cr } } \right), \left( {\matrix{0 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 2 & 0 \cr 0 & 0 & 0 & 3 \cr } } \right)} \}$

Then, set the variable $P$, solve the equation on the top.

P = x~Array~{n, n};
sol = Solve[(P~MatrixPower~-1).A.P == J, Flatten@P];


The original result is rather long as algebraic form. For simplicity, arbitrarily substituting some numerical value, here, the elements on the diagonal are all set to 1:

(P = P /. Flatten@sol /. Table[x[i, i] -> 1, {i, n}]) // MatrixForm


$\left(\begin{array}{cccc}1 & 0 & 0 & 0 \\-1 & 1 & 0 & 0 \\\frac{1}{2} & -2 & 1 & 0 \\-\frac{1}{6} & 2 & -3 & 1 \\\end{array}\right)$

Finally, show the numerical equation and verify the correctness.

Inactivate[MatrixPower[P, -1].A.P == J, Dot] // TraditionalForm
% // Activate


$\left(\begin{array}{cccc}1 & 0 & 0 & 0 \\1 & 1 & 0 & 0 \\\frac{3}{2} & 2 & 1 & 0 \\\frac{8}{3} & 4 & 3 & 1 \\\end{array}\right).\left(\begin{array}{cccc}0 & 0 & 0 & 0 \\1 & 1 & 0 & 0 \\1 & 2 & 2 & 0 \\1 & 2 & 3 & 3 \\\end{array}\right).\left(\begin{array}{cccc}1 & 0 & 0 & 0 \\-1 & 1 & 0 & 0 \\ \frac{1}{2} & -2 & 1 & 0 \\-\frac{1}{6} & 2 & -3 & 1 \\\end{array}\right)=\left(\begin{array}{cccc}0 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 2 & 0 \\0 & 0 & 0 & 3 \\\end{array}\right)$

True

All code runs completed in about 2.5 seconds, but when it comes to 5×5(n=5), the equation will take too long that there is no hope for waiting it finish out.

So where is the problem?

Thanks for helping!

P.S. All code is merged as follows:

n = 4;
rules = {{i_, i_} -> i - 1, {i_, j_} /; i >= j -> j};
lowT[n_, range__: All] := SparseArray[rules[[range]], {n, n}]

{A, J} = lowT[n, #] & /@ {All, 1} // Normal;
MatrixForm /@ {A, J}

P = x~Array~{n, n};
sol = Solve[(P~MatrixPower~-1).A.P == J, Flatten@P];
(P = P /. Flatten@sol /. Table[x[i, i] -> 1, {i, n}]) // MatrixForm

Inactivate[MatrixPower[P, -1].A.P == J, Dot] // TraditionalForm
% // Activate

• FWIW: Solve[] is not really necessary here. Since you mention already knowing about Jordan form, you can use JordanDecomposition[] to show that lowT[n, All] is diagonalizable, and the first matrix returned by JordanDecomposition[lowT[n, All]] is the similarity transformation you are looking for. Aug 24, 2020 at 12:41

MatrixPower is the bottleneck, and it's unnecessary to call it:

Remove[A, i, j, J, lowT, n, P, range, rules, sol, x]
n = 6;

rules = {{i_, i_} -> i - 1, {i_, j_} /; i >= j -> j};
lowT[n_, range__: All] := SparseArray[rules[[range]], {n, n}]

{A, J} = lowT[n, #] & /@ {All, 1} // Normal;
MatrixForm /@ {A, J}

P = x~Array~{n, n};
sol = Solve[Thread[Flatten[A.P] == Flatten[P.J]], Drop[Flatten@P, {1, -1, n + 1}]];

(P = P /. Flatten@sol /. Table[x[i, i] -> 1, {i, n}]) // MatrixForm


$\left( \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & -2 & 1 & 0 & 0 & 0 \\ -\frac{1}{6} & 2 & -3 & 1 & 0 & 0 \\ \frac{1}{24} & -\frac{4}{3} & \frac{9}{2} & -4 & 1 & 0 \\ -\frac{1}{120} & \frac{2}{3} & -\frac{9}{2} & 8 & -5 & 1 \\ \end{array} \right)$