In the course of my research (somewhat related to compressive sensing), I am trying to determine a complex, symmetric matrix $L$ (i.e. $L = L^T$) through the following optimization formulation:
$$ \hat{L} = \mathrm{argmin}_L W \circ ||AL||^2_2$$
where $\circ$ is the element-wise product. I have computed what I understand to be the Riemannian gradient as: $$ \nabla f_{pq} = \sum_i W_{iq}\bar{A}_{ip}(AL)_{iq}$$ By the argument presented in Srinivasan and Panda (https://arxiv.org/abs/1911.06491), I projected this gradient to the set of symmetric matrices as as: $$\mathrm{grad} f = \frac{1}{2}(\nabla f + \nabla f^T)$$
and update $L$ as: $$L_{k+1} = L_k - \alpha \cdot \mathrm{grad} f(L_k) $$
However, any initial point seems to be stationary for a few iterations before spontaneously becoming numerically unstable in MATLAB, so I feel that I am missing something here. I would appreciate any insight.
