**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 is $d \sigma$


-------------

### Verification in code  

Setting up the variables:

    SeedRandom[0];
    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](https://www.sciencedirect.com/science/article/pii/S0024379516000331) that works on the sympletic analogue cites this [paper](https://www.sciencedirect.com/science/article/pii/0024379593000566) 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.

-------------------------