Randomized Methods for Linear Constraints: Convergence Rates and Conditioning

Introduction

The condition number of a problem measures the sensitivity of a solution to small perturbations in its input data. For many problems that arise in numerical analysis, there is often a simple relationship between the condition number of a problem instance and the distance to the set of ill-posed problems—those problem instances whose condition numbers are infinite . For example, with respect to the problem of inverting a matrix $A$, it is known (see , for example) that if $A$ is perturbed to $A+E$ for sufficiently small $E$, then

Thus, a condition measure for this problem may be taken as $\|A^{-1}\|$. Associated with this is the classical Eckart-Young theorem found in , relating the above condition measure to the distance to ill-posedness.

We are typically concerned with relative condition numbers as introduced by Demmel in . For example, with respect to the problem of matrix inversion, the relative condition number is $k(A):=\|A\|\|A^{-1}\|$, the commonly used condition measure.

Condition numbers are also important from an algorithmic perspective. In the above example of matrix inversion, for instance, the sensitivity of a problem under perturbations could come into prominence regarding errors in either the initial problem data or accumulated computational error due to rounding. Hence, it seems natural that condition numbers should affect algorithm speed. For example, in the context of linear programming, Renegar defined a condition measure based on the distance to ill-posedness in –in a similar sense as the Eckart-Young result–and showed its effect on the convergence rate of interior point methods in .

For another example, consider the problem of finding a solution to the system $Ax=b$ where $A$ is a positive-definite matrix. It was shown in that the steepest descent method is linearly convergent with rate $(\frac{k(A)-1}{k(A)+1})^{2}$ and that this bound is asymptotically tight for almost all choices of initial iterates. Similarly, it is well known (see ) that the conjugate gradient method applied to the same problem is also linearly convergent with rate $\frac{\sqrt{k(A)}-1}{\sqrt{k(A)}+1}$.

From a computational perspective, a related and important area of study is that of error bounds. Given a subset of a Euclidean space, an error bound is an inequality that bounds the distance from a test vector to the specified subset in terms of some residual function that is typically easy to compute. In that sense, an error bound can be used both as part of a stopping rule during implementation of an algorithm as well as an aide in proving algorithmic convergence. A comprehensive survey of error bounds for a variety of problems arising in optimization can be found in .

With regards to the problem of solving a nonsingular linear system $Ax=b$, one connection between condition measures and error bounds is immediate. Let $x^{*}$ be a solution to the system and $x$ be any other vector. Then

so the distance to the solution set is bounded by a constant multiple of the residual vector, $\|Ax-b\|$, and this constant is the same one that appears in the context of conditioning and distance to infeasibility. As we discuss later, this result is not confined to systems of linear equations.

As a result, error bounds and the related condition numbers often make a prominent appearance in convergence proofs for a variety of algorithms. In this paper, motivated by a recent randomized iterated projection scheme for systems of linear equations due to Strohmer and Vershynin in , we revisit some classical algorithms and show that, with an appropriate randomization scheme, we can demonstrate convergence rates directly in terms of these natural condition measures. The rest of the paper is organized as follows. In the remainder of this section, we define some notation used throughout the rest of this paper. In Section 3, we consider the problem of solving a linear system $Ax=b$ and show that a randomized coordinate descent scheme, implemented according to a specific probability distribution, is linearly convergent with a rate expressible in terms of traditional conditioning measures. In Section 4, we build upon the work of Strohmer and Vershynin in by considering randomized iterated projection algorithms for linear inequality systems. In particular, we show how randomization can provide convergence rates in terms of the traditional Hoffman error bound in as well as in terms of Renegar’s distance to infeasibility from . In Section 5, we consider randomized iterated projection algorithms for general convex sets and, under appropriate metric regularity assumptions, obtain local convergence rates in terms of the modulus of regularity.

