Paved with Good Intentions: Analysis of a Randomized Block Kaczmarz Method

Introduction

The Kaczmarz method [Kac37] is an iterative algorithm for solving overdetermined least-squares problems. Because of its simplicity and performance, this scheme has found application in fields ranging from image reconstruction to digital signal processing [SS87, CFM+92, FS95, Nat01]. At each iteration, the basic Kaczmarz method makes progress by enforcing a single constraint, while the block Kaczmarz method [Elf80] enforces many constraints at once. This paper introduces a randomized version of the block Kaczmarz method that converges with an expected linear rate, and we characterize the performance of this algorithm using geometric properties of the blocks of equations. This analysis leads us to consider the concept of a row paving of a matrix, a partition of the rows into well-conditioned blocks. We summarize the literature on row pavings, and we explain how this theory interacts with the block Kaczmarz method. Together, these results yield an efficient block Kaczmarz scheme that applies to many overdetermined least-squares problems.

Let $\bm{A}$ be a real or complex $n\times d$ matrix with full column rank, and suppose that $\bm{b}$ is a vector with dimension $n$. Consider the overdetermined least-squares problem

We say that $\bm{A}$ is standardized when (1.2) is in force.

The spectral norm is denoted by $\lVert{\cdot}\rVert$, while $\lVert{\cdot}\rVert_{\rm F}$ represents the Frobenius norm. When applied to an Hermitian matrix, the maps $\lambda_{\min}$ and $\lambda_{\max}$ return the algebraic minimum and maximum eigenvalues. For a $p\times q$ matrix $\bm{W}$, we arrange the singular values as follows.

The minimum singular value is positive if and only if $\bm{WW}^{*}$ or $\bm{W}^{*}\bm{W}$ is nonsingular. We define the condition number $\kappa(\bm{W}):=\sigma_{\max}(\bm{W})/\sigma_{\min}(\bm{W})$. The dagger $\dagger$ denotes the Moore–Penrose pseudoinverse. When $\bm{W}$ has full row rank, its pseudoinverse is determined by the formula $\bm{W}^{\dagger}:=\bm{W}^{*}(\bm{WW}^{*})^{-1}$.

2. The Simple Kaczmarz Method

The Kaczmarz method is an iterative algorithm that produces an approximation to the minimizer $\bm{x}_{\star}$ of the least-squares problem (1.1). The method commences with an arbitrary guess $\bm{x}_{0}$ for the solution. At the $j$th iteration, we select a row index $t=t(j)$ of the matrix $\bm{A}$, and we project the current iterate $\bm{x}_{j-1}$ onto the solution space of the equation $\left\langle{\bm{a}_{t}},\,{\bm{x}}\right\rangle=b_{t}$. That is,

This process continues until it triggers an appropriate convergence criterion.

To develop a complete algorithm, we also need a control mechanism that specifies how to select rows. For example, the most classical approach cycles through the rows in order. Instead, we focus on a modern formulation that uses a randomized control mechanism. Randomization has several benefits: the resulting algorithm is easy to analyze, it is simple to implement, and it is often effective in practice.

Our primary reference is the randomized Kaczmarz algorithm recently proposed by Strohmer and Vershynin [SV09b]. When $\bm{A}$ is standardized, their method operates as follows. At iteration $j$, independently of all previous random choices, the algorithm draws the row index $t(j)$ uniformly at random from the set $\{1,\dots,n\}$ of all row indices. Then the current iterate is updated using the rule (1.3). The paper [SV09b] provides a short, elegant proof that this iteration converges at an expected linear rate to the solution $\bm{x}_{\star}$ of a consistent least-squares problem (i.e., where the residual $\bm{e}$ is zero).

Needell [Nee10] has extended the argument of [SV09b] to the case of an inconsistent least-squares problem. For a standardized matrix $\bm{A}$, Needell’s error estimate reads

In words, the randomized Kaczmarz method converges in expectation at a linear rateMathematicians often use the term exponential convergence for the concept numerical analysts call linear convergence. until it reaches a fixed ball about the true solution $\bm{x}_{\star}$, at which point the error may cease to decay. The radius of this ball roughly equals the second term in (1.4), while the convergence rate is controlled by the bracket. When the residual $\bm{e}$ is zero, the bound (1.4) reduces to the error estimate from [SV09b].

When $\bm{A}$ is an $n\times d$ standardized matrix whose columns are well conditioned, the minimum singular value $\sigma^{2}_{\min}(\bm{A})\geq\textrm{const}\cdot n/d$. In this case, the error bound (1.4) simplifies to

3. The Block Kaczmarz Method

In some situations [EHL81], practitioners prefer to use a block version of the Kaczmarz method to solve the least-squares problem (1.1). We consider a formulation due to Elfving [Elf80]. This procedure begins with an initial guess $\bm{x}_{0}$ for the solution. At each iteration $j$, we select a subset $\tau=\tau(j)$ of the row indices of $\bm{A}$, and we project the current iterate $\bm{x}_{j-1}$ onto the solution space of $\bm{A}_{\tau}\bm{x}=\bm{b}_{\tau}$, the set of equations listed in $\tau$. That is,

This process continues until it has converged. We have written $\bm{A}_{\tau}$ for the row submatrix of $\bm{A}$ indexed by $\tau$, while $\bm{b}_{\tau}$ is the subvector of $\bm{b}$ with components listed in $\tau$. We assume that the row submatrix $\bm{A}_{\tau}$ is fatA $p\times q$ matrix is fat when $p\leq q$., so the pseudoinverse (1.5) returns the solution to an underdetermined least-squares problem.

To specify a block Kaczmarz algorithm, one must decide what blocks of indices are permissible, as well as the mechanism for selecting a block at each iteration. In this paper, we study a version that is based on two design decisions. First, this algorithm requires a partition $T=\{\tau_{1},\dots,\tau_{m}\}$ of the row indices of $\bm{A}$. The method only considers blocks of indices that appear in the partition $T$. Second, we use a simple randomized control scheme to choose which block to enforce. At each iteration, independently of all previous choices, we draw a block $\tau$ uniformly at random from the partition $T$. These decisions lead to Algorithm 1.3. We postpone a detailed discussion on implementation to Section 4.

We make no claim that randomized selection provides the optimal sequence for choosing blocks. As with the simple Kaczmarz method, the scaling of the rows of the matrix can play a significant role in the behavior of the algorithm. See [CHJ09, SV09a] for a discussion of this issue.

