Preconditioning Kaczmarz method by sketching

TL;DR

This paper introduces a novel preconditioning technique for the Kaczmarz method using sketching, which aims to improve convergence speed while reducing computational complexity through random sketching and QR decomposition.

Contribution

The paper presents a new sketching-based preconditioning approach for the Kaczmarz method, offering a computationally efficient alternative to QR-based preconditioning.

Findings

Sketching-based preconditioning improves convergence speed.

The method is effective on both random and real systems.

Numerical experiments validate the approach.

Abstract

We propose a new method for preconditioning Kaczmarz method by sketching. Kaczmarz method is a stochastic method for solving overdetermined linear systems based on a sampling of matrix rows. The standard approach to speed up convergence of iterative methods is using preconditioner. As we show the best possible preconditioner for this method can be constructed from QR decomposition of the system matrix, but the complexity of this procedure is too high. Therefore, to reduce this complexity, we use random sketching and compare it with the Kaczmarz method without preconditioning. The developed method is applicable for different modifications of classical Kaczmarz method that were proposed recently. We provide numerical experiments to show the performance of the developed methods on solving both random and real overdetermined linear systems.

Equations

$$Ax=b,$$

$$\begin{array}{rl}& x_{k+1}=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert x_k-x{\Vert }_2^2,\\ \text{s.t. }& a_{i}x=b_i,\end{array}$$

$$x_{k+1}=x_k-\frac{a_i{x}_k-b_i}{\Vert a_i{\Vert }_2^2}{a}_i^{\top }.$$

$$\mathbb{E}\Vert x^{\ast }-x_k{\Vert }_2^2\le (1-{\kappa }_F(A{)}^{-2}{)}^k\Vert x^{\ast }-x_0{\Vert }_2^2,$$

$${\kappa }_F(A)\le \sqrt{n}\cdot k(A).$$

$$k(A{P}_R^{\ast })=1,$$

$$A=QR,$$

$$Q=A{R}^{-1},P_R^{\ast }=R^{-1}.$$

$$A{P}_{R}y=b,x=P_{R}y.$$

$$x_{k+1}=x_k-\frac{a_i{P}_R{x}_k-b_i}{\Vert P_R{a}_i{\Vert }_2^2}(a_i{P}_R{)}^{\top },$$

$$s_{ij}=\{\begin{array}{ll}1& j=k_i,\\ 0& \text{otherwise},\end{array}$$

$$\hat{A}=\hat{Q}\hat{R},$$

$$\mu (A)=\underset{i=1,...,m}{\max }\Vert q_i{\Vert }_2^2,$$

$$f(x)=a^{\top }x+c+\alpha \exp \left(-\frac{(\Vert x{\Vert }_2-\mu {)}^2}{{\sigma }_f}\right),$$

Full text

{AMS}

15A06, 15B52, 65F08, 65F20, 68W20

1 Introduction

This paper presents a new method of preconditioning of the Kaczmarz method [9] based on sketching technique [23]. The Kaczmarz method is the stochastic method for solving overdetermined linear systems $Ax=b$, where $A\in\mathbb{R}^{m\times n}$ and $m\gg n$. We assume that we can only sample rows of $A$, and can not use the full $A$. Existing convergence results [22] show that convergence is defined by $\kappa_{F}(A)=\|A\|_{F}\|A^{-1}\|_{2}$. We show by extensive numerical simulations that this value can be minimized by using preconditioner. To reduce the time required for a preconditioner computation we approximate the matrix $A$ by sketching and compute the preconditioner for a sketched matrix. According to our assumptions, we can sketch the matrix only by random sampling of a fixed number of rows. The number of sketched rows is a parameter of the proposed method, and we test different numbers proportional to the number of columns $n$.

Numerical experiments demonstrate the performance of the proposed method for both random consistent linear systems, test tomography problems and non-linear function approximation problem with random Fourier features technique. Also, we study random linear systems with noisy right-hand sides and show how the variance of the noise affects the convergence.

Main contributions of this paper are the following.

• •

We study the preconditioning of the Kaczmarz method by sketching, analyze its complexity and provide limitations of this approach. (Section 2, 3).

• •

We propose to use preconditioning as the fine-tuning technique to get highly accurate solutions of overdetermined linear systems. (Section 4)

• •

We compare the proposed approach on both synthetic data, test tomography data and function approximation problem using explicit feature map. (Section 4)

1.1 Related work

The Kaczmarz method for solving overdetermined linear systems is a well-known randomized iterative method, which has numerous interpretations and generalizations [6, 14]. For example, it can be re-written as a special case of the stochastic gradient descent or coordinate descent applied to the linear least-squares problem corresponding to the considered system [16, 18]. The idea of preconditioning of such kind of methods is mentioned in the paper [24] but in the context of using stochastic gradient method for solving $\ell_{p}$ regression problem. In contrast, we emphasize on the Kaczmarz method and its convergence speed up with preconditioning. The convergence of the Kaczmarz method for a noisy linear system is presented in [13]. Our study shows that in the case of small noise the preconditioned Kaczmarz method gives smaller relative error than without preconditioning. Also, there is an extension [25] of the Kaczmarz method for the inconsistent linear system that uses random orthogonal projection technique to address the inconsistency. One of the main issues in the Kaczmarz method is the correct sampling of rows: if two sequential rows are correlated than the progress in the convergence will be very small. The paper [15] addresses this issue and proposes a two-subspace projection method.

2 Kaczmarz method and its preconditioning

Let $A\in\mathbb{R}^{m\times n}$ be a rectangular matrix, where $m\gg n$ and $b\in\mathbb{R}^{m}$. We consider the problem of solving an overdetermined linear system

[TABLE]

with the following assumptions

the least-squares solution of the problem using a full matrix $A$ is not available 2. 2)

