The optimization problem is defined as
$\min_{\mathbf X\in \mathbb R^{m\times n}}\|\mathbf{Y}-\mathbf{A}_1\mathbf{XA}_2\|_F^2+\lambda\|\mathbf{X}\|^2_F$
where $\lambda>0$, and $\mathbf{A}_1\in \mathbb{R}^{m\times m}, \mathbf{A}_2\in \mathbb{R}^{n\times n}$ are circulant matrices.
What came to my mind is that let the derivative $\mathbf{A}_1^T(\mathbf{Y}-\mathbf{A}_1\mathbf{XA}_2)\mathbf{A}_2^T + \lambda\mathbf{X}=0$
I noticed that these circulant matrices can be decomposed using DFT as $\mathbf{A}_1=F_m\Lambda_1 F_m^H$, $\mathbf{A}_2=F_n\Lambda_1 F_n^H$,
but I still do not know how to simplify it further.
I assume that $F_m$ denotes the unitary $m \times m$ Fourier transform.
If that is the case, you can rewrite the first norm as $$ \begin{aligned} \|Y - F_m \Lambda_1 F_m^* X F_n \Lambda_2 F_n^* \|_{\mathsf{F}} &= \| F_m F_m^* Y F_n F_n^* - F_m \Lambda_1 F_m^* X F_n \Lambda_2 F_n^* \|_{\mathsf{F}} \\ &= \|F_m^* Y F_n - \Lambda_1 F_m^* X F_n \Lambda_2\|_{\mathsf{F}}, \end{aligned} $$ since $F_m$ and $F_n$ are unitary matrices. Similarly, the second norm is $$ \|X\|_{\mathsf{F}} = \|F_m F_m^* X F_n^* F_n\|_{\mathsf{F}} = \|F_m^* X F_n\|_{\mathsf{F}} $$
Now redefine $\tilde{X} := F_m^* X F_n$, $\tilde{Y} = F_m^* Y F_n$ and solve
$$ \min_{\tilde{X}} \|\tilde{Y} - \Lambda_1 \tilde{X} \Lambda_2\|_{\mathsf{F}}^2 + \lambda \| \tilde{X} \|_{\mathsf{F}}^2. $$
After you solve for $\tilde{X}$, you can obtain the original solution by $X = F_m \tilde{X} F_n^*$.