The classical, deterministic versions of the simple algorithms we consider here have been widely studied, in part due to the extreme simplicity of each iteration: their linear convergence is well-known. However, as remarked for linear systems of equations in , randomized versions are interesting for several reasons. The randomized iterated projection method for linear equations from which this work originated may have some practical promise, even compared with conjugate gradients, for example . Our emphasis here, however, is theoretical: randomization here provides a framework for simplifying the analysis of algorithms, allowing easy bounds on the rates of linear convergence in terms of natural linear-algebraic condition measures, such as relative condition numbers, Hoffman constants, and the modulus of metric regularity.

Notation

On the Euclidean space ${\bf R}^{n}$, we denote the Euclidean norm by $\|\cdot\|$. Let $e_{i}$ denote the column vector with a 1 in the $i^{th}$ position and zeros elsewhere.

We consider $m$-by-$n$ real matrices $A$. We denote the set of rows of $A$ by $\{a_{1}^{T},\ldots,a_{m}^{T}\}$ and the set of columns is denoted $\{A_{1},\ldots,A_{n}\}$. The spectral norm of $A$ is the quantity $\|A\|_{2}:=\max_{\|x\|=1}\|Ax\|$ and the Frobenius norm is $\|A\|_{F}:=\sum_{i,j}a_{ij}^{2}$. These norms satisfy the following inequality :

For an arbitrary matrix, $A$, let $\|A^{-1}\|_{2}$ be the smallest constant $M$ such that $\|Ax\|_{2}\geq\frac{1}{M}\|x\|_{2}$ for all vectors $x$. In the case $m\geq n$, if $A$ has singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{n}$, then $M$ can also be expressed as the reciprocal of the minimum singular value $\sigma_{n}$, and, if $A$ is invertible, this quantity equals the spectral norm of $A^{-1}$.

The relative condition number of $A$ is the quantity $k(A):=\|A\|_{2}\|A^{-1}\|_{2}$; related to this is the scaled condition number introduced by Demmel in , defined by $\kappa(A):=\|A\|_{F}\|A^{-1}\|_{2}$. From this, it is easy to verify (using the singular value decomposition, for example) the following relationship between condition numbers:

Now suppose the matrix $A$ is $n$-by-$n$ symmetric and positive definite. The energy norm (or $A$-norm), denoted $\|\cdot\|_{A}$, is defined by $\|x\|_{A}:=\sqrt{x^{T}Ax}$. The inequality

is useful later. Further, if $A$ is simply positive semi-definite, we can generalize Inequality 2.2:

where $\underline{\lambda}(A)$ is the smallest non-zero eigenvalue of $A$. We denote the trace of $A$ by $\mbox{\rm tr}\,A$: it satisfies the inequality

Given a nonempty closed convex set $S$, let $P_{S}(x)$ be the projection of $x$ onto $S$: that is, $P_{S}(x)$ is the vector $y$ that is the optimal solution to $\min_{z\in S}\|x-z\|_{2}$. Additionally, define the distance from $x$ to a set $S$ by

The following useful inequality is standard:

Randomized Coordinate Descent

Let $A$ be an $n$-by-$n$ positive-definite matrix. We consider a linear system of the form $Ax=b$, with solution $x^{*}=A^{-1}b$. We consider the equivalent problem of minimizing the strictly convex quadratic function

Suppose our current iterate is $x$ and we obtain a new iterate $x_{+}$ by performing an exact line search in the nonzero direction $d$: that is, $x_{+}$ is the solution to $\min_{t\in{\bf R}}f(x+td)$. This gives us

One natural choice of a set of easily-computable search directions is to choose $d$ from the set of coordinate directions, $\{e_{1},\ldots,e_{n}\}$. Note that, when using search direction $e_{i}$, we can compute the new point

using only $2n+2$ arithmetic operations. If the search direction is chosen at each iteration by successively cycling through the set of coordinate directions, then the algorithm is known to be linearly convergent but with a rate not easily expressible in terms of typical matrix quantities (see or ) . However, by choosing a coordinate direction as a search direction randomly according to an appropriate probability distribution, we can obtain a convergence rate in terms of the relative condition number. This is expressed in the following result.