we can sample rows of matrix $A$ according to some probability distribution.

To fulfill these assumptions, we consider the Kaczmarz method [9] for solving (1). This method is based on the iterative update of the current approximation $x_{k}$ using the solution of the following optimization problem:

[TABLE]

where $a_{i}$ is the $i$-th row of the matrix $A$ and $b_{i}$ is the $i$-th element of the vector $b$. Problem (2) has an analytical solution

[TABLE]

The rows of the matrix $A$ can be sampled in different ways. The original paper [9] proposed a sequential selection of rows one by one. Later, the study [22] proposed a Randomized Kaczmarz method (RKM), which samples rows with probabilities $\frac{\|a_{i}\|_{2}^{2}}{\|A\|_{F}^{2}}$. The convergence of the RKM is characterized by the following Theorem.

Theorem 2.1** ([22])**

Let $x^{*}$ be the solution of (1). Then RKM converges in expectation with the average error as

[TABLE]

*where $\kappa_{F}(A)=\|A\|_{F}\|A^{-1}\|_{2}=\frac{\|A\|_{F}}{\sigma_{\min}(A)}$.

From this Theorem, it follows that decreasing $\kappa_{F}(A)$ leads to faster convergence of RKM. Also, let $k(A)=\|A\|_{2}\|A^{-1}\|_{2}$ be a condition number of the matrix $A$. The following well-known inequality between $\kappa_{F}(A)$ and $k(A)$ holds

[TABLE]

Thus, $k(A)$ is an upper bound of $\kappa_{F}(A)$ up to the factor $\sqrt{n}$, and therefore instead of the minimization of $\kappa_{F}(A)$ we can minimize $k(A)$. The standard approach for the reduction of the condition number is to use preconditioners [20]. We will denote by $P_{L}\in\mathbb{R}^{m\times m}$ and by $P_{R}\in\mathbb{R}^{n\times n}$ the left and right preconditioners, respectively. In our case, the matrix $A$ is rectangular of size $m\times n$ and $m\gg n$. Therefore, the right preconditioner is preferable since it requires $\mathcal{O}(n^{2})$ memory in contrast to the left preconditioner, which requires $\mathcal{O}(m^{2})$. An ideal right preconditioner $P_{R}^{*}$ satisfies

[TABLE]