[tb] Block Kaczmarz Method with Uniform Random Control

\alginout • Matrix $\bm{A}$ with dimension $n\times d$ • Right-hand side $\bm{b}$ with dimension $n$ • Partition $T=\{\tau_{1},\dots,\tau_{m}\}$ of the row indices $\{1,\dots,n\}$ • Initial iterate $\bm{x}_{0}$ with dimension $d$ • Convergence tolerance $\varepsilon>0$ An estimate $\hat{\bm{x}}$ for the solution to $\min_{\bm{x}}\lVert{\bm{A}\bm{x}-\bm{b}}\rVert_{2}^{2}$ {algtab*} $j\leftarrow 0$\algrepeat $j\leftarrow j+1$ Choose a block $\tau$ uniformly at random from $T$ $\bm{x}_{j}\leftarrow\bm{x}_{j-1}+(\bm{A}_{\tau})^{\dagger}(\bm{b}_{\tau}-\bm{A}_{\tau}\bm{x}_{j-1})$ { Solve least-squares problem }\alguntil $\lVert{\bm{A}\bm{x}_{j}-\bm{b}}\rVert_{2}^{2}\leq\varepsilon^{2}$ $\hat{\bm{x}}\leftarrow\bm{x}_{j}$

4. Desiderata for the Partition

The implementation and behavior of the block Kaczmarz method depend heavily on the properties of the submatrices $\bm{A}_{\tau}$ indexed by the blocks $\tau$ in the partition $T$. Let us explain how the structure of the submatrices plays a role in the implementation; the claims about performance will emerge from the theoretical results in Section 1.5.

The most expensive (arithmetic) step in Algorithm 1.3 occurs when we apply the pseudoinverse $\bm{A}_{\tau}^{\dagger}$ to a vector. We can perform this calculation efficiently provided that each submatrix $\bm{A}_{\tau}$ has well-conditioned rows. Indeed, in this case, we can invoke an iterative least-squares solver [Bjö96], such as CGLS, to apply the pseudoinverse $\bm{A}_{\tau}^{\dagger}$ approximately using a small number of matrix–vector multiplies with $\bm{A}_{\tau}$ and $\bm{A}_{\tau}^{*}$. In particular, we never need to form the pseudoinverse.

This observation highlights how important it is to control the geometric properties of the submatrices $\bm{A}_{\tau}$ induced by $T$. Let us make a definition that encapsulates the information that we will need.

An $(m,\alpha,\beta)$ row paving of a matrix $\bm{A}$ is a partition $T=\{\tau_{1},\dots,\tau_{m}\}$ of the row indices that verifies

The number $m$ of blocks is called the size of the paving. The numbers $\alpha$ and $\beta$ are called lower and upper paving bounds. The ratio $\beta/\alpha$ gives a uniform bound on the squared condition number $\kappa^{2}(\bm{A}_{\tau})$ for each $\tau$. Note that $\alpha=0$ unless each submatrix $\bm{A}_{\tau}$ is fat.

Every partition $T$ of the rows of a matrix $\bm{A}$ has associated paving parameters $(m,\alpha,\beta)$. In a moment, we will see how these quantities play a role in the performance of the algorithm. Roughly speaking, it is best that the size $m$, the upper bound $\beta$, and the conditioning $\beta/\alpha$ of the paving are small. Later, in Sections 1.7 and 3, we will discuss what kind of bounds we can expect on the paving parameters, as well as computational methods for producing good pavings. Note that, for a row paving to be useful in our context, the cost of producing the paving must not exceed the cost of solving the least-squares problem by other means!

5. Convergence of Randomized Block Kaczmarz

The main result of this paper provides information about the convergence properties of the randomized block Kaczmarz method, Algorithm 1.3, in terms of the parameters of the row paving $T$.

Suppose $\bm{A}$ is a matrix with full column rank that admits an $(m,\alpha,\beta)$ row paving $T$. Consider the least-squares problem

Let $\bm{x}_{\star}$ be the unique minimizer, and define the residual $\bm{e}:=\bm{A}\bm{x}_{\star}-\bm{b}$. For any initial estimate $\bm{x}_{0}$, the randomized block Kaczmarz method, Algorithm 1.3, produces a sequence $\{\bm{x}_{j}:j\geq 0\}$ of iterates that satisfies

Turn to Section 2 for the proof of Theorem 1.2.

The expression (1.6) states that the block Kaczmarz method exhibits an expected linear rate of convergence until it reaches a ball about the true solution. The radius of this ball, which we call the convergence horizon, is comparable with the second term on the right-hand side of (1.6). The bracket controls the convergence rate. The minimum singular value of $\bm{A}$ affects both the rate of convergence and the convergence horizon. In each case, we prefer $\sigma_{\min}(\bm{A})$ to be as large as possible.

The properties of the row paving play an interesting role in Theorem 1.2. Curiously, the rate of convergence depends only on the upper paving bound $\beta$ and the number $m$ of blocks in the paving. On the other hand, the convergence horizon reflects the conditioning $\beta/\alpha$ of the paving. Thus, the conditioning of the paving only affects the error bound when the least-squares problem is inconsistent (i.e., $\bm{e}$ is nonzero). Nevertheless, as Section 1.4 suggests, we usually want the paving to be well conditioned to ensure that we can apply the block update rule (1.5) efficiently.

6. Simple Kaczmarz versus Block Kaczmarz

First, notice that Theorem 1.2 improves on the earlier result (1.4) for the simple Kaczmarz method. Indeed, the simple Kaczmarz method is equivalent to using a row paving with $n$ blocks, where each block contains exactly one of the $n$ rows. When $\bm{A}$ is standardized, the paving constants satisfy $\alpha=\beta=1$, and we reach the error bound

The convergence horizon $\lVert{\bm{e}}\rVert^{2}\leq n\lVert{\bm{e}}\rVert_{\infty}^{2}$, so the displayed bound beats (1.4) when $\bm{A}$ is standardized.

More generally, suppose $\bm{A}$ is a standardized matrix with an $(m,\alpha,\beta)$ row paving $T$. Let us compare the simple Kaczmarz algorithm with uniformly random control (Section 1.2) to the block Kaczmarz method, Algorithm 1.3. Both methods satisfy an error bound of the form

where the convergence rate $\rho$ and the convergence horizon $h$ depend on the choice of algorithm.

First, we compare the convergence rates of the two methods. The bounds (1.4) and (1.6) imply

