# Programming with Mathematica Syntax - Dunford decomposition

I am always amazed while I am reading some solutions of problems in mathematical problems of this site with Mathematica.

For instance, look at the problem 48.

Whereas I am writing 100 lines of code to solve this problem, a proposed 1-line solution is:

romDigits@Take[IntegerDigits[Plus@@Table[i^i,{i,1000}]],-10]


In this post, I would like to "factor" my code using the powerful syntax of Mathematica. Below you will find my code to decompose a complex matrix A as a sum D+N where $$D$$ is a diagonalisable matrix and $$N$$ is a nilpotent matrix (called Newton-Raphson method for Dunford decomposition of a matrix).

ClearAll;

Unprotect[MatrixPower]
MatrixPower[m_?SquareMatrixQ, 0] := IdentityMatrix[Length[m]]
Protect[MatrixPower]

DunfordEffective[A_] := Module[{Chi,P,Q,An,Anp1,coeffsP,coeffsQ,PAn,QAn,d,n},
Print[A // MatrixForm];
Chi=CharacteristicPolynomial[A,X];
P=Simplify[Chi/PolynomialGCD[Chi,D[Chi,X]]];
coeffsP=CoefficientList[P,X];
(* Derivative polynomial of P *)
Q=D[P,X];
coeffsQ=CoefficientList[Q,X];
An=A;
PAn=Sum[MatrixPower[An, j - 1] coeffsP[[j]], {j, 1, Length[coeffsP]}];
QAn=Sum[MatrixPower[An, j - 1] coeffsQ[[j]], {j, 1, Length[coeffsQ]}];
Anp1=An-PAn.Inverse[QAn];
While[An != Anp1,
An=Anp1;
PAn=Sum[MatrixPower[An, j - 1] coeffsP[[j]], {j, 1, Length[coeffsP]}];
QAn=Sum[MatrixPower[An, j - 1] coeffsQ[[j]], {j, 1, Length[coeffsQ]}];
Anp1=An-PAn.Inverse[QAn];
];
d=An;
n=A-An;
{d,n,d+n-A} // MatrixForm
]

A={{1,0,1},{1,0,0},{0,0,1}};
DunfordEffective[A]


Could you provide me some critics about this code ? How would it be nicer to read at ? Please, don't hesitate to provide me any critics, I want to improve my programming skills. The aim is not to treat big sized matrices but some programming remarks in Mathematica which could used syntax proper to Mathematica.

Thanks

• Notice a single ClearAll doesn't make sense. Please press F1 and check the document for its correct usage. – xzczd Nov 2 '18 at 16:09
• If the matrices involve approximate numbers, one can use SchurDecomposition for this. – Daniel Lichtblau Nov 2 '18 at 16:23

Here is a slightly modified version of your code.

Clear[DunfordEffective2]

DunfordEffective2[A_, TOL_: 1. 10^-8] :=
Module[{Chi, P, Q, Anp1, coeffsP, coeffsQ, An, PAn, QAn, X, powers, id, iter},
Chi = CharacteristicPolynomial[A, X];
P = Simplify[Chi/PolynomialGCD[Chi, D[Chi, X]]];
Q = D[P, X];
coeffsP = CoefficientList[P, X];
coeffsQ = CoefficientList[Q, X];
id = IdentityMatrix[Length[A]];
An = A;

powers = NestList[An.# &, id, Length[coeffsP] - 1];
PAn = coeffsP.powers;
QAn = coeffsQ.powers[[1 ;; Length[coeffsQ]]];
Anp1 = An - Transpose[LinearSolve[QAn][Transpose[PAn], "T"]];
iter = 1;
While[Max[Abs[An - Anp1]] > TOL,
iter++;
An = Anp1;
powers = NestList[An.# &, id, Length[coeffsP] - 1];
PAn = coeffsP.powers;
QAn = coeffsQ.Most[powers];
Anp1 = An - Transpose[LinearSolve[QAn][Transpose[PAn], "T"]];
];
<|"Diagonalizable" -> An, "Nilpotent" -> A - An, "Iterations" -> iter|>
]


Essential changes are:

• Getting rid of MatrixForm: really, it is only meant for displaying purposes.

• Usage of a tolerance in the while loop in case matrices of inexact precision are treated.

• Using LinearSolve instead of inverse; it is a bit faster for larger matrices and should be more accurate.

• Using NestList to compute matrix powers.

• Replacing the sums by coeffsP.powers and coeffsQ.Most[powers]: This way, the powers have to be computed only once per iteration.

• Using a somewhat more self-explaining way for the output.

For binary matrices, this is already much faster:

n = 30;
SeedRandom[1];
A = RandomInteger[{0, 1}, {n, n}];
a = DunfordEffective[A]; // AbsoluteTiming // First
b = DunfordEffective2[A]; // AbsoluteTiming // First
a[[1]] == b[["Diagonalizable"]]


0.736236

0.252682

True

But it is way faster for matrices in machine precision:

n = 30;
A = RandomReal[{-1, 1}, {n, n}];
b = DunfordEffective2[A]; // AbsoluteTiming // First


0.01196

But in cannot beat Carl Woll's suggestion (of course, JordanDecomposition is highly optimized):

c = "Diagonalizable" /. dunfordDecomposition[A]; // AbsoluteTiming
Max[Abs[b[["Diagonalizable"]] - c]]


0.001261

6.83786*10^-15

Moreover, the algorithm here is not very efficient. The computational complexity of finding the eigenvalues should be proportional to the complexity of matrix-matrix multiplication, where $$n$$ is the number of rows of $$A$$. With typical hardware and typical algorithms, this is $$O(n^3)$$. Due to this article by Beelen and Van Dooren, also the Jordan decomposition itself has computational complexity at most $$O(n^3)$$ (if an appropriate algorithm is used). This algorithm here needs at least $$n$$ matrix-matrix multiplications, so its complexity should be at least $$O(n^4)$$.

• That's great, thanks ! – Smilia Nov 2 '18 at 23:13

Why not just use JordanDecomposition?

dunfordDecomposition[m_] := Module[{s, j, id = IdentityMatrix @ Length @ m},
{s, j} = JordanDecomposition[m];
{
"Diagonalizable" -> s . (j id) . Inverse[s],
"Nilpotent" -> s . (j (1-id)) . Inverse[s]
}
]


dunfordDecomposition[{{1, 0, 1}, {1, 0, 0}, {0, 0, 1}}]