which means that $AP_{R}^{*}$ is unitary. Therefore, we can obtain $P_{R}^{*}$ from the QR decomposition of the matrix $A$. In particular, if

[TABLE]

where $Q\in\mathbb{R}^{m\times n}$ is an orthogonal matrix and $R\in\mathbb{R}^{n\times n}$ is an upper triangular matrix, then

[TABLE]

The complexity of every step in this procedure for deriving optimal right preconditioner $P^{*}_{R}$ is the following:

• •

QR decomposition of the matrix $A\in\mathbb{R}^{m\times n}$ costs $\mathcal{O}(mn^{2})$,

• •

inversion of upper triangular matrix $R$ costs $\mathcal{O}(n^{3})$.

Thus, the total complexity is $\mathcal{O}(mn^{2})+\mathcal{O}(n^{3})$, which is too high in the case of $m\gg n$. Section 3 shows how to reduce the complexity of preconditioner computation down to $\mathcal{O}(n^{3})$.

After introducing the right preconditioner $P_{R}$ we have the following modification of the problem (1):

[TABLE]

The corresponding modification of the Kaczmarz method is

[TABLE]

where $a_{i}$ is the $i$-th row of the matrix $A$ sampled according to some distribution.

Note that the cost of the one iteration of the method (8) is $\mathcal{O}(m^{2}+n)$, which is higher than corresponding cost $\mathcal{O}(n)$ of the original Kaczmarz method (3). Also, we need to perform preprocessing for computing $P_{R}$. However, the expected speedup asymptotically dominates the increased cost of every iteration and gives a smaller relative error for the same time. Examples of such behavior and runtime comparison are provided in Section 4.

Although we show in (5) and (6) how to compute the optimal right preconditioner, this computation is intractable because we can not use the full matrix $A$. To overcome this limitation we propose to use sketching [23, 5, 7] of the matrix $A$ to compute the approximation of $P^{*}_{R}$. The next section discusses some ways to sketch the matrix $A$ and resulting sketching-based preconditioners.

3 Sketching-based preconditioners

This section gives an overview of

sketching idea, how it helps in constructing the preconditioner for the Kaczmarz method and limitations of this approach.

The sketching is the procedure of replacing original matrix $A\in\mathbb{R}^{m\times n}$ by a sketched matrix $\hat{A}=SA\in\mathbb{R}^{r\times n}$, where $S\in\mathbb{R}^{r\times m}$ is some given matrix and

$r\ll m$ [23, 6]. Since the complexity of computing the product $SA$ for arbitrary dense matrices $S$ and $A$ is $\mathcal{O}(rmn)$, which is intractable, we consider matrix $S$ from a set $\mathcal{S}$, such that $S=[s_{ij}]\in\mathcal{S}\subset\mathbb{R}^{r\times n}$ if and only if for $i=1,\ldots,r$,

[TABLE]

where $k_{1}<k_{2}<\ldots<k_{r}$ and $\{k_{1},k_{2},\ldots,k_{r}\}\subset\{1,\ldots,m\}$. For such matrices $S\in\mathcal{S}$, product $SA$ is equivalent to the selection of $r$ rows of the matrix $A$ with indices $\{k_{1},\ldots,k_{r}\}$ such that $k_{1}<k_{2}<\ldots<k_{r}$. Details of rows selection we discuss later in this section. Now assume we have a sketched matrix $\hat{A}$ such that $r\geq n$. Then we compute QR decomposition of this matrix

[TABLE]

where $\hat{Q}\in\mathbb{R}^{r\times n}$ is orthogonal and $\hat{R}\in\mathbb{R}^{n\times n}$ is upper triangular matrix, and use $\hat{P}_{R}=\hat{R}^{-1}$ as an approximation of $P^{*}_{R}$. Comparing with (5) and (6), we reduce the complexity of the QR decomposition from $\mathcal{O}(mn^{2})$ to $\mathcal{O}(rn^{2})$, which can be rewritten in the form $\mathcal{O}(n^{3})$ if $r=\gamma n$, where $\gamma\geq 1$ is some given factor. Therefore, the total complexity of computing preconditioner $\hat{P}_{R}$ is $\mathcal{O}(n^{3})$ and does not depend on the number of rows $m$ of the matrix $A$.

