Consider the linear matrix equation where $A\in {\mathbb{R}}^{m\times p}$, $B\in {\mathbb{R}}^{q\times n}$ and $C\in {\mathbb{R}}^{m\times n}$. Such problems arise in linear control and filtering theory for continuous or discrete-time large-scale dynamical systems. If $AXB=C$ is consistent, $X^{*}=A^{†}C{B}^{†}$ is the minimum Frobenius norm solution. If $AXB=C$ is inconsistent, $X^{*}=A^{†}C{B}^{†}$ is the minimum Frobenius norm least squares solution. When the matrices A and B are small and dense, direct methods based on QR-fraction are attractive [1,2]. However, for large matrices A and B, iterative methods have attracted much attention [3,4,5,6,7]. Recently, Du et al. proposed the randomized block coordinate descent (RBCD) method for solving the matrix least-squares problem $\underset{X\in {\mathbb{R}}^{p\times q}}{min}{\parallel C-AXB\parallel }_F$ without strong convexity assumption in [8]. This method requires that matrix B is full row rank. Wu et al. [9] introduced two kinds of Kaczmarz-type methods to solve consistent matrix equation $AXB=C$: relaxed greedy randomized Kaczmarz (ME-RGRK) and maximal weighted residual Kaczmarz (ME-MWRK). Although the row and column index selection strategy is time-consuming, the ideas of these two methods are suitable for solving large-scale consistent matrix equations.

In this paper, randomized Kaczmarz method [10] and randomized extended Kaczmarz method [11] are used to solve consistent and inconsistent matrix equations (1) by the product of matrix and vector.

All the results in this paper hold in the complex field. But for the sake of simplicity, we only discuss it in the real number field.

In this paper, we denote $A^T$, $A^{†}$, $r\left(A\right)$, $R\left(A\right)$, ${\parallel A\parallel }_F=\sqrt{\mathrm{trace}\left(A^{T}A\right)}$ and ${\langle A,B\rangle }_F=\mathrm{trace}\left(A^{T}B\right)$ as the transpose, the Moore-Penrose generalized inverse, the rank of A, the column space of A, the Frobenius norm of A and the inner product of two matrices A and B, respectively. For an integer $n\ge 1$, let $\left[n\right]=\{1,2,\cdots ,n\}$. We use I to denote the identity matrix whose order is clear from the context. In addition, for a given matrix $G=\left(g_{ij}\right)\in {\mathbb{R}}^{m\times n}$, $G_{i,:}$, $G_{:,j}$, ${\sigma }_{max}\left(G\right)$ and ${\sigma }_{min}\left(G\right)$, are used to denote ith row, jth column, the maximum singular value and the smallest nonzero singular value of G respectively. Let $E_k$ denote the expected value conditional on the first k iterations, that is, ${\mathbb{E}}_k[·]=\mathbb{E}[·|i_0,j_0,i_1,j_1,...,i_{k-1},j_{k-1}],$ where $i_s$ and $j_s(s=0,1,...,k-1)$ are the row and the column chosen at the sth iteration. Let the conditional expectations with respect to the random row index be ${\mathbb{E}}_k^i[·]=\mathbb{E}[·|i_0,j_0,i_1,j_1,...,i_{k-1},j_{k-1},j_k]$ and with respect to the random column index be ${\mathbb{E}}_k^j[·]=\mathbb{E}[·|i_0,j_0,i_1,j_1,...,i_{k-1},j_{k-1},i_k].$ By the law of total expectation, it holds that ${\mathbb{E}}_k[·]={\mathbb{E}}_k^i\left[{\mathbb{E}}_k^j[·]\right]$.

The organization of this paper is as follows. In Section 2, we will discuss Kaczmarz method (ME-RBK and ME-PRBK) for finding the minimal F-norm solution ($A^{†}C{B}^{†}$) of the consistent matrix Eq. (1). In Section 3, we discuss the extended Kaczmarz method (IME-REBK) for finding the minimal F-norm least-squares solution of the matrix Eq. (1). In Section 4, some numerical examples are provided to illustrate the effectiveness of our new methods. Finally, some brief concluding remarks are described in Section 5.

