3
$\begingroup$

Motivated by this post and this paper is an attempt to implement this simple relation in Mathematica to diagonalize a $(2,2)$ matrix. However, the result is not diagonal after iterating discretized version 100k times, can anyone see the correct way to implement this?

brockett = 
  Compile[{{X0, _Real, 2}, {A, _Real, 
     2}, {lr, _Real}, {iters, _Integer}}, 
   Module[{i, comm, X}, X = X0;
    For[i = 1, i <= iters, i += 1, comm = X . A - A . X;
     X += lr  (X . comm - comm . X);];
    X]];

X0 = {{3, 1}, {1, 2}} // N;
A = {{5, 2}, {2, 7}} // N;
{"A", A // MatrixForm}
{"X0", X0 // MatrixForm}
{"Xt", brockett[X0, A, .0001, 100000] // MatrixForm}
$\endgroup$
2
  • 2
    $\begingroup$ It looks like A needs to be diagonal. $\endgroup$ Feb 27 at 9:20
  • 1
    $\begingroup$ Yes, as @DanielHuber notes, this will work fine with A = N@{{5, 0}, {0 7}} (for example). Also, for nomenclature, clearly one should have brockettBracket[m1_, m2_] := m1 . (m1 . m2 - m2 . m1) - (m1 . m2 - m2 . m1) . m1 (One should always have handy a pocket-packet of brockettBrackets) $\endgroup$ Feb 27 at 17:01

1 Answer 1

3
$\begingroup$

Looks like A needs to be diagonal.

Also, couldn't figure out why I couldn't reuse definition of commutator inside the compiled function until reading this post, I guess I need to compile expressions first in order to inline them.

ClearAll["Global`*"];
SeedRandom[1];

(* commutator *)
comm = Compile[{{X, _Real, 2}, {Y, _Real, 2}}, X . Y - Y . X];
(* symmetrizer *)
symm[X0_] := (X0 + X0\[Transpose])/2;

d = 10;
X0 = symm[RandomVariate[NormalDistribution[], {d, d}]/Sqrt[d]];
A = DiagonalMatrix[RandomVariate[NormalDistribution[], {d}]]/Sqrt[d];

brockett = 
  Compile[{{X0, _Real, 2}, {A, _Real, 
     2}, {lr, _Real}, {iters, _Integer}},
   Module[{i, X},
    X = X0;
    For[i = 1, i <= iters, i += 1,
     X += lr   comm[X, comm[X, A]];
     ];
    X
    ],
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

relError[a_, b_] := Norm[a - b, 1]/Min[Norm[a, 1], Norm[b, 1]];
Xt = brockett[X0, A, .2, 50000];
Print["eig approx error: ", 
  relError[Sort@Diagonal[Xt], Sort@Eigenvalues[X0]]/d];
Print["Xt offdiagonal: ", Total[Abs@UpperTriangularize[Xt, 1], 2]];
ListLinePlot[{Sort@Eigenvalues@X0, Sort@Diagonal@Xt}, 
 PlotLegends -> {"true eigs", "Xt diagonal"}]
ArrayPlot@Xt

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Could also use NDSolve, like so. In[350]:= x0 = {{3, 1}, {1, 2}} // N; amat = {{5, 0}, {0, 7}} // N; xmat = Map[#[t] &, Array[x, {2, 2}], {2}]; odes = Thread[ Flatten[D[xmat, t] - brockettBracket[xmat, amat]] == 0]; inits = Thread[Flatten[xmat /. t -> 0] == Flatten[x0]]; solns = NDSolveValue[Join[odes, inits], Flatten[xmat], {t, 0, 100}]; Chop[solns /. t -> 100] Out[356]= {1.38197, 0, 0, 3.61803} $\endgroup$ Feb 27 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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