Programming with Mathematica Syntax - Dunford decomposition

Question

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

Answer 1

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

About computational complexity

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

Answer 2

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]
    }
]

Then, for your example:

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

{"Diagonalizable" -> {{1, 0, 0}, {1, 0, -1}, {0, 0, 1}}, "Nilpotent" -> {{0, 0, 1}, {0, 0, 1}, {0, 0, 0}}}