because $\log(1-t)\leq-t$ when $t<1$. In other words, the simple method requires a factor $n/(\beta m)$ more iterations than the block method to achieve the same reduction in error.

Next, we examine the convergence horizons. The bounds (1.4) and (1.6) yield

We see that the convergence horizon $h_{\rm block}$ for the block method never exceeds $h_{\rm simp}$ by a factor larger than the conditioning $\beta/\alpha$ of the row paving. But $h_{\rm block}$ may be substantially smaller than $h_{\rm simp}$ if the components of the residual vector are highly nonuniform.

To make more detailed claims about the relative merits of the two algorithms, we describe two situations where we have additional information about the structure of the matrix $\bm{A}$, or a lack thereof.

First, we consider the case where the computational cost of the block update rule (1.5) is roughly comparable with the cost of the simple update rule (1.3). This situation can occur when

Each submatrix $\bm{A}_{\tau}$ admits a fast multiply; and

Each submatrix $\bm{A}_{\tau}$ is well conditioned.

See Section 4.3.1 for a numerical example where these properties hold.

In this setting, one iteration of the block method has roughly the same cost as one iteration of the standard method. As a consequence, the comparison (1.7) of the convergence rates provides a reasonable assessment of how much each algorithm reduces the error per unit of arithmetic. We see that the block algorithm really is about $n/(\beta m)$ times faster than the simple Kaczmarz method. When $\beta m$ is small in comparison with $n$, this represents a massive acceleration.

6.2. Example 2: Unstructured Submatrices

On the other hand, suppose that each submatrix $\bm{A}_{\tau}$ is unstructured and dense. Then the block method may involve much more arithmetic per iteration than the simple method. As a consequence, the comparison (1.7) is unfair to the simple method.

In this case, it is more appropriate to examine the convergence rate per epoch, the minimum number of iterations it takes the algorithm to touch each row of $\bm{A}$ once. For the simple Kaczmarz method, an epoch consists of $n$ iterations; for the block method, an epoch consists of $m$ iterations. In this setting, each algorithm requires about the same amount of arithmetic in one epoch, so we consider the per-epoch convergence rates:

We see that, in theory, the per-epoch convergence of the block method is worse, and the disadvantage increases with the upper bound $\beta$ on the paving. The best case for the block method occurs if $\beta=1$. This may happen, for example, when each block in the paving contains a single row of the matrix ($m=n$).

In practice, the block method displays better convergence behavior than this estimate suggests. The block method accrues further advantages because of subtle computational issues involving data transfer and basic linear algebra subroutines (BLASx). We discuss these points in Section 4.1, and we provide some numerical support in Section 4.3.2.

7. Existence of Good Pavings

So far, we have assumed that the matrix $\bm{A}$ comes packaged with a natural row paving $T$. For some of the applications we have in mind, this hypothesis is reasonable. Nevertheless, the block Kaczmarz method would be more versatile if we could construct row pavings for a broad class of matrices. To that end, Popa [Pop99] has developed an approach for producing a paving of a sparse matrix; see also [Pop01, Pop04]. But we can travel much farther down this road. It is an astonishing fact that every standardized matrix admits a good row paving.

Fix a number $\delta\in(0,1)$. Let $\bm{A}$ be a standardized matrix with $n$ rows. Then $\bm{A}$ admits a row paving whose parameters satisfy

Proposition 1.3 follows most directly from the recent results [Ver06, Cor. 1.5] and [Tro09, Thm. 1.2], whose provenance can be traced to the celebrated papers [BT87, BT91]. Although Proposition 1.3 is only an existential result, the literature describes several efficient algorithms for constructing row pavings. In particular, under some additional conditions, it is possible to pave a matrix by partitioning its rows at random. See Section 3 for more results and background on paving.

Before we continue, let us take a moment to explain some of the key aspects of Proposition 1.3. The main point is that the size of the paving depends only on the spectral norm of the matrix—not on the smallest singular value. As a consequence, it is possible to pave matrices with substantial null spaces!

A second point is that we can make the squared conditioning $\beta/\alpha$ as close to one as we desire. In particular, the choice $\delta=0.5$ yields $\beta/\alpha\leq 3$. This property licenses us to apply an iterative algorithm to solve least-squares problems involving the submatrices induced by the paving, as described in Section 1.4. We remark that the dependence of the paving size $m$ on the parameter $\delta$ is optimalTo verify this point, consider a large matrix $\bm{A}$ whose entries have equal magnitude and independent random signs. Use the Bai–Yin Law [BY93, Thm. 2] to estimate singular values and norms. as $\delta\to 0$.

Next, we develop some intuition about the role of the spectral norm in Proposition 1.3. Suppose that $\bm{A}$ is a standardized matrix with $n$ rows. If the lower paving bound $\alpha>0$, then a row paving $T$ of $\bm{A}$ must contain at least $n/{\rm rank}(\bm{A})$ blocks. Otherwise, $\bm{A}_{\tau}$ is rank deficient for some $\tau\in T$ simply because this submatrix has more than ${\rm rank}(\bm{A})$ rows. Now, using the fact $\lVert{\bm{A}}\rVert_{\rm F}^{2}\leq{\rm rank}(\bm{A})\lVert{\bm{A}}\rVert^{2}$, we easily verify that

This bound is sharp. (Consider the two extreme examples: a matrix with orthonormal rows and a matrix with identical rows.) Therefore, we can view the squared spectral norm as a proxy for the minimal number of blocks in a row paving whose lower bound $\alpha>0$.

We conclude that Proposition 1.3 delivers a row paving whose size falls within a logarithmic factor of optimal. One may wonder whether it is possible to remove the logarithm. For general matrices, this question remains open. It is known [And79, BHKW88, BT91] that an affirmative answer would imply the long-standing conjecture of Kadison and Singer [KS59].

8. Paved with Good Intentions

We conclude the Introduction by merging our theorem on the convergence of the block Kaczmarz method with the result on the existence of pavings.

Suppose that $\bm{A}$ is a standardized matrix with full column rank. Let $T$ be a good row paving of $\bm{A}$, as guaranteed by Proposition 1.3 with $\delta=1/2$. Under the notation of Theorem 1.2, the block Kaczmarz method, Algorithm 1.3 admits the convergence estimate

