On the Convergence of the Randomized Block Kaczmarz Algorithm for Solving a Matrix Equation

Abstract

A randomized block Kaczmarz method and a randomized extended block Kaczmarz method are proposed for solving the matrix equation $AXB=C$, where the matrices A and B may be full-rank or rank-deficient. These methods are iterative methods without matrix multiplication, and are especially suitable for solving large-scale matrix equations. It is theoretically proved that these methods converge to the solution or least-square solution of the matrix equation. The numerical results show that these methods are more efficient than the existing algorithms for high-dimensional matrix equations.

1. Introduction

Consider the linear matrix equation

$$AXB=C,$$

(1)

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 many practical applications such as surface fitting in computer-aided geometric design (CAGD), signal and image processing, photogrammetry, etc.; see, for example, [1,2,3,4] and the large body of literature therein. 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 fractions are attractive [5,6]. However, for large A and B matrices, iterative methods have attracted a lot of attention [7,8,9,10,11]. 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 [12]. This method requires that matrix B is a full row-rank matrix. Wu et al. [13] introduced two kinds of Kaczmarz-type methods to solve the 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, the randomized Kaczmarz method [14] and the randomized extended Kaczmarz method [15] are used to solve consistent and inconsistent matrix equation (1) with the product of the matrix and vector.

