# How to solve for $\mathbf{X}$ in a matrix equation $\mathbf{X A} \mathbf{X}^\top = \mathbf{B}$?

I need to find a matrix $$\mathbf{X}$$ such that it belongs to $$O(D,D)$$ and satisfies the equation:

$$\mathbf{X A X}^\top = \mathbf{B}$$

where both matrices are square, $$\mathbf{A}$$ is symmetric and $$\mathbf{B}$$ is diagonal. I have no idea how to easily do this in Mathematica. Any tips?

Thanks!

• What is O(D,D)? And please provide some code about A and B. Jul 27, 2021 at 22:51
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 27, 2021 at 23:24
• Do A and B have any kind of symmetry? Are they square and the same size? Jul 28, 2021 at 5:45
• From a mathematics point of view, see this: mathoverflow.net/questions/145225/… Jul 28, 2021 at 9:14
• I suppose you realize that neither answer so far has ${\bf X} \in O(D,D)$....I don't know a linear algebra solution off hand, brute force methods are not proving robust. Note $O(D,D) \ne O(2D)$ and generally its matrices are not orthogonal. Instead they satisfy $\sigma . X^t . \sigma = X^{-1}$, where $\sigma = diag(1,\dots,1,-1,\dots,-1)$ with $D$ ones and minus-ones. I'm not that familiar with, though. Jul 30, 2021 at 19:16

For real, square and symmetric matrices A, B, X, the solution of $$XAX^T=B,$$ is given by $$X=B (B^{-1}A^{-1})^{1/2}.$$

For example, in Mathematica this can be written as

A = {{1, 2}, {2, 3}};
X = {{3, 4}, {4, -1}};
B = X.A.Transpose[X];


and the solution can be verified to be

In:= X - B.MatrixPower[Inverse[B].Inverse[A], 1/2] // N // Chop
Out= {{0, 0}, {0, 0}}


as expected.

In the more general case, where the matrices are non-symmetric etc, the situation is more complicated, but still a numerical solution may be achieved using SchurDecomposition, see the classic paper "A Schur Method for Solving Algebraic Riccati Equations" by Laub, A.J., which can be found here.

• Thank you so incredibly much! My X matrix is a complete mess but it seems to have worked. I'm very much a Mathematica beginner, I wonder, is Chop needed here because after all the engine does in fact a numerical calculation? Jul 28, 2021 at 16:25
• I'm happy to see it worked! I used Chop as I took the difference of two nearly identically quantities which ends up being nearly zero (~10^{-15} etc). If you just need X itself, Chop is not needed. Jul 28, 2021 at 16:27
• @4815162342 By the way, if you are happy with the solution, please consider accepting the answer by clicking the tick symbol on the left :-) Jul 28, 2021 at 16:28
• I just did. Sorry, still learning how to use the platform in terms of asking questions after years and years of merely lurking around :-) Jul 28, 2021 at 16:32
• @4815162342 Often we ask people to wait ~24hrs so that people around the world in each timezone have a chance to contribute an answer. Marking a Q&A as answered early sometimes discourages people from even reading the question. (And in this case, the question was initially stated incompletely.) Jul 28, 2021 at 23:45

The SVD of a real symmetric matrix $$A$$ is simple, fast, accurate way to get $$X$$ and $$B$$ (and $$A$$ need not be invertible):

(* set up *)
SeedRandom;
amat = #\[Transpose] . # &@RandomReal[{-1, 1}, {5, 5}];

(* solution *)
{xmat, bmat} = Most@SingularValueDecomposition[amat];
xmat = xmat\[Transpose];


Check:

Chop[xmat . amat . xmat\[Transpose]] == bmat
(*  True  *)

xmat . amat . xmat\[Transpose] // Threshold // MatrixForm
bmat // MatrixForm It should be clear that the above may be adapted to the original problem statement, amended with the suitable hypotheses:

Given two real, symmetric, orthogonally-similar matrices $$A$$ and $$B$$, not necessarily invertible, find an orthogonal matrix $$X$$ such that $$X.A.X^t=B$$.

(* set up: 5x5, rank 4 *)
SeedRandom;
amat = #\[Transpose] . # &@RandomReal[{-1, 1}, {4, 5}];
conj = Orthogonalize@RandomReal[{-1, 1}, {5, 5}];
bmat = conj . amat . conj\[Transpose];

(* solution *)
u1 = First@SingularValueDecomposition[amat];
u2 = First@SingularValueDecomposition[bmat];
xmat = u2 . u1\[Transpose];


Check:

Norm[xmat . amat . xmat\[Transpose] - bmat]/Norm[amat]
(*  3.12068*10^-16  *)


