Abstract: In order to solve the least square solution of the matrix equation AXB + CXD = F and its best approximate solution, a multi-step iterative algorithm was designed. The minimum Frobenius norm solution of the matrix sequence generated by the multi-step iterative algorithm converges to the least squares problem of AXB + CXD = F is proved. By modifying the coefficient matrix F, the optimal approximation solution of the least square problem of the matrix sequence generated by the multi-step iterative algorithm converges to the matrix equation AXB + CXD = F is proved. At the same time, the numerical comparison of multi-step iterative algorithm with fixed point algorithm and conjugate gradient algorithm is given. The experimental results show that the multi-step iterative algorithm is effective.