Consider an $n$-by-$n$ positive semidefinite system $Ax=b$ and let $x_{0}\in{\bf R}^{n}$ be an arbitrary starting point. For $j=0,1,\ldots,$ compute

where, at each iteration $j$, the index $i$ is chosen independently at random from the set $\{1,\ldots,n\}$, with distribution

Notice in the algorithm that the matrix $A$ may be singular, but that nonetheless $a_{ii}>0$ almost surely. If $A$ is merely positive semidefinite, solutions of the system $Ax=b$ coincide with minimizers of the function $f$, and consistency of the system is equivalent to $f$ being bounded below. We now have the following result.

Consider a consistent positive-semidefinite system $Ax=b$, and define the corresponding objective and error by

Then Algorithm 3.3 is linearly convergent in expectation: indeed, for each iteration $j=0,1,2,\ldots$,

In particular, if $A$ is positive-definite and $x^{*}=A^{-1}b$, we have the equivalent property

Hence the expected reduction in the squared error $\|x_{j}-x^{*}\|_{A}^{2}$ is at least a factor

Proof Note that if coordinate direction $e_{i}$ is chosen during iteration $j$, then Equation 3.2 shows

Using Inequality 2.3 and Equation 3.1, we easily verify

and the first result follows. Applying Equation 3.1 provides the second result. The final result comes from applying Inequalities 2.1 and 2.4. $\Box$

The simple idea behind the proof of Theorem 3.4 is the main engine driving the remaining results in this paper. Fundamentally, the idea is to choose a probability distribution so that the expected distance to the solution from the new iterate is the distance to the solution from the old iterate minus some multiple of a residual. Then, using some type of error bound to bound the distance to a solution in terms of the residual, we obtain expected linear convergence of the algorithm.

Now let us consider the more general problem of finding a solution to a linear system $Ax=b$ where $A$ is an $m\times n$. More generally, since the system might be inconsistent, we seek a “least squares solution” by minimizing the function $\|Ax-b\|^{2}$. The minimizers are exactly the solutions of the positive-semidefinite system $A^{T}Ax=A^{T}b$, to which we could easily apply the previous algorithm; however, as usual, we wish to avoid computing the new matrix $A^{T}A$ explicitly. Instead, we can proceed as follows.

Consider a linear system $Ax=b$ for a nonzero $m$-by-$n$ matrix $A$. Let $x_{0}\in{\bf R}^{n}$ be an arbitrary initial point and let $r_{0}=b-Ax_{0}$ be the initial residual. For each $j=0,1,\ldots,$ compute

where, at each iteration $j$, the index $i$ is chosen independently at random from the set $\{1,\ldots,n\}$, with distribution

(In the formula for $\alpha_{j}$, notice by assumption that $A_{i}\neq 0$ almost surely.)

Note that the step size at each iteration can be obtained by directly minimizing the residual in the respective coordinate direction. However, the algorithm can also be viewed as the application of the algorithm for positive definite systems on the system of normal equations, $A^{T}Ax=A^{T}b$, without actually having to compute the matrix $A^{T}A$. Given the motivation of directly minimizing the residual, we would expect that Algorithm 3.5 would converge to a least squares solution, even in the case where the underlying system is inconsistent. The next result shows that this is, in fact, the case.

Consider any linear system $Ax=b$, where the matrix $A$ is nonzero. Define the least-squares residual and the error by

Then Algorithm 3.5 is linearly convergent in expectation to a least squares solution for the system: for each iteration $j=0,1,2\ldots$,

In particular, if $A$ has full column rank, we have the equivalent property

where $\hat{x}=(A^{T}A)^{-1}A^{T}b$ is the unique least-squares solution.

Proof It is easy to verify, by induction on $j$, that the iterates $x_{j}$ are exactly the same as the iterates generated by Algorithm 3.3, when applied to the positive semi-definite system $A^{T}Ax=A^{T}b$, and furthermore that the residuals satisfy $r_{j}=b-Ax_{j}$ for all $j=0,1,2,\ldots$. Hence, the results follow directly by Theorem 3.4. $\Box$

