On the Randomized Kaczmarz Algorithm

I Problem Statement

The KA can be described as follows. Let us define the hyperplane $H_{i}$ as:

where the $i$-th row of $A$ is denoted as $\mathbf{a}_{i}^{T}$ and the $i$-th element of $\mathbf{b}$ is denoted as $b_{i}$. Geometrically, the solution of (1) can be thought as the intersection of all hyperplanes $\{H_{i}\}_{i=1}^{m}$, and the KA seeks to find the solution by successively projecting to the hyperplanes from an initial approximation $\mathbf{x}_{0}$. The process is mathematically written as

where $i=mod(k,m)+1.$ Here we use the Matlab convention $mod(\cdot,\cdot)$ to denote the modulus after the division operation. Fig. 1 illustrates the algorithm in a low dimensional case.

The key difference between the RKA and the KA is that RKA chooses the rows following a specified probability distribution. More precisely, the probability for selecting $\mathbf{a}_{i}^{T}$ is given as $\frac{\|\mathbf{a}_{i}\|_{2}^{2}}{\|A\|_{F}^{2}}.$ Note that this probability is proportional to the row norms.

Although the KA is simple to state, its rate of convergence is still not completely explored. While for the RKA, with the predescribed choice of the probability distribution, the following convergence result is set up in :

However, it is argued in that ’Assigning probabilities corresponding to the row norms is in general certainly not optimal’. In the follows, we will try to find an optimized probability distribution for selecting the rows from $A$, so that a better performance can be obtained. The distribution vector is derived by minimizing an upper bound to the convergence rate which can be obtained by solving a convex optimization problem.

This note is organized as follows. The next section discusses the main results; In section 3, we discuss how to approximately solve the arising Semi-Definite-Programming (SDP) problem with smaller computational cost; In section 4, illustrative experiments will be conducted to verify the findings; Finally, we draw some conclusions in section 5.

II Optimized RKA

i.e. every row of the matrix $B$ is a normalized version of the corresponding row of matrix $A$.

Assume that currently we have $\mathbf{x}_{k-1}$, and based on $\mathbf{x}_{k-1}$, the next approximation $\mathbf{x}_{k}$ is given by (2), in which the index $i$ is chosen randomly according to $\mathbf{p}$. By the property of the projection operation, we have that

in which $\alpha_{i}$ denotes the angle between $\mathbf{x}_{k-1}-\mathbf{x}$ and the selected $\mathbf{b}_{i}$, i.e. the normal direction of the chosen hyperplane.

Based on the previous formula, we have that

in which $\beta_{i}$ denotes the angle between $\mathbf{y}$ and $\mathbf{b}_{i}$.

Based on the relations in (6), (7) and (8), we have that

By iterating the relations given in eq. (9) and eq. (10), the following results follow.

in which the expectations are taken with respect to all the random choices of the rows up to time $k$.

Note that $\Omega_{1}<1$ can be guaranteed if $\mathbf{p}$ is a strictly positive vector. This can be proven by a contradiction argument as follows. If $\Omega_{1}=1$, and since $\sin^{2}(\beta_{i})\leq 1$ for any $i$ and $\sum_{i=1}^{m}p_{i}=1$, we have that $\sin^{2}(\beta_{i})=1$, i.e. $\cos(\beta_{i})=0$ holds for all $i$. Considering that $rank(A)=n$, i.e. $rank(B)=n$, hence $\mathbf{x}_{k}-\mathbf{x}$ can not be orthogonal to the vectors $\{\mathbf{b}_{i}\}_{i=1}^{m}$, and the result follows. Based on this observation, we can see that exponential convergence in expectation can be obtained by a wide range of probability distribution vectors. This finding extends the result in , which only guarantees the exponential convergence for a given specific choice of the probability distribution vector. $\blacksquare$

According to Theorem 1, in order to get a better performance, we need to find a probability distribution vector, such that $\Omega_{1}$ can be made as small as possible. When the optimized $\Omega_{1}$ is obtained, we can also have a lower bound to the convergence speed of the RKA based on $\Omega_{2}$. In the following, we will first derive a closed form for $\Omega_{1}$ and $\Omega_{2}$, and then introduce a convex optimization problem to calculate the probability distribution vector $\hat{\mathbf{p}}$ which minimizes $\Omega_{1}$.

so in order to minimize $\Omega_{1}$, equivalently, we can maximize the following

If we restrict $\|\mathbf{y}\|_{2}=1$, then we have that

in which $\sigma_{n}(\cdot)$ denotes the smallest singular value of the matrix. The previous discussions can be summarized as:

in which $\sigma_{1}(\cdot)$ denotes the maximal singular value of the matrix.

Notice that minimizing $\Omega_{1}$ is equivalent to maximizing $\sigma_{n}(B^{T}\operatorname*{diag}(\mathbf{p})B)$, then we can solve the following problem instead:

This problem can be rewritten as the following SDP problem, in which $\hat{t}$ denotes the optimized $\sigma_{n}$ and $\hat{\mathbf{p}}$ denotes the corresponding probability distribution vector:

After solving the optimization problem of (16), $\hat{\mathbf{p}}$ is applied to the RKA to select the rows. Such a scheme will be abbreviated as ORKA in the following.

There exist cases such that $\Omega_{1}=\Omega_{2}$, i.e. there exists a vector $\mathbf{p}$, such that

