Nonlinear optimization: Optimizing a matrix to make its square is close to a given matrix. I'm trying to solve a  minimization problem whose purpose is to optimize a matrix whose square is close to another given matrix. But I can't find an effective tool to solve it.
Here is my problem:

Assume we have an unknown  Q with parameter $q11, q12,q14,q21,q22,q23,q32,q33,q34,q41,q43,q44$,  and a given matrix G, that is, 
  $Q=\begin{pmatrix}
 q11&q12 &0 &q14 \\q21&q22& q23&0\\ 0&q32& q33&q34\\ q41&0& q43&q44\\ 
 \end{pmatrix} $,   $G=\begin{pmatrix}
 0.48&0.24 &0.16 &0.12 \\ 0.48&0.24 &0.16 &0.12\\0.48&0.24 &0.16 &0.12\\0.48&0.24 &0.16 &0.12
 \end{pmatrix} $,
The problem is how to find the values of $q11, q12,q14,q21,q22,q23,q32,q33,q34,q41,q43,q44$ such that the square of $Q$ is very close to matrix $G$.

I choose to minimize the Frobenius norm of their difference, that is, 

$ Q* ={\arg\min}_{Q}   \| Q^2-G\|_F$
s.t.  $0\leq q11, q12,q14,q21,q22,q23,q32,q33,q34,q41,q43,q44 \leq 1$,$\quad$
     $\quad$ $q11+q12+q14=1$,
     $\quad$ $q21+q22+q23=1$,
     $\quad$ $q32+q33+q34=1$,
     $\quad$ $q41+q43+q44=1$.

During those days, I am frustrated to find a effective tool to execute the above optimization algorithm, can someone help me to realize it?

For the modified question, let me try to give an answer that can address situations with matrices having the structure that you have. 
Basically, your matrix $G$ has the following structure $$G=uv^T$$ where I have taken $u$ as the all $1$'s vector and $v$ a vector of positive coordinates such that $v^Tu=1$, i.e. $G$ is a stochastic matrix. It is also easy to see that $G$ is idempotent i.e. $G^2=G$. 
You need to find out $Q$ which is a stochastic such that $\|Q^2-G\|_F$ is minimized. Obviously the only solution is $Q=G$. But if you do not want to take $Q=G$, then you cannot find an optimal solution. Instead you can choose a matrix $\delta$ having the property $\delta u=0$ and then can form the matrix $Q=G+\delta$, then you need to ensure that, for a chosen $\epsilon>0$, $$\|Q^2-G\|_F<\epsilon\implies \|G\delta+\delta G+\delta^2\|_F<\epsilon\\\implies \|uv^T\delta+\delta^2\|_F<\epsilon$$ since $\delta u=0$.
Edit: regarding how to calculate the matrix $\delta$, I am not sure it is always possible to get closed form solutions of $\delta$ and you have to use some numerical technique. Basically, for a given $\epsilon$, you have to calculate the Frobenius norm which will yield,
$$Tr(N\delta^Tvv^T \delta+2\delta^T vu^T\delta^2+\delta^2(\delta^T)^2)<\epsilon$$ As I mentioned earlier, in general, this will ensure that the elements of $\delta$ remain inside a hyperellipse of degree $4$. 
