Find a matrix $X$ that block-diagonalizes a particular group of matrices

Essentially, I want to find a single matrix $$X$$ such that conjugation by $$X$$ sends:

$$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \mapsto \begin{bmatrix} -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix}$$ and

$$\frac{1}{2} \begin{bmatrix} -i & -i & 0 & i & i \\ \sqrt{\frac{3}{2}} & 0 & -1 & 0 & \sqrt{\frac{3}{2}} \\ i & -i & 0 & i & -i \\ -\frac{1}{2} & 1 & -\sqrt{\frac{3}{2}} & 1 & -\frac{1}{2} \\ \end{bmatrix} \mapsto \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & e^{\frac{2 i \pi }{3}} & 0 \\ 0 & 0 & 0 & 0 & e^{-\frac{1}{3} (2 i \pi )} \\ \end{bmatrix}$$

The matrices are unitary and I have a mathematical guarantee this is possible.

Some math details: The two matrices on the left generate the tetrahedral group, and the two matrices on the right also generate the tetrahedral group. The matrices on the left are in the Spin(2) representation of the rotation group (rotation by 180 degrees about the z axis and rotation by 120 degrees about the $$(1,1,1)/\sqrt{3}$$ axis). Through some fancy representation theory (looking at the traces of the matrices on the left), you can show that you have to be able to decompose it into the direct sum of irreducible representations, and those give the matrices on the right.

I was able to get close to solving the problem. If I simply worry about the more complicated 2nd matrix, I can use the eigensystem to ensure the 2nd map is satisfied. The first map is then a mess with the top-left 3x3 block not how I want it. Also, this generalizes poorly to higher dimensions (2 sets of (2n+1)x(2n+1) matrices that I'm guaranteed are similar to other (2n+1)x(2n+1) matrices).

Is there any magic command or simple solution to this problem? I can probably do some sort of iterative procedure that keeps detailed track of all the eigenvectors/values, first finding the correct X for the first matrix, and then the correct X for the second matrix compatible with the first one. But that seems like it would be very involved!

This question is similar to Is there a clean way to extract the subspaces invariant under a list of matrices?, but I think my problem is much simpler because I know the invariant subspaces and in detail what I want the output to look like.

{s1,s2}={{{1,0,0,0,0},{0,-1,0,0,0},{0,0,1,0,0},{0,0,0,-1,0},{0,0,0,0,1}},{{-(1/4),-(1/2),-(Sqrt[(3/2)]/2),-(1/2),-(1/4)},{-(I/2),-(I/2),0,I/2,I/2},{Sqrt[3/2]/2,0,-(1/2),0,Sqrt[3/2]/2},{I/2,-(I/2),0,I/2,-(I/2)},{-(1/4),1/2,-(Sqrt[(3/2)]/2),1/2,-(1/4)}}};
{t1,t2}={{{-1,0,0,0,0},{0,-1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}},{{0,1,0,0,0},{0,0,1,0,0},{1,0,0,0,0},{0,0,0,E^(2 Pi I/3),0},{0,0,0,0,E^(4 Pi I/3)}}};
ct=ConjugateTranspose;
(* "Arg", "SortBy" and "Transpose" are used to make sure that the eigenvalues and eigenvectors of s2 and t2 are in the same order. *)
{t2d,t2x}=N@Transpose[SortBy[Transpose[Eigensystem[t2]],Arg[N[#[]]]&]];
(* the spectrum is degenerate, so we have to manually orthonormalize. This also ensures t2x and s2x are both unitary matrices *)
t2x=Transpose[Chop[ Orthogonalize[t2x]]];
{s2d,s2x}=N@Transpose[SortBy[Transpose[Eigensystem[s2]],Arg[N[#[]]]&]];
s2x=Transpose[Chop[ Orthogonalize[s2x]]];
X=ct[s2x.ct[t2x]];
MatrixForm[Chop[X.s1.ct[X]]]
MatrixForm[Chop[X.s2.ct[X]]]

Maybe newtons method:

X = ToExpression[Table["r" <> ToString[i] <> ToString[j], {i, 5}, {j, 5}]] +
I ToExpression[Table["i" <> ToString[i] <> ToString[j], {i, 5}, {j, 5}]];

eqs = DeleteCases[Flatten[{
ComplexExpand[Re[X.s1.ct[X]]] - Re[t1],
ComplexExpand[Im[X.s1.ct[X]]] - Im[t1],
ComplexExpand[Re[X.s2.ct[X]]] - Re[t2],
ComplexExpand[Im[X.s2.ct[X]]] - Im[t2]}], 0];
vars = Flatten[ComplexExpand[ReIm[X]]];
cur = RandomReal[{-1, 1}, Length[vars]];

jac = D[eqs, {vars, 1}];
Do[cur -= PseudoInverse[jac /. #].(eqs /. #) &[Thread[vars -> cur]], {20}]

X = Chop[X /. Thread[vars -> cur]]

It gives a different solution everytime because the starting point is chosen randomly and the problem is underdetermined.

By adding three zero equalities:

eqs = DeleteCases[Flatten[{i12, r41, i51,
ComplexExpand[Re[X.s1.ct[X]]] - Re[t1],
ComplexExpand[Im[X.s1.ct[X]]] - Im[t1],
ComplexExpand[Re[X.s2.ct[X]]] - Re[t2],
ComplexExpand[Im[X.s2.ct[X]]] - Im[t2]}], 0];

The system of equations remain solvable and the numerical solution happens to coincide with the particularly simple matrix:

$$\left( \begin{array}{ccccc} 0 & -1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 & 0 \\ 1 & 0 & 0 & 0 & -1 \\ 1 & 0 & -1 & 0 & 1 \\ 1 & 0 & -1 & 0 & 1 \\ \end{array} \right)\sqrt{\left( \begin{array}{ccccc} 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 \\ \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\ -\frac{1}{4} & 0 & \frac{1}{2} & 0 & -\frac{1}{4} \\ \frac{1}{4} & 0 & -\frac{1}{2} & 0 & \frac{1}{4} \\ \end{array} \right)}$$

where the multiplication and the square root are entry-wise.

• Fantastic!! I think the matrix is unique up to phases (you can left-multiply X by any diagonal matrix with entries Exp[-I*{d1,d1,d1,d2,d3}] and still get another unitary that solves the problem). There must be some simple structure I've missed; each row represents a different superposition of states, so you can look at the first row and go "spin 1 plus spin -1", or the fifth and see "spin 2 minus spin 0 plus spin -2". A simple way to plot: Array[ ParametricPlot3D[(Abs[X[[#]].Table[SphericalHarmonicY[2,m,a,b],{m,2,-2,-1}]] ){Cos[b]Sin[a],Sin[b]Sin[a],Cos[a]},{a,0,Pi},{b,0,2 Pi}]&,5] – David Nov 28 '18 at 7:42