Projection onto the Set of Circulant Matrices Defining $ \mathcal{C}_{n} $ the set of Real Circulant Matrices.  
The orthogonal projection of a given matrix $ Y \in \mathbb{R}^{n \times n} $ onto the set is given by the following minimization problem:
$$ \arg \min_{X \in \mathcal{C}_{n} } \frac{1}{2} {\left\| X - Y \right\|}_{F}^{2} $$
Where $ {\left\| \cdot \right\|}_{F} $ is the Frobenius Norm.
Since the set of Circulant Matrices is Convex Set and the objective function is convex the whole problem is a convex optimization problem.
 A: Solution 001
Defining the Forward Shift Matrix:
$$ \Pi = \begin{bmatrix}
0 & 1 & 0 & \cdots & 0 \\ 
0 & 0 & 1 &  & 0 \\ 
\vdots &  & \ddots & \ddots & \vdots \\ 
0 &  &  &  & 1 \\ 
1 & 0 &  & \cdots & 0
\end{bmatrix} $$
One could see that the set $ \left\{ {\Pi}^{k} \mid k \in \left\{ 0, 1, \ldots, n - 1 \right\} \right\} $ is Orthogonal Basis (Can become Orthonormal with division by $ \sqrt{n} $) with respect to Frobenius Inner Product.
Since Circulant Matrix is composed by shifting forward the first row of the matrix one could see that any Circulant Matrix $ C $ with its first row given by $ c $ can be built as following:
$$ C = \sum_{k = 1}^{n} {c}_{k} {\Pi}^{k - 1} $$
Namely, this Shift Forward matrix is the basis of Circulant Matrices.
Hence, just like we build Fourier Series Coefficients by projection onto the basis, given a Matrix $ A $ we can create its projection by:
$$ {a}_{k} = \frac{1}{n} \left \langle A, {\Pi}^{k - 1} \right \rangle, \; k \in \left\{ 1, 2, \ldots, n \right\} $$
Where $ \left \langle \cdot, \cdot \right \rangle $ is the Frobenius Inner Product - $ \left \langle A, B \right \rangle = {A}^{H} B $.  
Then the Circulant Matrix closest to A is given by:
$$ {P}_{\mathcal{C}^{n}} \left( A \right) = \sum_{k = 1}^{n} {a}_{k} {\Pi}^{k - 1} $$
The funny thing is this gives me better result than my solution above even if I allow Complex Matrices (Which means something isn't right).
A: Solution 002
Since $ X $ is Circulant Matrix it is diagonalizable by the DFT Matrix.
Hence $ X = {F}^{H} \operatorname{diag} \left( x \right) F $ where $ x $ is the Fourier Transform of the first row of $ X $.
Now, Since $ F $ is Unitary Matrix and our objective function is Euclidean Distance which is invariant to multiplication by Unitary Matrix the problem can be rewritten as following:
$$ \arg \min_{X \in \mathcal{C}_{n} } \frac{1}{2} {\left\| X - Y \right\|}_{F}^{2} = \arg \min_{X \in \mathcal{C}_{n} } \frac{1}{2} {\left\| F X {F}^{H} - F Y {F}^{H} \right\|}_{F}^{2} = \arg \min_{x \in \mathbb{C}^{n} } \frac{1}{2} {\left\| \operatorname{diag} \left( x \right) - F Y {F}^{H} \right\|}_{F}^{2} $$
The above is clearly minimized by $ {x}_{i} = {\left( F Y {F}^{H} \right)}_{ii} $.
Namely $ x $ is the diagonal of $ F Y {F}^{H} $.  
To generate the corresponding Real Circulant Matrix one should use the same deocmposition:
$$ X = {F}^{H} \operatorname{diag} \left( x \right) F $$
In cases $ Y $ is real matrix the above will yield a real matrix as well (Up to numerical issues due to quantization).