To summarize, Proposition 1.3 yields a small paving of the matrix $\bm{A}$ with exceptional conditioning. With this choice of paving, we can perform the block Kaczmarz update (1.5) quickly using an iterative least-squares algorithm. (Indeed, to apply $\bm{A}_{\tau}^{\dagger}$ with a fixed level of precision, it suffices to perform a constant number of matrix–vector multiplies with $\bm{A}_{\tau}$ and $\bm{A}_{\tau}^{*}$.) Furthermore, we see that the block Kaczmarz method converges linearly with a rate that is controlled by the condition number of the matrix $\bm{A}$, and the convergence horizon is on the same order as the size of the residual. This is essentially the best outcome one might hope for.

9. Organization

The rest of the paper has the following structure. Section 2 contains a proof of the main result, Theorem 1.2. In Section 3, we give an overview of the literature on pavings. Section 4 discusses numerical aspects of the block Kaczmarz method. We discuss related work on Kaczmarz methods and future directions in Section 5. Finally, Appendix A offers a proof of a supplemental result.

Analysis of the Randomized Block Kaczmarz Algorithm

This section contains the proof of the main result, Theorem 1.2, on the convergence of Algorithm 1.3. We commence with two simple lemmas. The first step provides a deterministic bound on how much one iteration of the algorithm reduces the error.

Instate the hypotheses and notation of Theorem 1.2. Then the error at iteration $j$ satisfies the deterministic bound

where $\tau=\tau(j)$ is the block selected at iteration $j$.

According to the update rule (1.5), block Kaczmarz computes

where we have introduced the decomposition $\bm{b}=\bm{A}\bm{x}_{\star}-\bm{e}$, restricted to the coordinates listed in $\tau$. Subtract $\bm{x}_{\star}$ from both sides to obtain

The range of $\bm{A}_{\tau}^{\dagger}$ and the range of $\mathbf{I}-\bm{A}_{\tau}^{\dagger}\bm{A}_{\tau}$ are orthogonal, so we may invoke the Pythagorean Theorem to reach

The second term on the right-hand side satisfies

where $\alpha$ is the lower bound on the row paving $T$. Combine the last two displays to wrap up. ∎

The second lemma gives us a means to average the two quantities appearing in Lemma 2.1 over a random choice of the block $\tau=\tau(j)$.

Instate the hypotheses and notation of Theorem 1.2. Suppose that $\tau$ is chosen uniformly at random from the row paving $T$. For fixed vectors $\bm{u}$ and $\bm{v}$, it holds that

The second identity emerges from a very short calculation:

which depends on the fact that the blocks of $T$ partition the components of $\bm{v}$.

Since $\bm{A}_{\tau}^{\dagger}\bm{A}_{\tau}$ is an orthogonal projector, we may apply the Pythagorean Theorem to obtain the relation

We control the remaining expectation as follows.

The second inequality depends on the bound $\sigma_{\min}^{2}(\bm{A}_{\tau}^{\dagger})=\sigma_{\max}^{-2}(\bm{A}_{\tau})\geq\beta^{-1}$. The fourth relation holds because the blocks in a paving partition the row indices of $\bm{A}$. To complete the proof, we simply combine the last two displays. ∎

The main result follows quickly once we merge the two lemmas.

First, we bound the expected error at iteration $j$ in terms of the error at iteration $j-1$. Average the bound from Lemma 2.1 over the randomness in $\tau=\tau(j)$ to reach

The second inequality follows from Lemma 2.2, with $\bm{u}=\bm{x}_{j-1}-\bm{x}_{\star}$ and $\bm{v}=\bm{e}$.

By applying this result repeatedly, we can control the expected error after $j$ iterations in terms of the initial error. Abbreviating $\gamma:=1-\sigma_{\min}^{2}(\bm{A})/(\beta m)$, we obtain the estimate

Reintroduce the value of $\gamma$ in this expression, and simplify to complete the proof. ∎

A Conversation about Pavings

As we have seen, the properties of the row paving $T$ of the matrix $\bm{A}$ have a significant effect on the behavior of the block Kaczmarz method, Algorithm 1.3. It is natural to ask if every standardized matrix $\bm{A}$ admits a good paving, and—if so—how we can exhibit such a paving.

This section summarizes the main results from the literature on row pavings, with particular attention to algorithmic techniques for constructing pavings. We focus on two methods in particular. The first approach, which is more general, extracts a well-conditioned row submatrix from $\bm{A}$ and repeats this process until the paving is complete. The second approach, which is more automatic, simply forms a random partition of the rows with an appropriate number of blocks.

For clarity, we only discuss row paving theorems for standardized matrices. For general matrices, it is often more natural to consider an alternative definition of a paving where the rows of the matrix are reweighted. For the reader’s convenience, we include some citations that address the general case.

The operator theory literature uses the term paving to refer to a partition $T=\{\tau_{1},\dots,\tau_{m}\}$ of the coordinates of a square matrix $\bm{H}$ with a zero diagonal in which each diagonal block satisfies the bound

where the parameter $\delta\in(0,1)$. We can obtain row paving results for a standardized matrix $\bm{A}$ by applying paving results for a square matrix to the hollow Gram matrix $\bm{H}=\bm{AA}^{*}-\mathbf{I}$. This approach generally leads to an estimate for the size $m$ of the paving that has an excessive dependence on the spectral norm of $\bm{A}$.

The first approach to paving relies on a type of result called a subset selection theorem. This class of result asserts that, under appropriate conditions, a matrix contains a (large) set of rows with distinguished geometric properties. Proposition 1.3 depends on the following subset selection theorem.

Fix a number $\delta\in(0,1)$. Let $\bm{A}$ be a standardized matrix with $n$ rows. Then there exists a subset $\tau$ of row indices with the properties

Proposition 3.1 follows from [Tro09, Thm. 1.2], once we track the parameter $\delta$ through the proof. This result has been attributed to Bourgain and Tzafriri [BT91], but the earliest reference seems to be Vershynin’s paper [Ver06, Cor. 1.5]. See Section 3.1.1 for further background.

Proposition 3.1 ensures that each standardized matrix contains a large set of well-conditioned rows. To construct a paving of a standardized matrix $\bm{A}$ with row indices $\{1,\dots,n\}$, we apply this result to identify a large subset $\tau_{1}$ of the row indices. We apply the same result to the set of remaining rows $\{1,\dots,n\}\setminus\tau_{1}$ to bite off another subset $\tau_{2}$, and so forth. After

steps, we have exhausted the entire matrix. This argument yields Proposition 1.3.