Since u1 and u2 are not unique, xmat is not necessarily the same as conj, and the solution xmat is unique only up to the uniqueness of the SVD. Given another solution such as conj, one can find its relationship to u1 and u2 as follows:

cj = u2\[Transpose] . conj . u1;
conj - u2 . cj . u1\[Transpose] // Norm
(*  9.36499*10^-16  *)


Perhaps the efficiency of SingularValueDecomposition on numerical problems is not widely appreciated. It is much more accurate than the Inverse/MatrixPower power approach and quite a bit faster. Here is a code for comparing the two methods. (Takes about a minute to run for 100 data points.)

SeedRandom;
(data = Transpose[#, {3, 1, 2}] &@Table[
(* Set up: random symm. mat. of random sizes *)
n = RandomInteger[{2, 800}];
amat = #\[Transpose] . # &@RandomReal[{-1, 1}, {n, n}];
conj = Orthogonalize@RandomReal[{-1, 1}, {n, n}];
bmat = conj . amat . conj\[Transpose];
(* SVD *)
t1 = (v1 = Last@SingularValueDecomposition[amat];
v2 = Last@SingularValueDecomposition[bmat];
xmat1 = v2 . v1\[Transpose];) // AbsoluteTiming // First;
(* Inverse *)
t2 = (xmat2 =
bmat . MatrixPower[Inverse[bmat] . Inverse[amat], 1/2];) //
AbsoluteTiming // First;
{{{n,
Norm[xmat1 . amat . xmat1\[Transpose] - bmat]/
Norm[amat]}, {n,
Norm[xmat2 . amat . xmat2\[Transpose] - bmat]/Norm[amat]}},
{{n, t1}, {n, t2}}},
{100}]) //
Frame -> True,
FrameLabel -> {"Dimension", #2},
PlotLegends -> {"SVD", "Inverse"}] &,
{#, {"Relative Error (2-norm)", "Timing"}}] &


Inverse throws a few warnings:

...
Inverse::luc: Result for Inverse of badly conditioned matrix {<<1>>} may contain significant numerical errors.
...
General::stop: Further output of Inverse::luc will be suppressed during this calculation.

The first plot shows the relative error of the two methods on random matrices versus dimension. It illustrates the stability of the SVD method. The second shows the time each run took. • Maybe it should be v1 = Last@SingularValueDecomposition[amat]; and so forth. They're the same as u1, u2, since $A$, $B$ are symmetric, but $V$ operates on the domain... Jul 29, 2021 at 18:08

Note : As far as I understand, the question posed by OP has not been answered as it does not seem that the $$O(D,D)$$ symmetry has been addressed. That the question does not seem to be answered can be checked by looking at the last few comments below the question by OP.

### Discussion and mathematical derivation of the algorithm

algorithm section below

Consider first that such a relationship exists.

Then, denoting $$\sigma$$ as

Diag$$(1,1,...,-1,-1,..,-1)=\{ \{ 1_{n \times n},0_{n \times n}\},\{0_{n \times n},-1_{n \times n}\}\}$$

or in code

DiagonalMatrix[ConstantArray[1,n]~Join~ConstantArray[-1,n]],

as $$X\in O(n,n)$$ we find :

$$XAX^{T}=X A \sigma X^{-1} \sigma=B$$

Thus, using $$\sigma^2=1$$

$$X A \sigma X^{-1}=B\sigma .$$

As both $$B$$ and $$\sigma$$ are diagonal, that last equation represents a diagonalization of $$A \sigma$$. The questions are thus:

• Is $$A \sigma$$ diagonalizable ?

• If we were to diagonalize $$A\sigma$$, would the change of basis matrix $$X$$ really belong to $$O(n,n)$$ ?

I see no reason that $$A \sigma$$ should be diagonalizable but as diagonalizable matrices are dense in the space of matrices we will consider that, for a random symmetric matrix $$A$$, $$A \sigma$$ is likely diagonalizable.

For the second point, we will assume a hypothetical natural extension of results from a paper that is cited in the last discussion below. That discussion concerns results from an analogue scenario of sympletic matrices. Assuming that the paper can be extended in a natural way, there exists (I explain why in the last section, the reader can go there now if they wish to do so) an $$\tilde{X}$$ that verifies the eigen factorization mentioned before

$$\tilde{X} A \sigma \tilde{X}^{-1}=B\sigma$$

and $$\tilde{X}$$ is an element of $$O(n,n)$$ that is,

$$\tilde{X}^T=\sigma \tilde{X}^{-1} \sigma .$$

Notice that there is no constraint for $$X$$ to be equal to $$\tilde{X}$$. Indeed for any invertible matrix $$C$$ that commutes with $$B\sigma$$, $$C \tilde{X}$$ also verifies the last eigen factorization :

$$C \tilde{X} A \sigma \tilde{X}^{-1} C^{-1} =B\sigma$$

is equivalent to

$$\tilde{X} A \sigma \tilde{X}^{-1} = C^{-1} B\sigma C$$

since $$C$$ commutes with $$B\sigma$$ we obtain the same equation.

Two matrices that commute are simultaneously diagonalizable. For a generic random matrix, the eigenvalues are distinct and so up to normalization and ordering there is a unique basis on which $$C$$ is also diagonal. The objective is then to find a diagonal and invertible $$C$$ such that :

$$X A \sigma X^{-1}=B\sigma$$

and for $$\tilde{X}= C X$$

$$\tilde{X}^T=\sigma \tilde{X}^{-1} \sigma$$

that is,

$$X^T C = \sigma X^{-1} C^{-1} \sigma$$.

Or equivalently for the last equation using $$\sigma^2$$ :

$$X \sigma X^T = \sigma C^{-2}$$

That last equation allows us to find C by computing $$X \sigma X^T$$ and taking the inverse of the square root for the first $$n$$ components of the diagonal, then, for $$i^2=-1$$, taking $$i$$ times the inverse of the square root of the last $$n$$ components. That is, $$X \sigma X^T=\textrm{Diag}_1 \oplus \textrm{Diag}_2$$ and $$C=\left(\textrm{Diag}_1\right)^{-1/2} \oplus i\left(\textrm{Diag}_2\right)^{-1/2}$$.

### The algorithm

code section below

• Step 1: Find the eigen decomposition of $$A \sigma$$. The output is (X,d) where d is diagonal and X is the change of basis matrix.

• Step 2: Compute $$X \sigma X^T$$ which will be of the form $$\textrm{Diag}_1(n) \oplus \textrm{Diag}_2(n)$$ , compute $$C=(\textrm{Diag}_1)^{-1/2} \oplus i (\textrm{Diag}_2)^{-1/2}$$. The solution is then $$\tilde{X}=X C$$ and the diagonal matrix is $$d \sigma$$

### Verification in code

Setting up the variables:

SeedRandom;
n = 8;
A = (#\[Transpose] + #) &@RandomReal[{-1, 1}, {n, n}];
σ = DiagonalMatrix[ConstantArray[1, n/2]~Join
~ConstantArray[-1, n/2]];
Aσ = A . σ;


Implementing the algorithm:

{iX, d} = JordanDecomposition[Aσ];
Xaux = Inverse@iX;
aux = Xaux . σ . Transpose[Xaux] // Chop // Diagonal;
aux2 = DiagonalMatrix[
Sqrt[aux[[1 ;; n/2]]]~Join~(I*Sqrt[aux[[n/2 + 1 ;; n]]])];
X = Inverse[aux2] . Xaux;


Verification:

d // DiagonalMatrixQ
X . A . Transpose[X] - d . σ // Norm (* check that
the factorization works *)
Inverse[X] - σ . Transpose[X] . σ // Norm (* check that
X is an element of O(n,n) *)

