Given a symmetric positive-definite matrix $M$, find all $A$ such that $A^\top M A=M$

Given $M$ a real symmetric positive-definite matrix, I would like to characterise all matrices $A$ such that $A^\top M A=M$. Note that the question of finding $A$ solutions to $A^\top M A=M$ for all positive-definite matrices $M$ has already been answered here.

Necessarily $\det(A)^2=1$.

Reducing the problem to diagonal p.q. matrices $M$ is symmetric, so diagonalisable in an orthonormal basis: $M=Q^\top D Q$ with $D$ diagonal and $QQ^\top=I$. Then,

\begin{align} &A^\top Q^\top D Q A=Q^\top D Q \\ &(QA)^\top D(QA)=Q^\top D Q \\ &(QAQ^\top)^\top D (QAQ^\top)=D \end{align}

so I think the problem can be simplified in

$$\text{Find the matrices B such that B^\top D B=D}$$ where $D=\operatorname{diag}(d_1,\dots,d_n)$, $d_i>0$.

and then recover $A$ from $A=Q^\top B Q$ ($Q$ and $D$ are known since $M$ is given). The $d_i>0$ stems for $M$ positive-definite.

Trying to solve this problem The $2^n$ matrices $B=\operatorname{diag}(\pm 1,\dots, \pm 1)$ are obvious solutions. What about the other ones? Writing $B^\top D B=D$ in components, the term $(i,j)$ is given by $$\sum_k B_{k,i}d_k B_{k,j}=d_i\delta_{ij}$$ In particular, $$\sum_k B_{k,i}^2d_k =d_i.$$ However, I don't see any clear conclusions using these equalities.

Question: How can to find/characterise the $B$ such that $B^\top D B=D$?

The matrix $D$ has a positive square root, name it $E$, which is a diagonal matrix with entries $\sqrt{d_i}>0$. Hence $E^{\top}=E$ and $(E^{-1})^{\top}=E^{-1}$.

From $B^{\top}DB=D$ we obtain $E^{-1}B^{\top}E=EB^{-1}E^{-1}$, hence $$(EBE^{-1})^{\top}=(E^{-1})^{\top}B^{\top}E^{\top}=EB^{-1}E^{-1}=(EBE^{-1})^{-1}.$$ So $EBE^{-1}$ is an orthogonal (unitary) matrix.

Starting from any orthogonal (unitary) matrix $C$ we set $B:=E^{-1}CE$. Then $B^{\top}=E^{\top}C^{\top}(E^{-1})^{\top}=EC^{\top}E^{-1}$ and so $$B^{\top}DB=B^{\top}EEB=EC^{\top}E^{-1}EEE^{-1}CE=EC^{\top}CE=E^2=D.$$

So, it turns out that the set of $n \times n$ matrices is isomorphic to the set of bilinear forms (bilinear mapping $B:V\times V \rightarrow \mathbb{F}$, $\mathbb{F}$ the underlying field of $V$), and bilinear forms are able to be expressed with matrix representations. From there, it's not difficult to show that a bilinear form is non-degenerate if and only if its matrix representation is non-singular, and that it is symmetric if and only if its matrix representation is symmetric.

In the case where the bilinear form is symmetric, there are a lot of really nice properties that rise out of it. Given a symmetric or antisymmetric bilinear form, there exists an associated automorphism group, $Aut(B)$:

$$Aut(B) = \{ T \in GL(V): B(T(v),T(w)) = B(v,w), \forall v, w \in V \}$$

From here, it is very intuitive to think of the matrix representation:

$$B(T(v),T(w)) = B(v,w), \forall v, w$$ $$\Leftrightarrow [x]^T[B][w] = {([T][x])}^T[B]([T][w])$$ $$\Leftrightarrow [x]^T[T]^T[B][T][w]$$

Which shows that the automorphism group of $B$ is precisely the set of matrices $[A]^T[B][A] = [B]$, called $G_{B,\mathcal{X}}$ for a basis $\mathcal{X}$.

The reason why this is relevant is because that, there aren't really any special formulas here that'll show you much more about the automorphism group of a bilinear form. They have some special properties but a big reason why those properties are meaningful and interesting, is because of their relationship to bilinear forms.

One subset of bilinear forms you might remember from linear algebra is the inner product, and how it was used to determine information concerning magnitude and angle. In the case of symmetric, positive definite matrices $[B]$, by examining them as bilinear forms, they are precisely the set of inner products on a real vector space. Then, each automorphism group is the associated set of linear transformations corresponding to rotations plus reflections, with the restriction to $A \in Aut(B)$ such that $Det(A) = +1$ being the proper rotation group.

When examining symmetric matrices as representations of bilinear forms, they break up into congruence classes typified by special automorphism groups, orthogonal group over $\mathbb{R}$ of signature $(p,q)$, or the indefinite orthogonal group $(p,q)$.

Edit: Also, an important fact I forgot to mention--the set of linear operators $Aut(B)$ actually loses some its meaning in its matrix representation. Specifically, for distinct bases $\mathcal{X}, \mathcal{Y}$,

$$G_{B,\mathcal{X}} = S^{-1} G_{B,\mathcal{Y}} S$$

Where $S = [Id]_{\mathcal{X}\mathcal{X}}$. In other words, the matrix representation of the automorphism group of a bilinear form $B$ varies under similarity transformation, so examining $[B]$ as the representation of a linear operator removes the connection between $G_{B,\mathcal{X}}$ and $Aut(B)$, and also shows that it is actually not the case that matrices $A$ such that $A^TMA = M$ for symmetric, positive definite $M$ are the same as matrices $C$ such that $C^TDC$, where $D$ is the diagonal matrix similar to $M$.