How to minimize quartic function for regularizing commutativity?

Question

My objective is regularizing a function over $X\in \mathbb{R}^{I\times J}$ such that its covariance is jointly diagonalizable with a positive semi-definite matrix $A \in \mathbb{R}^{I\times I}$. For this, I am trying to minimize an objective with commutativity between the covariance of $X$ and $A$. The objective roughly looks like: \begin{gather} f(X)+\|XX^\top A- AXX^\top\|_F^2, \end{gather} where $X$ is assumed to be zero mean. My problem is finding an optimizer rather than strictly solving the form I propose. So I am open to change the Frobenius norm of the regularization to some other norm if helpful.

I am planning on using an augmented Lagrangian to introduce auxiliary variables to separate regularization from $f(X)$. This changes my function for the regularization to something like: \begin{gather} \delta\|X-Y-\Lambda\|_F^2+\|YY^\top A- AYY^\top\|_F^2. \end{gather}

I have came upon polynomial sum of squares in my search but couldn't represent my function in that form.

My questions (any of them is sufficient for me) are:

1. Could a gradient descent approach minimize this function as it is quartic in $Y$?
2. If this was quadratic, I would have solved this by setting the derivative to zero but in this case the derivative is a cubic function, so it would have multiple roots. Would it suffice to select any root?
3. Would it make it easier to use $\ell_1$-norm rather than Frobenius norm to find an optimizer?

Notes: I am an electrical engineer and this is my first question on StackExchange. Please let me know if I missed anything.