3
$\begingroup$

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

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

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

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

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

$\endgroup$
  • $\begingroup$ I am using effective algorithms that is I don't compute eigenvalues of the matrix A to find its decomposition. It is required because I would like to compare some methods. $\endgroup$ – Smilia Nov 2 '18 at 16:15
  • $\begingroup$ @Smilia I guess you meant to say "efficient" instead of "effective"? I would not say that your algorithm is very efficient. See also my remark on computational complexity under my post. $\endgroup$ – Henrik Schumacher Nov 2 '18 at 19:57
  • $\begingroup$ ok, I mean effective in the "mathematical" sense: en.wikipedia.org/wiki/Effective_method $\endgroup$ – Smilia Nov 2 '18 at 23:12

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.