By the coordinate descent nature of this algorithm, once we have computed the initial residual $r_{0}$ and column norms $\{\|A_{i}\|^{2}\}_{i=1}^{n}$, we can perform each iteration in $O(n)$ time, just as in the positive-definite case. Specifically, this new iteration takes $4n+1$ arithmetic operations, compared with $2n+2$ for the positive-definite case.

For a computational example, we apply Algorithm 3.5 to random $500\times n$ matrices where each element of $A$ and $b$ is an independent Gaussian random variable and we let $n$ take values 50, 100, 150 and 200.

Note that in the above examples, the theoretical bound provided by Theorem 3.6 predicts the actual behavior of the algorithm reasonably well.

Randomized Iterated Projections

Iterated projection algorithms share some important characteristics with coordinate descent algorithms. Both are well studied and much convergence theory exists; a comprehensive overview on iterated projections can be found in . However, even for linear systems of equations, standard developments do not provide bounds on convergence rates in terms of usual condition numbers. By contrast, in the recent paper , Strohmer and Vershynin obtained such bounds via the following randomized iterated projection algorithm, which also provided the motivation for our work in the previous section.

Consider a linear system $Ax=b$ for a nonzero $m$-by-$n$ matrix $A$. Let $x_{0}\in{\bf R}^{n}$ be an arbitrary initial point. For each $j=0,1,\ldots,$ compute

where, at each iteration $j$, the index $i$ is chosen independently at random from the set $\{1,\ldots,m\}$, with distribution

Notice that the new iterate $x_{j+1}$ is simply the orthogonal projection of the old iterate $x_{j}$ onto the hyperplane $\{x:a_{i}^{T}x=b_{i}\}$. At first sight, the choice of probability distribution may seem curious, since we could rescale the equations arbitrarily without having any impact on the projection operations. However, following , we emphasize that the aim is to understand linear convergence rates in terms of linear-algebraic condition measures associated with the original system, rather than in terms of geometric notions associated with the hyperplanes. This randomized algorithm has the following behavior.

Given any matrix $A$ with full column rank, suppose the linear system $Ax=b$ has solution $x^{*}$. Then Algorithm 4.1 converges linearly in expectation: for each iteration $j=0,1,2,\ldots$,

We seek a way of generalizing the above algorithm and convergence result to more general systems of linear inequalities, of the form

where the disjoint index sets $I_{\leq}$ and $I_{=}$ partition the set $\{1,2,\ldots,m\}$. To do so, staying with the techniques of the previous section, we need a corresponding error bound for a system of linear inequalities. First, given a vector $x\in{\bf R}^{n}$, define the vector $x^{+}$ by $(x^{+})_{i}=\max\{x_{i},0\}$. Then a starting point for this subject is a result by Hoffman in .

For any right-hand side vector $b\in{\bf R}^{m}$, let $S_{b}$ be the set of feasible solutions of the linear system (4.3). Then there exists a constant $L$, independent of $b$, with the following property:

where the function $e\colon{\bf R}^{m}\to{\bf R}^{m}$ is defined by

In the above result, each component of the vector $e(Ax-b)$ indicates the error in the corresponding inequality or equation. In particular $e(Ax-b)=0$ if and only if $x\in S_{b}$. Thus Hoffman’s result provides a linear bound for the distance from a trial point $x$ to the feasible region in terms of the size of the “a posteriori error” associated with $x$.

We call the minimum constant $L$ such that property (4.5) holds the Hoffman constant for the system (4.3). Several authors give geometric or algebraic meaning to this constant, or exact expressions for it, including , , ; for a more thorough treatment of the subject, see . In the case of linear equations (that is, $I_{\leq}=\emptyset$), an easy calculation using the singular value decomposition shows that the Hoffman constant is just the reciprocal of the smallest nonzero singular value of the matrix $A$, and hence equals $\|A^{-1}\|_{2}$ when $A$ has full column rank.