3.1 Selection of rows for sketched matrix $\hat{A}$

The remaining question in the procedure of construction $\hat{P}_{R}$ is how to obtain a sketched matrix $\hat{A}$ or in particular how to select rows in the matrix $A$ in such a way that the matrix $\hat{P}_{R}$ approximates the matrix $P^{*}_{R}$.

Since we can only sample rows of matrix $A$, but can not use full matrix, we use a simple uniform random sampling strategy. Unfortunately, this strategy is not universal and can give a very poor approximation of $P^{*}_{R}$ or even singular matrix $\hat{R}$. This phenomenon relates to the notion of coherence of the matrix.

Definition 3.1** ([3])**

Let $\mu(A)$ be the coherence of matrix $A$:

[TABLE]

*where $q_{i}$ is the $i$-th row of the matrix $Q$ from the QR decomposition of the matrix $A$. *

According to [3], the coherence $\mu(A)$ is between $n/m$ and $1$ and the larger coherence of the matrix, the poorer approximation $\hat{R}$ is given by a uniform random sampled sketched matrix $\hat{A}$. To address this issue, authors of [3] propose several row mixing methods that reduce the coherence of the matrix with unitary transformations, e.g, Discrete Cosine Transform. Unfortunately, such kind of methods require $\mathcal{O}(m\log m)$ operations and are too costly, since in our case $m\gg n$. Thus, the used approach of constructing sketched matrix $\hat{A}$ is applicable only to matrices that have low coherence.

4 Numerical experiments

In this section, we provide numerical experiments to demonstrate the performance of the proposed approach to speed-up the convergence of the Kaczmarz method. The source code can be found at GitHub***https://github.com/amkatrutsa/preckacz. We use random matrices, standard tomography data from AIRTool II [8] and matrices taken from non-linear function approximation problem. We reduce the latter problem to the problem of solving the overdetermined linear system with random Fourier features technique [1]. In all experiments, we initialize $x_{0}$ with zero vector. Considered algorithms are implemented in Python with NumPy library [17]. Since we assume that the full matrix $A$ is not available, we use uniform random sampling of rows in the Kaczmarz method and its preconditioning modifications. Also, we test the proposed approach for a range of factors $\gamma$ required to compute $\hat{P}_{R}$. We take into account the preprocessing time required to compute $\hat{P}_{R}$ as a starting point for all convergence plots shown below.

4.1 Synthetic data

We generate a random matrix of size $m=10000,n=100$ with condition number $k\sim 10^{6}$. Entries of the matrix are generated from the standard normal distribution $\mathcal{N}(0,1)$. After that, we generate a ground truth solution $x^{*}\in\mathbb{R}^{n}$ from $\mathcal{N}(0,I)$. We compare the proposed preconditioning method in two settings: consistent system, where $b=Ax^{*}$ and consistent system with noisy right-hand side: $b=Ax^{*}+\varepsilon$, $\varepsilon\sim\mathcal{N}(0,\sigma^{2}I)$. In both settings, we test different values of $\gamma$. On the one hand, the larger $\gamma$, the better approximation quality of $\hat{P}_{R}$. But on the other hand, the larger $\gamma$, the more time we need to compute $\hat{P}_{R}$ in pre-processing step. Therefore, we identify the minimal value of $\gamma$ among the tested values that give faster convergence than the Kaczmarz method without preconditioning (“No preconditioning” label in figures below).

4.1.1 Consistent system setting

In this experiment setting we generate the right-hand side $b=Ax^{*}$. Figure 1 shows the dependence of the relative error on the convergence time for the compared methods. In the case of the Kaczmarz method, which is labeled as “No preconditioning”, $\hat{P}_{R}=I$. This plot demonstrates that if $\gamma$ is not large enough, e.g. $\gamma=1$, the convergence of the preconditioned Kaczmarz method can be significantly slower than the convergence without preconditioning. In the same time, already for $\gamma=2$ we observe a significant speedup compared with the Kaczmarz method without preconditioning. Further increase of $\gamma$ shows that $\gamma=3$ is already large enough and a further increase is not reasonable.

