In this paper, an iterative method is proposed for solving matrix equation Ps j=1 AjXjBj = E. This method is based on the global least
squares (GL-LSQR) method for solving the linear system of equations with
the multiple right hand sides. For applying the GL-LSQR algorithm to
solve the above matrix equation, a new linear operator, its adjoint and a
new inner product are defined. It is proved that the new iterative method
obtains the least norm solution of the mentioned matrix equation within finite iteration steps in the exact arithmetic, when the above matrix equation
is consistent. Moreover, the optimal approximate solution (X1 ?, X2 ?, . . . , Xs ?)
to a given multiple matrices (X ¯1, X ¯2, . . . , X ¯s) can be derived by finding the
least norm solution of a new matrix equation. Finally, some numerical experiments are given to illustrate the efficiency of the new method