At the kth iteration, Kaczmarz method selects randomized a row $i\in \left[m\right]$ of A and does an orthogonal projection of the current estimate matrix $X^{\left(k\right)}$ onto the corresponding hyperplane $H_i=\{X\in {\mathbb{R}}^{p\times q}|A_{i,:}XB=C_{i,:}\}$, that is

The Lagrangian function of the conditional optimization problem (2) is where $Y\in {\mathbb{R}}^{1\times n}$ is Lagrangian multiplier. By the matrix differentiation, we get the gradient of $L(X,Y)$ and set $\nabla L(X,Y)=0$ for finding the stationary matrix:

By the first equation of (4), we have $X^{(k+1)}=X^{\left(k\right)}-A_{i,:}^{T}Y{B}^T$. Substituting this into the second equation of (4), we can get $Y=-\frac{1}{\parallel A_{i,:}{\parallel }_2^2}(C_{i,:}-A_{i,:}{X}^{\left(k\right)}){\left(B^{T}B\right)}^{†}$. So the projected randomized block Kacmarz (ME-PRBK) for solving $AXB=C$ iterates as

However, in practice, it is very expensive to calculate the pseudoinverse of large-scale matrices. Next, we generalize the average block Kaczmarz method [12] for solving linear equations to matrix equations.

At the kth step, we get the approximate solution $X^{(k+1)}$ by projecting the current estimate $X^{\left(k\right)}$ onto the hyperplane $H_{i,j}=\{X\in {\mathbb{R}}^{p\times q}|A_{i,:}{X}^{\left(k\right)}{B}_{:,j}=C_{i,j}\}$. Using the Lagrangian multiplier method, we can obtain the following Kaczmarz method for $AXB=C$:

Inspired by the idea of average block Kaczmaz algorithm for $Ax=b$, we consider the average block Kaczmaz method for $AXB=C$ respect to B. where $\lambda$ is stepsize and $v_j^k$ be the weights that satisfy $v_j^k\ge 0$ and $\sum _{j=1}^n{v}_j^k=1$. If $v_j^k=\frac{\parallel B_{:,j}{\parallel }_2^2}{{\parallel B\parallel }_F^2}$, then

Set $\alpha =\frac{\lambda }{{\parallel B\parallel }_F^2}>0$, we obtain the following randomized block Kaczmarz iteration where i is selected with probability $p_i=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$. The cost of each iteration of this method is $4q(n+p)+p+1-2q$ if the square of the row norm of A has been calculated in advance. We describe this method as Algorithm 1, which is called the ME-RBK algorithm.

First, we give the following lemma whose proof can be found in [8].

In the following theorem, with the idea of the RK method [10], we will prove that Algorithm 1 converges to the the least F-norm solution of $AXB=C$.

In [11,15,16], the authors proved that the Kaczmarz method does not converge to the least squares solution of $AX=b$ when $AX=b$ is inconsistent. Analogically, if the matrix Eq. (1) is inconsistent, the above ME-PRBK method dose not converge to $A^{†}C{B}^{†}$. The following theorem give the error bound of the for the inconsistent matrix equation.

Next, we use the idea of randomized extended Kaczmarz method (see [16,17,18] for details) for solving the least squares solution of the inconsistent Eq. (1). At each iteration, $Z^{\left(k\right)}$ is the kth iterate of ME-RBK applied to $A^{T}Z{B}^T=0$ with the initial guess $Z^{\left(0\right)}$, and $X^{\left(k\right)}$ is one-step ME-RBK update for $AXB=C-Z^{\left(k\right)}$. we can get the following randomized extended block Kaczmarz iteration where $\alpha >0$ is the step size, and i and j are selected with probability $p_i=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$ and ${\widehat{p}}_j\left(A\right)=\frac{\parallel A_{:,j}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$, respectively. The cost of each iteration of this method is $4n(q+m)+m+1-2n-q$ for updating $Z^{\left(k\right)}$ and $(4q+1)(n+p)+1-2q$ for updating $X^{\left(k\right)}$ if the square of the row norm and the column norm of A have been calculated in advance. We describe this method as Algorithm 2, which is called the ME-REBK algorithm.