(Edit: the original claim was much more ambitious).

It is possible to transform the representation into a representation which is slightly less general, without extending the field.

Let $\mu:G\to GL_d(\mathbb{R})$ s.t. $\rho(G)=M\mu(G)M^{-1}$ for $G\in GL_d(\mathbb{C})$. W.l.o.g. $\mu$ is orthogonal, i.e. $hh^\top=I$ for any $h\in\mu(G)$. Thus $$|G|I=\sum_{h\in\mu(G)} hh^\top=\sum_{h\in\mu(G)} hh^*= \sum_{g\in\rho(G)} M^{-1}gM (M^{-1}gM)^*=\sum_{g\in\rho(G)} M^{-1}gM M^*g^*(M^*)^{-1}.$$ Multiplying on the right by $M^*$ and on the left by $M$ one obtains $$|G|MM^*=\sum_{g\in\rho(G)} gM M^*g^*.$$

Up to scaling by a positive real, there is unique, by irredicibility of $\rho$, Hermitean p.s.d. matrix $H$ providing $G$-invariant sesquilinear scalar product, i.e. $gHg^*=H$ for any $g\in\rho(G)$. Therefore $H=MM^*$. Note that $H$ only records the info about $M$ up to a unitary factor. Namely, $M=PU$, with $U$ unitary and $P$ Hermitean positive definite ($PU$ is called left polar decomposition of $M$). Thus $H=MM^*=PUU^*P^*=PP^*=P^2$.

Note that $H$ affords a (unique, in fact), $LDL^*$ decomposition, with $D=diag(D_1,\dots,D_d)$, and $D_k>0$ for all $1\leq k\leq d$.

It is well-known that $L$ and $D$ can be computed by solving a system of linear equations, without extending the field (e.g. as on p.43 of our preprint (=MSc thesis of my student K.Hymabaccus)).

Set $M=LD^{1/2}U$. Then $$\mu(G)=M^{-1}\rho(G)M=U^*D^{-1/2}L^{-1}\rho(G)LD^{1/2}U.$$ Thus we can perform the conjugation by $L$, and so reduce to the case of $M=D^{1/2}U$.

(Edit: the original claim was much more ambitious).

It is possible to transform the representation into a representation which is slightly less general, without extending the field. Namely, this can be done if the transformation $M=PU$ which was used to make $\rho$ non-real, with $U$ Hermitean p.s.d. and $U$ unitary, satisfies $U=I$ (the decomposition $M=PU$, a polar decomposition, always exists).

Let $\mu:G\to GL_d(\mathbb{R})$ s.t. $\rho(G)=M\mu(G)M^{-1}$ for $G\in GL_d(\mathbb{C})$. W.l.o.g. $\mu$ is orthogonal, i.e. $hh^\top=I$ for any $h\in\mu(G)$. Thus $$|G|I=\sum_{h\in\mu(G)} hh^\top=\sum_{h\in\mu(G)} hh^*= \sum_{g\in\rho(G)} M^{-1}gM (M^{-1}gM)^*=\sum_{g\in\rho(G)} M^{-1}gM M^*g^*(M^*)^{-1}.$$ Multiplying on the right by $M^*$ and on the left by $M$ one obtains $$|G|MM^*=\sum_{g\in\rho(G)} gM M^*g^*.$$

Up to scaling by a positive real, there is unique, by irredicibility of $\rho$, Hermitean p.s.d. matrix $H$ providing $G$-invariant sesquilinear scalar product, i.e. $gHg^*=H$ for any $g\in\rho(G)$. Therefore $H=MM^*$. Note that $H$ only records the info about $M$ up to a unitary factor. Namely, $M=PU$, with $U$ unitary and $P$ Hermitean positive definite ($PU$ is called left polar decomposition of $M$). Thus $H=MM^*=PUU^*P^*=PP^*=P^2$.

Note that $H$ affords a (unique, in fact), $LDL^*$ decomposition, with $D=diag(D_1,\dots,D_d)$, and $D_k>0$ for all $1\leq k\leq d$.

It is well-known that $L$ and $D$ can be computed by solving a system of linear equations, without extending the field (e.g. as on p.43 of our preprint (=MSc thesis of my student K.Hymabaccus)).

Set $M=LD^{1/2}U$. Then $$\mu(G)=M^{-1}\rho(G)M=U^*D^{-1/2}L^{-1}\rho(G)LD^{1/2}U.$$ Thus we can perform the conjugation by $L$, and so reduce to the case of $M=D^{1/2}U$.

