First, a warm welcome to StackExchange! Consider this as a bit of brainstorming.
Given that the diagonal elements are positive, I think there is a particular geometry of your parameters: They live in the Generalized Stiefel Manifold, where each point is a scaled orthonormal frame of dimension $T$ embedded in $\mathbb{R}^N$.
Let me rewrite :
$$
\begin{align*}
SS^\top &= \Lambda\\
\frac{1}{\sqrt{\Lambda}}SS^\top\frac{1}{\sqrt{\Lambda}^\top} &= \frac{1}{\sqrt{\Lambda}}\Lambda\frac{1}{\sqrt{\Lambda}^\top}\\
B SS^\top B^\top &= B\Lambda B^\top\\
OO^\top &= I
\end{align*}
$$
where $\Lambda$ is the $N\times N$ diagonal matrix, $B=\frac{1}{\sqrt{\Lambda}}\succ0$ and $O=BS$. Similarly, $S=\sqrt{\Lambda}O$. Division and square roots are element-wise. Now $O$ lives on the usual Stiefel manifold. When $B$ is known, the problem can simply be solved by a Riemannian optimization on the Generalized Stiefel Manifold (see manopt). However, now, this is not the case.
So let's take a second look at the problem at hand:
$$\min_{S} f(S) \mbox{ s.t. } SS^T \mbox{is diagonal}$$
This can be formulated in terms of the new variables:
$$\min_{\Lambda,O} f(\sqrt{\Lambda}O) \mbox{ s.t. } OO^T =I \,\wedge\, \Lambda>0$$
where $O$ is Stiefel, and $\Lambda>0$. The latter can also be denoted by a vector $\mathbf{\lambda}\equiv \text{diag}(\Lambda)$ and $\lambda>0$. The true exponential and logarithm maps for the Stiefel manifold are available (for instance in Manopt) and $\lambda$ has a simpler constraint (for projection). So from this point on, I am being more heuristic and would suggest to alternatively optimize these two:
- Fix $\lambda$, optimize $O$ on the Stiefel manifold using the chain rule:
$$
\nabla_O f(\sqrt{\Lambda}O) = \sqrt{\Lambda} \nabla_{S} f(S)
$$
- Fix $O$, optimize for $\sqrt{\lambda}$ (can also use positivity constraint):
$$
\nabla_{\sqrt{\Lambda}} f(\sqrt{\Lambda}O) = O \nabla_{S} f(S)
$$
You can of course use gradient descent but also other sophisticate algorithms like Riemannian-LBFGS / trust-region.