For the problem of finding a solution to a system of linear inequalities, we consider a randomized algorithm generalizing Algorithm 4.1.

Consider the system of inequalities (4.3). Let $x_{0}$ be an arbitrary initial point. For each $j=0,1,\ldots,$ compute

where, at each iteration $j$, the index $i$ is chosen independently at random from the set $\{1,\ldots,m\}$, with distribution

In the above algorithm, notice $\beta_{j}=e(Ax_{j}-b)_{i}$. We can now generalize Theorem 4.2 as follows.

Suppose the system (4.3) has nonempty feasible region $S$. Then Algorithm 4.6 converges linearly in expectation: for each iteration $j=0,1,2,\ldots$,

Proof Note that if the index $i$ is chosen during iteration $j$, then it follows that

Note $P_{S}(x_{j})\in S$. Hence if $i\in I_{\leq}$, then $a_{i}^{T}P_{S}(x_{j})\leq b_{i}$, and $e(Ax_{j}-b)_{i}\geq 0$, so

On the other hand, if $i\in I_{=}$, then $a_{i}^{T}P_{S}(x_{j})=b_{i}$, so

Putting these two cases together with the previous inequality shows

Taking the expectation with respect to the specified probability distribution, it follows that

and the result now follows by the Hoffman bound. $\Box$

Since Hoffman’s bound is not independent of the scaling of the matrix $A$, it is not surprising that a normalizing constant like $\|A\|^{2}_{F}$ term appears in the result.

For a computational example, we consider linear inequality systems $Ax\leq b$ where the elements of $A$ are independent standard Gaussian random variables and $b$ is chosen so that the resulting system has a non-empty interior. We consider matrices $A$ which are $500\times n$, letting $n$ take values 50, 100, 150 and 200. We then apply Algorithm 4.6 to these problems and observe the following computational results.

Another natural conditioning measure for linear inequality systems is the distance to infeasibility, defined by Renegar in , and shown in to govern the convergence rate of interior point methods for linear programming. It is interesting, therefore, from a theoretical perspective, to obtain a linear convergence rate for iterated projection algorithms in terms of this condition measure as well. For simplicity, we concentrate on the inequality case, $Ax\leq b$. To begin, let us recall the following results.

The distance to infeasibility for the system $Ax\leq b$ is the number

Suppose $\mu>0$. Then there exists a point $\hat{x}\in S$ satisfying $\|\hat{x}\|\leq\|b\|/\mu$. Furthermore, any point $x\in{\bf R}^{n}$ satisfies the inequality

Using this, we can bound the linear convergence rate for the Algorithm 4.6 in terms of the distance to infeasibility, as follows. Notice first that $\|x_{j}-\hat{x}\|$ is nonincreasing in $j$, by Inequality 2.5. Suppose we start Algorithm 4.6 at the initial point $x_{0}=0$. Applying Theorem 4.8, we see that for all $j=1,2,\ldots,$

Using this inequality in place of Hoffman’s bound in the proof of Theorem 4.7 gives

Although this bound may not be the best possible (and, in fact, it may not be as good as the bound provided in Theorem 4.7), this result simply emphasizes a relationship between algorithm speed and conditioning measures that appears naturally in other contexts. In the next section, we proceed with these ideas in a more general framework.

Metric Regularity and Local Convergence

The previous section concerned global rates of linear convergence. If instead we are interested in local rates, we can re-examine a generalization of our problem through an alternative perspective of set-valued mappings. Consider a set-valued mapping \Phi:{\bf R}^{n}\;{\lower 1.0pt\hbox{\rightarrow}}\kern-12.0pt\hbox{\raise 2.8pt\hbox{\rightarrow}}\;{\bf R}^{m} and the problem of solving the associated constraint system of the form $b\in\Phi(x)$ for the unknown vector $x$. For example, finding a feasible solution to $Ax\leq b$ is equivalent to finding an $x$ such that

Related to this is the idea of metric regularity of set-valued mappings. We say the set-valued mapping $\Phi$ is metrically regular at $\bar{x}$ for $\bar{b}\in\Phi(\bar{x})$ if there exists $\gamma>0$ such that