The paper [Tro09] contains an efficient computational method for identifying the subset $\tau$ promised by Proposition 3.1. This algorithm chooses a random set $\omega$ of rows from the matrix $\bm{A}$ with twice the cardinality of the desired subset $\tau$. Then it computes a matrix factorization of the submatrix $\bm{A}_{\omega}$ that exposes a well-conditioned subset $\tau$ of rows inside $\omega$.

The last few years have witnessed some striking advances in this area. Indeed, Spielman and Srivastava [SS12] have recently invented an elementary proof of the Restricted Invertibility Principle. Their method only involves linear algebra, and it leads to sharp constants. Youssef [You12b, Thm. 4.2] has adapted these ideas to obtain an elementary proof of the Kashin–Tzafriri theorem [KT94]. These results are appealing because they construct the required subsets using an algorithmic procedure that admits a polynomial-time implementation. See [Sri10, Nao11] for further exposition.

Vershynin [Ver01] has obtained a theory of subset selection for matrices that are not necessarily standardized. In his results, the squared Frobenius norm $\lVert{\bm{A}}\rVert_{\rm F}^{2}$ plays the role of the number $n$ of rows. Srivastava’s dissertation [Sri10, Chap. 3] contains an algorithm for (weighted) subset selection that applies to general matrices. See Youssef’s works [You12b, You12a] for the latest developments in this direction.

Finally, let us mention that subset selection theorems have applications throughout mathematics and engineering. See the paper [CT06] for a discussion and references.

2. Randomized Methods for Paving

The second approach to paving has the benefit of utmost simplicity, but it is more limited in scope. The idea is to divide the rows of the matrix into random blocks of approximately equal size. Under additional assumptions, each submatrix induced by this partition is likely to be well conditioned. To describe this idea in more detail, we first introduce the concept of a random partition.

Suppose that $\pi$ is a permutation on $\{1,2,\ldots,n\}$, chosen uniformly at random. For each $i=1,2,\ldots,m$, define the set

It is clear that $T=\{\tau_{1},\tau_{2},\ldots,\tau_{m}\}$ is a partition of $\{1,2,\ldots,n\}$ into $m$ blocks of approximately equal size. We say that $T$ is a random partition of $\{1,2,\dots,n\}$ into $m$ blocks.

For every standardized matrix, we can use a random partition to construct a paving whose upper bound $\beta$ is relatively small.

Let $\bm{A}$ be a tall A $p\times q$ matrix is tall when $p\geq q$., standardized matrix with $n$ rows. Consider a randomized partition $T$ of the row indices with $m\geq\lVert{\bm{A}}\rVert^{2}$ blocks. Then $T$ is a row paving with upper bound $\beta\leq 6\log(1+n)$, with probability at least $1-n^{-1}$.

Proposition 3.3 results from an argument based on the matrix Chernoff inequality [Tro12, Thm. 1.1] and a union bound. A model for this type of proof appears in the paper [Tro11]. We omit the details.

In contrast, if we wish to construct a paving with a nontrivial lower bound $\alpha$, we must place additional assumptions on the matrix. An example [BT91, Ex. 2.2] of Bourgain and Tzafriri implies that we must assume the rows of the matrix $\bm{A}$ are weakly correlated to obtain a random paving with $\alpha>0$. It is natural to carve out a class of matrices that meet this requirement.

Suppose $\bm{A}$ is a matrix with rows $\bm{a}_{1},\bm{a}_{2},\ldots,\bm{a}_{n}$. We say that $\bm{A}$ is incoherent when

Incoherent matrices arise, for example, in signal processing problems [DH01]. Every incoherent, standardized matrix admits a random paving with controlled lower and upper bounds.

Proposition 3.5 follows from [Tro08a, Cor. 5.2], along with some standard arguments [BT91, Tro08b]. The paper [CD12] contains superior estimates for the constants in this analysis. See Section 3.2.3 for some further background.

Random paving is a striking idea because it is almost completely automatic. Given a guarantee that the matrix $\bm{A}$ is incoherent and an estimate for the spectral norm, we can obtain a good paving of the matrix without any further computation.

Prima facie, random paving has limited applicability because the incoherence hypothesis in Proposition 3.5 is rather stringent. Fortunately, there is a fast linear transformation that can be used to convert any standardized matrix into an incoherent matrix that is nearly standardized.

The fast incoherence transform is the $n\times n$ random matrix $\bm{S}=\mathbf{F}\bm{E}$ where $\mathbf{F}$ is the $n\times n$ unitary discrete Fourier transform (DFT) and $\bm{E}$ is an $n\times n$ diagonal matrix whose entries are independent RademacherA Rademacher random variable takes the values $\pm 1$ with equal probability. random variables.

The critical fact is that the fast incoherence transform converts a large class of standardized matrices into incoherent matrices that are nearly standardized.

Suppose that $\bm{A}$ is a standardized matrix with $n$ rows whose norm satisfies

Let $\bm{S}$ be an $n\times n$ fast incoherence transform. Then the matrix $\bm{W}=\bm{SA}$ satisfies the probability bound

We do not have a reference for Proposition 3.7, but the result is probably not new. Appendix A offers a short proof based on the Hanson–Wright inequality [HW71] for Rademacher chaos.

Let us pause to examine the hypothesis (3.1). Recall that, for a standardized matrix $\bm{A}$ with $n$ rows, the squared spectral norm $\lVert{\bm{A}}\rVert^{2}$ attains its maximal value $n$ when the rows are identical. Therefore, the bound (3.1) stipulates that the rows of $\bm{A}$ must exhibit a small amount of diversity. In this case, the fast incoherence transform spreads out the diversity evenly so that no pair of columns is correlated too strongly.

2.2. Incoherence, Paving, Kaczmarz

This discussion suggests the following approach for solving the overdetermined least-squares problem (1.1) with a standardized matrix $\bm{A}$:

Apply the fast incoherence transform $\bm{S}$ to the objective function:

2.3. Related Results

The aforementioned theorems all yield the wrong dependence on the spectral norm in the size of a random row paving. The précis [Tro08a] shows how to obtain the correct quadratic dependence that is quoted in Proposition 3.5. If we were to strengthen the incoherence requirement in Proposition 3.5, it is likely that we could remove the logarithmic factor from the size of the paving by adapting an argument [BT91, Prop. 2.7] of Bourgain and Tzafriri. On the other hand, the logarithmic factor is necessary at the incoherence level we have imposed [BT91, Ex. 2.2].

