In this paper, we discuss least squares symmetrizable solutions of matrix equations (AX = B, XC = D) and its optimal approximation solution. With the matrix row stacking, Kronecker product and special relations between two linear subspaces are topological isomorphism, and we derive the general solutions of least squares problem. With the invariance of the Frobenius norm under orthogonal transformations, we obtain the unique solution of optimal approximation problem. In addition, we present an algorithm and numerical experiment to obtain the optimal approximation solution.