4.1.2 Consistent system with noise setting

In this setting, we generate the right-hand side $b=Ax+\varepsilon$, where $\varepsilon\in\mathcal{N}(0,\sigma^{2}I)$. We test the dependence of the relative error on the required time for different values of $\gamma$ and $\sigma$. Figure 2 shows the same result with respect to the considered values of $\gamma$ as for the consistent setting, see Figure 1. In particular, $\gamma=1$ is not large enough to give a good approximation of $P^{*}_{R}$, but $\gamma=2$ or larger is already large enough to get faster convergence. Also, this plot shows the dependence of the performance of the considered methods for different noise levels. The higher the level of noise, the less gain from the preconditioned Kaczmarz method we get. That is because preconditioning gives a very small relative error, but in the high noise case, we can not get very accurate $x_{k}$ due to the noisy right-hand side.

Thus, we experimentally show the effectiveness of the preconditioning Kaczmarz method in the cases of randomly generated consistent linear system and randomly generated linear system with a small level of noise in the right-hand side. In the next section, we use matrices coming from different types of tomography to show the performance of the proposed method in the real application.

4.2 Tomography data

To generate the tomography data we use the package AIR Tool II [8]. In particular, to generate matrices $A$ we use functions paralleltomo, fancurvedtomo and fanlineartomo with default parameters and number of discretization intervals $q=50$. The number of intervals $q$ corresponds to a $q\times q$ square image that we reconstruct through solving an overdetermined system. With default parameters and $q=50$ these functions give matrices $12780\times 2500$. Denote by $A_{p},A_{c}$ and $A_{l}$ matrices given by functions paralleltomo, fancurvedtomo and fanlineartomo respectively. As the ground truth solution we use a standard test image shown in Figure 3 and denote it by $X^{*}\in\mathbb{R}^{q\times q}$. This image is represented as $q\times q$ matrix and is reconstructed by the considered methods in the reshaped form $q^{2}\times 1$.

Figure 4a shows the convergence of the same methods that were studied in Section 4.1. Here we again observe that $\gamma=1$ is not enough for faster convergence. Moreover, for the matrix $A_{p}$, $\gamma=1$ gives a singular matrix $\hat{R}$ and we use its pseudoinverse to compute $\hat{P}_{R}$. For larger values of $\gamma$ we get nonsingular matrices $\hat{R}$ and observe asymptotically faster convergence compared with no preconditioning case. However, we see two phases in the convergence of preconditioned methods: the first phase gives less accurate approximation than the Kaczmarz method without preconditioning, but at some moment the Kaczmarz method gives higher relative error than the preconditioned Kaczmarz methods. For example, Figure 4a shows that after 300 seconds the preconditioned Kaczmarz method, where $\gamma=3$, gives lower relative error than the Kaczmarz method without preconditioning. We observe the same convergence in the case of matrices $A_{c}$ and $A_{l}$, see Figure 4b, 4c. To make the convergence plot more clear, we skip the case $\gamma=1$ since it gives non-singular but very ill-conditioned matrix $\hat{R}$.

In addition to the convergence plots in Figure 4, we provide the normalized elementwise error between ground truth image and images reconstructed with the considered methods, see Figure 5. These plots show regions of the reconstructed images with high and low reconstruction quality, which is measured as $\frac{|X^{*}-\hat{X}|}{\|X^{*}\|_{F}}$, where $X^{*}$ is ground truth solution and $\hat{X}$ is approximation given by the corresponding method. The nominator is the elementwise absolute value of the difference $X^{*}$ and $\hat{X}$.

Thus, the proposed preconditioned Kaczmarz method can be used as a fine-tuning technique that will improve the accuracy given by the Kaczmarz method without preconditioning. Details of this approach are in the next section.

4.2.1 Fine-tuning with the preconditioned Kaczmarz method