The fast incoherence transform is based on ideas of Ailon and Chazelle [AC09], who use the random matrix $\bm{S}$ to perform dimension reduction. We believe that the first application of $\bm{S}$ for randomized linear algebra appears in the paper Woolfe et al. [WLRT08], where they use this transform to aid in computing matrix decompositions. See the works [AMT10, HMT11, BG12] for further results in this direction. Liberty’s dissertation [Lib09] describes other randomized maps that can play a similar role.

Numerical Aspects of Block Kaczmarz

The main goal of this paper is to study the theoretical properties of the block Kaczmarz method, but we believe that a short discussion about numerics can also provide some useful insights. In Section 4.1, we outline some of the situations where we expect the block Kaczmarz method to outperform the simple Kaczmarz algorithm. Afterward, in Section 4.2, we consider some of the questions that arise when implementing the block Kaczmarz method. Finally, in Section 4.3, we offer some simple computational examples to illustrate our main points.

It is natural to ask when it might be better to use the block Kaczmarz scheme instead of the simple Kaczmarz scheme. The answer to this question depends heavily on the specific application at hand, as well as the architecture of the computers on which the algorithm is implemented.

First, least-squares problems sometimes involve a matrix that admits a natural row paving. In this case, it may be advantageous to exploit the paving algorithmically. For example, certain signal processing applications involve multi-sampling schemes, where we collect several batches of uniform time samples of a signal. Each sample set produces a set of equations that is easy to solve. The block Kaczmarz method provides an effective way to use this structure [FS95]. See Section 4.3.1 for a related numerical example.

Second, the block Kaczmarz algorithm can be implemented more efficiently than the simple Kaczmarz algorithm in many computer architectures. This claim rests on two facts. First, data transfer now plays a major role in the cost of numerical algorithms. The block Kaczmarz algorithm is efficient in this regard, because it moves a large block of equations into working memory and operates with it for some time. In contrast, the simple Kaczmarz method repeatedly transfers new equations into working memory. Second, the block Kaczmarz algorithm can exploit high-level basic linear algebra subroutines (BLASx). Indeed, inner products dominate the arithmetic in the simple Kaczmarz method, while it is possible to implement the block Kaczmarz algorithm using matrix–vector products. As a result, block Kaczmarz relies on BLAS2, rather than BLAS1. A more detailed discussion of these issues falls outside the scope of this paper.

2. Implementing Block Kaczmarz Methods

The randomized block Kaczmarz method, Algorithm 1.3, is easy to describe, but there remain several implementation issues that require attention.

Most of the arithmetic in the algorithm occurs when we apply the pseudoinverse $\bm{A}_{\tau}^{\dagger}$ to a vector in the update (1.5). Equivalently, we must solve an (underdetermined) least-squares problem at each iteration. The appropriate numerical method depends on a wide variety of issues [Bjö96], so we cannot give a universal prescription. Here are some factors worth considering.

When the blocks in the paving are very small, a direct method based on QR decomposition or the SVD is likely to be fastest. A direct method may also be appropriate when the conditioning $\beta/\alpha$ of the paving is high and the matrix is dense.

In case the conditioning $\beta/\alpha$ of the paving is small, we recommend using an iterative method, such as CGLS, LSQR, or the Chebyshev semi-iterative method. These techniques may also be appropriate for very sparse problems. It is not necessary to run these iterations until they have converged fully, and the algorithms will benefit from warm starts provided by the convergence of the outer iteration.

2.2. Setting the Convergence Tolerance

Theorem 1.2 implies that Algorithm 1.3 can reduce the error $\lVert{\hat{\bm{x}}-\bm{x}_{\star}}\rVert_{2}^{2}$ to a level comparable with the convergence horizon. Let us examine how this fact affects our choice of the convergence tolerance.

Combine the convergence criterion $\lVert{\bm{A}\hat{\bm{x}}-\bm{b}}\rVert_{2}^{2}\leq\varepsilon^{2}$ with the decomposition $\bm{b}=\bm{A}\bm{x}_{\star}-\bm{e}$ to obtain

It follows that we must set the convergence tolerance $\varepsilon$ so that

Otherwise, we have no guarantee that the algorithm will terminate.

2.3. Checking for Convergence

The convergence criterion $\lVert{\bm{A}\hat{\bm{x}}-\bm{b}}\rVert^{2}\leq\varepsilon^{2}$ may be somewhat expensive to verify because it involves a multiply with the full matrix $\bm{A}$. As a consequence, it is better to check for convergence only on occasion. For instance, if the paving consists of $m$ blocks, we might only evaluate the convergence criterion every $m$ iterations.

2.4. Randomized Cyclic Control

Finally, in practice, the block Kaczmarz algorithm is more effective if we use a randomized control scheme different from the one we have analyzed. Recall that Algorithm 1.3 samples a uniformly random block at each iteration, independent of all previous choices. Instead, we recommend using an alternative scheme that samples blocks without replacement. We can express this method formally as follows.

For each $q=0,1,2,\dots$, draw an independent, uniformly random permutation $\pi_{q}$ on the indices $\{1,\dots,m\}$ of the blocks.

For each iteration $j\geq 1$, decompose $j=qm+r$, where $q$ and $r$ are nonnegative integers with $0\leq r<m$. Select the block $\tau(j)=\tau_{\pi_{q}(r+1)}$.

In other words, at each epoch, we cycle through all the blocks in random order.

In Section 4.3.1, we offer some numerical evidence that this alternative control scheme is more effective than the approach used in Algorithm 1.3. At present, compelling explanations for this phenomenon are lacking. See [RR12] for some discussion and conjectures.

3. Numerical Experiments

In this section, we present some numerical experiments to complement our discussions about the implementation and theoretical performance of the randomized block Kaczmarz method, Algorithm 1.3. In Section 4.3.1, we consider an example where the block method has a clear advantage over the simple method. In Section 4.3.2, we examine some other cases where the benefits of the block method are more subtle.

First, we study the situation where the matrix $\bm{A}$ comes with a natural partition of the rows into well-conditioned blocks, each admitting a fast matrix–vector multiply. As we have discussed in Section 1.6.1, the block method has a significant advantage over the simple method in this setting. Our experiments bear out this point.

We build a $300\times 100$ matrix $\bm{A}_{\rm circ}$ by stacking $15$ partial random circulant matrices:

$\bm{R}$ is the restriction to the first $20$ coordinates,

$\mathbf{F}$ is the $100\times 100$ unitary DFT, and

