# Find solution for specific similarity transformation

I have a situation where I need to find a similarity transformation of a specific form. I have some Hermitian matrix $$X$$, and I need to find a matrix $$T$$ that satisfies $$X = T J T^\dagger$$ where,

$$J = \text{diag}(1, -1; \dots; 1, -1).$$

I've tried using solve. For example:

Module[{X, j, T, n = 2},
X = {{-2, 0}, {0, 2}};
j = DiagonalMatrix[{1, -1}];
T = Array[t, {n, n}];
Solve[Simplify[
X - T.j.T\[ConjugateTranspose] == ConstantArray[0, {n, n}] //
ComplexExpand], Flatten[T]]]


In this case, a solution would be, $$T = \sqrt{2} \begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}$$

However I obtain a warning/error and lots of solutions that depend on other solutions

Solve::svars: Equations may not give solutions for all "solve" variables.


and it gives me a list of solutions like this:

{{t[1, 1] -> -Sqrt[-2 + t[1, 2]^2], t[2, 1] -> -t[1, 2],
t[2, 2] -> Sqrt[-2 + t[1, 2]^2]}, {t[1, 1] -> -Sqrt[-2 + t[1, 2]^2],
t[2, 1] -> t[1, 2],
t[2, 2] -> -Sqrt[-2 + t[1, 2]^2]}, {t[1, 1] -> Sqrt[-2 + t[1, 2]^2],
t[2, 1] -> -t[1, 2],
t[2, 2] -> -Sqrt[-2 + t[1, 2]^2]}, {t[1, 1] -> Sqrt[-2 + t[1, 2]^2],
t[2, 1] -> t[1, 2],
t[2, 2] -> Sqrt[-2 + t[1, 2]^2]}, {t[1, 1] -> -I Sqrt[2],
t[1, 2] -> 0, t[2, 1] -> 0,
t[2, 2] -> -I Sqrt[2]}, {t[1, 1] -> -I Sqrt[2], t[1, 2] -> 0,
t[2, 1] -> 0, t[2, 2] -> I Sqrt[2]}, {t[1, 1] -> I Sqrt[2],
t[1, 2] -> 0, t[2, 1] -> 0,
t[2, 2] -> -I Sqrt[2]}, {t[1, 1] -> I Sqrt[2], t[1, 2] -> 0,
t[2, 1] -> 0, t[2, 2] -> I Sqrt[2]}}


Many of which don't seem to even be a solution (for example the last one is wrong by a minus sign). What is the best way that I can compute a single solution to this?

Thanks

EDIT:

Usually my matrix X contains symbols

• (1) Is T restricted to be unitary? If so, this can be taken to a linear algebra problem by multiplying both sides on the right by T. (2) ComplexExpand will assume all variables (in particular, the t[i,j]) are real-valued unless told otherwise. Commented Jun 26, 2020 at 16:01
• @DanielLichtblau Thanks for the tip on ComplexExpand, and no T is not unitary (I think in that case J would be the diagonal matrix of eigenvalues rather than the form I specify) Commented Jun 28, 2020 at 8:15
• Ooof, that's right. Let me retry that: are the rows of T orthogonal to one another? Commented Jun 28, 2020 at 14:55

Here is a method that usually works, and, when it does, delivers a "nicer" result than the Eigensystem-based method in the case of integer or rational matrices.

The idea is to use a "symmetrized" LU decomposition, but in a way that does not give permutations. This requires knowledge of the pivot strategy (lucky for me, I wrote it). This strategy chooses the smallest nonzero pivot, so to reduce the likelihood of permuting we weight successive rows/columns by increasing powers of 100.

So the example goes like so.

SeedRandom[1234]
dim = 4;
mat = RandomInteger[{-5, 5}, {dim, dim}];
symmat = mat + Transpose[mat];

mult = DiagonalMatrix[Table[100^j, {j, 0, dim - 1}]];
invmult = Inverse[mult];
h2 = mult.symmat.mult;
{lu, perm, cnum} = LUDecomposition[h2];
diag = DiagonalMatrix[Diagonal[lu]];
lower = LowerTriangularize[lu] - diag + IdentityMatrix[dim];


Check that no permuting took place.

In[1506]:= perm

(* Out[1506]= {1, 2, 3, 4} *)


Now form the signed diagonal and the transformation matrix. We can use the "lower" part of the LU decomposition since the input is symmetric and no permuting was done.

newdiag = Sign[diag]
sqrroots = Sqrt[Abs[diag]];
tmatrix = invmult.lower.sqrroots

