# 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. – Daniel Lichtblau Jun 26 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) – Joe Bentley Jun 28 at 8:15
• Ooof, that's right. Let me retry that: are the rows of T orthogonal to one another? – Daniel Lichtblau Jun 28 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? – Joe Bentley Jul 3 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. – Daniel Lichtblau Jul 3 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! – Joe Bentley Jun 26 at 15:01
• Hmm, what if X contains symbols? – Joe Bentley Jun 26 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. – Suba Thomas Jun 26 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 – Joe Bentley Jun 29 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 – Joe Bentley Jun 29 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. – Daniel Lichtblau Jun 29 at 14:01
• @SubaThomas As for Normalize, I simply didn't think of it. Feel free to make edits. – Daniel Lichtblau Jun 30 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. – Daniel Lichtblau Jun 30 at 16:45