(* True

4.85014*10^-15

5.88943*10^-15 *)


### Discussion on the mathematical validity of the algorithm

This paper that works on the sympletic analogue cites this paper in proposition 5 and states that (the following does not need to be understood), if two matrices are auto adjoint with respect to the symplectic scalar product and they are similar then the change of basis matrix is unitarily symplectic.

If I take the hypothesis that the same applies here then basically from the eigen decomposition:

$$X A \sigma X^{-1}=B\sigma$$

we see that $$A \sigma$$ and $$B\sigma$$ are similar that is there is a change of basis that relates the two. Extending the results from the papers cited to this case, if $$A \sigma$$ and $$B \sigma$$ are similar and are auto adjoint (a generalization of being symmetric) with respect to the pseudo scalar product whose matrix is $$\sigma$$=Diag$$\left(1,1,...-1,-1...-1\right)$$, that is they verify the equation :

$$\sigma Y^{T} \sigma= Y$$

where the left hand side is the adjoint with respect to $$\sigma$$, then there exists a change of matrix basis $$\tilde{X}$$ that is an element of $$O(n,n)$$. Basically, if both matrices respect the $$O(n,n)$$ symmetry and are similar then their is a change of basis matrix that respects the symmetry and that is the analogue of a unitary matrix.

It suffices then to show that:

$$\sigma (A\sigma)^{T} \sigma= A\sigma$$

and

$$\sigma B^{T} \sigma= B$$

which can be checked using $$\sigma^T=\sigma$$, the fact that diagonal matrices commute and that $$\sigma^2$$.

Moreover, taking the hypothesis that one may extend the results from the first paper, the minimal constraint for the decomposition in the original question to work is that A is symmetric and $$A \sigma A \sigma$$ is diagonalizable.