where $\Phi^{-1}(b)=\{x:b\in\Phi(x)\}$. Further, the modulus of regularity is the infimum of all constants $\gamma$ such that Equation 5.2 holds. Metric regularity is strongly connected with a variety of ideas from variational analysis: a good background reference is .

Metric regularity generalizes the error bounds discussed in previous sections at the expense of only guaranteeing a bound in local terms. For example, if $\Phi$ is a single-valued linear map, then the modulus of regularity (at any $\bar{x}$ for any $\bar{b}$) corresponds to the typical conditioning measure $\|\Phi^{-1}\|$ (with $\|\Phi^{-1}\|=\infty$ implying the map is not metrically regular) and if $\Phi$ is a smooth single-valued mapping, then the modulus of regularity is the reciprocal of the minimum singular value of the Jacobian, $\nabla\Phi(x)$. From an alternative perspective, metric regularity provides a framework for generalizing the Eckart-Young result on the distance to ill-posedness of linear mappings cited in Theorem 1.1. Specifically, if we define the radius of metric regularity at $\bar{x}$ for $\bar{b}$ for a set-valued mapping $\Phi$ between finite dimensional spaces by

where the infimum is over all linear functions $E$, then one obtains the strikingly simple relationship (see )

We will not be directly using the above result. Here, we simply use the fundamental idea of metric regularity which says that the distance from a point to the solution set, $d(x,\Phi^{-1}(b))$, is locally bounded by some constant times a ”residual”. For example, in the case where $\Phi$ corresponds to the linear inequality system (5.1), we have that $d(b,\Phi(x))=\|(Ax-b)^{+}\|$ implies that the modulus of regularity is in fact a global bound and equals the Hoffman bound. More generally, we wish to emphasize that metric regularity ties together several of the ideas from previous sections at the expense of those results now only holding locally instead of globally.

In what follows, assume all distances are Euclidean distances. We wish to consider how the modulus of regularity of $\Phi$ affects the convergence rate of iterated projection algorithms. We remark that linear convergence for iterated projection methods on convex sets has been very widely studied: see , for example. Our aim here is to observe, by analogy with previous sections, how randomization makes the linear convergence rate easy to interpret in terms of metric regularity.

Let $S_{1},S_{2},\ldots,S_{m}$ be closed convex sets in a Euclidean space ${\bf E}$ such that $\cap_{i}S_{i}\neq\emptyset$. Then, in a manner similar to , we can endow the product space ${\bf E}^{m}$ with the inner product

and consider the set-valued mapping $\Phi:{\bf E}\rightarrow{\bf E}^{m}$ given by

Then it clearly follows that $\bar{x}\in\cap_{i}S_{i}\Leftrightarrow 0\in\Phi(\bar{x})$. Under appropriate regularity assumptions, we obtain the following local convergence result.

Suppose the set-valued mapping $\Phi$ given by Equation 5.3 is metrically regular at $\bar{x}$ for 0 with regularity modulus $\gamma$. Let $\bar{\gamma}$ be any constant strictly larger than $\gamma$ and let $x_{0}$ be any initial point sufficiently close to $\bar{x}$. Further, suppose that $x_{j+1}=P_{S_{i}}(x_{j})$ with probability $\frac{1}{m}$ for $i=1,\ldots,m$. Then

Proof First, note that by Inequality 2.5, the distance $\|x_{j}-\bar{x}\|$ is nonincreasing in $j$. Hence if $x_{0}$ is sufficiently close to $\bar{x}$, then $x_{j}$ is as well for all $j\geq 0$. Then, again using Inequality 2.5 (applied to the set $S_{i}$), we have, for all points $x\in S\subset S_{i}$,

Taking the minimum over $x\in S$, we deduce

using the definition of metric regularity. $\Box$