$\bm{E}_{i}$ is a random diagonal matrix whose entries are independent Rademacher random variables. Each $\bm{E}_{i}$ is drawn independently from the others.

This construction ensures that each $\bm{C}_{i}$ is a section of a circulant matrix with orthonormal rows. As a consequence, the pseudoinverse equals the adjoint: $\bm{C}_{i}^{\dagger}=\bm{C}_{i}^{*}$. Furthermore, we can apply either $\bm{C}_{i}$ or $\bm{C}_{i}^{*}$ to a vector quicklyRecall that an FFT or inverse FFT of length $d$ requires roughly $d\log_{2}(d)$ complex floating-point operations (flops), although the precise count depends on arithmetic properties of the number $d$. using one FFT and one inverse FFT, each of length $d$.

Using this matrix, we perform a small Matlab experiment to compare the behavior of randomized block Kaczmarz and randomized simple Kaczmarz.

Draw a random matrix $\bm{A}_{\rm circ}$ according to the recipe above, and fix it. Let $T$ be the row partition with 15 blocks induced by the circulant structure.

Let $\bm{x}_{\star}=[1,\dots,1]^{*}$, and form $\bm{b}=\bm{A}\bm{x}_{\star}$. Set the initial iterate $\bm{x}_{0}=\bm{0}$.

Apply the simple Kaczmarz method from Section 1.2 to produce iterates $\{\bm{x}_{j}^{\rm simp}:j\geq 0\}$.

Apply Algorithm 1.3 to produce iterates $\{\bm{x}_{j}^{\rm block}:j\geq 0\}$.

For each algorithm, at each iteration $j$, compute the minimum, median, and maximum value of the approximation error $\lVert{\bm{x}^{\rm alg}_{j}-\bm{x}_{\star}}\rVert_{2}$ over the 100 trials.

In this experiment, we form the matrix $\bm{A}_{\rm circ}$ in the first step and store it. To perform the update (1.3), the simple Kaczmarz method loads the required row from memory without doing any extra computation, so the cost of the simple update rule is just $4d$ complex flops, where $d=100$. The block Kaczmarz method uses the structure of the circulant blocks, as described in the previous paragraph, to perform the update (1.5) with about $4d\log_{2}(d)+4d$ complex flops.

In Figure 1 [left panel], we plot the approximation error as a function of the number of flops expended by each algorithm. The heavy lines indicate the median error over 100 trials, while the minimum and maximum errors describe the boundaries of the shaded region. The block algorithm reduces the error to $10^{-11}$ in about $1.6\times 10^{6}$ complex flops, while the simple algorithm requires about $3.2\times 10^{7}$ complex flops to achieve the same result. This amounts to a 20-fold reduction in the amount of arithmetic!

In Section 4.2.4, we claim that the block Kaczmarz algorithm is more effective when we use an alternative control scheme that samples blocks without replacement. To illustrate this point, let us perform an experiment with the block circulant matrix using the same methodology as above. Figure 1 [right panel] compares the performance of the block Kaczmarz method when we sample blocks with and without replacement. This chart shows that sampling without replacement reduces the amount of computation by about 15%. In other experiments, we have seen even more substantial improvements.

3.2. Unstructured Submatrices

Next, we consider some examples where we cannot exploit the structure of the submatrices to accelerate the block Kaczmarz method. Although our theory does not yield higher rates of convergence, the numerical evidence still suggests that the block Kaczmarz method offers a decisive improvement over the simple Kaczmarz method.

We perform an experiment with this matrix that follows the same methodology as in Section 4.3.2. To implement the block Kaczmarz update rule (1.5), we use the CGLS algorithm to solve the underdetermined least-squares problem. We use warm starts, and we halt the least-squares iteration before it has fully converged to control the computational cost. Our code counts the actual number of (real) flops during each iteration, and it measures the actual CPU time that is expended.

Figure 2 shows three different views of the data we collected during this experiment. The left panel shows the approximation error for the simple and block Kaczmarz method as a function of the number of epochs elapsed. As we expect from the discussion in Section 1.6.2, the two algorithms show a similar rate of convergence in this view. In fact, the block method does somewhat better than the theory predicts. In the center panel, we chart the approximation error as a function of the number of flops performed. The simple algorithm requires about $5\times 10^{6}$ flops to achieve an error of $10^{-11}$, while the block algorithm takes $2\times 10^{7}$ flops to reach the same point. Thus, the block method needs four times as much arithmetic. But the story is not over. The right panel displays the approximation error as a function of CPU time. We discover that our implementation of the block method is 10 times faster than the simple method! It is dangerous to draw conclusions about the efficiency of an algorithm from the behavior of a Matlab script. Nevertheless, we regard this experiment as limited evidence of the computational advantages of the block method that we outlined in Section 4.1.

We perform the same experiment for $\bm{A}_{\rm coh}$ as we did with $\bm{A}_{\rm rdm}$ to compare the behavior of the simple Kaczmarz method and the randomized Kaczmarz method. Figure 3 shows the results of this trial. We see that the simple Kaczmarz method scarcely reduces the error at all, while the block method achieves a healthy rate of convergence. The paper [NW12] provides an analysis of this example in the case where each block contains two rows, but we do not yet have a complete explanation for the performance of Algorithm 1.3 when the blocks are larger.

Related Work and Future Directions

We conclude with a short discussion about recent research on randomized Kaczmarz methods and block variants of the simple Kaczmarz method. Finally, we offer some musings about opportunities for further research.

The Kaczmarz method was originally introduced in the paper [Kac37]. It was reinvented by researchers in tomography [GBH70] under the appellation “algebraic reconstruction technique” (ART). See Byrne’s book [Byr08] for a contemporary summary of this literature.

The classical variants of the Kaczmarz method rely on deterministic mechanisms for selecting a row at each iteration. Indeed, the simplest version just cycles through the rows in order. It has long been known that the cyclic control scheme performs badly when the rows are arranged in an unhappy order [HS78]. The literature contains empirical evidence [HM93] that randomized control mechanisms may be more effective, but until recently there was no compelling theoretical analysis to support this observation.

The paper [SV09b] of Strohmer and Vershynin is significant because it provides the first explicit convergence proof for a randomized variant of the Kaczmarz algorithm. This work establishes that a randomized control scheme leads to an expected linear convergence rate, which can be written in terms of geometric properties of the matrix. In contrast, deterministic convergence analyses that appear in the literature often lead to expressions, e.g., [XZ02, Eqn. (1.2)], whose geometric meaning is not evident.

