So this means opting for approach (2) and registering a function to compute C(x) := -A^{-1}(x) B and its derivative? Otherwise I only see how to encode a full n^2 constraint Y as
\begin{align}
Y = C D C^{\top} \iff A Y A^{\top} = B D B^{\top}
\end{align}
Or is there a way to rearrange the computation to enforce only y = {\rm diag}(Y) that I’m missing?