Note that metric regularity at $\bar{x}$ for 0 is a slightly stronger assumption than actually necessary for this result. Specifically, the above result holds as long as Equation 5.2 holds for all $x$ near $\bar{x}$ with $\bar{b}=0$ fixed, as opposed to the above definition requiring it to hold for all $b$ near $\bar{b}$ as well.

For a moment, let $m=2$ and consider the sequence of iterates $\{x_{j}\}_{j\geq 0}$ generated by the randomized iterated projection algorithm. By idempotency of the projection operator, there’s no benefit to projecting onto the same set in two consecutive iterations, so the subsequence consisting of different iterates corresponds exactly to that of the non-randomized iterated projection algorithm. In particular, if $x_{j}\in S_{1}$, then

since $d(x_{j},S_{1})=0$. This gives us the following corollary, which also follows through more standard deterministic arguments.

If $\Phi$ is metrically regular at $\bar{x}$ for 0 with regularity modulus $\gamma$ and $\bar{\gamma}$ is larger than $\gamma$, then for $x_{0}$ sufficiently close to $\bar{x}$, the 2-set iterated projection algorithm is linearly convergent and

Further, consider the following refined version of the $m$-set randomized algorithm. Suppose $x_{0}\in S_{1}$ and $i_{0}=1$. Then for $j=1,2,\ldots,$ let $i_{j}$ be chosen uniformly at random from $\{1,\ldots,m\}\backslash\{i_{j-1}\}$ and $x_{j+1}=P_{S_{i_{j}}}(x_{j})$. Then we obtain the following similar result.

If $\Phi$ is metrically regular at $\bar{x}$ for 0 with regularity modulus $\gamma$ and $\bar{\gamma}$ is larger than $\gamma$, then for $x_{0}$ sufficiently close to $\bar{x}$, the refined $m$-set randomized iterated projection algorithm is linearly convergent in expectation and

A simple but effective product space formulation by Pierra in has the benefit of reducing the problem of finding a point in the intersection of finitely many sets to the problem of finding a point in the intersection of 2 sets. Using the notation above, we consider the closed set in the product space given by

where the linear mapping $A:{\bf E}\rightarrow{\bf E}^{m}$ is defined by $Ax=(x,x,\ldots,x)$. Again, notice that $\bar{x}\in\cap_{i}S_{i}\Leftrightarrow(\bar{x},\ldots,\bar{x})\in T\cap L$. One interesting aspect of this formulation is that projections in the product space ${\bf E}^{m}$ relate back to projections in the original space ${\bf E}$ by

This formulation provides a nice analytical framework: we can use the above equivalence of projections to consider the method of averaged projections directly, defined as follows.

Let $S_{1},\ldots,S_{m}\subseteq E$ be nonempty closed convex sets. Let $x_{0}$ be an initial point. For $j=1,2,\ldots$, let

Simply put, at each iteration, the algorithm projects the current iterate onto each set individually and takes the average of those projections as the next iterate. In the product space formulation, this is equivalent to $x_{j+1}=P_{L}(P_{T}(x_{j}))$. Expanding on the work of Pierra in , additional convergence theory for this algorithm has been examined by Bauschke and Borwein in . Under appropriate regularity conditions, the general idea is that convergence of the iterated projection algorithm for two sets implies convergence of the averaged projection algorithm for $m$ sets. In a similar sense, we prove the following result in terms of randomized projections.

Suppose $S=\cap_{i=1}^{m}S_{i}$ is non-empty. If the randomized projection algorithm of Theorem 5.4 is linearly convergent in expectation with rate $\alpha$, then so is Algorithm 5.7.

Proof Let $x_{j}$ be the current iterate, $x_{j+1}^{AP}$ be the new iterate in the method of averaged projections and $x_{j+1}^{RP}$ be the new iterate in the method of uniformly randomized projections. Then note that:

By convexity of the $S_{i}$’s, it follows that

Hence, the method of averaged projections converges no more slowly than the method of uniformly random projections. In particular, under the assumptions of Theorem 5.4, the method of averaged projections converges with rate no larger than $1-\frac{1}{m\bar{\gamma}^{2}}$.

References