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

Question

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

Answer 1

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)$