# Implementing Brockett's Flow in Mathematica

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}

• It looks like A needs to be diagonal. Feb 27 at 9:20
• 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) Feb 27 at 17:01

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


• 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}` Feb 27 at 19:42