(* Out[1511]= {{-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}

Out[1513]= {{Sqrt[10], 0, 0, 0}, {-3 Sqrt[2/5], 4 Sqrt[2/5], 0,
0}, {1/Sqrt[10], -(11/(2 Sqrt[10])), 3/(2 Sqrt[2]),
0}, {-2 Sqrt[2/5], 13/(4 Sqrt[10]), 13/(12 Sqrt[2]), Sqrt[437/2]/6}} *)


Check the result.

In[1514]:= tmatrix.newdiag.Transpose[tmatrix] == symmat

(* Out[1514]= True *)


This method also scales considerably better than the one that uses an eigensystem since an exact LU is typically much faster than an exact eigensystem computation, when working over an integer matrix. For example, scaling to 20 x 20 the eigensystem method takes 51 seconds on my desktop, whereas the LU approach takes 0.03 seconds. For 200 x 200 the LU approach was 37 seconds.

• This worked, thanks! Is this method documented anywhere? Commented Jul 3, 2020 at 8:54
• No idea if it is documented. I coded a variant a long time ago to answer a question raised in the Usenet group sci.math.symbolic. Now I cannot locate that thread. It may have had a reference.What would be nice (but I have not worked out) would be a way to avoid having to force there to be no permutations. That would, among other things, allow for handling of matrices that have zeros as eigenvalues. Commented Jul 3, 2020 at 12:42

FindInstance gives the solution you expect.

Module[{X, j, T, n = 2}, X = {{-2, 0}, {0, 2}};
j = DiagonalMatrix[{1, -1}];
T = Array[t, {n, n}];
T /. FindInstance[
Simplify[
X - T.j.T\[ConjugateTranspose] == ConstantArray[0, {n, n}] //
ComplexExpand], Flatten[T]]]

{{{0, -Sqrt[2]}, {-Sqrt[2], 0}}}


Use Reduce for symbolic matrices.

Module[{X, j, T, n = 2},
X = {{-m, 0}, {0, m}};
j = DiagonalMatrix[{1, -1}];
T = Array[t, {n, n}];
Reduce[Thread[Flatten[X - T.j.T\[Transpose]] == 0], Flatten[T], Reals][[1]]
]


m < 0 && ((t[1, 1] == -Sqrt[-m] && t[1, 2] == 0 && t[2, 1] == 0 && (t[2, 2] == -Sqrt[-m] || t[2, 2] == Sqrt[-m])) || (t[1, 1] == Sqrt[-m] && t[1, 2] == 0 && t[2, 1] == 0 && (t[2, 2] == -Sqrt[-m] || t[2, 2] == Sqrt[-m])))

• Indeed, that is wonderful, thank you! Commented Jun 26, 2020 at 15:01
• Hmm, what if X contains symbols? Commented Jun 26, 2020 at 15:07
• You can use Reduce or use Solve and Select those solutions that are free of T. Also, it may be better to add the condition Det[T] != 0 to all of these approaches. Commented Jun 26, 2020 at 15:41

Turns out this appeared in the forum sci.math.symbolic back in the 90's. But I can no longer locate the thread. Anyway, here goes. First create a symmetric example matrix.

SeedRandom[1234]
mat = RandomInteger[{-5, 5}, {4, 4}];
symmat = mat + Transpose[mat]

(* Out[474]= {{-10, 6, -1, 4}, {6, -10, 5, -5},
{-1, 5, -2, 3}, {4, -5, 3, 4}} *)


Compute the eigensystem and normalize the eigenvectors.

{vals, vecs} = Eigensystem[symmat];
newvecs = Map[#/Sqrt[#.#] &, vecs];


We can form the signature diagonal matrix and the transformation matrix like so:

sqrroots = DiagonalMatrix[Sqrt[Abs[vals]]];
diagmat = DiagonalMatrix[Sign[vals]];
tmat = Transpose[newvecs].sqrroots;


The diagonal matrix is what we might expect.

In[508]:= diagmat

(* Out[508]= {{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, 1}} *)


I will show the numeric values of the transformation matrix since the exact form is large.

In[507]:= N[tmat]

(* Out[507]= {{-2.63202, 0.397454,
1.80326, -0.145953}, {3.04025, -0.407664,
1.05579, -0.43772}, {-1.24813,
0.565715, -1.08974, -0.652167}, {1.2636, 2.36751, 0.139461,
0.104966}} *)


Check the result.

tmat.diagmat.Transpose[tmat] - symmat // Expand // N // Chop

(* Out[503]= {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}} *)


Appendix

For a Hermitian matrix $$X$$, with eigenvectors arranged as columns of a matrix $$V$$ and eigenvectors arranged as a diagonal matrix $$\Lambda$$, we have

$$X\ V\ = \ V \ \Lambda$$

$$X\ V\ V^H = \ V \ \Lambda\ V^H$$

$$X\ = \ V \ \sqrt{|\Lambda|}\ J \sqrt{|\Lambda|}\ V^H$$

Here $$J$$ is a diagonal matrix with entries $$-1$$ or $$1$$, $$|\Lambda|$$ is the matrix of absolute element values (not Det), and the transformation turns out as $$T=V\sqrt{|\Lambda| }$$.

• Thank you, this is similar to the method I had previously implemented in Python, but much cleaner! I will end up using your answer since the result is easier to use Commented Jun 29, 2020 at 8:31
• Actually, in this case I notice that the answer is not reproduced, since diagmat is the opposite sign to my J matrix Commented Jun 29, 2020 at 8:33
• There are two negatives and two positives, so reversing signs will (fortunately) not change that.I did not try to canonicalize e.g. by having positives before negatives. Doing so, if that is wanted, would involve reordering eigenvalues and corresponding eigenvectors. Commented Jun 29, 2020 at 14:01
• @SubaThomas As for Normalize, I simply didn't think of it. Feel free to make edits. Commented Jun 30, 2020 at 16:27
• Also I realized I have a somewhat better method, one that delivers a more concise result. I'll post it once I brush off the 25 year old cyberdust. Commented Jun 30, 2020 at 16:45