i.e. $B^{T}\operatorname*{diag}(\mathbf{p})B=\frac{1}{n}I_{n}.$ In such cases, $\Omega_{1}=\Omega_{2}=1-\frac{1}{n}$, and the optimized probability distribution obtained by solving eq. (16) is the same as suggested in . It can be verified that when the columns of $A$ are orthogonal and of equal norm, then such property will hold. $\blacksquare$

The optimization problem (16) can also be formulated as

in the sense that $\hat{t}=\frac{1}{\mathbf{1}^{T}\hat{\mathbf{q}}}$ and $\hat{\mathbf{p}}=\hat{t}\hat{\mathbf{q}}.$

Since $\mathbf{q}$ in (17) is nonnegative, one has that $\mathbf{1}^{T}\mathbf{q}=\|\mathbf{q}\|_{1}$. It is known that the $l_{1}$ norm minimization problem is likely to return sparse solutions, which gives that $\hat{\mathbf{q}}$ is likely to be sparse. In the experiment section, we will also illustrate this phenomena. $\blacksquare$

Next, we discuss the relation between the ORKA and the RKA. It is obvious that the projection operations in (2) depend only on the corresponding normal vectors of the hyperplanes $\{H_{i}\}_{i=1}^{m}$, so we can optimize $\kappa(A)=\|A\|_{F}\|A^{\dagger}\|_{2}$ subject to the norms of the rows of matrix $A$. The optimization problem is given as

Set $\mathbf{1}^{T}\mathbf{q}=1$ and notice the fact that $A^{T}A=B^{T}\operatorname*{diag}(\mathbf{q})B$, then we can rewrite the previous problem as follows

It can be observed that this optimization is equivalent to the problem given by (16).

We conclude this observation in the following theorem.

The ORKA can do at least as good as the RKA, in the sense that if we optimize $\kappa(A)$ over the norms of rows of $A$, we obtain the same probability distribution vector as the one obtained by the ORKA.

III Further Discussions

Note that although the formulation in (16) is convex, it is still time consuming to solve this SDP optimization problem. In this section, we will discuss two possibilities to solve it approximately , which can alleviate some of the computational cost. One approximation of (16) is obtained by relaxing the constraint $B^{T}\operatorname*{diag}(\mathbf{p})B-tI_{n}\succeq 0$ by the following linear constraints:

In order to get a better relaxation, we introduce another approximation method which relates to the research of Optimal Input Design . Notice that $tr(B^{T}\operatorname*{diag}(\mathbf{p})B)=1$, i.e. the summation of all the singular values of $B^{T}\operatorname*{diag}(\mathbf{p})B$ is fixed, then maximizing $\sigma_{n}(B^{T}\operatorname*{diag}(\mathbf{p})B)$ means that we want all the singular values of $B^{T}\operatorname*{diag}(\mathbf{p})B$ to be close. This leads us to consider maximizing the product of the singular values of $B^{T}\operatorname*{diag}(\mathbf{p})B$, or maximizing the determinant of $B^{T}\operatorname*{diag}(\mathbf{p})B$. As the $\log$ function is monotonically increasing, we can optimize the following

in which $|\cdot|$ denotes the matrix determinant. Optimizing this quantity subject to the same constraints of (15) boils down to solve the so-called D-Optimal Design problem. One simple iterative algorithm to solve such problem has been suggested in , which is given as

Here, $\mathbf{p}^{t}$ denotes the estimation at time $t$, and $p_{i}^{t}$ denotes its $i$-th element. It has been proven in that for this algorithm, $\log|B^{T}\operatorname*{diag}(\mathbf{p}^{t})B|$ decreases monotonically w.r.t. $t$. We will make use of such property to approximately solve (15) when the objective function is replace by (20). More discussions will be given in next section.

IV Experiments

In this section, we will conduct experiments to illustrate the efficacy of the presented methods. The setup of our experiment is given as follows. The matrix $A$ is first generated by randn(m,n) in Matlab with $m=200$ and $n=20$, after that, each row is normalized, and then scaled with a random number which is uniformly distributed in $$. The reason for generating$A$as such is that in the first stage, the generated rows of$A$will have different directions which are uniformly distributed on the sphere$S^{n-1}$; and in the second stage, different rows of$A$with be assigned with different norms, which is directly related to the probability distribution vector chosen in .$\mathbf{x}$is generated by randn(n,1), and$\mathbf{b}$is generated as$\mathbf{b}=A\mathbf{x}$. We will compare the Mean Square Error (MSE) along the projection path obtained by all these methods, the first is the one suggested in (abbreviated as RKA), the second is the one obtained by the SDP optimization given by (16) (abbreviated as ORKA) and the third is the one obtained by the LP approximations given by (19) (abbreviated as LPORKA), the last is the one obtained by the iterative method to solve the D-Optimal Design criteria (abbreviated as ITEORKA). We iterate (21) for$10$times in this experiment. For each method, we run the experiment 2000 times to get the averaged performance. The CVX toolboxhttp://cvxr.com/ is used to solve the SDP and LP optimization problems. From the experiment, we can observe that the time for solving the LP problem in LPORKA is close to the time needed for the$10$ iterations of (21), and the time needed for solving (16) in ORKA is approximately 7 times as them.

V Conclusion

This note discusses the possibility and methodology to find a probability distribution vector for selecting the rows of $A$ to result in a better convergence speed of the Randomize Kaczmarz Algorithm. The lower bound and upper bound for the convergence speed is derived first. Then an optimized probability distribution vector is obtained by minimizing the upper bound, which turns to be given by solving a convex optimization problem. Properties of the approach are also discussed along the note.

References