As we show in the previous section the preconditioned Kaczmarz asymptotically gives a smaller relative error. However, the Kaczmarz method without preconditioning gives a smaller relative error during the first iterations. Therefore, we propose to start solving a given system without using preconditioner and then switch to the preconditioned modification of the Kaczmarz method. This strategy gives a lower error rate for the same runtime, see Figure 6, where we denote by $\tau$ the time to switch to the preconditioned Kaczmarz method. We use $\tau=120$ seconds in all experiments. In the case of tomography given by $A_{p}$ we do not consider fine-tuning for $\gamma=1$, since it gives poor convergence that we see in Figure 4a.

4.3 Kernel regression problem

One more source of the consistent overdetermined linear systems is the kernel regression problem with explicit feature maps, e.g. random Fourier features [2, 12]. These maps reduce the non-linear regression problem to the linear regression problem with the modified matrix. We use this approach in the problem of non-linear function approximation such that the parametric function model is not known. Also, it is known that introducing Gaussian kernel reduces linear regression problem to non-linear and can give a very accurate approximation of smooth functions [11].

Consider the following function $f:\mathbb{R}^{2}\to\mathbb{R}$ that is used as the test function in the package George [1]:

[TABLE]

and use the following parameters: $a=(1,1),c=1,\alpha=0.1,\mu=0,\sigma_{f}=1$. To compose an overdetermined linear system we sample $m=10000$ points uniform in $D=[0,1]\times[0,1]$ and compute the function values in these points. So, we have a matrix $Z\in\mathbb{R}^{m\times 2}$, where every row is a sampled point from the set $D$ and the right-hand side $b$ is the vector of the function values in the corresponding points. Since the test function (10) is non-linear, the composed linear system $Z\tilde{x}=b$ is not consistent, but the implicit kernel trick with random Fourier features helps in constructing a consistent linear system $Ax=b$. This explicit feature map gives an approximation of introducing the Gaussian kernel and the accuracy of the approximation exponentially grows with the dimension $d$ [19]. To construct the matrix $A$, we generate the matrix $M\in\mathbb{R}^{2\times d}$ with entries from the normal distribution $\mathcal{N}(0,\sigma^{2})$, where $d\geq 2$ and $\sigma$ are hyper-parameters that have to be tuned. After that, we compute $B=ZM\in\mathbb{R}^{m\times d}$, $B_{c}=\cos(B)$ and $B_{s}=\sin(B)$, where sine and cosine are computed elementwise. Now we are ready to construct the matrix $A\in\mathbb{R}^{m\times 2d}$ that gives a consistent linear system. We use the columns of the matrices $B_{c}$ and $B_{s}$ in the following way:

the $i$-th column of the matrix $B_{c}$ is the $(2i-1)$-th column of the matrix $A$ for $i=1,\ldots,d$ 2. 2)

the $i$-th column of the matrix $B_{s}$ is the $2i$-th column of the matrix $A$ for $i=1,\ldots,d$.

Thus, for sufficiently large $d$ and small $\sigma$ we have a consistent linear system $Ax=b$ [2].

Figure 7 shows the convergence of the considered methods for different values of $d$. Experiments show that $\sigma=1$ gives a sufficiently accurate solution, therefore we fix it and show the dependence only on the dimension $d$. Figure 7b shows that already for $d=5$ the preconditioned Kaczmarz method converges much faster for all considered values of $\gamma$.

5 Conclusion

This paper presents the preconditioned Kaczmarz method by sketching technique and demonstrates its performance in different applications. However, the proposed approach requires a number of sketched rows that is an unknown hyperparameter, therefore theoretical investigation on the minimal, sufficient number of the sketched rows is the interesting further study. Another ingredient of the Kaczmarz method is the distribution for a sampling of rows. To sample rows in the right way, the methods for sampling from unnormalized distributions can be used. Also, the preconditioning approach can be generalized to non-linear systems and inspire modifications of the corresponding non-linear stochastic least-squares solvers. One more interesting extension of this work is to consider state-of-the-art adaptive stochastic gradient methods [4, 10] from the preconditioning perspective like in [21]. This point of view can lead to more efficient and theoretically motivated optimization methods.