Problem:
In the context of decision making in some optimization problems, I found that it is meaningful to compute $AGA^T$ with $A\in\mathbb R^{m\times n}$ and $G\in\mathbb R^{n\times n}$ a PSD symmetric matrix.
If we just multiply the matrices in order, we get a runtime of $O(mn^2+n^2m)$, I am wondering if this can be improved and if yes (most probably) by how much.
Attempt:
I tried to make a divide and conquer algorithm like Strassen's but it feels hard to do something better than computing $AGB^T$ as follow \begin{align*} &AGB^T\\ =&\begin{bmatrix} A_{11}G_{11}B_{11}^T+A_{11}G_{12}B_{12}^T+A_{12}G_{21}B_{11}^T+A_{12}G_{22}B_{12}^T& A_{11}G_{11}B_{21}^T+A_{11}G_{12}B_{22}^T+A_{12}G_{21}B_{21}^T+A_{12}G_{22}B_{22}^T\\ A_{21}G_{11}B_{11}^T+A_{21}G_{12}B_{12}^T+A_{22}G_{21}B_{11}^T+A_{22}G_{22}B_{12}^T& A_{21}G_{11}B_{21}^T+A_{21}G_{12}B_{22}^T+A_{22}G_{21}B_{21}^T+A_{22}G_{22}B_{22}^T\\ \end{bmatrix} \end{align*} Therefore we need to compute \begin{align*} A_{11}G_{11}B_{11}^T&& A_{11}G_{12}B_{12}^T&& A_{12}G_{21}B_{11}^T&& A_{12}G_{22}B_{12}^T\\ A_{11}G_{11}B_{21}^T&& A_{11}G_{12}B_{22}^T&& A_{12}G_{21}B_{21}^T&& A_{12}G_{22}B_{22}^T\\ A_{21}G_{11}B_{11}^T&& A_{21}G_{12}B_{12}^T&& A_{22}G_{21}B_{11}^T&& A_{22}G_{22}B_{12}^T\\ A_{21}G_{11}B_{21}^T&& A_{21}G_{12}B_{22}^T&& A_{22}G_{21}B_{21}^T&& A_{22}G_{22}B_{22}^T\\ \end{align*}
Now there are two cases, either we don't have much information about $A$, $G$ and $B^T$, or $G$ is symmetric and $B=A$, the third raw of the above is redundant and given by the second raw transposed. The gain is little, and probably even less than doing normal multiplication without computing the bottom left part of the matrix.
If $G$ is PSD, then there exists a matrix $B$ such that $G = BB^\top$. Then $$AGA^\top = (AB) (AB)^\top.$$ The latter is not necessarily faster (you would also have to compute $B$) but at least it is more compact. Depending on your purposes with the matrix $AGA^\top$ it might offer some improvements. For example, $$x^\top A GA^\top x = \| B^\top A^\top x \|^2.$$ Also, are both matrices $G$ and $A$ new every time you want to compute this? If only one of those matrices is changing, you might get some additional efficiency by exploiting certain matrix decompositions.