5
$\begingroup$

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!

$\endgroup$
9
  • 5
    $\begingroup$ What is O(D,D)? And please provide some code about A and B. $\endgroup$
    – cvgmt
    Jul 27 at 22:51
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Jul 27 at 23:24
  • 1
    $\begingroup$ Do A and B have any kind of symmetry? Are they square and the same size? $\endgroup$
    – mikado
    Jul 28 at 5:45
  • $\begingroup$ From a mathematics point of view, see this: mathoverflow.net/questions/145225/… $\endgroup$
    – Hans Olo
    Jul 28 at 9:14
  • 2
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Jul 30 at 19:16
8
$\begingroup$

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[4]:= X - B.MatrixPower[Inverse[B].Inverse[A], 1/2] // N // Chop
Out[4]= {{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.

$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – 4815162342
    Jul 28 at 16:25
  • $\begingroup$ 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. $\endgroup$
    – Hans Olo
    Jul 28 at 16:27
  • $\begingroup$ @4815162342 By the way, if you are happy with the solution, please consider accepting the answer by clicking the tick symbol on the left :-) $\endgroup$
    – Hans Olo
    Jul 28 at 16:28
  • 1
    $\begingroup$ I just did. Sorry, still learning how to use the platform in terms of asking questions after years and years of merely lurking around :-) $\endgroup$
    – 4815162342
    Jul 28 at 16:32
  • $\begingroup$ @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.) $\endgroup$
    – Michael E2
    Jul 28 at 23:45
6
$\begingroup$

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[0];
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

Addendum

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[2];
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  *)

Addendum 2

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[3];
(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}]) //
  MapThread[ListLogPlot[#,
      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.

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ 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... $\endgroup$
    – Michael E2
    Jul 29 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.