All the results in this paper hold in the complex field. But for the sake of simplicity, we only discuss them in terms of 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,\dots ,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 the ith row, the jth column, the maximum singular value and the smallest nonzero singular value of G, respectively. Let ${\mathbb{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,\dots ,i_{k-1},j_{k-1}],$ where $i_s$ and $j_s(s=0,1,\dots ,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,\dots ,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\dots ,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 the block Kaczmarz method (ME-RBK) for finding the minimal F-norm solution ($A^{†}C{B}^{†}$) of consistent matrix equation (1). In Section 3, we discuss the extended block Kaczmarz method (IME-REBK) for finding the minimal F-norm least-squares solution of matrix equation (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.

2. The Randomized Block Kaczmarz Method for Consistent Equation

At the kth iteration, the Kaczmarz method selects randomized a row $i\in \left[m\right]$ of A and performs 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,

$$\underset{X\in {\mathbb{R}}^{p\times q}}{min}\frac{1}{2}{\parallel X-X^{\left(k\right)}\parallel }_F^{2}s.t.A_{i,:}XB=C_{i,:}.$$

(2)

The Lagrangian function of the conditional optimization problem (2) is

$$L(X,Y)=\frac{1}{2}{\parallel X-X^{\left(k\right)}\parallel }_F^2+\langle Y,A_{i,:}XB-C_{i,:}\rangle ,$$

(3)

where $Y\in {\mathbb{R}}^{1\times n}$ is a Lagrangian multiplier. Via the matrix differentiation, we obtain the gradient of $L(X,Y)$ and set $\nabla L(X,Y)=0$ to find the stationary matrix:

$$\left\{\begin{array}{c}{\nabla }_{X}L(X,Y){|}_{X^{(k+1)}}=X^{(k+1)}-X^{\left(k\right)}+A_{i,:}^{T}Y{B}^T=0,\hfill \\ \\ {\nabla }_{Y}L(X,Y){|}_{X^{(k+1)}}=A_{i,:}{X}^{(k+1)}B-C_{i,:}=0.\hfill \end{array}\right.$$

(4)

Using 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 obtain $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

$$X^{(k+1)}=X^{\left(k\right)}+\frac{A_{i,:}^T}{\parallel A_{i,:}{\parallel }_2^2}(C_{i,:}-A_{i,:}{X}^{\left(k\right)}B)B^{†}.$$

(5)

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

At the kth step, we obtain 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$:

$$X^{(k+1)}=X^{\left(k\right)}+\frac{A_{i,:}^T(C_{i,:}-A_{i,:}{X}^{\left(k\right)}{B}_{:,j})B_{:,j}^T}{{\parallel A_{i,:}\parallel }_2^2{\parallel B_{:,j}\parallel }_2^2}.$$

Inspired by the idea of the average block Kaczmaz algorithm for $Ax=b$, we consider the average block Kaczmaz method for $AXB=C$ with respect to B.

$$X^{(k+1)}=X^{\left(k\right)}+\lambda \frac{A_{i,:}^T}{\parallel A_{i,:}{\parallel }_2^2}\sum _{j=1}^n{v}_j^k\frac{(C_{i,:}-A_{i,:}{X}^{\left(k\right)}{B}_{:,j})B_{:,j}^T}{\parallel B_{:,j}{\parallel }_2^2},$$

where $\lambda$ is stepsize and $v_j^k$ are the weights that satisfy $v_j^k\ge 0$ and ${\displaystyle \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

$$X^{(k+1)}=X^{\left(k\right)}+\frac{\lambda }{{\parallel B\parallel }_F^2}\frac{A_{i,:}^T}{\parallel A_{i,:}{\parallel }_2^2}(C_{i,:}-A_{i,:}{X}^{\left(k\right)}B)B^T.$$

Setting $\alpha =\frac{\lambda }{{\parallel B\parallel }_F^2}>0$, we obtain the following randomized block Kaczmarz iteration:

$$X^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T\left((C_{i,:}-\left(A_{i,:}{X}^{\left(k\right)}\right)B)B^T\right),k=0,1,2,\dots ,$$

(6)

where i is selected with probability $p_i=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$. We describe this method as Algorithm 1, which is called the ME-RBK algorithm.

Algorithm 1 Randomized Block Kaczmarz Method for $AXB=C$ (ME-RBK)

Input:  $A\in R^{m\times p}$, $B\in R^{q\times n}$, $C\in R^{m\times n}$, $X^{\left(0\right)}=0\in R^{p\times q}$ 1: for $k=0,1,2,\dots ,$do 2:    Pick i with probability $p_i\left(A\right)=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$ 3:    Compute $X^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T\left((C_{i,:}-\left(A_{i,:}{X}^{\left(k\right)}\right)B)B^T\right)$ 4: end for

We arrange the computational process of calculating $X^{(k+1)}$ in

Table 1. The complexities of computing $X^{(k+1)}$ in ME-RBK.

$y_1=A_{i,:}{X}^{\left(k\right)}$	$y_2=C_{i,:}-y_{1}B$	$y_3=y_2{B}^T$	$y_4^T=\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T$	$Y_1=y_4^T{y}_3$	$X^{\left(k\right)}+Y_1$
$(2p-1)q$	$(2q-1)n+n$	$(2n-1)q$	$1+p$	$pq$	$pq$

, which only costs $4q(n+p)+p+1-2q$ flopping operations (flops) if the square of the row norm of A has been calculated in advance.

Remark 1. 

Note that the problem of finding a solution of $AXB=C$ can be posed as the following linear least-squares problem:

$$\underset{X\in {\mathbb{R}}^{p\times q}}{min}\frac{1}{2}{\parallel AXB-C\parallel }_F^2=\underset{X\in {\mathbb{R}}^{p\times q}}{min}\frac{1}{2}\sum _{i=1}^m{\parallel A_{i,:}XB-C_{i,:}\parallel }_2^2.$$

(7)

Define the component function

$$f_i\left(X\right)=\frac{1}{2}{\parallel A_{i,:}XB-C_{i,:}\parallel }_2^2,$$

then differentiate with X to obtain its gradient

$$\nabla f_i\left(X\right)=A_{i,:}^T(A_{i,:}XB-C_{i,:})B^T.$$

Therefore, the randomized block Kaczmarz method (6) is equivalent to one step of the stochastic gradient descent method [17] applied to (7) with stepsize $\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}$.

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

Lemma 1 

([12]). Let $A\in {\mathbb{R}}^{m\times p}$ and $B\in {\mathbb{R}}^{q\times n}$ be any nonzero matrix. Let

$$\mathcal{M}=\{M\in R^{p\times q}|\exists Y\in {\mathbb{R}}^{m\times n}s.t.M=A^{T}Y{B}^T\}.$$

For any matrix $M\in \mathcal{M}$, it holds that

$${\parallel AMB\parallel }_F^2\ge {\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right){\parallel M\parallel }_F^2.$$

Remark 2. 

$M\in \mathcal{M}$ means that $M_{:,j}\in R\left(A^T\right),j=1,2,\dots ,q$ and ${\left(M_{i,:}\right)}^T\in R\left(B\right),$$i=1,2,\dots ,p$. In fact, $\mathcal{M}$ is well defined because $0\in \mathcal{M}$ and $A^{†}C{B}^{†}\in \mathcal{M}$.

In the following theorem, with the idea of the RK method [14], we will prove that $X^{\left(k\right)}$ generated by Algorithm 1 converges to the least F-norm solution of $AXB=C$.

Theorem 1. 

Assume $0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}$. If matrix equation (1) is consistent, the sequence $\left\{X^{\left(k\right)}\right\}$ generated by the ME-RBK method starting from the initial matrix $X^{\left(0\right)}\in {\mathbb{R}}^{p\times q}$, in which $X_{:,j}^{\left(0\right)}\in R\left(A^T\right),$$j=1,2,\dots ,q$ and ${\left(X_{i,:}^{\left(0\right)}\right)}^T\in R\left(B\right),i=1,2,\dots ,p$, converges linearly to $A^{†}C{B}^{†}$ in mean square form. Moreover, the solution error in expectation for the iteration sequence $X^{\left(k\right)}$ obeys

$$E\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]\le {\rho }^k{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2,$$

(8)

where $\rho =1-\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right)$, and $i\in \left[m\right]$ picked with probability $p_i\left(A\right)=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$.

Proof. 

For $k=0,1,2,\dots$, by (6) and $A_{i,:}{A}^{†}C{B}^{†}B=C_{i,:}$ (consistency), we have

$$\begin{array}{cc}\hfill X^{(k+1)}-A^{†}C{B}^{†}& =X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T(C_{i,:}-A_{i,:}{X}^{\left(k\right)}B)B^T-A^{†}C{B}^{†}\hfill \\ \hfill & =(X^{\left(k\right)}-A^{†}C{B}^{†})-\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T,\hfill \end{array}$$

then

$$\begin{array}{cc}\hfill \parallel X^{(k+1)}-A^{†}C{B}^{†}{\parallel }_F^2& ={∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2+\frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^4}{\parallel A_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\parallel }_F^2\hfill \\ \hfill & -\frac{2\alpha }{\parallel A_{i,:}{\parallel }_2^2}{\langle X^{\left(k\right)}-A^{†}C{B}^{†},A_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\rangle }_F.\hfill \end{array}$$

It follows from

$$\begin{array}{cc}\hfill & \frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^4}{\parallel A_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\parallel }_F^2\hfill \\ \hfill & =\frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^4}\mathrm{trace}\left(B{B}^T{(X^{\left(k\right)}-A^{†}C{B}^{†})}^T{A}_{i,:}^T{A}_{i,:}{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\right)\hfill \\ \hfill & =\frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^2}\mathrm{trace}\left(B{B}^T{(X^{\left(k\right)}-A^{†}C{B}^{†})}^T{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\right)\hfill \\ \hfill & =\frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^2}\parallel A_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T{\parallel }_2^2(\mathrm{by}\mathrm{trace}\left(u{u}^T\right)={\parallel u\parallel }_2^2\mathrm{for}\mathrm{any}\mathrm{vector}u)\hfill \\ \hfill & \le \frac{{\alpha }^2{\parallel B\parallel }_2^2}{\parallel A_{i,:}{\parallel }_2^2}\parallel A_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†}){B\parallel }_2^2(\mathrm{by}\parallel u^T{B}^T{\parallel }_2={\parallel Bu\parallel }_2\le {\parallel B\parallel }_2{\parallel u\parallel }_2),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & \frac{2\alpha }{\parallel A_{i,:}{\parallel }_2^2}{\langle X^{\left(k\right)}-A^{†}C{B}^{†},A_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\rangle }_F\hfill \\ \hfill & =\frac{2\alpha }{\parallel A_{i,:}{\parallel }_2^2}\mathrm{trace}\left(B^T{(X^{\left(k\right)}-A^{†}C{B}^{†})}^T{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B\right)\hfill \\ \hfill & =\frac{2\alpha }{\parallel A_{i,:}{\parallel }_2^2}{\parallel A_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B\parallel }_2^2\hfill \end{array}$$

that

$${\parallel X^{(k+1)}-A^{†}C{B}^{†}\parallel }_F^2\le {∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2-\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{\parallel A_{i,:}{\parallel }_2^2}{\parallel A_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B\parallel }_2^2.$$

By taking the conditional expectation, we have

$$\begin{array}{cc}\hfill {\mathbb{E}}_k\left[{∥X^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]\le & {∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\hfill \\ \hfill & -{\mathbb{E}}_k\left[\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{\parallel A_{i,:}{\parallel }_2^2}{\parallel A_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B\parallel }_2^2\right]\hfill \\ \hfill =& {∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2-\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}{\parallel A(X^{\left(k\right)}-A^{†}C{B}^{†})B\parallel }_F^2.\hfill \end{array}$$

From $X^{\left(0\right)}\in \mathcal{M}$ and $A^{†}C{B}^{†}\in \mathcal{M}$, we have $X^{\left(0\right)}-A^{†}C{B}^{†}\in \mathcal{M}$. Noting $A_{i,:}^T=A^T{I}_{:,i}$, it is easy to show that $X^{\left(k\right)}-A^{†}C{B}^{†}\in \mathcal{M}$ through induction. Then, from Lemma 1 and $0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}$, we can obtain

$$\begin{array}{cc}\hfill {\mathbb{E}}_k\left[{∥X^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]& \le {∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2-\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}\hfill \\ \hfill & ·{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right){∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\hfill \\ \hfill & =\left(1-\frac{2\alpha -{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right)\right){∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2.\hfill \end{array}$$

(9)

Finally, from (9) and induction on the iteration index k, we obtain the estimate (8).    □

Remark 3. 

Using a similar approach to that used in the proof of Theorem 1, we can prove that the iterate $X^{\left(k\right)}$ generated by ME-PRBK (5) satisfies the following estimate:

$$\mathbb{E}\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]\le {\widehat{\rho }}^k{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2,$$

where $\widehat{\rho }=1-\frac{{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right)}{{\parallel A\parallel }_F^2{\parallel B\parallel }_2^2}$. The convergence factor of GRK in [18] is ${\rho }_{GRK}=1-\frac{{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right)}{{\parallel A\parallel }_F^2{\parallel B\parallel }_F^2}$. It is obvious that

$${\rho }_{GRK}>\underset{0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}}{min}\rho =1-\frac{{\sigma }_{min}^2\left(A\right){\sigma }_{min}^2\left(B\right)}{{\parallel A\parallel }_F^2{\parallel B\parallel }_2^2}=\widehat{\rho }$$

and $\rho <{\rho }_{GRK}$ when $\frac{1-\sqrt{1-\frac{{\parallel B\parallel }_2^2}{{\parallel B\parallel }_F^2}}}{{\parallel B\parallel }_2^2}<\alpha <\frac{1+\sqrt{1-\frac{{\parallel B\parallel }_2^2}{{\parallel B\parallel }_F^2}}}{{\parallel B\parallel }_2^2}$. This means that the convergence factor of ME-PRBK is the smallest and the factor of ME-RBK can be smaller than that of GRK when α is properly selected.

3. The Randomized Extended Block Kaczmarz Method for Inconsistent Equation

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

Theorem 2. 

Assume that the consistent equation $AXB=C$ has a solution $X^{*}=A^{†}C{B}^{†}$. Let ${\widehat{X}}^{\left(k\right)}$ denote the kth iterate of the ME-PRBK method applied to the inconsistent equation $AXB=C+W$ for any $W\in {\mathbb{R}}^{m\times n}$ starting from the initial matrix $X^{\left(0\right)}\in R^{p\times q}$, in which $X_{:,j}^{\left(0\right)}\in R\left(A^T\right),j=1,2,\dots ,q$ and ${\left(X_{i,:}^{\left(0\right)}\right)}^T\in R\left(B\right),i=1,2,\dots ,p$. In exact arithmetic, it follows that

$$\mathbb{E}\left[{∥{\widehat{X}}^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]\le {\widehat{\rho }}^k{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2+\frac{1-{\widehat{\rho }}^k}{1-\widehat{\rho }}\frac{\parallel W{B}^{†}{\parallel }_F^2}{{\parallel A\parallel }_F^2},$$

(10)

Proof. 

Set $H_i=\left\{X\right|A_{i}XB=C_i\}$, ${\widehat{H}}_i=\left\{X\right|A_{i}XB=C_i+W_i\}$. Let Y denote the iterate of the PRBK method applied to the consistent equation $AXB=C$ at the kth step, that is,

$$Y={\widehat{X}}^{\left(k\right)}+\frac{A_{i,:}^T}{\parallel A_{i,:}{\parallel }_2^2}(C_{i,:}-A_{i,:}{\widehat{X}}^{\left(k\right)}B)B^{†}.$$

It follows from

$$\begin{array}{cc}\hfill {\langle Y-A^{†}C{B}^{†},{\widehat{X}}^{(k+1)}-A^{†}C{B}^{†}\rangle }_F& ={\langle Y-A^{†}C{B}^{†},\frac{A_{i,:}^T}{\parallel A_{i,:}{\parallel }_2^2}{W}_i{B}^{†}\rangle }_F\hfill \\ \hfill & =\mathrm{trace}\left({\left(B^{†}\right)}^T{W}_i^T\frac{A_{i,:}}{\parallel A_{i,:}{\parallel }_2^2}(Y-A^{†}C{B}^{†})\right)\hfill \\ \hfill & =\mathrm{trace}\left(W_i^T\frac{A_{i,:}Y-A_{i,:}{A}^{†}C{B}^{†}}{\parallel A_{i,:}{\parallel }_2^2}{\left(B^{†}\right)}^T\right)\hfill \\ \hfill & =\mathrm{trace}\left(W_i^T\frac{A_{i,:}Y-A_{i,:}{A}^{†}C{B}^{†}}{\parallel A_{i,:}{\parallel }_2^2}B{\left(B^{T}B\right)}^{†}\right)\hfill \\ \hfill & (\mathrm{by}{\left(B^{†}\right)}^T=B{\left(B^{T}B\right)}^{†})\hfill \\ \hfill & =\mathrm{trace}\left(W_i^T\frac{A_{i,:}YB-A_{i,:}{A}^{†}C{B}^{†}B}{\parallel A_{i,:}{\parallel }_2^2}{\left(B^{T}B\right)}^{†}\right)=0\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {∥{\widehat{X}}^{(k+1)}-Y∥}_F^2=\frac{{∥A_{i,:}^T{W}_i{B}^{†}∥}_F^2}{\parallel A_{i,:}{\parallel }_2^4}=\frac{1}{\parallel A_{i,:}{\parallel }_2^4}\mathrm{trace}\left({\left(B^{†}\right)}^T{W}_i^T{A}_{i,:}{A}_{i,:}^T{W}_i{B}^{†}\right)=\frac{{∥W_i{B}^{†}∥}_2^2}{\parallel A_{i,:}{\parallel }_2^2}\end{array}$$

that

$$\begin{array}{cc}\hfill {∥{\widehat{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2& ={∥Y-A^{†}C{B}^{†}∥}_F^2+{∥{\widehat{X}}^{(k+1)}-Y∥}_F^2\hfill \\ \hfill & ={∥Y-A^{†}C{B}^{†}∥}_F^2+\frac{{∥W_i{B}^{†}∥}_2^2}{\parallel A_{i,:}{\parallel }_2^2}.\hfill \end{array}$$

(11)

By taking the conditional expectation on both sides of (11), we can obtain

$$\begin{array}{cc}\hfill {\mathbb{E}}_k\left[{∥{\widehat{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]& ={\mathbb{E}}_k\left[{∥Y-A^{†}C{B}^{†}∥}_F^2\right]+{\mathbb{E}}_k\left[\frac{{∥W_i{B}^{†}∥}_2^2}{\parallel A_{i,:}{\parallel }_2^2}\right]\hfill \\ \hfill & \le \widehat{\rho }{∥{\widehat{X}}^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2+\frac{{∥W{B}^{†}∥}_F^2}{{\parallel A\parallel }_F^2}\hfill \end{array}$$

The inequality is obtained using Remark 3. Applying this recursive relation iteratively, we have

$$\begin{array}{cc}\hfill \mathbb{E}\left[{∥{\widehat{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]& \le \widehat{\rho }\mathbb{E}\left[{∥{\widehat{X}}^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]+\frac{{∥W{B}^{†}∥}_F^2}{{\parallel A\parallel }_F^2}\hfill \\ \hfill & \le {\widehat{\rho }}^2\mathbb{E}\left[{∥{\widehat{X}}^{(k-1)}-A^{†}C{B}^{†}∥}_F^2\right]+(\widehat{\rho }+1)\frac{{∥W{B}^{†}∥}_F^2}{{\parallel A\parallel }_F^2}\hfill \\ \hfill & \le \dots \hfill \\ \hfill & \le {\widehat{\rho }}^{k+1}{∥{\widehat{X}}^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2+({\widehat{\rho }}^k+\dots +\widehat{\rho }+1)\frac{{∥W{B}^{†}∥}_F^2}{{\parallel A\parallel }_F^2}\hfill \\ \hfill & ={\widehat{\rho }}^{k+1}{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2+\frac{1-{\widehat{\rho }}^{k+1}}{1-\widehat{\rho }}\frac{\parallel W{B}^{†}{\parallel }_F^2}{{\parallel A\parallel }_F^2}.\hfill \end{array}$$

This completes the proof.    □

Next, we use the idea of the randomized extended Kaczmarz method (see [20,21,22] for details) to solve the least-squares solution of the inconsistent Equation (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 the one-step ME-RBK update for $AXB=C-Z^{\left(k\right)}$. We can obtain the following randomized extended block Kaczmarz iteration:

$$\left\{\begin{array}{c}{Z}^{(k+1)}=Z^{\left(k\right)}-\frac{\alpha }{\parallel A_{:,j}{\parallel }_2^2}{A}_{:,j}\left(\left(\left(A_{:,j}^T{Z}^{\left(k\right)}\right)B^T\right)B\right),\hfill \\ \\ X^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T\left((C_{i,:}-Z_{i,:}^{(k+1)}-\left(A_{i,:}{X}^{\left(k\right)}\right)B)B^T\right),\hfill \end{array}\right.$$

(12)

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.

Algorithm 2 Randomized Extended Block Kaczmarz Method for $AXB=C$ (ME-REBK)

Input:  $A\in R^{m\times p}$, $B\in R^{q\times n}$, $C\in R^{m\times n}$, $X^{\left(0\right)}=0\in R^{p\times q}$, $Z^{\left(0\right)}=C$ 1: for $k=0,1,2,\dots ,$do 2:    Pick j with probability ${\widehat{p}}_j\left(A\right)=\frac{\parallel A_{:,j}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$ 3:    Compute $Z^{(k+1)}=Z^{\left(k\right)}-\frac{\alpha }{\parallel A_{:,j}{\parallel }_2^2}{A}_{:,j}\left(\left(\left(A_{:,j}^T{Z}^{\left(k\right)}\right)B^T\right)B\right)$ 4:    Pick i with probability $p_i\left(A\right)=\frac{\parallel A_{i,:}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$ 5:    Compute $X^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T\left((C_{i,:}-Z_{i,:}^{(k+1)}-\left(A_{i,:}{X}^{\left(k\right)}\right)B)B^T\right)$ 6: end for

Theorem 3. 

Assume $0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}$. Let $\left\{Z^{\left(k\right)}\right\}$ denote the kth iteration of ME-RBK applied to $A^{T}Z{B}^T=0$ starting from the initial matrix $Z^{\left(0\right)}\in R^{m\times n}$, in which $Z_{:,j}^{\left(0\right)}\in C_{:,j}+R\left(A\right),$$j=1,2,\dots ,n$ and ${\left(Z_{i,:}^{\left(0\right)}\right)}^T\in {\left(C_{i,:}\right)}^T+R\left(B^T\right),i=1,2,\dots ,m$. Then, $\left\{Z^{\left(k\right)}\right\}$ converges linearly to $C-A{A}^{†}C{B}^{†}B$ in mean square form, and the solution error in expectation for the iteration sequence $X^{\left(k\right)}$ obeys

$$E\left[{∥Z^{\left(k\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\right]\le {\rho }^k{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2,$$

(13)

where the jth column of A is selected with probability ${\widehat{p}}_j\left(A\right)=\frac{\parallel A_{:,j}{\parallel }_2^2}{{\parallel A\parallel }_F^2}$.

Proof. 

In Theorem 1, replacing A with $A^T$, B with $B^T$ and C with 0, we can prove Theorem 3 based on the result of Theorem 1. For the sake of conciseness, we omit the proof process. □

Theorem 4. 

Assume $0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}$. The sequence $\left\{X^{\left(k\right)}\right\}$ is generated using the ME-REBK method for $AXB=C$, starting from the initial matrix $X^{\left(0\right)}\in R^{p\times n}$ and $Z^{\left(0\right)}\in R^{m\times n}$, where $X_{:,j}^{\left(0\right)}\in R\left(A^T\right),$$j=1,2,\dots ,q$, ${\left(X_{i,:}^{\left(0\right)}\right)}^T\in R\left(B\right),i=1,2,\dots ,p$$Z_{:,j}^{\left(0\right)}\in C_{:,j}+R\left(A\right),$$j=1,2,\dots ,n$ and ${\left(Z_{i,:}^{\left(0\right)}\right)}^T\in {\left(C_{i,:}\right)}^T+R\left(B^T\right),i=1,2,\dots ,m$. For any $\epsilon >0$, it holds that

$$\begin{array}{cc}\hfill E\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]\le & \frac{{(1+\epsilon )}^{k+1}-(1+\epsilon )}{{\epsilon }^2}\frac{{\alpha }^2{\parallel B\parallel }_2^2{\rho }^k}{{\parallel A\parallel }_F^2}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\hfill \\ \hfill & +{(1+\epsilon )}^k{\rho }^k{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2,\hfill \end{array}$$

(14)

where $i\in \left[m\right]$, $j\in \left[p\right]$ are picked with probability $p_i\left(A\right)=\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.

Proof. 

Let $X^{\left(k\right)}$ denote the kth iteration of the ME-REBK method for $AXB=C$, and ${\tilde{X}}^{(k+1)}$ be the one-step Kaczmarz update for the matrix equation $AXB=A{A}^{†}C{B}^{†}B$ from $X^{\left(k\right)}$, i.e.,

$${\tilde{X}}^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T(A_{i,:}{A}^{†}C{B}^{†}B-A_{i,:}{X}^{\left(k\right)}B)B^T.$$

We have

$$\begin{array}{cc}\hfill {\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}& =X^{\left(k\right)}-A^{†}C{B}^{†}-\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T\hfill \end{array}$$

and

$$X^{(k+1)}-{\tilde{X}}^{(k+1)}=\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T(C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B)B^T.$$

For any $\epsilon >0$, via triangle inequality and Young’s inequality, we can obtain

$$\begin{array}{cc}\hfill {∥X^{(k+1)}-A^{†}C{B}^{†}∥}_F^2=& {∥(X^{(k+1)}-{\tilde{X}}^{(k+1)})+({\tilde{X}}^{(k+1)}-A^{†}C{B}^{†})∥}_F^2\hfill \\ \hfill \le & {({∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F+{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F)}^2\hfill \\ \hfill \le & {∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2+{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\hfill \\ \hfill & +2{∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F\hfill \\ \hfill \le & (1+\frac{1}{\epsilon }){∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2+(1+\epsilon ){∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2.\hfill \end{array}$$

(15)

By taking the conditional expectation on the both sides of (15), we have

$$\begin{array}{cc}\hfill {\mathbb{E}}_k\left[{∥X^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]\le & (1+\frac{1}{\epsilon }){\mathbb{E}}_k\left[{∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2\right]\hfill \\ \hfill & +(1+\epsilon ){\mathbb{E}}_k\left[{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right].\hfill \end{array}$$

(16)

It follows from

$$\begin{array}{cc}\hfill & {∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2={∥\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T(C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B)B^T∥}_F^2\hfill \\ \hfill & =\frac{{\alpha }^2}{\parallel A_{i,:}{\parallel }_2^2}\mathrm{trace}\left(B{(C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B)}^T(C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B)B^T\right)\hfill \\ \hfill & \le \frac{{\alpha }^2{\parallel B\parallel }_2^2}{\parallel A_{i,:}{\parallel }_2^2}{∥C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B∥}_2^2\hfill \end{array}$$

that

$$\begin{array}{cc}\hfill {\mathbb{E}}_k\left[{∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2\right]& \le {\alpha }^2{\parallel B\parallel }_2^2{\mathbb{E}}_k^j{\mathbb{E}}_k^i\left[\frac{{∥C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B∥}_2^2}{\parallel A_{i,:}{\parallel }_2^2}\right]\hfill \\ \hfill & ={\alpha }^2{\parallel B\parallel }_2^2{\mathbb{E}}_k^j\left[\frac{1}{{\parallel A\parallel }_F^2}\sum _{i=1}^m{∥C_{i,:}-Z_{i,:}^{(k+1)}-A_{i,:}{A}^{†}C{B}^{†}B∥}_2^2\right]\hfill \\ \hfill & =\frac{{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}{\mathbb{E}}_k\left[{∥Z^{(k+1)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\right].\hfill \end{array}$$

By Theorem 3, it yields

$$\begin{array}{cc}\hfill \mathbb{E}\left[{∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2\right]& \le \frac{{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}\mathbb{E}\left[{∥Z^{(k+1)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\right]\hfill \\ \hfill & \le \frac{{\alpha }^2{\parallel B\parallel }_2^2}{{\parallel A\parallel }_F^2}{\rho }^{k+1}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2.\hfill \end{array}$$

(17)

From $X^{\left(0\right)}-A^{†}C{B}^{†}\in \mathcal{M}$, we have $X^{\left(k\right)}-A^{†}C{B}^{†}\in \mathcal{M}$. Then, by using Theorem 1, we can obtain

$$\begin{array}{cc}\hfill {\mathbb{E}}_k[\parallel {\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}{\parallel }_F^2]& ={\mathbb{E}}_k\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}-\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T{A}_{i,:}(X^{\left(k\right)}-A^{†}C{B}^{†})B{B}^T∥}_F^2\right]\hfill \\ \hfill & \le \rho {∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2,\hfill \end{array}$$

then

$$\begin{array}{c}\hfill \mathbb{E}\left[{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]\le \rho \mathbb{E}\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right].\end{array}$$

(18)

Combining (16)–(18) yields

$$\begin{array}{cc}\hfill \mathbb{E}\left[{∥X^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]\le & (1+\frac{1}{\epsilon })\mathbb{E}\left[{∥X^{(k+1)}-{\tilde{X}}^{(k+1)}∥}_F^2\right]\hfill \\ \hfill & +(1+\epsilon )\mathbb{E}\left[{∥{\tilde{X}}^{(k+1)}-A^{†}C{B}^{†}∥}_F^2\right]\hfill \\ \hfill \le & (1+\frac{1}{\epsilon })\frac{{\alpha }^2{\parallel B\parallel }_2^2{\rho }^{k+1}}{{\parallel A\parallel }_F^2}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\hfill \\ \hfill & +(1+\epsilon )\rho \mathbb{E}\left[{∥X^{\left(k\right)}-A^{†}C{B}^{†}∥}_F^2\right]\hfill \\ \hfill \le & (1+\frac{1}{\epsilon })\frac{{\alpha }^2{\parallel B\parallel }_2^2{\rho }^{k+1}}{{\parallel A\parallel }_F^2}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2[1+(1+\epsilon )]\hfill \\ \hfill & +{(1+\epsilon )}^2{\rho }^2\mathbb{E}\left[{∥X^{(k-1)}-A^{†}C{B}^{†}∥}_F^2\right]\hfill \\ \hfill \le & \dots \hfill \\ \hfill \le & (1+\frac{1}{\epsilon })\frac{{\alpha }^2{\parallel B\parallel }_2^2{\rho }^{k+1}}{{\parallel A\parallel }_F^2}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\sum _{i=0}^k{(1+\epsilon )}^i\hfill \\ \hfill & +{(1+\epsilon )}^{k+1}{\rho }^{k+1}{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2\hfill \\ \hfill =& \frac{{(1+\epsilon )}^{k+2}-(1+\epsilon )}{{\epsilon }^2}\frac{{\alpha }^2{\parallel B\parallel }_2^2{\rho }^{k+1}}{{\parallel A\parallel }_F^2}{∥Z^{\left(0\right)}-(C-A{A}^{†}C{B}^{†}B)∥}_F^2\hfill \\ \hfill & +{(1+\epsilon )}^{k+1}{\rho }^{k+1}{∥X^{\left(0\right)}-A^{†}C{B}^{†}∥}_F^2.\hfill \end{array}$$

This completes the proof. □

Remark 4. 

Replacing $B^T$ in (12) with $B^{†}$, we obtain the following projection-based randomized extended block Kaczmarz mathod (ME-PREBK) iteration:

$$\left\{\begin{array}{c}{Z}^{(k+1)}=Z^{\left(k\right)}-\frac{\alpha }{\parallel A_{:,j}{\parallel }_2^2}{A}_{:,j}\left(\left(\left(A_{:,j}^T{Z}^{\left(k\right)}\right)B^T\right){\left(B^{†}\right)}^T\right),\hfill \\ \\ X^{(k+1)}=X^{\left(k\right)}+\frac{\alpha }{\parallel A_{i,:}{\parallel }_2^2}{A}_{i,:}^T\left((C_{i,:}-Z_{i,:}^{(k+1)}-\left(A_{i,:}{X}^{\left(k\right)}\right)B)B^{†}\right),\hfill \end{array}\right.$$

(19)

4. Numerical Experiments

In this section, we will present some experimental results of the proposed algorithms for solving various matrix equations, and compare them with ME-RGRK and ME-MWRK in [13] for consistent matrix equations and RBCD in [12] for inconsistent matrix equations. All experiments were carried out using MATLAB (version R2020a) on a DESKTOP-8CBRR86 with Intel(R) Core(TM) i7-4712MQ CPU @2.30GHz 2.29GHz, RAM 8GB and Windows 10.

All computations start from the initial guess $X^{\left(0\right)}=0$, and are terminated once the relative error (RE) of the solution, defined by

$$RE=\frac{\parallel X^{\left(k\right)}-X^{*}{\parallel }_F^2}{\parallel X^{*}{\parallel }_F^2}$$

at the current iteration $X^{\left(k\right)}$, satisfies $RE<{10}^{-6}$ or exceeds the maximum iteration K = 50,000, where $X^{*}=A^{†}C{B}^{†}$. We report the average number of iterations (denoted as “IT”) and the average computing time in seconds (denoted as“CPU”) for 20 repeated trial runs of the corresponding method. Three examples are tested, and A and B are generated as follows.

• Type I: For given $m,p,q,n$, the entries of A and B are generated from a standard normal distribution, i.e., $A=randn(m,p),B=randn(q,n).$
• Type II: Like [18], for given $m,p$, and $r_1=rank\left(A\right)$, we construct a matrix A by $A=U_1{D}_1{V}_1^T$, where $U_1\in R^{m\times r_1}$ and $V_1\in R^{p\times r_1}$ are orthogonal column matrices, $D\in R^{r_1\times r_1}$ is a diagonal matrix whose first $r-2$ diagonal entries are uniformly distributed numbers in $[{\sigma }_{min}\left(A\right),{\sigma }_{max}\left(A\right)]$, and the last two diagonal entries are ${\sigma }_{max}\left(A\right),{\sigma }_{min}\left(A\right)$. The entries of B are generated using a similar method with parameters $q,n,r_2=rank\left(B\right)$.
• Type III: The real-world sparse data come from the Florida sparse matrix collection [23]. Table 2 lists the features of these sparse matrices.

Table 2. The detailed features of sparse matrices from [23].

Name	Size	Rank	Sparsity
ash219	$219\times 85$	85	$2.3529\%$
ash958	$958\times 292$	292	$0.68493\%$
divorce	$50\times 9$	9	$50\%$

Table 2. The detailed features of sparse matrices from [23].

Name	Size	Rank	Sparsity
ash219	$219\times 85$	85	$2.3529\%$
ash958	$958\times 292$	292	$0.68493\%$
divorce	$50\times 9$	9	$50\%$

4.1. Consistent Matrix Equation

Given $A,B$, we set $C=A{X}^{*}B$ with $X^{*}=randn(p,q)$ to construct a consistent matrix equation. First, we test the impact of $\alpha$ in the ME-RBK method on the experimental results. Figure 1 plots the IT and CPU versus different “para” with different matrices in

Table 3. IT and CPU of ME-RGRK, ME-MWRK, ME-RBK and ME-PRBK for the consistent matrix equations with Type I.

No.	m	p	q	n		ME-RGRK	ME-MWRK	ME-RBK	ME-PRBK
1	100	40	40	100	IT	49,707	27,579	7834.5	1152.8
					CPU	0.78	2.01	0.20	0.04
2	40	100	100	40	IT	>	49,332.7	6334.7	1507.2
					CPU	>	2.15	0.33	0.08
3	500	100	100	500	IT	>	32,109	4021.8	1866.1
					CPU	>	57.61	0.52	0.19
4	1000	200	300	2000	IT	>	>	6429.6	4450.4
					CPU	>	>	3.95	0.72

, where $para=$ 0.1:0.1:1.9 so that $\alpha =\frac{para}{{\parallel B\parallel }_2^2}$ satisfies $0<\alpha <\frac{2}{{\parallel B\parallel }_2^2}$ in Theorem 1. From Figure 1, it can be seen that the number of iteration steps and the running time decrease with the increase in parameters. However, when $para=1.9$, both IT and CPU begin to increase. The same situations occur when solving consistent or inconsistent equations with different matrices in

Table 4. IT and CPU of ME-RGRK, ME-MWRK, ME-RBK and ME-PRBK for the consistent matrix equations with Type II.

m	p	${\mathit{r}}_1$	$[{\mathit{\sigma }}_{min}\left(\mathit{A}\right),{\mathit{\sigma }}_{max}\left(\mathit{A}\right)]$	q	n	${\mathit{r}}_2$	$[{\mathit{\sigma }}_{min}\left(\mathit{B}\right),{\mathit{\sigma }}_{max}\left(\mathit{B}\right)]$		ME-RGRK	ME-MWRK	ME-RBK	ME-PRBK
100	40	20	[1, 2]	40	100	40	[1, 2]	IT	4361.0	1987.1	503.1	332.4
								CPU	0.06	0.17	0.01	0.008
100	40	20	[1, 5]	40	100	20	[1, 5]	IT	22,423.2	6439.8	9307.6	1056.3
								CPU	0.33	0.57	0.20	0.02
1000	200	100	[1, 2]	100	1000	50	[1, 2]	IT	20,055.5	7047.8	2587.4	1674.3
								CPU	78.45	70.56	0.42	0.23
1000	100	50	[1, 5]	200	1000	200	[1, 5]	IT	>	>	18,898.3	2833.6
								CPU	>	>	3.61	0.42

and

Table 5. IT and CPU of ME-RGRK, ME-MWRK, ME-RBK and ME-PRBK for the consistent matrix equations with Type III.

A	B		ME-RGRK	ME-MWRK	ME-RBK	ME-PRBK
divorce	ash219${}^T$	IT	43,927.8	14,164.4	10,993.4	3873.5
		CPU	1.15	1.35	0.36	0.13
divorce	ash219	IT	40,198.7	17,251.4	11,557.4	3124.7
		CPU	0.63	0.80	0.43	0.11
ash219	ash958${}^T$	IT	>	>	6042.3	2267.0
		CPU	>	>	1.04	0.36
ash219	ash958	IT	>	>	5745.4	2114.2
		CPU	>	>	1.22	0.42

. Therefore, we set $\alpha =\frac{1.8}{{\parallel B\parallel }_2^2}$ in all experiments.

In Table 3, Table 4 and Table 5, we report the average IT and CPU of the ME-RGRK, ME-MWRK, ME-RBK and ME-PRBK methods for solving consistent eqautions. In the following tables, the item “>” represents that the number of iteration steps exceeds the maximum iteration (50,000), and the item “-” represents that the method does not converge.

From these tables, we can see that the ME-RBK and ME-PRBK methods vastly outperform the ME-RGRK and ME-MWRK methods in terms of both IT and CPU times regardless of whether the matrices A and B are full column/row rank or not. As the matrix dimension increases, the CPU time of the ME-RBK and ME-PRBK methods increases slowly, while the running time of ME-RGRK and ME-MWRK increases dramatically.

In addition, when the matrix size is small, the ME-PRBK method is competitive, because the pseudoinverse is less expensive and the number of iteration steps is small. When the matrix size is large, the matrix is large, and the ME-RBK method is more challenging because it does not need to calculate the pseudoinverse (see the last line in Table 3).

4.2. Inconsistent Matrix Equation

To construct an inconsistent matrix equation, we set $C=A{X}^{*}B+R$, where $X^{*}$ and R are random matrices which are generated by $X^{*}=randn(p,q)$ and $R=\delta *randn(p,q),$$\delta \in (0,1)$. Numerical results of the RBCD, IME-REBK and IME-PREBK methods are listed in

Table 6. IT and CPU of RBCD, IME-REBK and IME-PREBK for the inconsistent matrix equations with Type I.

m	p	q	n		RBCD	IME-REBK	IME-PREBK
100	40	40	100	IT	17,212.2	21,270.5	2469.4
				CPU	0.69	1.57	0.19
40	100	100	40	IT	-	24,708.3	2174.6
				CPU	-	2.38	0.21
500	100	100	500	IT	6059.1	7352.8	2512.3
				CPU	4.97	7.28	2.70
1000	200	300	2000	IT	14,209.2	12,490.4	5183.5
				CPU	152.32	246.56	99.68

,

Table 7. IT and CPU of RBCD, IME-REBK and IME-PREBK for the inconsistent matrix equations with Type II.

m	p	${\mathit{r}}_1$	$[{\mathit{\sigma }}_{min}\left(\mathit{A}\right),{\mathit{\sigma }}_{max}\left(\mathit{A}\right)]$	q	n	${\mathit{r}}_2$	$[{\mathit{\sigma }}_{min}\left(\mathit{B}\right),{\mathit{\sigma }}_{max}\left(\mathit{B}\right)]$		RBCD	IME-REBK	IME-PREBK
100	40	20	[1, 2]	40	100	40	[1, 2]	IT	1035.9	762.2	384.5
								CPU	0.04	0.05	0.03
100	40	20	[1, 5]	40	100	20	[1, 5]	IT	-	16,224.7	1507.8
								CPU	-	1.23	0.12
1000	200	100	[1, 2]	100	1000	50	[1, 2]	IT	-	4067.5	2217.8
								CPU	-	24.23	13.67
1000	200	100	[1, 5]	100	1000	100	[1, 5]	IT	>	48,269.6	4328.0
								CPU	>	365.78	30.25

and

Table 8. IT and CPU of RBCD, IME-REBK and IME-PREBK for the inconsistent matrix equations with Type III.

A	B		RBCD	IME-REBK	IME-PREBK
divorce	ash219${}^T$	IT	>	20,026.3	4308.4
		CPU	>	2.26	0.48
divorce	ash219	IT	-	19,199.1	4026.2
		CPU	-	1.49	0.31
ash219	ash958${}^T$	IT	22,313.4	10,823.5	2561.8
		CPU	15.69	19.03	5.89
ash219	ash958	IT	-	10,020.7	2363.5
		CPU	-	15.83	4.72

. From these tables, we can see that the IME-PREBK method is better than the RBCD method in terms of IT and CPU time, especially when the $\frac{{\sigma }_{max}}{{\sigma }_{min}}$ is large (see the last line in Table 7). The IME-REBK method is not competitive for B with full row rank because it needs to solve two equations. However, when B does not have full row rank, the RBCD method does not converge, while the IME-REBK and IME-PREBK methods do.

5. Conclusions

In this paper, we have proposed a randomized block Kaczmarz algorithm for solving the consistent matrix equation and its extended version for the inconsistent case. Theoretically, we have proved that the proposed algorithms converge linearly to the unique minimal F-norm solution or least-squares solution (i.e., $A^{†}C{B}^{†}$) without requirements on A and B having full column/row rank. The numerical results show the effectiveness of the algorithms. Since the proposed algorithms only require one row or one column of A at each iteration without a matrix–matrix product, they are suitable for the scenarios where the matrix A is too large to fit in the memory or matrix multiplication is considerably expensive.