# How to find a unitary transform between two sets of matrices?

In an earlier post, I asked a question about unitary transformation between two given matrices. Some really nice solutions were provided. However, what I'm really looking for is a unitary matrix U between two sets of matrices (that is, the same U is a unitary transform for all 5 corresponding pairs simultaneously). I verified that the ones provided do not achieve the goal for the following two sets,

G1={{{0, 0, 0, I}, {0, 0, I, 0}, {0, -I, 0, 0}, {-I, 0, 0, 0}}, {{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}, {{0, 0, I, 0}, {0, 0, 0, -I}, {-I, 0, 0, 0}, {0, I, 0, 0}}, {{0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}, {0, 1, 0, 0}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}};

G2={{{0, 0, 0, -I}, {0, 0, -I, 0}, {0, I, 0, 0}, {I, 0, 0, 0}}, {{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}, {{0, 0, -I, 0}, {0, 0, 0, I}, {I, 0, 0, 0}, {0, -I, 0, 0}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}, {{0, 0, -1, 0}, {0, 0, 0, -1}, {-1, 0, 0, 0}, {0, -1, 0, 0}}};


This problem arises from the need to switch between two gamma matrix conventions in physics. I know a solution exists up to a constant phase factor.

Just solve the equations $$g_i U=U g'_i,\qquad i=1,2,...,5$$ for $$g\in G_1$$ and $$g'\in G_2$$, and then impose unitarity on $$U$$:

Module[{G1, G2},
G1 = {{{0, 0, 0, I}, {0, 0, I, 0}, {0, -I, 0, 0}, {-I, 0, 0, 0}}, {{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}, {{0, 0, I, 0}, {0, 0, 0, -I}, {-I, 0, 0, 0}, {0, I, 0, 0}}, {{0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}, {0, 1, 0, 0}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}};
G2 = {{{0, 0, 0, -I}, {0, 0, -I, 0}, {0, I, 0, 0}, {I, 0, 0, 0}}, {{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}, {{0, 0, -I, 0}, {0, 0, 0, I}, {I, 0, 0, 0}, {0, -I, 0, 0}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}, {{0, 0, -1, 0}, {0, 0, 0, -1}, {-1, 0, 0, 0}, {0, -1, 0, 0}}};
Array[u, {4, 4}] /. (Table[G1[[i]].Array[u, {4, 4}] == Array[u, {4, 4}].G2[[i]], {i, 1, 5}] // Solve)[[1]]
]
% /. FindInstance[%.ConjugateTranspose[%] == IdentityMatrix[4], Variables[%]][[1]] // MatrixForm


• Wow, this is a beautiful solution. Thanks. I didn't know the power of FindInstance. Accepted. Aug 1, 2021 at 16:49
• @Frank, just a general point of order that we typically like to wait ~24 hours to mark an answer as accepted to allow the possibility for those in other time zones to submit their own solutions which may have their own unique elegance to them. Aug 1, 2021 at 18:17
• @AccidentalFourierTransform, actually, it has a small bug. The G1 and G2 should be switched in the code. It should be G2.U=U.G1. Aug 1, 2021 at 18:43
• @Frank I don't know what you mean by your last comment. Unitary equivalence is $Ua=bU$ or $aU=Ub$, both are quite literally the exact same condition. An equivalence relation is by its very definition symmetric. Aug 2, 2021 at 8:04
• @AccidentalFourierTransform, you formulated it as G1=U.G2.U^dagger. I wanted G2=U.G1.U^dagger for the given G1 and G2, which gives a different solution. Aug 2, 2021 at 12:34

First you make a matrix U containing all the variables.

U = Array[
ToExpression[StringJoin["u", Map[ToString, List[##]]]] &, {4, 4}];


This gives us the matrix:

{{u11, u12, u13, u14},
{u21, u22, u23, u24},
{u31, u32, u33, u34},
{u41, u42, u43, u44}}


where u11 ... u44 are simply free variables (they will be our 16 unknowns).

Then you need to express all the equations following from your constraints on U, G1 and G2. For example $$U$$ is a unitary matrix. That means

U.Transpose[U]- IdentityMatrix[4]


gives 0. But this is actually a conjunction of 16 distinct equations (one for each index of the resulting matrix). We write:

eq1 =  Apply[And,
Map[# == 0 &, Flatten[U.Transpose[U] - IdentityMatrix[4]]]];


This is our first constraint. U is unitary. Our second one relates to $$G1$$ and $$G2$$:

eq2 = Apply[And,
Map[# == 0 &, Flatten[Map[U.#.Transpose[U] &, G1] - G2]]];


Now we solve for our unknowns:

sol = Solve[eq1 && eq2, Variables[U]];


This gives us the solutions as a list of possible values for our unknowns:

{{u11 -> 0, u12 -> -(1/Sqrt[2]), u13 -> 0, u14 -> -(1/Sqrt[2]), ... }}


But we're not interested in those values per say, only in how our final matrix might look. With:

values = Map[U /. # &, sol]


we get a list of two solutions for $$U$$, which are:

{{{0, -(1/Sqrt[2]), 0, -(1/Sqrt[2])},
{1/Sqrt[2], 0, 1/Sqrt[2], 0},
{0, 1/Sqrt[2], 0, -(1/Sqrt[2])},
{-(1/Sqrt[2]), 0, 1/Sqrt[2],  0}},

{{0, 1/Sqrt[2], 0, 1/Sqrt[2]},
{-(1/Sqrt[2]), 0,-(1/Sqrt[2]),  0},
{0, -(1/Sqrt[2]), 0, 1/Sqrt[2]},
{1/Sqrt[2], 0, -(1/Sqrt[2]), 0}}}


Tough this only gives the real solutions. Explicitly computing all solutions for the complex domain might be too complicated for Mathematica to do directly.

• This works! Appreciate the step-by-step explanations. The two solutions differ by an overall sign, which is ok. Aug 1, 2021 at 19:44