Yes, it is possible to transform the representation into a real form without extending the field. Let $\mu:G\to GL_d(\mathbb{R})$ s.t. $\rho(G)=M\mu(G)M^{-1}$ for $G\in GL_d(\mathbb{C})$. W.l.o.g. $\mu$ is orthogonal, i.e. $hh^\top=I$ for any $h\in\mu(G)$. Thus $$|G|I=\sum_{h\in\mu(G)} hh^\top=\sum_{h\in\mu(G)} hh^*= \sum_{g\in\rho(G)} M^{-1}gM (M^{-1}gM)^*=\sum_{g\in\rho(G)} M^{-1}gM M^*g^*(M^*)^{-1}.$$ Multiplying on the right by $M^*$ and on the left by $M$ one obtains $$|G|MM^*=\sum_{g\in\rho(G)} gM M^*g^*.$$

Up to scaling by a positive real, there is unique, by irredicibility of $\rho$, Hermitean p.s.d. matrix $H$ providing $G$-invariant sesquilinear scalar product, i.e. $gHg^*=H$ for any $g\in\rho(G)$. Therefore $H=MM^*$. Note that $H$ affords a (unique, in fact), $LDL^*$ decomposition, with $D=diag(D_1,\dots,D_d)$, and $D_k>0$ for all $1\leq k\leq d$.

It is well-known that $L$ and $D$ can be computed by solving a system of linear equations, without extending the field (e.g. as on p.43 of our preprint (=MSc thesis of my student K.Hymabaccus)).

Set $M=LD^{1/2}$. Then $$\mu(G)=M^{-1}\rho(G)M=D^{-1/2}L^{-1}\rho(G)LD^{1/2}.$$ Finally, note that $D^{1/2}\mu(G)D^{-1/2}=L^{-1}\rho(G)L\subset GL_d(\mathbb{R})$, as $D$ is real.

Computationally, $H$ can only be found of the form $\alpha MM^*$, for some $\alpha>0$. Then $\alpha$ gets absorbed in $D$, so the same computation can actually be carried out. I suppose we should add this functionality to our GAP package RepnDecomp.

(Edit: the original claim was much more ambitious).

It is possible to transform the representation into a representation which is slightly less general, without extending the field.

Let $\mu:G\to GL_d(\mathbb{R})$ s.t. $\rho(G)=M\mu(G)M^{-1}$ for $G\in GL_d(\mathbb{C})$. W.l.o.g. $\mu$ is orthogonal, i.e. $hh^\top=I$ for any $h\in\mu(G)$. Thus $$|G|I=\sum_{h\in\mu(G)} hh^\top=\sum_{h\in\mu(G)} hh^*= \sum_{g\in\rho(G)} M^{-1}gM (M^{-1}gM)^*=\sum_{g\in\rho(G)} M^{-1}gM M^*g^*(M^*)^{-1}.$$ Multiplying on the right by $M^*$ and on the left by $M$ one obtains $$|G|MM^*=\sum_{g\in\rho(G)} gM M^*g^*.$$

Up to scaling by a positive real, there is unique, by irredicibility of $\rho$, Hermitean p.s.d. matrix $H$ providing $G$-invariant sesquilinear scalar product, i.e. $gHg^*=H$ for any $g\in\rho(G)$. Therefore $H=MM^*$. Note that $H$ only records the info about $M$ up to a unitary factor. Namely, $M=PU$, with $U$ unitary and $P$ Hermitean positive definite ($PU$ is called left polar decomposition of $M$). Thus $H=MM^*=PUU^*P^*=PP^*=P^2$.

Note that $H$ affords a (unique, in fact), $LDL^*$ decomposition, with $D=diag(D_1,\dots,D_d)$, and $D_k>0$ for all $1\leq k\leq d$.

It is well-known that $L$ and $D$ can be computed by solving a system of linear equations, without extending the field (e.g. as on p.43 of our preprint (=MSc thesis of my student K.Hymabaccus)).

Set $M=LD^{1/2}U$. Then $$\mu(G)=M^{-1}\rho(G)M=U^*D^{-1/2}L^{-1}\rho(G)LD^{1/2}U.$$ Thus we can perform the conjugation by $L$, and so reduce to the case of $M=D^{1/2}U$.