In the wake of Strohmer and Vershynin’s work [SV09b], several other researchers have written about randomized versions of the Kaczmarz scheme and related topics. In particular, Needell demonstrates that the randomized Kaczmarz method converges, even when the linear system is inconsistent [Nee10]. Zouzias and Freris [ZF12] exhibit a randomized procedure, based on ideas from [Pop98], that can reduce the size of the residual $\bm{e}$. Leventhal and Lewis [LL10] provide an analysis of a randomized iteration for solving least-squares problems with polyhedral constraints, while Richtárik and Takáč have extended these ideas to more general optimization problems [RT11]. Some other references include [EN11, RT11, CP12, NW12].

2. Block Kaczmarz Methods

The block Kaczmarz update rule (1.5) we are studying is originally due to Elfving [Elf80, Eqn. (2.2)]. This update is a special case of a general framework due to Eggermont et al. [EHL81]. Byrne describes a number of other block Kaczmarz methods in his book [Byr08, Chap. 9].

By now, there is an extensive literature on the convergence behavior of block projection methods, a class of algorithms that includes block Kaczmarz schemes. In particular, we call out the work of Xu and Zikatanov [XZ02], which contains a refined convergence analysis that applies to a wide range of algorithms. Nevertheless, to our knowledge, Algorithm 1.3 is the only block Kaczmarz method that offers an (expected) linear rate of convergence that depends explicitly on geometric properties of the system matrix $\bm{A}$ and its submatrices $\bm{A}_{\tau}$.

Most of the literature on block Kaczmarz methods assumes that a partition $T$ of the rows of the matrix is provided as part of the problem data. We are aware of some research on the prospects for partitioning a matrix in a manner that is favorable for block Kaczmarz methods. In particular, Popa [Pop99] has introduced an algorithm for partitioning a sparse matrix so that each block contains mutually orthogonal rows. Popa has pursued this idea in a sequence of papers, including [Pop01, Pop04]. We believe that our work is the first to recognize the natural connection between the paving literature and the block Kaczmarz method.

3. Some Future Directions

There are a number of interesting open questions connected with Algorithm 1.3. First, empirical experiments make it clear that a control strategy based on sampling without replacement is far more effective than our strategy based on sampling with replacement. At present, there is no compelling explanation for this phenomenon (but see [RR12]). Second, we think that other types of block Kaczmarz update rules [Byr08, Chap. 9] would also submit to a convergence analysis similar to the proof of Theorem 1.2. Third, it might be worthwhile to extend the argument of Zouzias and Freris [ZF12] to develop a method with a smaller convergence horizon.

Block algorithms for numerical linear algebra are almost as old as the field itself. Gauss himself suggested a block version of his algorithm for solving linear systems [Gau95, Ben09]. Many other algorithms for numerical linear algebra [GVL96] and least-squares problems [Bjö96] admit natural block variants. The field of optimization also contains a wide variety of block methods, such as [AC89, Tse93].

We believe that row pavings can also play a role in the development and analysis of other types of block algorithms. For instance, it is straightforward to develop a block version of the randomized Gauss–Seidel algorithm introduced in [LL10]. Using the techniques in this paper, we can easily bound the rate of convergence for this algorithm in terms of the properties of a column paving of the system matrix. It would be interesting to pursue this example and others.

Appendix A The Fast Incoherence Transform

The goal of this Appendix is to prove Proposition 3.7, which states that the fast incoherence transform makes a standardized matrix incoherent with high probability. To that end, suppose that $\bm{A}$ is a matrix with $n$ unit-norm rows, denoted $\bm{a}_{1},\dots,\bm{a}_{n}$. Recall that the fast incoherence transform is the matrix $\bm{S}=\mathbf{F}\bm{E}$, where $\mathbf{F}$ is the $n\times n$ unitary DFT and $\bm{E}$ is a diagonal matrix whose entries $\xi_{1},\dots,\xi_{n}$ are independent Rademacher random variables.

where we write ${\rm f}_{ij}$ for the $(i,j)$ entry of the DFT matrix. (We use the convention that the inner product is antilinear in the second coordinate.) This expression allows us to calculate the expectation of the Gram matrix with ease. Since the rows $\mathbf{f}_{1},\dots,\mathbf{f}_{n}$ of the DFT form an orthonormal family,

We can rewrite $Y$ is a more symmetric manner:

We see that the random variable $Y$ is a symmetric, homogeneous, second-order Rademacher chaos. We can bound the probability that $Y$ is large by invoking a result of Hanson and Wright [HW71]; see [FR12, Chap. 8] for a modern proof.

Consider the chaos variable $Y$ defined in (A.1). Let $\bm{Y}$ be the Hermitian matrix whose $(j,k)$ entry is $y_{jk}$ and whose diagonal entries are zero. Then

To apply this result, we need to obtain bounds for the Frobenius norm and spectral norm of the matrix $\bm{Y}$. Note that

and $\operatorname{diag}(\cdot)$ converts a vector into a diagonal matrix. We may now calculate the required norms of $\bm{Y}$. First,

The strict inequality holds because $\bm{AA}^{*}$ has a unit diagonal, which the identity matrix cancels off. The final bound follows from the interpolation $\lVert{\bm{P}}\rVert_{\rm F}^{2}\leq\operatorname*{trace}(\bm{P})\lVert{\bm{P}}\rVert$, valid when $\bm{P}$ is positive semidefinite.

Next, let us instate the hypothesis (3.1) from Proposition 3.7:

Introduce this assumption into the bounds (A.2) and (A.3), and apply the Hanson–Wright inequality, Proposition A.1, to reach

In other words, it is unlikely that any single pair of rows from $\bm{W}=\bm{SA}$ has an inner product much different from its expectation.

Therefore, with high probability $\bm{W}$ is an incoherent matrix that is nearly standardized.

Acknowledgments

We would like to thank Michael Mahoney, Ben Recht, Thomas Strohmer, Steve Wright for helpful discussions about randomized linear algebra and numerical experiments. Roman Vershynin provided insight on the random paving literature. Michael McCoy explained advanced plotting techniques in Matlab, and Margot Stokol shared her expertise on color theory. JAT was supported in part by ONR awards N00014-08-1-0883 and N00014-11-1002, AFOSR award FA9550-09-1-0643, DARPA award N66001-08-1-2065, and a Sloan Research Fellowship.

References