1
$\begingroup$

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
$\endgroup$
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.