Paved with Good Intentions: \\ Analysis of a Randomized Block Kaczmarz Method
Deanna Needell and Joel A. Tropp
The block Kaczmarz method is an iterative scheme for solving overdetermined least-squares problems. At each step, the algorithm projects the current iterate onto the solution space of a subset of the constraints. This paper describes a block Kaczmarz algorithm that uses a randomized control scheme to choose the subset at each step. This algorithm is the first block Kaczmarz method with an (expected) linear rate of convergence that can be expressed in terms of the geometric properties of the matrix and its submatrices. The analysis reveals that the algorithm is most effective when it is given a good row paving of the matrix, a partition of the rows into well-conditioned blocks. The operator theory literature provides detailed information about the existence and construction of good row pavings. Together, these results yield an efficient block Kaczmarz scheme that applies to many overdetermined least-squares problem. %The Kaczmarz method is an iterative method for solving overcomplete linear systems of equations $\mtx{A}\vct{x}=\vct{b}$. The randomized version of the Kaczmarz method put forth by Strohmer and Vershynin iteratively projects onto a randomly chosen solution space given by a single row of the matrix $\mtx{A}$. In this paper we analyze a block version of the method which projects the current iterate onto a subspace given by a block $\mtx{A}_\tau$ of the matrix rather than a single row. Our approach utilizes a paving of the matrix $\mtx{A}$ to guarantee exponential convergence, and suggests that paving yields a significant improvement in performance in certain regimes. We show that the paving can be constructed randomly in many cases which leads to a completely randomized block Kaczmarz algorithm. Numerical experiments suggest that the block version with paving can indeed lead to advantages in practice.
NT-rev.tex
\documentclass[11pt,reqno,final]{amsart}
%\documentclass[11pt,twocolumn]{article}
\usepackage{amsmath, amssymb, epsfig}
\usepackage{mathrsfs}
\usepackage{color}
\usepackage{alg}
\usepackage{url}

\usepackage[font=small,margin=10pt,labelfont={bf},labelsep={space}]{caption}
\usepackage{subfig}
\usepackage{wrapfig}

\usepackage[scaled]{helvet} % I like Helvetica for sf
\usepackage{fourier} % Euler and Utopia font; bolder than your
%average bear
%\usepackage{MinionPro}
\usepackage{bm} % boldmath must be called after the package

\renewcommand{\thealgorithmfloat}{\thesection.\arabic{algorithmfloat}}

\makeatletter
\def\algspacing{\alg@unmargin}
\makeatother

\newcommand{\algfor}[1]{\textbf{for} #1\\\algbegin}
%\newcommand{\alguntill}[1]{\textbf{until} {#1}\\\algbegin}

%%%%% Page setup

\usepackage[margin=1in]{geometry}

%
%\evensidemargin=0in
%\oddsidemargin=0in
%\textwidth=6.5in
%\topmargin=-0.33in
%\headheight=0.25in
%%\headsep=0.25in
%\textheight=9in

\newlength{\algorithmwidth}
\algorithmwidth=0.95\textwidth

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{fact}[theorem]{Fact}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{question}[theorem]{Question}

\theoremstyle{remark}
\newtheorem*{remark}{Remark}

\numberwithin{theorem}{section}
\numberwithin{equation}{section}

\newcommand{\<}{\left\langle}
\renewcommand{\>}{\right\rangle}
\newcommand{\bigO}{\mathrm{O}}
\def \remarks {\noindent {\bf Remarks.}\ \ }
\def \example {\vspace{0.5cm} \noindent {\bf  Example.}\ \ }

%\newcommand{\enorm}[1]{\norm{#1}_2}

\DeclareMathOperator*{\supp}{supp}
\DeclareMathOperator*{\Dspan}{span}
\DeclareMathOperator*{\trace}{trace}
\DeclareMathOperator*{\range}{range}
\DeclareMathOperator*{\codim}{codim}
\DeclareMathOperator*{\Null}{Null}
\DeclareMathOperator*{\mn}{mn}
\DeclareMathOperator*{\Idg}{Id_\Gamma}
\DeclareMathOperator*{\argmin}{arg min} 
\DeclareMathOperator*{\argmax}{arg max} 

\newcommand{\notate}[1]{\textcolor{red}{\textbf{[#1]}}}
\newcommand{\defby}{\overset{\mathrm{\scriptscriptstyle{def}}}{=}}
\def \Re {\mathbb{R}}
\def \C {\mathbb{C}}
\def \P {\mathbb{P}}
\def \E {\mathbb{E}}
\def \A {\mathcal{A}}
\def \M {\mathcal{M}}
\def \G {\mathcal{G}}
\def \S {\mathcal{S}}
\def \H {\mathcal{H}}
\def \R {\mathcal{R}}
\def \V {\mathcal{V}}
\def \U {\mathcal{U}}
\def \W {\mathcal{W}}
\newcommand{\e}{\mathrm{e}} % exponential
\def \eps {\varepsilon}
%\def \< {\langle}
%\def \> {\rangle}
\def \sign {{\rm sign}}
\def \Span {{\rm span}}
\def \rank {{\rm rank }}

\newcommand{\vct}[1]{\bm{#1}}
\newcommand{\mtx}[1]{\bm{#1}}

%\newcommand{\pr}[2]{\langle {#1} , {#2} \rangle}


\newcommand{\ip}[2]{\left\langle {#1}, \, {#2} \right\rangle}
\newcommand{\norm}[1]{\lVert {#1} \rVert}
%\newcommand{\smnorm}[1]{\bigg\Vert {#1} \bigg\Vert}
\newcommand{\normsq}[1]{\norm{#1}^2}
\newcommand{\enorm}[1]{\norm{#1}_2}
%\newcommand{\smenorm}[1]{\smnorm{#1}_2}
\newcommand{\enormsq}[1]{\enorm{#1}^2}
%\newcommand{\smenormsq}[1]{\smenorm{#1}^2}
\newcommand{\fnorm}[1]{\norm{#1}_{\rm F}}
\newcommand{\fnormsq}[1]{\fnorm{#1}^2}
\newcommand{\infnorm}[1]{\norm{#1}_{\infty}}
\newcommand{\infnormsq}[1]{\infnorm{#1}^2}
\newcommand{\Expect}{\E}
\newcommand{\abs}[1]{\left\vert {#1} \right\vert}
\newcommand{\abssq}[1]{\abs{#1}^2}
\newcommand{\pinv}{\dagger}
\newcommand{\cnst}[1]{\mathrm{#1}}
\newcommand{\Id}{\mathbf{I}}
\newcommand{\Prob}[1]{\mathbb{P}\left\{ #1 \right\}}
\newcommand{\diag}{\operatorname{diag}}

\renewcommand{\baselinestretch}{1.05}

\begin{document}

%\setlength{\parindent}{0in}
%\parskip 7.2pt


\title[Paved with Good Intentions]{Paved with Good Intentions: \\ Analysis of a Randomized Block Kaczmarz Method}
\author{Deanna Needell and Joel A. Tropp}
\date{17 August 2012}




\begin{abstract}
The block Kaczmarz method is an iterative scheme for solving overdetermined least-squares problems.  At each step, the algorithm projects the current iterate onto the solution space of a subset of the constraints.  This paper describes a block Kaczmarz algorithm that uses a randomized control scheme to choose the subset at each step.  This algorithm is the first block Kaczmarz method with an (expected) linear rate of convergence that can be expressed in terms of the geometric properties of the matrix and its submatrices.  The analysis reveals that the algorithm is most effective when it is given a good \emph{row paving} of the matrix, a partition of the rows into well-conditioned blocks.  The operator theory literature provides detailed information about the existence and construction of good row pavings.  Together, these results yield an efficient block Kaczmarz scheme that applies to many overdetermined least-squares problem.
%The Kaczmarz method is an iterative method for solving overcomplete linear systems of equations $\mtx{A}\vct{x}=\vct{b}$.  The randomized version of the Kaczmarz method put forth by Strohmer and Vershynin iteratively projects onto a randomly chosen solution space given by a single row of the matrix $\mtx{A}$.  In this paper we analyze a block version of the method which projects the current iterate onto a subspace given by a block $\mtx{A}_\tau$ of the matrix rather than a single row.  Our approach utilizes a paving of the matrix $\mtx{A}$ to guarantee exponential convergence, and suggests that paving yields a significant improvement in performance in certain regimes.  We show that the paving can be constructed randomly in many cases which leads to a completely randomized block Kaczmarz algorithm.  Numerical experiments suggest that the block version with paving can indeed lead to advantages in practice.
\end{abstract}

\maketitle

\section{Introduction} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\notate{Rewrite}
The Kaczmarz method~\cite{Kac37:Angenaeherte-Aufloesung} 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~\cite{SS87:Incorporation-A-Priori,CFMSS92:New-Variants,FS95:Kaczmarz-Based-Approach,Nat01:Mathematics-Computerized}.  At each iteration, the basic Kaczmarz method makes progress by enforcing a single constraint, while the block Kaczmarz method~\cite{Elf80:Block-Iterative-Methods} 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 \emph{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. % (Section~\ref{sec:inc-pave-kacz}).

%The Kaczmarz method~\cite{Kac37:Angenaeherte-Aufloesung} is a popular iterative solver for overdetermined least-squares problems.  At each step, the Kaczmarz method enforces a single equation.  Many practitioners use block Kaczmarz methods, which enforce many equations at once.  Because of their simplicity and performance, these schemes arise in applications ranging from image reconstruction to digital signal processing~\cite{Nat01:Mathematics-Computerized,SS87:Incorporation-A-Priori,CFMSS92:New-Variants}.
%The goal of this paper is to describe how the performance of a certain block Kaczmarz
%method depends on the metric properties of the data matrix.  We also argue that
%the blocking technique we consider, called a \emph{paving},
%has applications for solving other computational problems.
%\notate{First paragraph is still very weak.}


%The Kaczmarz method~\cite{K37:Angena} is a popular iterative solver of consistent overdetermined systems of linear equations $\mtx{A}\vct{x}=\vct{b}$.  Because of its simplicity and performance, the method and its derivatives are used in a range of applications from 
%image reconstruction to digital signal processing~\cite{CFMSS92,N86:TheMath,SS87:Applications}.  The method performs a series of orthogonal projections and iteratively converges to the solution of the system of equations.  It is therefore computationally feasible even for very large and overdetermined systems.  

\subsection{Standing Assumptions and Notation} \label{sec:assumptions}

Let $\mtx{A}$ be a real or complex $n \times d$ matrix with full column rank, and suppose that $\vct{b}$ is a vector with dimension $n$.
Consider the overdetermined least-squares problem
\begin{equation} \label{eqn:basic-ls}
\mathrm{minimize} \quad \enormsq{ \mtx{A} \vct{x} - \vct{b} }.
\end{equation}
The symbol $\norm{\cdot}_p$ refers to the $\ell_p$ vector norm for $p \in [1, \infty]$.
We write $\vct{x}_{\star}$ for the unique minimizer of~\eqref{eqn:basic-ls}, and we introduce the residual vector $\vct{e} := \mtx{A} \vct{x}_{\star} - \vct{b}$.

To streamline our discussion, we will often assume that each \emph{row} $\vct{a}_i$ of the matrix $\mtx{A}$ shares the same $\ell_2$ norm:
\begin{equation} \label{eqn:std-matrix}
\enorm{ \vct{a}_i } = 1
\quad\text{for each $i = 1, \dots, n$.}
\end{equation}
We say that $\mtx{A}$ is \emph{standardized} when~\eqref{eqn:std-matrix} is in force.
%Some of the results in this work hold for general matrices, but the details are significantly more complicated.

%\notate{Cut eigenvalue def if it's not used.}
The spectral norm is denoted by $\norm{\cdot}$, while $\fnorm{\cdot}$ 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 $\mtx{W}$, we arrange the singular values as follows.
$$
\sigma_{\max}(\mtx{W}) := \sigma_1(\mtx{W}) \geq \sigma_2(\mtx{W})
	\geq \dots \geq \sigma_{\min\{p, q\}}(\mtx{W}) =: \sigma_{\min}(\mtx{W}).
$$
The minimum singular value is positive if and only if $\mtx{WW}^*$ or $\mtx{W}^* \mtx{W}$
is nonsingular.  We define the \emph{condition number} $\kappa(\mtx{W}) := \sigma_{\max}(\mtx{W}) / \sigma_{\min}(\mtx{W})$.  The dagger $\pinv$ denotes the Moore--Penrose pseudoinverse.  When $\mtx{W}$ has full row rank, its pseudoinverse is determined by the formula $\mtx{W}^\pinv := \mtx{W}^* (\mtx{WW}^*)^{-1}$.



% are both singular.

%define the minimum singular value of a matrix to be the smallest \emph{nontrivial} singular value.  More precisely, when $\mtx{B}$ has dimension $p \times q$, we set
%$$
%\sigma_{\min}(\mtx{B}) := \sigma_{\min\{p, q\}}(\mtx{B}).
%$$
%where the singular values $\sigma_1 \geq \sigma_2 \geq \dots$ of a matrix are arranged in weakly decreasing order.


\subsection{The Simple Kaczmarz Method} \label{sec:simple-kacz}

%We begin with a short discussion of the Kaczmarz method and some recent results on its performance.

The Kaczmarz method is an iterative algorithm that produces an approximation to the minimizer $\vct{x}_{\star}$ of the least-squares problem~\eqref{eqn:basic-ls}.  The method commences with an arbitrary guess $\vct{x}_0$ for the solution.  At the $j$th iteration, we select a row index $t = t(j)$ of the matrix $\mtx{A}$, and we project the current iterate $\vct{x}_{j-1}$ onto the solution space of the equation $\ip{ \vct{a}_t }{ \vct{x} } = b_t$.  That is,
\begin{equation} \label{eqn:simple-kacz}
\vct{x}_j = \vct{x}_{j-1} + \frac{b_t - \ip{ \vct{a}_t }{ \vct{x}_{j-1} }}
{\enormsq{\vct{a}_t}} \, \vct{a}_t.
\end{equation}
This process continues until it triggers an appropriate convergence criterion.

%Given that $\mtx{A}$ is standardized, the update takes the form
%\begin{equation} \label{eqn:simple-kacz}
%\vct{x}_k = \vct{x}_{k-1} + \left( b_i - \ip{ \vct{a}_i }{ \vct{x}_{k-1} } \right) \vct{a}_i.
%\end{equation}

%%%% General form of update
%
%\begin{equation} \label{eqn:simple-kacz}
%\vct{x}_k = \vct{x}_{k-1} + \frac{b_i - \ip{ \vct{a}_i }{ \vct{x}_{k-1} }}
%{\enormsq{\vct{a}_i}} \, \vct{a}_i.
%\end{equation}

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 \emph{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~\cite{SV09:Randomized-Kaczmarz}.
%Recently, Strohmer and Vershynin have proposed and analyzed a randomized Kaczmarz algorithm~.
When $\mtx{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~\eqref{eqn:simple-kacz}.
%They provide a short, elegant proof that
The paper~\cite{SV09:Randomized-Kaczmarz} provides a short, elegant proof that
this iteration converges at an expected linear rate to the solution $\vct{x}_{\star}$ of a consistent least-squares problem (i.e., where the residual $\vct{e}$ is zero).

%In the paper~\cite{SV09:Randomized-Kaczmarz}, Strohmer and Vershynin propose and analyze a Kaczmarz algorithm based on a randomized control scheme.  At step $k$, independently of all previous random choices, their method selects a random row index $i(k)$ with probability proportional to the squared $\ell_2$ norm of the row.  That is,
%$$
%i(k) = j \quad\text{with probability}\quad \frac{\enormsq{\vct{a}_j}}{\fnormsq{\mtx{A}}}
%\quad\text{for $j = 1, \dots, n$}.
%$$
%Here, $\fnorm{\cdot}$ denotes the Frobenius norm.
%Strohmer and Vershynin provide a short, elegant proof that this iteration converges at an expected linear rate to the solution $\vct{x}_{\star}$ of a consistent least-squares problem (i.e., the case where the residual $\vct{e}$ is zero).

Needell~\cite{Nee10:Randomized-Kaczmarz} has extended the argument of~\cite{SV09:Randomized-Kaczmarz} to the case of an \emph{inconsistent} least-squares problem.
For a standardized matrix $\mtx{A}$, Needell's error estimate reads
\begin{equation}\label{eqn:Needell-rate}
\Expect \enormsq{ \vct{x}_j - \vct{x}_\star }
	\ \leq \ \left[ 1 - \frac{\sigma^2_{\min}(\mtx{A})}{n} \right]^j \enormsq{ \vct{x}_0 - \vct{x}_\star }
	\ + \ \frac{n \infnormsq{\vct{e}}}{\sigma_{\min}^2(\mtx{A})}.
\end{equation}
In words, the randomized Kaczmarz method converges in expectation at a linear rate%
\footnote{Mathematicians often use the term \emph{exponential convergence} for the concept numerical analysts call \emph{linear convergence}.}
until it reaches a fixed ball about the true solution $\vct{x}_{\star}$, at which point the error may cease to decay.  The radius of this ball roughly equals the second term in~\eqref{eqn:Needell-rate}, while
the convergence rate is controlled by the bracket.
When the residual $\vct{e}$ is zero, the bound~\eqref{eqn:Needell-rate} reduces to the error estimate from~\cite{SV09:Randomized-Kaczmarz}.

When $\mtx{A}$ is an $n\times d$ standardized matrix whose columns are well conditioned, the minimum singular value $\sigma^2_{\min}(\mtx{A}) \geq \textrm{const} \cdot n / d$.  In this case, the error bound~\eqref{eqn:Needell-rate} simplifies to
$$
\Expect \enormsq{ \vct{x}_j - \vct{x}_\star }
	\ \lesssim \ \left[ 1 - \frac{\textrm{const}}{d} \right]^j \enormsq{ \vct{x}_0 - \vct{x}_\star }
	\ + \ \frac{d}{\textrm{const}} \infnormsq{\vct{e}}.
$$
It follows that $\bigO(d)$ iterations of the Kaczmarz method suffice to reduce the error by a constant fraction, provided that the squared error is substantially larger than $d \infnormsq{\vct{e}}$.

%Needell~\cite{N10:rknoisy} extends the argument of~\cite{SV09:Randomized-Kaczmarz} to the case of an \emph{inconsistent} least-squares problem.
%For a standardized matrix $\mtx{A}$, Needell's error estimate reads
%\begin{equation}\label{eqn:Needell-rate}
%\Expect \enormsq{ \vct{x}_j - \vct{x}_\star }
%	\leq \left[ 1 - \frac{1}{R^2} \right]^j \enormsq{ \vct{x}_0 - \vct{x}_\star }
%	\ + \ R^2 \cdot \infnormsq{\vct{e}}
%\end{equation}
%where $R$ denotes the \emph{scaled condition number}
%\begin{equation} \label{eqn:scaled-cond}
%R^2 := \frac{n}{\sigma_{\min}^2(\mtx{A})}.
%\end{equation}
%When the residual $\vct{e}$ is zero, the bound~\eqref{eqn:Needell-rate} reduces to the error estimate from~\cite{SV09:Randomized-Kaczmarz}.

%%% Deanna's full estimate.
%
%\begin{equation}\label{eqn:Needell-rate}
%\mathbb{E} \enormsq{ \vct{x}_k - \vct{x}_\star }
%	\leq \left[ 1 - \frac{1}{R^2} \right]^k \enormsq{ \vct{x}_0 - \vct{x}_\star }
%	\ + \ R^2 \cdot \max\nolimits_i \frac{\abssq{e_i}}{\enormsq{\vct{a}_i}}
%\end{equation}

%In words, the randomized Kaczmarz method converges at a linear rate until it reaches a fixed ball about the true solution $\vct{x}_{\star}$, at which point the error may cease to decay.  The radius of this ball roughly equals the second term in~\eqref{eqn:Needell-rate}.
%, which reflects the size of the residual $\vct{e}$.

%When the scaled condition number $R$ is small, the algorithm converges quickly and reaches a small radius about the true solution; when $R$ is larger, the rate tends to slow and the final error may increase.

%The bound~\eqref{eqn:Needell-rate} yields a reasonably accurate assessment of how the Strohmer--Vershynin algorithm performs in practice.  



%Given a vector $\vct{b}$ and an $n \times p$ matrix $\mtx{A}$ with rows $\vct{a}_1, \vct{a}_2, \ldots \vct{a}_n$, the algorithm begins with an initial estimate $x_0$ and cyclically projects the estimation onto each of the solution spaces.  This process can be described as follows:
%
%$$
%\vct{x}_j = \vct{x}_{j-1} + \frac{\vct{b}[i] - \langle \vct{a}_i, \vct{x}_{j-1}\rangle}{\|\vct{a}_i\|_2^2}\vct{a}_i,
%$$
%
%where $\vct{b}[i]$ denotes the $i$th coordinate of $\vct{b}$, $\vct{x}_j$ is the estimation in the $j$th iteration, $\|\cdot\|_2$ denotes the usual $\ell_2$ vector norm, and $i = (j \mod n) + 1$ cycles through the rows of $\mtx{A}$. 
%
%Since the method cycles through the rows of $\mtx{A}$, the performance of the algorithm may depend heavily on the ordering of these rows.  A poor ordering may lead to very slow convergence.  To overcome this obstacle, one can select the row $\vct{a}_i$ at random to improve the convergence rate~\cite{HM93:Algebraic,N86:TheMath}.  Strohmer and Vershynin proposed and analyzed a method which selects a given row with probability proportional to its Euclidean norm~\cite{SV09:Randomized-Kaczmarz,SV06:Arandom}.  They show that with this selection strategy, the randomized Kaczmarz method has an expected exponential convergence rate to the unique solution $\vct{x}_\star$:
%
%\begin{equation}\label{eqn:RVrate}
%\mathbb{E}\|\vct{x}_j - \vct{x}_\star\|_2^2 \leq \left(1 - \frac{1}{R}\right)^j\|\vct{x}_0 - \vct{x}_\star\|_2^2,
%\end{equation}
%
%where $R$ is the scaled condition number, $R = \|\mtx{A}^{-1}\|^2\|\mtx{A}\|_F^2$, $\|\cdot\|_F$ denotes the Frobenius norm, and $\|\mtx{A}^{-1}\| \defby \inf\{M : M\|\mtx{A}\vct{x}\|_2 \geq \|\vct{x}\|_2 \text{ for all }x\}$ is well-defined since $\mtx{A}$ has full column rank.

%When the system $\mtx{A}\vct{x}_\star = \vct{b}$ is perturbed by noise and no longer consistent, the randomized Kaczmarz method still provides expected exponential convergence down to the error threshold~\cite{N10:rknoisy}.  This convergence rate~\eqref{eqn:RVrate} is essentially independent of the number of rows of $\mtx{A}$ and shows that for well-conditioned matrices, the randomized Kaczmarz method converges to the solution in just $\bigO(p)$ iterations~\cite{SV09:Randomized-Kaczmarz}.  This yields an overall runtime of $\bigO(p^2)$ which is much superior to others such as $\bigO(np^2)$ for Gaussian elimination.  There are also cases where randomized Kaczmarz even outperforms the conjugate gradient method, see the discussion in ~\cite{SV09:Randomized-Kaczmarz} for details.

\subsection{The Block Kaczmarz Method} \label{sec:block-kacz}

%\cite[Chap.~9]{Byr08:Applied-Iterative}

In some situations~\cite{EHL81:Iterative-Algorithms}, practitioners prefer to use a block version of the Kaczmarz method to solve the least-squares problem~\eqref{eqn:basic-ls}.  We consider a formulation due to Elfving~\cite{Elf80:Block-Iterative-Methods}.  This procedure begins with an initial guess $\vct{x}_0$ for the solution.  At each iteration $j$, we select a subset $\tau = \tau(j)$ of the row indices of $\mtx{A}$, and we project the current iterate $\vct{x}_{j-1}$ onto the solution space of $\mtx{A}_{\tau} \vct{x} = \vct{b}_{\tau}$, the set of equations listed in $\tau$.  That is,
\begin{equation}\label{eqn:block-kacz}
\vct{x}_j = \vct{x}_{j-1} + (\mtx{A}_\tau)^\pinv(\vct{b}_\tau - \mtx{A}_\tau\vct{x}_{j-1}).
\end{equation}
This process continues until it has converged.
We have written $\mtx{A}_{\tau}$ for the row submatrix of $\mtx{A}$ indexed by $\tau$, while $\vct{b}_{\tau}$ is the subvector of $\vct{b}$ with components listed in $\tau$.  We assume that the row submatrix $\mtx{A}_{\tau}$ is fat%
\footnote{A $p \times q$ matrix is \emph{fat} when $p \leq q$.},
so the pseudoinverse~\eqref{eqn:block-kacz} returns the solution to an underdetermined least-squares problem.

%In some situations~\cite{FS95:Kaczmarz-Based-Approach,Pop98:Extensions-Block-Projections,Byr08:Applied-Iterative}, practitioners prefer to use a block version of the Kaczmarz method to solve the least-squares problem~\eqref{eqn:basic-ls}.

%This type of algorithm has the following general form.

%\notate{At this point, still ambiguous whether blocks are fat or tall!}

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 $\mtx{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~\ref{alg:randomized-block-kaczmarz}.
%, which appears on page~\pageref{alg:randomized-block-kaczmarz}.
We postpone a detailed discussion on implementation to Section~\ref{sec:implementation}.

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~\cite{CHJ09:Note-Behavior,SV09:Comments-Randomized} for a discussion of this issue.


\begin{center}
\begin{algorithm}[tb]
\caption{Block Kaczmarz Method with Uniform Random Control}
	\label{alg:randomized-block-kaczmarz}
\begin{center} \fbox{
\begin{minipage}{.95\textwidth} 
\vspace{4pt}
\alginout{\begin{itemize}
\item	Matrix $\mtx{A}$ with dimension $n \times d$
\item	Right-hand side $\vct{b}$ with dimension $n$
\item	Partition $T = \{\tau_1, \dots, \tau_m\}$ of the row indices $\{1, \dots, n\}$
\item	Initial iterate $\vct{x}_0$ with dimension $d$
\item	Convergence tolerance $\varepsilon > 0$
\end{itemize}}
%A matrix $\mtx{A} \in \Re^{n \times d}$, a vector $\vct{b} \in \Re^n$, a partition $T$ of the row indices $\{1, \dots, n\}$,
%an initial iterate $\vct{x}_0 \in \Re^d$, a convergence tolerance $\varepsilon$}
{An estimate $\hat{\vct{x}}$ for the solution to $\min_{\vct{x}} \enormsq{ \mtx{A}\vct{x} - \vct{b} }$
}
\vspace{8pt}\hrule\vspace{8pt}

\begin{algtab*}
$j \leftarrow 0$

\algrepeat %\\
	$j \leftarrow j + 1$ \\
	Choose a block $\tau$ uniformly at random from $T$ \\
	$\vct{x}_j \leftarrow \vct{x}_{j-1} + (\mtx{A}_\tau)^\pinv (\vct{b}_\tau - \mtx{A}_\tau\vct{x}_{j-1})$ \hfill \{ Solve least-squares problem \} \\ %\\
\alguntil{$\enormsq{ \mtx{A}\vct{x}_j - \vct{b} } \leq \eps^2$}
$\hat{\vct{x}} \leftarrow \vct{x}_j$
\end{algtab*}
%\vspace{-4pt}
\end{minipage}}
\end{center}
%\vspace{0in}
\end{algorithm}
\end{center}


 

%\notate{Check this paragraph + rewrite.}
%It is appropriate to choose $\eps^2 \approx \enormsq{\vct{e}}$, assuming that a bound for the size of the residual is available.
%According to standard results from perturbation theory~\cite[location]{GvL96:}, the convergence criterion ensures that
%\begin{equation} \label{eqn:perturbation}
%\enormsq{ \hat{\vct{x}} - \vct{x}_{\star} }
%	\leq \kappa^2(\mtx{A}) \, \eps^2.
%\end{equation}
%We postpone a more complete discussion about the implementation of Algorithm~\ref{alg:randomized-block-kaczmarz} to Section~\notate{xxxxx}.


%In this paper, we analyze a natural variant of the block Kaczmarz method that uses a randomized control mechanism.  We consider the case where the allowable blocks of indices are fixed in advance.  More precisely, the input to the algorithm includes a partition $T = \{ \tau_1, \dots, \tau_m \}$ of the row indices of $\mtx{A}$.  At each iteration, independently of all previous choices, we select a block $\tau$ of indices uniformly at random from the partition $T$.  These decisions lead to Algorithm~\ref{alg:randomized-block-kaczmarz}.



%Can apply Pinv using direct method (e.g., QR decomposition + triangular solve),
%but in many settings it is more efficient to use an iterative algorithm
%such as CGLS, LSQR, or even Chebyshev semi-iterative.  Efficiency of these
%schemes will depend on the properties of the blocks.
%
%randomized cyclic control
%
%Only want to check convergence criterion once per epoch


%We next introduce the block version of the Kaczmarz method.  For simplicty and to avoid degenerate cases, we will assume that the rows of the matrix $\mtx{A}$ all have unit $\ell_2$ norm.  We refer to such matrices as \textit{standardized}.  Like the traditional method, this version iteratively projects the current estimation onto the solution spaces.  However, rather than using the solution space of a \textit{single} equation, the block method projects onto the solution space of many equations simultaneously by selecting a block of rows rather than a single row.  The important question involves the means for selecting this block of rows.  To make sufficient progress in each step of the algorithm, we desire that the block of rows be \textit{well-conditioned}.  We define the condition number $\kappa(\mtx{A})$ of $\mtx{A}$ by
%
%\begin{equation}\label{eqn:kappa}
%\kappa(\mtx{A}) \defby \sup\left\{\frac{\|\mtx{A}\vct{x}\|_2}{\|\mtx{A}\vct{y}\|_2} : \|\vct{x}\|_2 =1, \|\vct{y}\|_2 = 1\right\} = \|\mtx{A}\|\cdot\|\mtx{A}^{-1}\|_2.
%\end{equation}
%Here and elsewhere, the symbol $\|\cdot\|$ indicates
%the spectral norm of a matrix.  Note that if the matrix does not have full column rank, the condition number is assigned the value $+\infty$.
%
%To guarantee well-conditioned blocks, we will utilize a theory on specially designed parititions $T$ of the rows of $\mtx{A}$, called \textit{pavings}.  Then we will select a block $\tau\in T$ at random in each iteration. 
%
%\begin{definition}\label{def:paving}
%A $(m, K)$ \textit{row paving} of the matrix $\mtx{A}$ is a partition $T = \{\tau_1, \tau_2, \ldots, \tau_m\}$ of the row indices such that
%$$
%\kappa(\mtx{A}^*_\tau) \leq K \quad\text{for each }\tau\in T.
%$$
%The numbers $m$ and $K$ are respectively called the \textit{size} and \textit{conditioning} of the paving.
%\end{definition}
%
%We will see in Section~\ref{sec:pavings} that we can easily construct such a paving.  In fact, we can do so by choosing blocks of $\mtx{A}$ \textit{at random}.  
%
%For a subset $\tau \subset \{1, 2, \ldots, n\}$, denote by $\mtx{A}_\tau$ the submatrix of $\mtx{A}$ whose rows are indexed by $\tau$.  Equipped with a row paving $T$ of $\mtx{A}$, we again begin with an arbitrary guess $x_0$ for the solution of the system.  Then for each iteration $j\geq 1$, select a block $\tau = \tau_j$ of rows at random from the paving $T$.  To obtain the next iterate, we project the current estimation onto the solution space of the equations listed in $\tau$:
%
%\begin{equation}\label{eqn:iterates}
%\vct{x}_j = \vct{x}_{j-1} + (\mtx{A_\tau})^\dagger(\vct{b}_\tau - \mtx{A}_\tau\vct{x}_{j-1}),
%\end{equation} 
%
%where $\mtx{A}_\tau^\dagger$ denotes the pseudoinverse of $\mtx{A}_\tau$.  When $\mtx{A}_\tau$ has full row rank, its pseudoinverse is determined by the formula $\mtx{A}_\tau^\dagger = \mtx{A}_\tau^*(\mtx{A}_\tau^* \mtx{A}_\tau)^{-1}$.  
%
%The method iterates until the residual error falls below a specified tolerance:
%$$
%\|\mtx{A}\vct{x}_j - \vct{b}\|_2 \leq \varepsilon.
%$$
%
%Then if the algorithm terminates at iteration $J$, standard perturbation theory for linear systems (see e.g.~\cite{HJ85:Matrix-Analysis}) allows us to conclude that the error satisfies
%$$
%\|\vct{x}_J - \vct{x}_\star\|_2 \leq \varepsilon\kappa(\mtx{A}).
%$$
%
%Since the paving will guarantee that each submatrix $\mtx{A}_\tau$ is well conditioned, we propose using the conjugate gradient algorithm to apply the pseudoinverse in~\eqref{eqn:iterates}.  Our proposed block Kaczmarz algorithm is then described by the following.
%
%\begin{algorithm}[thb]
%\caption{Block Kaczmarz}
%	\label{alg:randomized-block-kaczmarz}
%\centering \fbox{
%\begin{minipage}{.99\textwidth} 
%\vspace{4pt}
%\alginout{Matrix $\mtx{A}$, vector $\vct{b}$, paving $T$, random control map $r : \mathbb{N} \rightarrow T$, tolerance $\varepsilon$}
%{An estimation $\vct{x}_j$ of the unique solution $\vct{x}_\star$ to $\mtx{A}\vct{x} = \vct{b}$
%}
%\vspace{8pt}\hrule\vspace{8pt}
%
%\begin{algtab*}
%Set $\vct{x}_0$.
%	 	\hfill \{ Trivial initial approximation \} \\		 
%$j \leftarrow 0$ \\
%
%\algrepeat
%	$j \leftarrow j + 1$ \\
%	$\tau \leftarrow \tau_{r(j)}$
%		\hfill \{ Randomly select a block $\tau\in T$ \} \\
%	$\vct{x}_j \leftarrow \vct{x}_{j-1} + (\mtx{A_\tau})^\dagger(\vct{b}_\tau - \mtx{A}_\tau\vct{x}_{j-1})$
%		\hfill \{ Project using conjugate gradient\} \\	
%\alguntil{$\|\mtx{A}\vct{x}_j - \vct{b}\|_2 \leq \varepsilon$} \\
%	
%%\alguntil{$\enorm{\vct{v}} \leq \eps$}
%% \{ Halting criterion \} \\
%	%end {\bf end} {\it until}
%
%%$\vct{a} \leftarrow \vct{a}^{k}$
%%	\hfill \{ Return final approximation \}
%\end{algtab*}
%\vspace{-14pt}
%\end{minipage}}
%%\vspace{-0.5in}
%\end{algorithm}
%
%\vfill


\subsection{Desiderata for the Partition} \label{sec:desiderata}

%As you may imagine,
The implementation and behavior of the block Kaczmarz method
depend heavily on the properties of the submatrices $\mtx{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~\ref{sec:intro-main}.

%\notate{Simplify...}
The most expensive (arithmetic) step in Algorithm~\ref{alg:randomized-block-kaczmarz}
occurs when we apply the pseudoinverse $\mtx{A}_{\tau}^\pinv$ to a vector.
We can perform this calculation efficiently \emph{provided that}
each submatrix $\mtx{A}_{\tau}$ has well-conditioned rows.
%In the extreme case where $\mtx{A}_{\tau}$ has orthonormal rows,
%the pseudoinverse equals the conjugate transpose: $\mtx{A}^\pinv = \mtx{A}^*$.
%More generally,
Indeed, in this case, we can invoke an iterative least-squares solver~\cite{Bjo96:Numerical-Methods},
such as {\tt CGLS},
to apply the pseudoinverse $\mtx{A}_\tau^\pinv$
approximately using a small number of matrix--vector
multiplies with $\mtx{A}_{\tau}$ and $\mtx{A}_{\tau}^*$.
In particular, we never need to form the pseudoinverse.

%It emerges that block Kaczmarz can offer significant
%advantages over simple Kaczmarz when each submatrix admits
%a fast multiply.


%  In particular,
%note that these matrix--vector multiplies can be ruthlessly efficient
%when each submatrix $\mtx{A}_{\tau}$ is sparse or structured (e.g., Toeplitz).


%  There
%are several situations where this computation can be performed efficiently:
%
%\begin{description} \setlength{\itemsep}{8pt}
%\item	[Well-Conditioned Blocks]
%When the rows of each submatrix $\mtx{A}_{\tau}$ form a well conditioned set,
%we can apply the pseudoinverse quickly using an iterative algorithm, such as
%CGLS, LSQR, or the Chebyshev semi-iterative method~\cite{Bjo96:Least-Squares}.
%In particular, when the condition number of $\mtx{A}_{\tau}$ is constant, we need
%$\bigO(\log(1/\eps))$ matrix--vector multiplies with $\mtx{A}_{\tau}$ to reduce
%the initial error by a factor of $\eps$.
%
%\item	[Orthogonal Blocks]
%When each submatrix $\mtx{A}_\tau$ has orthonormal rows, the pseudoinverse
%$\mtx{A}^\pinv = \mtx{A}^*$.  The cost is a single matrix--vector multiply.
%
%\item	[Structured Blocks]
%When each submatrix $\mtx{A}_\tau$ has algebraic structure (e.g., Toeplitz
%or partial Fourier transform), we can perform matrix--vector multiplies
%efficiently.
%\end{description}

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

%The behavior of the block Kaczmarz scheme,
%Algorithm~\ref{alg:randomized-block-kaczmarz},
%depends heavily on the metric properties of the blocks
%of equations that the algorithm can enforce.
%The following important definition summarizes
%the information that we will need.

%\begin{definition}[Row Paving]
%%Let $\mtx{A}$ be a matrix.
%An $(m, \alpha, \beta)$ \emph{row paving} of a matrix $\mtx{A}$ is a partition $T = \{ \tau_1, \dots, \tau_m \}$ of the row indices that verifies
%$$
%\alpha \leq \sigma_{\min}^2( \mtx{A}_\tau ) \leq
%%\quad\text{and}\quad
%\sigma_{\max}^2( \mtx{A}_{\tau} ) \leq \beta
%\quad\text{for each $\tau \in T$.}
%$$
%The number $m$ of blocks is called the \emph{size} of the paving.  The numbers $\alpha$ and $\beta$ are called \emph{lower} and \emph{upper bounds} on the paving.  The ratio $\beta / \alpha$ gives a uniform bound on the squared condition number of %of the rows of
%$\mtx{A}_{\tau}$ for each $\tau$.
%\end{definition}


%\notate{Still need to emphasize that we're thinking about fat blocks!}

\begin{definition}[Row Paving]
%Let $\mtx{A}$ be a matrix.
An $(m, \alpha, \beta)$ \emph{row paving} of a matrix $\mtx{A}$ is a partition $T = \{ \tau_1, \dots, \tau_m \}$ of the row indices that verifies
$$
\alpha \leq \lambda_{\min}( \mtx{A}_\tau \mtx{A}_{\tau}^* )
\quad\text{and}\quad
\lambda_{\max}( \mtx{A}_{\tau} \mtx{A}_{\tau}^* ) \leq \beta
\quad\text{for each $\tau \in T$.}
$$
The number $m$ of blocks is called the \emph{size} of the paving.  The numbers $\alpha$ and $\beta$ are called \emph{lower} and \emph{upper paving bounds}.  The ratio $\beta / \alpha$ gives a uniform bound on the squared condition number $\kappa^2(\mtx{A}_{\tau})$ for each $\tau$.  Note that $\alpha = 0$ unless each submatrix $\mtx{A}_{\tau}$ is fat.
\end{definition}

Every partition $T$ of the rows of a matrix $\mtx{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~\ref{sec:intro-paving} and~\ref{sec:pavings}, 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!

%The paving parameters have a few simple properties that are worth mentioning.
%\begin{itemize}
%\item	The size $m$ of the paving ranges between one and the number of rows of $\mtx{A}$.
%
%\item	The paving bounds satisfy $0 \leq \alpha \leq \beta \leq \norm{\mtx{A}}^2$.
%
%\item	Each block of the paving has full row rank if and only if $\alpha > 0$.
%
%\item	The squared condition number $\kappa^2$ provides a uniform bound for the conditioning
%of the matrices $\mtx{A}_\tau$ for $\tau \in T$.
%\end{itemize}


\subsection{Convergence of Randomized Block Kaczmarz} \label{sec:intro-main}

The main result of this paper provides information about the convergence properties of the randomized block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}, in terms of the parameters of the row paving $T$.

\begin{theorem}[Convergence] \label{thm:convergence}
Suppose $\mtx{A}$ is a matrix with full column rank that admits an $(m, \alpha, \beta)$ row paving $T$.  Consider the least-squares problem
$$
\mathrm{minimize} \quad \enormsq{ \mtx{A} \vct{x} - \vct{b} }.
$$
Let $\vct{x}_{\star}$ be the unique minimizer, and define the residual $\vct{e} := \mtx{A} \vct{x}_{\star} - \vct{b}$.
For any initial estimate $\vct{x}_0$, the randomized block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}, produces a sequence $\{\vct{x}_j : j \geq 0 \}$ of iterates that satisfies
\begin{equation} \label{eqn:block-rate}
\Expect \enormsq{ \vct{x}_j - \vct{x}_{\star} }
	\ \leq \ \left[ 1 - \frac{\sigma_{\min}^2(\mtx{A})}{\beta m} \right]^{j}
	\enormsq{ \vct{x}_0 - \vct{x}_{\star} }
	\ + \ \frac{\beta}{\alpha} \cdot \frac{\enormsq{\vct{e}}}{ \sigma_{\min}^2(\mtx{A}) }.
\end{equation}
%The symbol $\sigma_{\min}^2(\mtx{A})$ denotes the smallest positive singular value of $\mtx{A}$.
\end{theorem}

\noindent
Turn to Section~\ref{sec:results} for the proof of Theorem~\ref{thm:convergence}.

The expression~\eqref{eqn:block-rate} 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 \emph{convergence horizon}, is comparable with the second term on the right-hand side of~\eqref{eqn:block-rate}.  The bracket controls the convergence rate.  The minimum singular value of $\mtx{A}$ affects both the rate of convergence and the convergence horizon.  In each case, we prefer $\sigma_{\min}(\mtx{A})$ to be as large as possible.

The properties of the row paving play an interesting role in Theorem~\ref{thm:convergence}.  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., $\vct{e}$ is nonzero).
Nevertheless, as Section~\ref{sec:desiderata} suggests, we usually want the paving to be well conditioned to ensure that we can apply the block update rule~\eqref{eqn:block-kacz} efficiently.

%Referring back to the error bound~\eqref{eqn:Needell-rate}, we see that the simple Kaczmarz algorithm and the block Kaczmarz algorithm have similar qualitative behavior.


%\notate{Interpretation:}
%Similar to the simple
%randomized Kaczmarz method from Section~\ref{sec:simple-kacz}.
%Linear convergence in expectation until we reach a ball about the true solution.
%Radius of the ball is determined by the second term.

%\notate{Role of sigma min:}
%Controls both rate of convergence and the ultimate proximity to the solution.
%In both cases, small sigma min = bad news.  Roughly parallel with simple case.
%
%\notate{Paving parameters:}
%Curiously, rate of convergence is determined only by the upper bound and
%the number of blocks.  If we can control these, we get fast convergence.
%The radius of the error horizon, however, depends on the ratio
%$\kappa^2 := \beta/\alpha$.  For inconsistent problems, need control here.
%
%And don't forget that small $\kappa^2$ is important for implementation.
%We assume ls problems are solved exactly...
%
%To make a full comparison with earlier results, need to get hold of
%paving parameters...  In some examples, there's a natural paving
%where parameters are easily determined.  But in other situations,
%we may need to produce the paving ourselves...

\subsection{Simple Kaczmarz versus Block Kaczmarz}

First, notice that Theorem~\ref{thm:convergence} improves on the earlier result~\eqref{eqn:Needell-rate} 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 $\mtx{A}$ is standardized, the paving constants satisfy $\alpha = \beta = 1$, and we reach the error bound
$$
\Expect \enormsq{ \vct{x}_j - \vct{x}_{\star} }
	\ \leq \ \left[ 1 - \frac{\sigma_{\min}^2(\mtx{A})}{n} \right]^{j}
	\enormsq{ \vct{x}_0 - \vct{x}_{\star} }
	\ + \ \frac{\enormsq{\vct{e}}}{ \sigma_{\min}^2(\mtx{A}) }.
$$
The convergence horizon $\normsq{\vct{e}} \leq n \infnormsq{\vct{e}}$, so the displayed bound beats~\eqref{eqn:Needell-rate} when $\mtx{A}$ is standardized.

More generally,
suppose $\mtx{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~\ref{sec:simple-kacz}) to the block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}.  Both methods satisfy an error bound of the form
$$
\Expect \enormsq{\vct{x}_j - \vct{x}_{\star}}
	\ \leq \ \cnst{e}^{-j \rho} \cdot \enormsq{\vct{x}_0 - \vct{x}_{\star}} \ + \ h
$$
where the convergence rate $\rho$ and the convergence horizon $h$ depend on the choice of algorithm.

%cf.~\eqref{eqn:Needell-rate} and~\eqref{eqn:block-rate}.

First, we compare the convergence rates of the two methods.  The bounds~\eqref{eqn:Needell-rate}
and~\eqref{eqn:block-rate} imply
\begin{equation} \label{eqn:compare-rates}
\rho_{\rm simp} \geq \frac{\sigma_{\min}^2(\mtx{A})}{n}
\quad\text{and}\quad
\rho_{\rm block} \geq \frac{\sigma_{\min}^2(\mtx{A})}{\beta m}
\end{equation}
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.

%Comparing the two rates, we see that the number $J$ of iterations required
%to reduce the initial error by a factor of $\eps$ is approximately
%$$
%J_{\rm block} \approx \frac{\beta m}{\sigma_{\min}^2(\mtx{A})} \log(1/\eps)
%\quad\text{and}\quad
%J_{\rm kacz} \approx \frac{n}{\sigma_{\min}^2(\mtx{A})} \log(1/\eps),
%$$
%provided that $\eps$ is not too small.

Next, we examine the convergence horizons.  The bounds~\eqref{eqn:Needell-rate}
and~\eqref{eqn:block-rate} yield
$$
h_{\rm simp} = \frac{n \infnormsq{\vct{e}}}{\sigma_{\min}^2(\mtx{A})}
\quad\text{and}\quad
h_{\rm block} = \frac{\beta}{\alpha} \cdot \frac{\enormsq{\vct{e}}}{\sigma_{\min}^2(\mtx{A})}.
$$
Their ratio satisfies
$$
\frac{h_{\rm block}}{h_{\rm simp}}
	= \frac{\beta}{\alpha} \cdot \frac{\enormsq{\vct{e}}}{n \infnormsq{\vct{e}}}
	\leq \frac{\beta}{\alpha}.
$$
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 $\mtx{A}$, or a lack thereof.

%For further discussion of the simple randomized Kaczmarz algorithm, see the papers~\cite{SV09:Randomized-Kaczmarz,N10:rknoisy,LL10:Randomized-Methods}.


\subsubsection{Example 1: Fast Updates} \label{sec:fast-update}

First, we consider the case where the computational cost of the block update
rule~\eqref{eqn:block-kacz} is roughly comparable with the cost of the
simple update rule~\eqref{eqn:simple-kacz}.  This situation can occur when
\begin{itemize}
\item	Each submatrix $\mtx{A}_{\tau}$ admits a fast multiply; and

\item	Each submatrix $\mtx{A}_{\tau}$ is well conditioned.
\end{itemize}
See Section~\ref{sec:experiments-circ} 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~\eqref{eqn:compare-rates} 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.


\subsubsection{Example 2: Unstructured Submatrices} \label{sec:unstruct}

On the other hand, suppose that each submatrix $\mtx{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~\eqref{eqn:compare-rates}
is unfair to the simple method.

In this case, it is more appropriate to examine the convergence rate
per \emph{epoch}, the minimum number of iterations it takes the algorithm
to touch each row of $\mtx{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:
$$
m \cdot \rho_{\rm block} \geq \frac{\sigma^2_{\min}(\mtx{A})}{\beta}
\quad\text{and}\quad
n \cdot \rho_{\rm simp} \geq \sigma^2_{\min}(\mtx{A}).
$$
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 ({\tt BLASx}).
We discuss these points in Section~\ref{sec:why}, and we provide 
some numerical support in Section~\ref{sec:experiments-dense}.

%Furthermore, because of data transfer costs,
%the block method is more computationally efficient in many computer
%architectures.  Indeed, the block Kaczmarz algorithm transfers a
%large submatrix into working memory and operates with it for some time,
%whereas the simple Kaczmarz method repeatedly transfers new equations
%into working memory.  See Section~\ref{sec:xxx} for a numerical example
%that provides some evidence for this claim.


%\subsubsection{Comparison with Conjugate Gradient}
%
%\notate{Cut this?  If not, revise.}
%
%One epoch of block Kaczmarz has a comparable cost with a single iteration
%of a conjugate gradient algorithm.  For consistent systems, there is
%a CG method that minimizes the error in the
%unknown~\cite[Sec.~7.4.1]{Bjo96:Numerical-Methods}.
%It achieves the convergence rate
%$$
%\normsq{\vct{x}_j - \vct{x}_{\star}}
%	\leq 4 \left[ \frac{\kappa(\mtx{A}) - 1}{\kappa(\mtx{A}) + 1} \right]^{2j}
%	\normsq{\vct{x}_0 - \vct{x}_{\star}}
%$$
%Therefore, the convergence rate for this flavor of conjugate gradient satisfies
%$$
%\rho_{\rm cg} \geq \frac{4}{\kappa(\mtx{A}) + 1}.
%$$
%Remark: This is different from CGLS, which is what people actually use.
%But the comparison is clearer for this example.  Result suggests that
%block Kaczmarz can be more efficient in the case where $\sigma_{\min}(\mtx{A})$
%is not so small, while $\sigma_{\max}(\mtx{A})$ is pretty big.  We can test this.



\subsection{Existence of Good Pavings} \label{sec:intro-paving}

So far, we have assumed that the matrix $\mtx{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~\cite{Pop99:Block-Projections-Algorithms}
has developed an approach for producing a paving of a sparse matrix;
see also~\cite{Pop01:Fast-Kaczmarz-Kovarik,Pop04:Kaczmarz-Kovarik-Algorithm}.
But we can travel much farther down this road.
It is an astonishing fact that \emph{every}
standardized matrix admits a good row paving.

\begin{proposition}[Existence of Good Row Pavings] \label{prop:intro-paving}
Fix a number $\delta \in (0, 1)$.  Let $\mtx{A}$ be a standardized matrix with $n$ rows.
Then $\mtx{A}$ admits a row paving whose parameters satisfy
$$
m \leq \cnst{C_{pave}} \cdot \delta^{-2} \normsq{\mtx{A}} \log(1+n)
\quad\text{and}\quad
1 - \delta \leq \alpha \leq \beta \leq 1 + \delta.
$$
The number $\cnst{C_{pave}}$ is a positive, universal\,%
\footnote{A \emph{universal} constant has no dependence on any parameter.}
constant.
\end{proposition}

%Proposition~\ref{prop:intro-paving} is a folkloric result.

Proposition~\ref{prop:intro-paving} follows most directly
from the recent results~\cite[Cor.~1.5]{Ver06:Random-Sets}
and~\cite[Thm.~1.2]{Tro09:Column-Subset},
whose provenance can be traced to the celebrated
papers~\cite{BT87:Invertibility-Large,BT91:Problem-Kadison-Singer}.
%Section~\ref{sec:pavings} contains more information.
Although Proposition~\ref{prop:intro-paving} 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 \emph{at random}.  See Section~\ref{sec:pavings} for more results and background on paving.

%We have been unable to locate the first analysis that
%leads to Proposition~\ref{prop:intro-paving}, although
%the result is anchored firmly in the literature.
%For example, see~\cite[Thm.~1.2]{Tro09:Column-Subset},
%which mistakenly attributes Bourgain and
%Tzafriri~\cite[Cor.~1.2]{BT91:Problem-Kadison-Singer}.  See
%Section~\ref{sec:pavings} for further discussion.

%%%
%%% The basic problem is that the analysis of BT91 leads to
%%% a dependence on ||A||^4.  Combining restricted invertibility (BT87)
%%% and Kashin--Tzafriri (1994) gives a similar result,
%%% but there's not accurate control on the upper/lower paving bound. 
%%% Vershynin's analysis (Subsequences of frames 1999;
%%% John's decompositions, 2001) doesn't yield adequate
%%%% conditioning bounds either.


% but this claim follows
%as an easy consequence of the restricted invertibility
%principle of Bourgain--Tzafriri~\cite[Thm.~1.2]{BT87:Invertibility-Large}
%and a complementary result due to Kashin--Tzafriri~\cite{KT94:Some-Remarks}.
%Vershynin's work~\cite{Ver99:Subsequences-Frames,Ver01:Johns-Decompositions}
%contains some related analysis.  See Section~\ref{sec:pavings} for more
%information.

%See Section~\ref{sec:pavings} for further discussion.

%It is an unexpected surprise that every matrix admits a good paving.
%This is the content of a major result due to Bourgain and Tzafriri...
%\notate{Wrong attribution???}
%
%Additional references: Careful here... Not clear that BT get the right dependence on the norm... My paper makes this clear?  Roman 2001, 2006...  Maybe this belongs later...

%For the moment, we restrict ourselves to some short remarks about the incredible content of Proposition~\ref{prop:intro-paving}.

\subsubsection{Understanding the Paving Theorem}

Before we continue, let us take a moment to explain some of the key aspects of Proposition~\ref{prop:intro-paving}.
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~\ref{sec:desiderata}.  We remark that the dependence of the paving size $m$ on the parameter $\delta$ is optimal%
\footnote{To verify this point, consider a large matrix $\mtx{A}$ whose entries have equal magnitude and independent random signs.  Use the Bai--Yin Law~\cite[Thm.~2]{BY93:Limit-Smallest} to estimate singular values and norms.}
as $\delta \to 0$.

%Observe that Proposition~\ref{prop:intro-paving} describes a row paving whose squared conditioning $\beta/\alpha \leq 3$, which is essentially optimal.  This property licenses the use of an iterative algorithm to solve least-squares problems involving the submatrices induced by the paving, as described in Section~\ref{sec:desiderata}.  The numbers 0.5 and 1.5 have no particular significance---we could make them arbitrarily close to one---but the value of $\cnst{C}_{pave}$ does depend on these specific choices.

Next, we develop some intuition about the role of the spectral norm in Proposition~\ref{prop:intro-paving}.  Suppose that $\mtx{A}$ is a standardized matrix with $n$ rows.  If the lower paving bound $\alpha > 0$, then a row paving $T$ of $\mtx{A}$ must contain at least $n / \rank(\mtx{A})$ blocks.  Otherwise, $\mtx{A}_{\tau}$ is rank deficient for some $\tau \in T$ simply because this submatrix has more than $\rank(\mtx{A})$ rows.  Now, using the fact $\fnormsq{\mtx{A}} \leq \rank(\mtx{A}) \normsq{\mtx{A}}$, we easily verify that 
$$
\normsq{\mtx{A}} \geq \frac{n}{\rank(\mtx{A})}.
$$
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~\ref{prop:intro-paving} 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~\cite{And79:Extensions-Restrictions,BHKW88:Matrix-Norm,BT91:Problem-Kadison-Singer} that an affirmative answer would imply the long-standing conjecture of Kadison and Singer~\cite{KS59:Extensions-Pure}.

%Proposition~\ref{prop:intro-paving} guarantees that the size of the paving lies within a logarithmic factor of its optimal value.  To see this,
%let us consider the two extreme examples of standardized matrices:
%\begin{itemize} \setlength{\itemsep}{4pt}
%\item	When the rows of $\mtx{A}$ are orthonormal, the smallest paving contains one block.  Since $\normsq{\mtx{A}} = 1$, Proposition~\ref{prop:intro-paving} asserts a row paving with at most $\cnst{C_{pave}} \log(n)$ blocks.
%
%\item	When the $n$ rows of $\mtx{A}$ are identical, the smallest paving requires $n$ blocks.  Since $\normsq{\mtx{A}} = n$, Proposition~\ref{prop:intro-paving} is a trivial result: it offers us a row paving with at most $\cnst{C_{pave}} n \log(n)$ blocks.
%\end{itemize}
%\notate{say a bit more...}

%It is known that Proposition~\ref{prop:intro-paving} is optimal up to the logarithmic factor in the size $m$ of the paving~\cite[Thm.~3.4]{BHKW88:Matrix-Norm}.  Based on the examples above, one may wonder whether it is possible to remove the logarithm.  This question remains open, but it is known~\cite{And79:Extensions-Restrictions} that an affirmative answer would imply the conjecture of Kadison and Singer~\cite{KS59:Extensions-Pure}.



\subsection{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.

\begin{corollary}[Block Kaczmarz with a Good Row Paving] \label{cor:good-intentions}
Suppose that $\mtx{A}$ is a standardized matrix with full column rank.  Let $T$ be a good row paving of $\mtx{A}$, as guaranteed by Proposition~\ref{prop:intro-paving} with $\delta = 1/2$.  Under the notation of Theorem~\ref{thm:convergence}, the block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz} admits the convergence estimate
$$
\Expect \normsq{ \vct{x}_j - \vct{x}_{\star} }
	\leq \left[ 1 - \frac{1}{6 \cnst{C_{pave}} \kappa^2(\mtx{A}) \log(1+n)} \right]^j
	\enormsq{ \vct{x}_0 - \vct{x}_\star }
	\ + \
	\frac{3\enormsq{\vct{e}}}{\sigma_{\min}^2(\mtx{A})}
$$
where $\cnst{C_{pave}}$ is the constant from Proposition~\ref{prop:intro-paving}.
\end{corollary}

To summarize, Proposition~\ref{prop:intro-paving} yields a small paving of the matrix $\mtx{A}$ with exceptional conditioning.  With this choice of paving, we can perform the block Kaczmarz update~\eqref{eqn:block-kacz} quickly using an iterative least-squares algorithm.  (Indeed, to apply $\mtx{A}_\tau^\pinv$ with a fixed level of precision, it suffices to perform a constant number of matrix--vector multiplies with $\mtx{A}_{\tau}$ and $\mtx{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 $\mtx{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.


 % in the least-squares problem~\eqref{eqn:basic-ls}.
% In view of~\eqref{eqn:conv-tol}, we can set the convergence tolerance $\eps = 2 \enorm{\vct{e}}$.


\subsection{Organization}

The rest of the paper has the following structure.  Section~\ref{sec:results} contains a proof of the main result, Theorem~\ref{thm:convergence}.  In Section~\ref{sec:pavings}, we give an overview of the literature on pavings.  Section~\ref{sec:implementation} discusses numerical aspects of the block Kaczmarz method.  We discuss related work on Kaczmarz methods and future directions in Section~\ref{sec:future}.  Finally, Appendix~\ref{app:fit} offers a proof of a supplemental result.

%The rest of the paper has the following structure.  In Section~\ref{sec:pavings} we state results about random paving constructions which we will utilize in the block Kaczmarz method.  Section~\ref{sec:results} states the main results and their proofs.  We prove Theorem~\ref{thm:big} and a result about systems perturbed by noise.  We discuss these results and comparison with other methods in Section~\ref{sec:discuss}.  Finally, we provide numerical evidence for the convergence of the block Kaczmarz method in Section~\ref{sec:numerics}.  The appendix provides the proofs of the paving theorems presented in Section~\ref{sec:pavings}.


\section{Analysis of the Randomized Block Kaczmarz Algorithm}
\label{sec:results} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This section contains the proof of the main result, Theorem~\ref{thm:convergence},
on the convergence of Algorithm~\ref{alg:randomized-block-kaczmarz}.  We commence
with two simple lemmas.  The first step provides a deterministic bound on how
much one iteration of the algorithm reduces the error.

%\subsection{Convergence Analysis}

\begin{lemma} \label{lem:err-bd}
Instate the hypotheses and notation of Theorem~\ref{thm:convergence}.
Then the error at iteration $j$ satisfies the deterministic bound
$$
\enormsq{ \vct{x}_{j} - \vct{x}_{\star} }
	\leq \enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau})(\vct{x}_{j-1} - \vct{x}_{\star}) }
	\ + \
	\frac{1}{\alpha} \enormsq{ \vct{e}_{\tau} },
$$
where $\tau = \tau(j)$ is the block selected at iteration $j$.
\end{lemma}

\begin{proof}
%Abbreviate $\tau = \tau(j+1)$.
According to the update rule~\eqref{eqn:block-kacz}, block Kaczmarz computes
$$
\vct{x}_{j}
	= \vct{x}_{j-1} + \mtx{A}_{\tau}^\pinv (\vct{b}_{\tau} - \mtx{A}_{\tau} \vct{x}_{j-1})
	= \vct{x}_{j-1} + \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} (\vct{x}_\star - \vct{x}_{j-1})
	- \mtx{A}_{\tau}^\pinv \vct{e}_{\tau},
$$
where we have introduced the decomposition $\vct{b} = \mtx{A} \vct{x}_{\star} - \vct{e}$, restricted to the coordinates listed in $\tau$.
Subtract $\vct{x}_{\star}$ from both sides to obtain
$$
\vct{x}_{j} - \vct{x}_{\star}
	= (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j-1} - \vct{x}_\star)
	- \mtx{A}_{\tau}^\pinv \vct{e}_{\tau}.
$$
The range of $\mtx{A}_{\tau}^\pinv$ and the range of $\Id - \mtx{A}_\tau^\pinv \mtx{A}_\tau$ are orthogonal, so we may invoke the Pythagorean Theorem to reach
$$
\enormsq{ \vct{x}_{j} - \vct{x}_{\star} }
	= \enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j-1} - \vct{x}_\star) }
	+ \enormsq{ \mtx{A}_{\tau}^\pinv \vct{e}_{\tau} }.
$$
The second term on the right-hand side satisfies
$$
\enormsq{ \mtx{A}_{\tau}^\pinv \vct{e}_{\tau} }
	\leq \sigma_{\max}^2(\mtx{A}_{\tau}^\pinv) \enormsq{\vct{e}_{\tau}}
	\leq \frac{1}{\sigma_{\min}^2(\mtx{A}_{\tau})} \enormsq{\vct{e}_{\tau}}
	\leq \frac{1}{\alpha} \enormsq{\vct{e}_{\tau}},
$$
where $\alpha$ is the lower bound on the row paving $T$.
Combine the last two displays to wrap up.
\end{proof}

%\begin{lemma}
%Instate the hypotheses and notation of Theorem~\ref{thm:convergence}.
%Then the error at iteration $j+1$ satisfies the bound
%$$
%\Expect_{\tau(j+1)} \enormsq{ \vct{x}_{j+1} - \vct{x}_{\star} }
%	\leq \left[ 1 - \frac{\sigma_{\min}^2(\mtx{A})}{\beta m} \right]
%	\enormsq{ \vct{x}_j - \vct{x}_{\star} }
%	\ + \
%	\frac{1}{\alpha} \enormsq{ \vct{e}_{\star} }
%$$
%where the expectation is computed with respect the random variable $\tau = \tau(j+1)$.
%\end{lemma}
%
%\begin{proof}
%Abbreviate $\tau = \tau(j+1)$.
%According to the rule~\eqref{eqn:block-kacz}, block Kaczmarz makes the update
%$$
%\vct{x}_{j+1}
%	= \vct{x}_j + \mtx{A}_{\tau}^\pinv (\vct{b}_{\tau} - \mtx{A}_{\tau} \vct{x}_{j})
%	= \vct{x}_j + \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} (\vct{x}_\star - \vct{x}_{j})
%	+ \mtx{A}_{\tau}^\pinv \vct{e}_{\star},
%$$
%where we have used the decomposition $\vct{b} = \mtx{A} \vct{x}_{\star} - \vct{e}_{\star}$.
%Subtract $\vct{x}_{\star}$ from both sides to obtain
%$$
%\vct{x}_{j+1} - \vct{x}_{\star}
%	= (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j} - \vct{x}_\star)
%	+ \mtx{A}_{\tau}^\pinv \vct{e}_{\star}.
%$$
%Compute the squared $\ell_2$ norm using the Pythagorean theorem:
%$$
%\enormsq{ \vct{x}_{j+1} - \vct{x}_{\star} }
%	= \enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j} - \vct{x}_\star) }
%	+ \enormsq{ \mtx{A}_{\tau}^\pinv \vct{e}_{\star} }.
%$$
%The second term on the right-hand side satisfies
%$$
%\enormsq{ \mtx{A}_{\tau}^\pinv \vct{e}_{\star} }
%	\leq \sigma_{\max}^2(\mtx{A}_{\tau}^\pinv) \enormsq{\vct{e}_{\star}}
%	\leq \frac{1}{\alpha} \enormsq{\vct{e}_{\star}},
%$$
%where $\alpha$ is the lower bound on the row paving $T$.
%Combine the last two displays and average over the choice of the block $\tau$
%to reach
%\begin{align*}
%\Expect_{\tau} \enormsq{ \vct{x}_{j+1} - \vct{x}_{\star} }
%	&= \Expect_{\tau}
%	\enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j} - \vct{x}_\star) }
%	+ \frac{1}{\alpha} \enormsq{\vct{e}_{\star}} \\
%	&\leq \left[ 1 - \frac{\sigma^2_{\min}(\mtx{A})}{\beta m} \right]
%	\enormsq{ \vct{x}_j - \vct{x}_{\star} }
%	+ \frac{1}{\alpha} \enormsq{\vct{e}_{\star}}.
%\end{align*}
%The inequality follows from Lemma~\ref{lem:project-rdm-block}.
%\end{proof}

The second lemma gives us a means to average the two quantities appearing
in Lemma~\ref{lem:err-bd} over a random choice of the block $\tau = \tau(j)$.

\begin{lemma} \label{lem:project-rdm-block}
Instate the hypotheses and notation of Theorem~\ref{thm:convergence}.
Suppose that $\tau$ is chosen uniformly at random from the row paving $T$.
For fixed vectors $\vct{u}$ and $\vct{v}$, %that are statistically independent from $\tau$,
it holds that
%%\begin{align*}
$$
\Expect \enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau}) \vct{u} }
	\leq \left[ 1 - \frac{\sigma^2_{\min}(\mtx{A})}{\beta m} \right] \enormsq{\vct{u}}
\quad\text{and}\quad
\Expect \enormsq{ \vct{v}_{\tau} } = \frac{1}{m} \enormsq{\vct{v}}.
$$
\end{lemma}

\begin{proof}
The second identity emerges from a very short calculation:
$$
\Expect \enormsq{\vct{v}_\tau}
	= \frac{1}{m} \sum\nolimits_{\omega \in T} \enormsq{\vct{v}_{\omega}}
	= \frac{1}{m} \enormsq{\vct{v}},
$$
which depends on the fact that the blocks of $T$ partition the components of $\vct{v}$.

Since $\mtx{A}_{\tau}^\pinv \mtx{A}_{\tau}$ is an orthogonal projector, we may apply the Pythagorean Theorem to obtain the relation
$$
\Expect \enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau}) \vct{u} }
	= \enormsq{\vct{u}} - \Expect \enormsq{ \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} \vct{u} }.
$$
We control the remaining expectation as follows.
\begin{align*}
\Expect \enormsq{ \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} \vct{u} }
	&\geq \Expect \left[ \sigma_{\min}^2(\mtx{A}_{\tau}^\pinv) \enormsq{ \mtx{A}_{\tau} \vct{u} } \right]
	\geq \frac{1}{\beta} \cdot \Expect \enormsq{ \mtx{A}_{\tau} \vct{u} } \\
	&= \frac{1}{\beta m} \sum_{\omega \in T} \enormsq{ \mtx{A}_{\omega} \vct{u} }
	= \frac{1}{\beta m} \enormsq{ \mtx{A} \vct{u} }
	\geq \frac{\sigma^2_{\min}(\mtx{A})}{\beta m} \enormsq{\vct{u}}.
\end{align*}
The second inequality depends on the bound $\sigma_{\min}^2(\mtx{A}_{\tau}^\pinv) = \sigma_{\max}^{-2}(\mtx{A}_{\tau}) \geq \beta^{-1}$.  The fourth relation holds because the blocks in a paving partition the row indices of $\mtx{A}$.  To complete the proof, we simply combine the last two displays.
\end{proof}

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

\begin{proof}[Proof of Theorem~\ref{thm:convergence}]
First, we bound the expected error at iteration $j$ in terms of the error at iteration $j - 1$.
Average the bound from Lemma~\ref{lem:err-bd} over the randomness in $\tau = \tau(j)$
to reach
%we combine Lemmas~\ref{lem:err-bd} and~\ref{lem:project-rdm-block}.  Averaging
%over the randomness in $\tau = \tau(j)$, we reach
\begin{align*}
\Expect_{\tau} \enormsq{ \vct{x}_{j} - \vct{x}_{\star} }
	&\leq \Expect_{\tau}
	\enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} )(\vct{x}_{j-1} - \vct{x}_\star) }
	+ \frac{1}{\alpha} \cdot \Expect_{\tau} \enormsq{\vct{e}_{\tau}} \\
	&\leq \left[ 1 - \frac{\sigma^2_{\min}(\mtx{A})}{\beta m} \right]
	\enormsq{ \vct{x}_{j-1} - \vct{x}_{\star} }
	+ \frac{1}{\alpha m} \enormsq{\vct{e}}.
\end{align*}
The second inequality follows from Lemma~\ref{lem:project-rdm-block},
with $\vct{u} = \vct{x}_{j-1} - \vct{x}_{\star}$ and $\vct{v} = \vct{e}$.

%To control the second term on the right-hand side, we have used the
%observation $\Expect_{\tau} \enormsq{\vct{e}_{\tau}}
%= m^{-1} \sum_{\rho \in T} \enormsq{\vct{e}_{\rho}}
%= m^{-1} \enormsq{\vct{e}}$, which follows from the fact that
%$T$ is a partition.

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(\mtx{A})/(\beta m)$,
we obtain the estimate
\begin{align*}
\Expect \enormsq{ \vct{x}_{j} - \vct{x}_{\star} }
	&= \Expect_{\tau(1)} \Expect_{\tau(2)} \cdots \Expect_{\tau(j)}
	\enormsq{ \vct{x}_{j} - \vct{x}_{\star} } \\
	&\leq \gamma^{j} \enormsq{ \vct{x}_0 - \vct{x}_{\star} }
	+ \frac{1}{\alpha m} \enormsq{\vct{e}}
	\left( \sum\nolimits_{i=0}^{j-1} \gamma^i \right) \\
	&\leq \gamma^{j} \enormsq{ \vct{x}_0 - \vct{x}_{\star} }
	+ \frac{1}{\alpha m(1 - \gamma)} \enormsq{\vct{e}}.
\end{align*}
Reintroduce the value of $\gamma$ in this expression, and simplify to complete the proof.
\end{proof}

%\subsection{Reducing the Residual}
%
%Zousias and Freris proposed a method for reducing the size of the residual.
%There's a block version.  Oh.  But they use columns, not rows.
%
%Recall $\vct{e}_{\star} = \vct{b} - \mtx{A} \vct{x}_{\star}$.
%
%Consider the iteration.  $\vct{c}_0 = \vct{b}$.  At step $j$,
%choose $\tau$ uniformly at random from $T$.  Update
%$$
%\vct{c}_j = (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau}) \vct{c}_{j-1}.
%$$
%Repeat.
%
%Analysis: Let $\vct{b}_{\star} = (\Id - \mtx{P}) \vct{b}$,
%where $\mtx{P}$ is the orthogonal projector onto the range of $\mtx{A}$.
%So $\vct{b}_{\star}$ is orthogonal to the range of $\mtx{A}$.
%
%Define $\vct{r}_j = \vct{c}_j - \vct{b}_{\star}$.
%
%$$
%\vct{c}_j - \vct{b}_{\star}
%	= (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau}) (\vct{c}_{j-1} - \vct{b}_\star})
%$$
%
%
%$$
%\enormsq{ (\Id - \mtx{A}_{\tau}^\pinv \mtx{A}_\tau) \vct{u} }
%	= \enormsq{\vct{r}_{j-1}} - \enormsq{ \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} \vct{r}_{j-1} }.
%$$
%Look at second term...
%$$
%\enormsq{ \mtx{A}_{\tau}^\pinv \mtx{A}_{\tau} \vct{r}_{j-1} }
%	\geq \sigma_{\min}^2( \mtx{A}_{\tau}^\pinv ) \enormsq{ \mtx{A}_{\tau} \vct{r}_{j-1} }
%	\geq \frac{1}{\beta} 
%$$



\section{A Conversation about Pavings}\label{sec:pavings} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

As we have seen, the properties of the row paving $T$ of the matrix $\mtx{A}$ have a significant effect on the behavior of the block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}.  It is natural to ask if every standardized matrix $\mtx{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 $\mtx{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.

\begin{remark}[Paving a Square Matrix]
The operator theory literature uses the term \emph{paving} to refer to a partition $T = \{\tau_1, \dots, \tau_m \}$ of the coordinates of a \emph{square} matrix $\mtx{H}$ with a zero diagonal in which each diagonal block satisfies the bound
$$
\norm{ \mtx{H}_{\tau\tau} } \leq (1 - \delta) \norm{ \mtx{H} }
\quad\text{for each $\tau \in T$,}
$$
where the parameter $\delta \in (0, 1)$.  We can obtain row paving results for a standardized matrix $\mtx{A}$ by applying paving results for a square matrix to the hollow Gram matrix $\mtx{H} = \mtx{AA}^* - \Id$.  This approach generally leads to an estimate for the size $m$ of the paving that has an excessive dependence on the spectral norm of $\mtx{A}$.
\end{remark}

\subsection{Subset Selection Theorems}

The first approach to paving relies on a type of result called a \emph{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~\ref{prop:intro-paving} depends on the following subset selection theorem.

%The first approach to paving is based on a \emph{subset selection} theorem.  This type of result asserts that, under appropriate conditions, a matrix contains a (large) set of rows with distinguished metric properties.  Proposition~\ref{prop:intro-paving} depends on the following subset selection theorem.

\begin{proposition}[Subset Selection] \label{prop:subset-select}
Fix a number $\delta \in (0, 1)$.  Let $\mtx{A}$ be a standardized matrix with $n$ rows.  Then there exists a subset $\tau$ of row indices with the properties
$$
\abs{\tau} \geq \frac{\cnst{c_{ss}} \cdot \delta^2 n}{\normsq{\mtx{A}}}
\quad\text{and}\quad
1 - \delta \leq \lambda_{\min}(\mtx{A}_\tau \mtx{A}_{\tau}^*)
\leq \lambda_{\max}(\mtx{A}_\tau \mtx{A}_{\tau}^*) \leq 1 + \delta.
$$
The number $\cnst{c_{ss}}$ is a positive, universal constant.
\end{proposition}

Proposition~\ref{prop:subset-select} follows from~\cite[Thm.~1.2]{Tro09:Column-Subset}, once we track the parameter $\delta$ through the proof.  This result has been attributed to Bourgain and Tzafriri~\cite{BT91:Problem-Kadison-Singer}, but the earliest reference seems to be Vershynin's paper~\cite[Cor.~1.5]{Ver06:Random-Sets}.  See Section~\ref{sec:ss-related} for further background.

Proposition~\ref{prop:subset-select} ensures that each standardized matrix contains a large set of well-conditioned rows.  To construct a paving of a standardized matrix $\mtx{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
$$
m \leq \cnst{C_{pave}} \cdot \delta^{-2} \normsq{\mtx{A}} \log(1+n)
$$
steps, we have exhausted the entire matrix.  This argument yields Proposition~\ref{prop:intro-paving}.

The paper~\cite{Tro09:Column-Subset} contains an efficient computational method for identifying the subset $\tau$ promised by Proposition~\ref{prop:subset-select}.   This algorithm chooses a random set $\omega$ of rows from the matrix $\mtx{A}$ with twice the cardinality of the desired subset $\tau$.  Then it computes a matrix factorization of the submatrix $\mtx{A}_{\omega}$ that exposes a well-conditioned subset $\tau$ of rows inside $\omega$.

\subsubsection{Related Results} \label{sec:ss-related}

The literature contains a variety of subset selection theorems that can be used to construct pavings.  The first major result in this direction is the Restricted Invertibility Principle of Bourgain and Tzafriri~\cite[Thm.~1.2]{BT87:Invertibility-Large}, which guarantees that every standardized matrix contains a set of $\cnst{const} / \norm{\mtx{A}}^2$ rows whose minimal singular value is bounded away from zero.  In a similar vein, Kashin and Tzafriri~\cite{KT94:Some-Remarks} establish that every standardized matrix contains a set of $\cnst{const} / \norm{\mtx{A}}^{2}$ rows whose norm is constant.  Both results are based on random selection combined with matrix factorization.  Neither result yields a submatrix whose singular values lie arbitrarily close to one.

In another notable paper, Bourgain and Tzafriri~\cite[Cor.~1.2]{BT91:Problem-Kadison-Singer} prove that every standardized matrix contains a set of $\cnst{const} \cdot \delta^{2} /\norm{\mtx{A}}^{4}$ rows whose squared singular values fall in the range $[1 - \delta, 1 + \delta]$.  This theorem offers quantitative control on the conditioning of the submatrix, but it gives the wrong dependence on the spectral norm of the matrix.  Proposition~\ref{prop:subset-select} is based on the same type of argument as this result of Bourgain and Tzafriri.  The iterative argument we use to draw the paving result, Proposition~\ref{prop:intro-paving}, as a corollary of Proposition~\ref{prop:subset-select} also appears in the paper~\cite{BT91:Problem-Kadison-Singer}.

%, but it uses more modern tools to streamline the argument.
 
%Literature describes other subset selection theorems allow us to control the lower bound only (Restricted Invertibility Principle)~\cite{BT87:Invertibility-Large} or the upper bound only (Kashin--Tzafriri)~\cite{KT94:Some-Remarks}.  These results lead to pavings in the same fashion.

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

Vershynin~\cite{Ver01:Johns-Decompositions} has obtained a theory of subset selection for matrices that are not necessarily standardized.  In his results, the squared Frobenius norm $\fnormsq{\mtx{A}}$ plays the role of the number $n$ of rows.  Srivastava's dissertation~\cite[Chap.~3]{Sri10:Spectral-Sparsification} contains an algorithm for (weighted) subset selection that applies to general matrices.  See Youssef's works~\cite{You12:Restricted-Invertibility,You12:Note-Column} for the latest developments in this direction.

Finally, let us mention that subset selection theorems have applications throughout mathematics and engineering.  See the paper~\cite{CT06:Kadison-Singer-Problem} for a discussion and references.

\subsection{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.

\begin{definition}[Random Partition]\label{def:rp}
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
$$
\tau_i = \{\pi(k) : k = \lfloor (i-1)n/m\rfloor + 1, \lfloor (i-1)n/m\rfloor + 2, \ldots, \lfloor in/m\rfloor\}.
$$
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 \textit{random partition} of $\{1, 2, \dots, n\}$ into $m$ blocks.
%It is clear that $T = \{\tau_0, \tau_1, \ldots, \tau_{m-1}\}$ is a partition of $\{1, 2, \ldots, n\}$ into $m$ blocks of approximately equal size.  In a natural sense, each block contains random indices.  Thus, we say that $T$ is a \textit{random partition}.  
\end{definition}

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

\begin{proposition}[Random Paving: Upper Bound] \label{prop:randpave-upper}
Let $\mtx{A}$ be a tall\,%
\footnote{A $p \times q$ matrix is \emph{tall} when $p \geq q$.},
standardized matrix with $n$ rows.  Consider a randomized partition $T$ of the row indices with $m \geq \normsq{\mtx{A}}$ blocks.  Then $T$ is a row paving with upper bound $\beta \leq 6\log(1 + n)$, with probability at least $1-n^{-1}$.
\end{proposition}

Proposition~\ref{prop:randpave-upper} results from an argument based on the matrix Chernoff inequality~\cite[Thm.~1.1]{Tro12:User-Friendly} and a union bound.  A model for this type of proof appears in the paper~\cite{Tro11:Improved-Analysis}.  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~\cite[Ex.~2.2]{BT91:Problem-Kadison-Singer} of Bourgain and Tzafriri implies that we must assume the rows of the matrix $\mtx{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.

\begin{definition}[Incoherence] \label{def:incoherence}
Suppose $\mtx{A}$ is a matrix with rows $\vct{a}_1, \vct{a}_2, \ldots, \vct{a}_n$.  We say that $\mtx{A}$ is \emph{incoherent} when
$$
\max_{i \neq \ell} \ \abs{\ip{ \vct{a}_i }{ \vct{a}_\ell }} \leq\frac{\cnst{c_{inc}}}{\log(1+n)}.
$$
The number $\cnst{c_{inc}}$ is a positive, universal constant.
\end{definition}

Incoherent matrices arise, for example, in signal processing problems~\cite{DH00:Uncertainty-Principles}.  Every incoherent, standardized matrix admits a random paving with controlled lower and upper bounds.

\begin{proposition}[Random Paving]\label{prop:randpave}
Suppose that $\mtx{A}$ is an incoherent, standardized matrix with $n$ rows.  Let $T$ be a random partition of the row indices into $m$ blocks where $m \geq \cnst{C_{rand}} \cdot \delta^{-2} \normsq{\mtx{A}} \log(1+n)$.  Then $T$ is a row paving of $\mtx{A}$ whose paving bounds satisfy $1-\delta \leq \alpha \leq \beta \leq 1 + \delta$, with probability at least $1-n^{-1}$.
\end{proposition}

Proposition~\ref{prop:randpave} follows from~\cite[Cor.~5.2]{Tro08:Norms-Random}, along with some standard arguments~\cite{BT91:Problem-Kadison-Singer,Tro08:Random-Paving}.  The paper~\cite{CD12:Invertibility-Submatrices} contains superior estimates for the constants in this analysis.  See Section~\ref{sec:randpave-related} for some further background.

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

%This observation leads to a natural framework for solving incoherent, overdetermined least-squares problems:
%\begin{enumerate}
%\item	Construct a random partition $T$ using the guidelines from Proposition~\ref{prop:randpave}; and
%\item	Apply the randomized block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}, with the paving $T$.
%\end{enumerate}
%With this approach, we obtain the same convergence guarantees outlined in Corollary~\ref{cor:good-intentions}.


\subsubsection{The Fast Incoherence Transform}

Prima facie, random paving has limited applicability because the incoherence hypothesis in Proposition~\ref{prop:randpave} 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.

\begin{definition}[Fast Incoherence Transform]
The \emph{fast incoherence transform} is the $n \times n$ random matrix $\mtx{S} = \mathbf{F}\mtx{E}$ where $\mathbf{F}$ is the $n \times n$ unitary discrete Fourier transform (DFT) and $\mtx{E}$ is an $n \times n$ diagonal matrix whose entries are independent Rademacher%
\footnote{A \emph{Rademacher} random variable takes the values $\pm 1$ with equal probability.}
random variables.
\end{definition}

The matrix $\mtx{S}$ is unitary, so it preserves Euclidean structure in a problem.  Furthermore, $\mtx{S}$ can be applied to a vector in $\bigO(n \log n)$ arithmetic operations by means of the FFT algorithm.  Thus, for an $n \times d$ matrix $\mtx{A}$, we can form the product $\mtx{SA}$ with $\bigO(nd \log n)$ arithmetic operations.  When the matrix $\mtx{A}$ is sparse or admits a fast matrix--vector multiply, it may be more appropriate to work directly with the factorized representation $\mtx{SA}$, because it inherits a fast multiply from its two factors.  To emphasize, in many contexts, it is unnecessary to form $\mtx{SA}$ explicitly.

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

\begin{proposition}[Fast Incoherence Transform] \label{prop:fit}
Suppose that $\mtx{A}$ is a standardized matrix with $n$ rows whose norm satisfies
\begin{equation} \label{eqn:fit-hyp}
\normsq{\mtx{A}} \leq \frac{\cnst{c_{fit}} \cdot n}{\log^3(1+n)}.
\end{equation}
Let $\mtx{S}$ be an $n \times n$ fast incoherence transform.  Then the matrix $\mtx{W} = \mtx{SA}$ satisfies the probability bound
$$
\Prob{ \max_{i, \ell}\abs{ \ip{\vct{w}_i}{ \vct{w}_\ell} - \delta_{i\ell}} \geq \frac{\cnst{c_{inc}}}{\log(1+ n)} }
	\leq \frac{1}{n},
$$
where $\vct{w}_i$ is the $i$th row of $\mtx{W}$ and $\delta_{i\ell}$ is the Kronecker delta.
\end{proposition}

We do not have a reference for Proposition~\ref{prop:fit}, but the result is probably not new.  Appendix~\ref{app:fit} offers a short proof based on the Hanson--Wright inequality~\cite{HW71:Bound-Tail} for Rademacher chaos.

Let us pause to examine the hypothesis~\eqref{eqn:fit-hyp}.  Recall that, for a standardized matrix $\mtx{A}$ with $n$ rows, the squared spectral norm $\normsq{\mtx{A}}$ attains its maximal value $n$ when the rows are identical.  Therefore, the bound~\eqref{eqn:fit-hyp} stipulates that the rows of $\mtx{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.


\subsubsection{Incoherence, Paving, Kaczmarz} \label{sec:inc-pave-kacz}

This discussion suggests the following approach for solving the overdetermined least-squares problem~\eqref{eqn:basic-ls} with a standardized matrix $\mtx{A}$:
\begin{enumerate}
\item	Apply the fast incoherence transform $\mtx{S}$ to the objective function:
$$
\enormsq{ \mtx{A} \vct{x} - \vct{b} }
	= \enormsq{ \mtx{SA} \vct{x} - \mtx{S}\vct{b} }
	=: \enormsq{ \tilde{\mtx{A}} \vct{x} - \tilde{\vct{b}} }.
$$
\item	Draw a random partition $T$ of the rows of $\tilde{\mtx{A}}$ with the number of blocks determined according to Proposition~\ref{prop:randpave}.
\item	Apply the randomized block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}, with the matrix $\tilde{\mtx{A}}$, the right-hand side~$\tilde{\vct{b}}$, and the paving $T$ to solve the transformed problem.
\end{enumerate}
Provided that the hypothesis~\eqref{eqn:fit-hyp} is in force, Proposition~\ref{prop:fit} shows that the fast incoherence transform is likely to convert the matrix $\mtx{A}$ into an incoherent matrix $\tilde{\mtx{A}}$.  Proposition~\ref{prop:randpave} shows that we can obtain a good paving of $\tilde{\mtx{A}}$ from a random partition of the rows.  If these randomized procedures are successful, when we apply the block Kaczmarz algorithm to solve the transformed least-squares problem, we obtain the same convergence guarantees outlined in Corollary~\ref{cor:good-intentions}.

For a well-conditioned $n \times d$ matrix $\mtx{A}$ with $d \gg \log n$, we can use this approach to solve the least-squares problem to a fixed level of precision after $\bigO(nd \log n)$ arithmetic operations.  In theory, this bound is almost comparable with the conjugate gradient method, which requires $\bigO(nd)$ operations to achieve fixed precision in this setting.


\subsubsection{Related Results} \label{sec:randpave-related}

The idea that randomness might help us to construct a paving is already inherent in the Bourgain--Tzafriri paper~\cite{BT87:Invertibility-Large} on the Restricted Invertibility Principle, where they use randomized row selection and matrix factorization to perform subset selection.  A result~\cite[Thm.~2.3]{BT91:Problem-Kadison-Singer} in their subsequent paper demonstrates that, in the presence of incoherence assumptions, a random partition induces a row paving with $\cnst{const} \cdot \norm{\mtx{A}}^4 \log(1+n)$ blocks; the extra factorization step is not necessary.  Furthermore, under an incoherence assumption slightly stricter than Definition~\ref{def:incoherence}, they prove that a random partition induces a paving with $\cnst{const} \cdot \norm{\mtx{A}}^4$ blocks; no logarithmic factor is necessary.  See the paper~\cite{Tro08:Random-Paving} for a modern proof of the latter result.

The aforementioned theorems all yield the wrong dependence on the spectral norm in the size of a random row paving.  The pr{\'e}cis~\cite{Tro08:Norms-Random} shows how to obtain the correct quadratic dependence that is quoted in Proposition~\ref{prop:randpave}.  If we were to strengthen the incoherence requirement in Proposition~\ref{prop:randpave}, it is likely that we could remove the logarithmic factor from the size of the paving by adapting an argument~\cite[Prop.~2.7]{BT91:Problem-Kadison-Singer} of Bourgain and Tzafriri.  On the other hand, the logarithmic factor is necessary at the incoherence level we have imposed~\cite[Ex.~2.2]{BT91:Problem-Kadison-Singer}.

The fast incoherence transform is based on ideas of Ailon and Chazelle~\cite{AC09:Fast-Johnson-Lindenstrauss}, who use the random matrix $\mtx{S}$ to perform dimension reduction.  We believe that the first application of $\mtx{S}$ for randomized linear algebra appears in the paper Woolfe et al.~\cite{WLRT08:Fast-Randomized}, where they use this transform to aid in computing matrix decompositions.  See the works~\cite{AMT11:Blendenpik-Supercharging,HMT11:Finding-Structure,BG12:Improved-Matrix} for further results in this direction.  Liberty's dissertation~\cite{Lib09:Accelerated-Dense} describes other randomized maps that can play a similar role.


\section{Numerical Aspects of Block Kaczmarz} \label{sec:implementation}

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~\ref{sec:why}, we outline some of the situations where we expect the block Kaczmarz method to outperform the simple Kaczmarz algorithm.  Afterward, in Section~\ref{sec:how}, we consider some of the questions that arise when implementing the block Kaczmarz method.  Finally, in Section~\ref{sec:experiments}, we offer some simple computational examples to illustrate our main points.

\subsection{Why Use Block Kaczmarz?} \label{sec:why}

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~\cite{FS95:Kaczmarz-Based-Approach}.  See Section~\ref{sec:experiments-circ} for a related numerical example.

%There are some circumstances where there is a natural row partition, e.g.,
%multi-sampling schemes, domain decomposition, red-black ordering for Poisson equation.

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 ({\tt 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 {\tt BLAS2}, rather than {\tt BLAS1}.  A more detailed discussion of these issues falls outside the scope of this paper.

\subsection{Implementing Block Kaczmarz Methods} \label{sec:how}

The randomized block Kaczmarz method, Algorithm~\ref{alg:randomized-block-kaczmarz}, is easy to describe, but there remain several implementation issues that require attention.

\subsubsection{The Block Update Rule}

Most of the arithmetic in the algorithm occurs when we apply the pseudoinverse $\mtx{A}_\tau^\pinv$ to a vector in the update~\eqref{eqn:block-kacz}.  Equivalently, we must solve an (underdetermined) least-squares problem at each iteration.  The appropriate numerical method depends on a wide variety of issues~\cite{Bjo96:Numerical-Methods}, so we cannot give a universal prescription.
Here are some factors worth considering.
\begin{itemize}
\item	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.

\item	In case the conditioning $\beta/\alpha$ of the paving is small, we recommend using an iterative method, such as {\tt CGLS}, {\tt 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.
\end{itemize}

%\notate{More precise...}
%First, if the matrix $\mtx{A}_\tau$ has condition number $K$, each iteration of the conjugate gradient algorithm reduces the error by a factor no smaller than $(K - 1)/(K + 1)$.  Thus when $\mtx{A}_\tau$ is well conditioned, the conjugate gradient offers
%exceptional performance.  Second, as the Kaczmarz algorithm proceeds, the residual $(\vct{b}_\tau - \mtx{A}_\tau
%\vct{x}_j)$ approaches $\range(\mtx{A}_\tau)$, so conjugate gradient benefits from having an initial iterate with a small error.  Third, conjugate gradient uses only matrix--vector multiplications, so it becomes even more efficient for a structured or sparse matrix $\mtx{A}$.

\subsubsection{Setting the Convergence Tolerance}


Theorem~\ref{thm:convergence} implies that Algorithm~\ref{alg:randomized-block-kaczmarz} can reduce the error $\enormsq{\hat{\vct{x}} - \vct{x}_{\star}}$ 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 $\enormsq{ \mtx{A}\hat{\vct{x}} - \vct{b} } \leq \eps^2$ with the decomposition $\vct{b} = \mtx{A} \vct{x}_{\star} - \vct{e}$ to obtain
$$
\enormsq{ \mtx{A}(\hat{\vct{x}} - \vct{x}_{\star}) + \vct{e} } \leq \eps^2.
$$
The residual $\vct{e}$ is orthogonal to the range of $\mtx{A}$, so we can apply the Pythagorean theorem and bound the $\ell_2$ norm below to reach
\begin{equation*} %\label{eqn:perturbation}
\sigma_{\min}^2(\mtx{A}) \enormsq{ \hat{\vct{x}} - \vct{x}_{\star} }
	+ \enormsq{\vct{e}}
	\leq \eps^2.
\end{equation*}
It follows that we must set the convergence tolerance $\eps$ so that
\begin{equation} \label{eqn:conv-tol}
\eps^2 > \left[ 1 + \frac{\beta}{\alpha} \right] \enormsq{\vct{e}}.
\end{equation}
Otherwise, we have no guarantee that the algorithm will terminate.


\subsubsection{Checking for Convergence}

The convergence criterion $\normsq{\mtx{A}\hat{\vct{x}} - \vct{b}} \leq \eps^2$ may be somewhat expensive to verify because it involves a multiply with the full matrix $\mtx{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.


\subsubsection{Randomized Cyclic Control} \label{sec:rand-cyclic}

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~\ref{alg:randomized-block-kaczmarz} samples a uniformly random block at each iteration, independent of all previous choices.  Instead, we recommend using an alternative scheme that samples blocks \emph{without} replacement.  We can express this method formally as follows.

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

\item	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)}$.
\end{enumerate}

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

In Section~\ref{sec:experiments-circ}, we offer some numerical evidence that this alternative control scheme is more effective than the approach used in Algorithm~\ref{alg:randomized-block-kaczmarz}.  At present, compelling explanations for this phenomenon are lacking.  See~\cite{RR12:Beneath-Valley} for some discussion and conjectures.

%That is, we draw a random permutation $\pi$ of the blocks, and we cycle through the blocks in order: $\tau(j) = \tau_{\pi(j \bmod m)}$.  See Section~\ref{sec:experiments-circ} for some evidence.  

%\notate{This section should go in the discussion.}



\subsection{Numerical Experiments}\label{sec: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~\ref{alg:randomized-block-kaczmarz}.  In Section~\ref{sec:experiments-circ},
we consider an example where the block method has a clear advantage over the simple method.
In Section~\ref{sec:experiments-dense}, we examine some other cases where the benefits
of the block method are more subtle.

\subsubsection{Submatrices with Fast Multiplies} \label{sec:experiments-circ}

First, we study the situation where the matrix $\mtx{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~\ref{sec:fast-update},
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 $\mtx{A}_{\rm circ}$ by stacking $15$ partial random circulant
matrices:
$$
\mtx{A}_{\rm circ} = \begin{bmatrix} \mtx{C}_1 \\ \vdots \\ \mtx{C}_{15} \end{bmatrix}
\quad\text{where}\quad
\mtx{C}_{i} = \mtx{R} \mathbf{F}^* \mtx{E}_i \mathbf{F}
\quad\text{for $i = 1, \dots, 15$.}
$$
In this expression,
\begin{itemize}
\item	$\mtx{R}$ is the restriction to the first $20$ coordinates,
\item	$\mathbf{F}$ is the $100 \times 100$ unitary DFT, and
\item	$\mtx{E}_{i}$ is a random diagonal matrix whose entries are independent Rademacher
random variables.  Each $\mtx{E}_i$ is drawn independently from the others.
\end{itemize}
This construction ensures that each $\mtx{C}_i$ is a section of a circulant matrix
with orthonormal rows.
As a consequence, the pseudoinverse equals the adjoint: $\mtx{C}_i^\pinv = \mtx{C}_i^*$.
Furthermore, we can apply either $\mtx{C}_i$ or $\mtx{C}_i^*$ to a vector quickly%
\footnote{Recall 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$.

%Thus, we can apply the block update
%rule~\eqref{eqn:block-kacz} with $\mtx{A}_{\tau} = \mtx{C}_i$ in about
%$8 d \log_2 d + 4d$ floating-point operations (flops), where $d = 100$.


Using this matrix, we perform a small {\sc Matlab} experiment to compare the behavior of randomized block Kaczmarz and randomized simple Kaczmarz.
\begin{enumerate}
\item	Draw a random matrix $\mtx{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. 
\item	Let $\vct{x}_{\star} = [1, \dots, 1]^*$, and form $\vct{b} = \mtx{A} \vct{x}_{\star}$.
Set the initial iterate $\vct{x}_0 = \vct{0}$.
\item	For each of 100 trials,
\begin{enumerate}
\item	Apply the simple Kaczmarz method from Section~\ref{sec:simple-kacz} to produce iterates $\{ \vct{x}_j^{\rm simp} : j \geq 0 \}$.
\item	Apply Algorithm~\ref{alg:randomized-block-kaczmarz} to produce iterates $\{ \vct{x}_j^{\rm block} : j \geq 0 \}$.
%\item	Halt each algorithm with the (actual) approximation error $\enorm{\vct{x}_j^{\rm alg} - \vct{x}_{\star}}$ reaches $10^{-11}$.
\end{enumerate}
\item	For each algorithm, at each iteration $j$, compute the minimum, median, and maximum value of the approximation error $\enorm{\vct{x}^{\rm alg}_j - \vct{x}_{\star}}$ over the 100 trials.
\end{enumerate}
In this experiment, we form the matrix $\mtx{A}_{\rm circ}$ in the first step and store it.  To perform the update~\eqref{eqn:simple-kacz}, 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~\eqref{eqn:block-kacz} with about $4 d \log_2(d) + 4d$ complex flops.

In Figure~\ref{fig:circ} [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~\ref{sec:rand-cyclic}, 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~\ref{fig:circ} [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.


\begin{figure}[t]
\begin{center}
\includegraphics[height=2.5in]{CircErrPerFlop-fixed2.pdf}
\includegraphics[height=2.5in]{CircWithVersusWithout-fixed.pdf}
\caption{({Convergence with a Block Circulant Matrix}){\bf .}
The matrix $\mtx{A}_{\rm circ}$ is a fixed $300 \times 100$ matrix consisting
of $15$ partial circulant blocks.
For each algorithm, we chart the approximation
error $\enorm{ \vct{x}_j - \vct{x}_{\star} }$ as a function of the
number of complex flops performed.
The heavy line denotes the median error over 100 independent trials;
the minimum and maximum errors envelop the shaded region.
See Section~\ref{sec:experiments-circ} for details.
{\bf [{Left}]} Convergence of the randomized block Kaczmarz method,
Algorithm~\ref{alg:randomized-block-kaczmarz}, versus the
randomized simple Kaczmarz method (Section~\ref{sec:simple-kacz}).
{\bf [{Right}]}  Convergence of two variants  (Section~\ref{sec:rand-cyclic})
of the randomized block Kaczmarz method.
One samples blocks independently with replacement;
the other samples blocks without replacement in each epoch.
} \label{fig:circ}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[height=1.75in]{GaussErrPerEpoch-final.pdf}
\raisebox{0.025in}{\includegraphics[height=1.7in]{GaussErrPerFlop-final.pdf}
\includegraphics[height=1.7in]{GaussErrPerCPU-final.pdf}}
\caption{({Convergence with a Random Matrix}){\bf .}
The matrix $\mtx{A}_{\rm rdm}$ is a fixed $300 \times 100$ matrix whose rows are drawn
uniformly at random from the Euclidean unit sphere and partitioned into 10
blocks of equal size.
For the iterates generated by each algorithm, we plot the decay of the approximation error
$\enorm{\vct{x}_j - \vct{x}_{\star}}$
as a function of three computational resources.
The heavy line denotes the median error over 100 independent trials;
the minimum and maximum errors envelop the shaded region.
{\bf [{Left}]} Approximation error as a function of the number of epochs
elapsed.
{\bf [{Center}]} Approximation error as a function of flops.
{\bf [{Right}]} Approximation error as a function of CPU time.
See Section~\ref{sec:experiments-dense} for details.}
\label{fig:dense}
\end{center}
\end{figure}


\subsubsection{Unstructured Submatrices} \label{sec:experiments-dense}

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.

First, we consider a $300 \times 100$ real matrix $\mtx{A}_{\rm rdm}$ whose rows are drawn
independently and uniformly at random from the $\ell_2$ unit sphere.  We partition
the rows of this matrix into $10$ equal blocks of $30$ consecutive rows each.
A heuristic application of the Bai--Yin Law~\cite[Thm.~2]{BY93:Limit-Smallest}
shows that the singular values of each $30 \times 100$ submatrix fall in the
range $1 \pm \sqrt{30/100}$.  Thus, the paving parameters $\alpha \approx 0.2$
and $\beta \approx 2.4$, and the condition number of each block is bounded by
$\beta / \alpha \approx 4.8$.

We perform an experiment with this matrix that follows the same methodology as
in Section~\ref{sec:experiments-dense}.
%In this case, the simple
%update rule~\eqref{eqn:simple-kacz} costs $4d$ (real) flops, where $d = 100$.
To implement the block Kaczmarz update rule~\eqref{eqn:block-kacz}, we use
the {\tt 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.
% to ensure that the plots are accurate.

Figure~\ref{fig:dense} 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~\ref{sec:unstruct},
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 {\sc Matlab} script.
Nevertheless, we regard this experiment as limited evidence
of the computational advantages of the block method that we outlined in Section~\ref{sec:why}.

\begin{wrapfigure}{r}{0.5\textwidth}
\begin{center}
\includegraphics[width=.5\textwidth]{UnifPoint5to1ErrPerFlop-final.pdf}
\caption{({Convergence with a Coherent Matrix}){\bf .}
The matrix $\mtx{A}_{\rm coh}$ is a fixed $300 \times 100$ matrix with strongly correlated
rows that are partitioned into 10 blocks of equal size.
We compare the approximation error $\enorm{\vct{x}_j - \vct{x}_{\star}}$ for
the iterates generated by each algorithm as a function of the number of flops
performed.  The heavy line denotes the median error over 100 trials;
the minimum and maximum error describe the boundaries of the shaded region.
See Section~\ref{sec:experiments-dense} for details.}
\label{fig:coh}
\end{center}
%\vspace{-0.5pc}
\end{wrapfigure}


Finally, we consider one other type of unstructured matrix, where our analysis
offers little guidance.  Let $\tilde{\mtx{A}}_{\rm coh}$ be a $300 \times 100$ matrix whose
entries are drawn independently at random from the uniform distribution on $[0.5, 1]$.
We form a standardized matrix $\mtx{A}_{\rm coh}$ by normalizing the rows of
$\tilde{\mtx{A}}_{\rm coh}$.
Since this matrix is dense and positive, the rows of this matrix tend to be strongly
correlated with each other; the maximum inner product between rows is typically $0.98$
or higher.  Furthermore, the norm of this matrix is usually close to the maximum
possible value $\sqrt{300}$.
  Our theory provides no map for this
\emph{terra incognita}.
Yet we brazenly partition the rows into $10$ equal blocks, each with $30$ contiguous rows,
as in the previous example.

We perform the same experiment for $\mtx{A}_{\rm coh}$ as we did with
$\mtx{A}_{\rm rdm}$ to compare the behavior of the simple
Kaczmarz method and the randomized Kaczmarz method.
Figure~\ref{fig:coh} 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~\cite{NW12:Two-Subspace-Projection} 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~\ref{alg:randomized-block-kaczmarz} when the blocks are larger.
%This is a subject for further investigation.



%\subsubsection{Sampling without Replacement} \label{sec:experiments-randperm}


\section{Related Work and Future Directions} \label{sec:future}

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.


\subsection{Kaczmarz Methods with Randomized Control}

%\notate{deterministic vs. random control}
The Kaczmarz method was originally introduced in the paper~\cite{Kac37:Angenaeherte-Aufloesung}.  It was reinvented by researchers in tomography~\cite{GBH70:Algebraic-Reconstruction} under the appellation ``algebraic reconstruction technique'' (ART).  See Byrne's book~\cite{Byr08:Applied-Iterative} 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~\cite{HS78:Angles-Null}.
% and the analysis of deterministic Kaczmarz algorithms tends to be rather involved.
The literature contains empirical evidence~\cite{HM93:Algebraic-Reconstruction} that \emph{randomized} control mechanisms may be more effective, but until recently there was no compelling theoretical analysis to support this observation.
%\notate{Add cites to this paragraph.}

The paper~\cite{SV09:Randomized-Kaczmarz} 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.,~\cite[Eqn.~(1.2)]{XZ02:Method-Alternating}, whose geometric meaning is not evident.

%This work contains clean bounds for the convergence rate in terms of the metric properties of the matrix $\mtx{A}$. The classical literature does not seem to contain any results of this type for a sequential projection method.

%\notate{classical analysis of kaczmarz vs. new randomized approach}
    % $\mtx{A}$.

%\notate{Other work on randomized Kaczmarz}
In the wake of Strohmer and Vershynin's work~\cite{SV09:Randomized-Kaczmarz}, 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~\cite{Nee10:Randomized-Kaczmarz}.
%Chen and Powell~\cite{CP12:Almost-Sure-Convergence} show that the algorithm converges almost surely under some additional assumptions.
Zouzias and Freris~\cite{ZF12:Randomized-Extended} exhibit a randomized procedure, based on ideas from~\cite{Pop98:Extensions-Block-Projections}, that can reduce the size of the residual $\vct{e}$.  Leventhal and Lewis~\cite{LL10:Randomized-Methods} provide an analysis of a randomized iteration for solving least-squares problems with polyhedral constraints, while Richt{\'a}rik and Tak{\'a}{\v c} have extended these ideas to more general optimization problems~\cite{RT11:Iteration-Complexity}.  Some other references include~\cite{EN11:Acceleration-Randomized,RT11:Iteration-Complexity,CP12:Almost-Sure-Convergence,NW12:Two-Subspace-Projection}.


%Leventhal and Lewis extend the work of Strohmer and Vershynin by showing exponential convergence in expectation of the randomized Kaczmarz algorithm for systems of linear inequalities~\cite{LL10:Randomized-Methods}.  They prove that for feasible systems of linear inequalities, the randomized Kaczmarz method exhibits a convergence rate of $1 - \frac{1}{L\|\mtx{A}\|_F^2}$ where $L$ denotes the Hoffman constant~\cite{Hof52:Approximate-Solutions}.  Note that in the case where the feasible set consists of a single point (i.e. the system is equivalent to a consistent system of linear \textit{equalities}), that $L = \|\mtx{A}^{-1}\|$ when $\mtx{A}$ has full column rank.  This work thus recovers the same convergence rate as that of Strohmer and Vershynin in this case.

\subsection{Block Kaczmarz Methods}

The block Kaczmarz update rule~\eqref{eqn:block-kacz} we are studying is originally due to Elfving~\cite[Eqn.~(2.2)]{Elf80:Block-Iterative-Methods}.  This update is a special case of a general framework due to Eggermont et al.~\cite{EHL81:Iterative-Algorithms}.  Byrne describes a number of other block Kaczmarz methods in his book~\cite[Chap.~9]{Byr08:Applied-Iterative}.

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~\cite{XZ02:Method-Alternating}, which contains a refined convergence analysis that applies to a wide range of algorithms.  Nevertheless, to our knowledge, Algorithm~\ref{alg:randomized-block-kaczmarz} is the only block Kaczmarz method that offers an (expected) linear rate of convergence that depends explicitly on geometric properties of the system matrix $\mtx{A}$ and its submatrices $\mtx{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~\cite{Pop99:Block-Projections-Algorithms} 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~\cite{Pop01:Fast-Kaczmarz-Kovarik,Pop04:Kaczmarz-Kovarik-Algorithm}.  We believe that our work is the first to recognize the natural connection between the paving literature and the block Kaczmarz method.

%Talk about Popa's stuff...
%
%In the paper~\cite{Pop98:Extensions-Block-Projections}, Popa proves that, when the matrix $\mtx{A}$ is rank-deficient and the least-squares problem~\eqref{eqn:basic-ls} is inconsistent, that the block Kaczmarz algorithm with cyclic control converges to the least-squares solution whenever the blocks are nonsingular and the initial estimate $\vct{x}_0$ is chosen to be a row of the matrix (or zero).  Popa has also developed a method for partitioning the matrix $\mtx{A}$ into blocks of mutually orthogonal rows, and he uses these decompositions to implement a block projection method~\cite{Pop99:Block-Projections}.  His work includes an empirical demonstration that block Kaczmarz methods often require significantly fewer iterations than simple Kaczmarz methods.

% ( standard counterparts (with relaxation parameters).

%Byrne has introduced an iterative block method for finding the least-squares solution to an inconsistent system of equations that arises in image reconstruction~\cite{Byr97:Convergent-Block-Iterative}.  In subsequent work~\cite{Byr08:Applied-Iterative,Byr09:Block-Iterative-Algorithms}, he describes a block variant of the algebraic reconstruction technique (BI-ART), similar to the block Kaczmarz method, that uses a relaxation parameter to accelerate convergence.  He proves that, for a consistent linear system, when the relaxation parameter is bounded by the reciprocal of the spectral radius of each matrix $\mtx{A}_{\tau}^* \mtx{A}_{\tau}$, that the method converges to the solution.  He also argues heuristically that, with a proper choice of relaxation parameters, the convergence rate of the block algorithm should improve on the convergence rate of the standard version.

%
%also provides an analysis of iterative block methods for finding the least squares solution to an inconsistent system of equations for image reconstruction~\cite{byrne1997convergent}.  He describes a block variant of the algebraic reconstruction technique (BI-ART), analogous to the block Kaczmarz method, using a relaxation parameter~\cite{byrne2008applied,byrne2009block}.  He proves that when the relaxation parameter is bounded by the inverse of the spectral radius of $\mtx{A}^\dagger_{\tau}\mtx{A}_{\tau}$ for each block $\mtx{A}_{\tau}$, that the method converges to the solution of the consistent linear system.  He also argues heuristically that with proper choice of relaxation parameters, the convergence rate of the block version of ART should improve upon the standard version by a constant factor.  However, these parameters may in general be difficult to compute.

%There are certain modern applications that lead to large systems of linear equations with a natural block structure.  For example, in wireless communication, it is common to consider signals with a bandwidth up to one Gigahertz.  In this application, it is necessary to use analog-to-digital converters (ADCs) that sample at very high rates.  It can be more cost-effective to use \emph{time-interleaved} ADCs, which sample 
%
%  In modern applications, it may be quite natural to block rows for large systems of equations.  Advanced wireless communication, for example, exhibits a heavily increasing demand for data, with transmission signals having a bandwidth up to one gigahertz.  In applications like these, it is necessary for analog to digital converters (ADC) to sample at very high rates.  Rather than sampling at very high frequencies with one expensive ADC, one uses many \textit{time-interleaved} ADCs which sample at less frequent rates in parallel.  In practice, each ADC experiences its own time-shift so that the collection of sample times is a nonuniformly shifted union of uniform sampling points.  This application called \textit{periodic non-uniform sampling} lends itself well to block methods, each block corresponding to an individual ADC (see e.g.~\cite{ST06:Fast-Reconstruction}).  In this setting, Feichtinger and Strohmer demonstrate that through the use of efficient matrix operations, the block Kaczmarz method significantly outperforms the standard method in terms of convergence rate and flop count~\cite{FS95:Kaczmarz-Based-Approach}.

\subsection{Some Future Directions}

There are a number of interesting open questions connected with Algorithm~\ref{alg:randomized-block-kaczmarz}.  First, empirical experiments make it clear that a control strategy based on sampling \emph{without} replacement is far more effective than our strategy based on sampling \emph{with} replacement.  At present, there is no compelling explanation for this phenomenon (but see~\cite{RR12:Beneath-Valley}).  Second, we think that other types of block Kaczmarz update rules~\cite[Chap.~9]{Byr08:Applied-Iterative} would also submit to a convergence analysis similar to the proof of Theorem~\ref{thm:convergence}.
%Second, it would be valuable to understand the effect of relaxation parameters on this method.
Third, it might be worthwhile to extend the argument of Zouzias and Freris~\cite{ZF12:Randomized-Extended} 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~\cite{Gau95:Theoria-Combinationis,Ben09:Key-Moments}.  Many other algorithms for numerical linear algebra~\cite{GVL96:Matrix-Computations} and least-squares problems~\cite{Bjo96:Numerical-Methods} admit natural block variants.  The field of optimization also contains a wide variety of block methods, such as~\cite{AC89:Block-Iterative-Projection,Tse93:Dual-Coordinate}.

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~\cite{LL10:Randomized-Methods}.   Using the techniques in this paper, we can easily bound the rate of convergence for this algorithm in terms of the properties of a \emph{column paving} of the system matrix.  It would be interesting to pursue this example and others.



%\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\appendix

\section{The Fast Incoherence Transform} \label{app:fit}

The goal of this Appendix is to prove Proposition~\ref{prop:fit}, which states that the fast incoherence transform makes a standardized matrix incoherent with high probability.  To that end, suppose that $\mtx{A}$ is a matrix with $n$ unit-norm rows, denoted $\vct{a}_1, \dots, \vct{a}_n$.  Recall that the fast incoherence transform is the matrix $\mtx{S} = \mathbf{F}\mtx{E}$, where $\mathbf{F}$ is the $n \times n$ unitary DFT and $\mtx{E}$ is a diagonal matrix whose entries $\xi_1, \dots, \xi_n$ are independent Rademacher random variables.

Consider the matrix $\mtx{W} := \mtx{SA}$, and introduce the Gram matrix of its rows $\mtx{G} := \mtx{WW}^* = \mtx{SAA}^* \mtx{S}^*$. Expanding this product, we find that the $(i, \ell)$ entry of $\mtx{G}$ takes the form
$$
g_{i\ell} = \sum\nolimits_{jk} \xi_j \xi_k \cdot
	{\rm f}_{ij} \bar{\rm f}_{\ell k} \cdot \ip{ \vct{a}_j }{ \vct{a}_k },
$$
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,
$$
\Expect g_{i\ell} = \sum\nolimits_j {\rm f}_{ij} \bar{\rm f}_{\ell k} \enormsq{\vct{a}_j}
	= \ip{ \mathbf{f}_{i} }{ \mathbf{f}_\ell } = \delta_{i\ell},
$$
where $\delta_{i\ell}$ is the Kronecker delta.  Note that we have invoked the standardization assumption here.

The real content of the argument is to obtain a bound on \emph{how much} the entries of the Gram matrix deviate from their expected values.  Fix a pair $(i, \ell)$ of indices, not necessarily distinct, and define the random variable
$$
Y := \abs{ g_{i\ell} - \delta_{i\ell} }.
$$
We can rewrite $Y$ is a more symmetric manner:
\begin{equation} \label{eqn:Y-chaos}
Y = \abs{ \sum\nolimits_{j \neq k} \xi_j \xi_k \cdot y_{jk} }
\quad\text{where}\quad
y_{jk} := \frac{1}{2} \left( {\rm f}_{ij} \bar{\rm f}_{\ell k}
	+ \bar{\rm f}_{ik} {\rm f}_{\ell j}  \right) \ip{ \vct{a}_j }{ \vct{a}_k }.
\end{equation}
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~\cite{HW71:Bound-Tail}; see~\cite[Chap.~8]{FR12:Mathematical-Introduction} for a modern proof.

\begin{proposition}[Hanson--Wright] \label{prop:hw}
Consider the chaos variable $Y$ defined in~\eqref{eqn:Y-chaos}.  Let $\mtx{Y}$ be the Hermitian matrix whose $(j, k)$ entry is $y_{jk}$ and whose diagonal entries are zero.  Then
$$
\Prob{ Y \geq t }
	\leq 2 \exp\left\{ - \cnst{c_{hw}} \cdot \min\left\{ \frac{t}{\norm{\mtx{Y}}}, \frac{t^2}{\fnormsq{\mtx{Y}}} \right\} \right\}
\quad\text{for $t \geq 0$.}
$$
The number $\cnst{c_{hw}}$ is a positive, universal constant.
\end{proposition}

To apply this result, we need to obtain bounds for the Frobenius norm and spectral norm of the matrix $\mtx{Y}$.  Note that
$$
\mtx{Y} = \frac{1}{2} (\mtx{Z} + \mtx{Z}^*)
\quad\text{where}\quad
\mtx{Z} := \diag(\mathbf{f}_i) \cdot ( \mtx{AA}^* - \Id ) \cdot \diag(\mathbf{f}_\ell)^*,
$$
and $\diag(\cdot)$ converts a vector into a diagonal matrix.  We may now calculate the required norms of $\mtx{Y}$.  First,
\begin{equation} \label{eqn:Y-specnorm}
\norm{\mtx{Y}} \leq \norm{\mtx{Z}}
	\leq \norm{\diag(\mathbf{f}_i)} \cdot \norm{\mtx{AA}^* - \Id} \cdot \norm{\diag( \mathbf{f}_\ell )}
	= \frac{1}{n} \norm{\mtx{AA}^* - \Id}
	\leq \frac{1}{n} \normsq{\mtx{A}}.
\end{equation}
The first inequality depends on the convexity of the spectral norm and its invariance under the conjugate transpose.  The third relation holds because the entries of the vectors $\mathbf{f}_i$ and $\mathbf{f}_\ell$ all have magnitude $n^{-1/2}$.  The last inequality follows because $\norm{\mtx{AA}^* - \Id} \leq \max\{ \normsq{\mtx{A}} - 1, 1 \} \leq \normsq{\mtx{A}}$.  For similar reasons,
\begin{equation} \label{eqn:Y-frobnorm}
\fnormsq{\mtx{Y}} \leq \fnormsq{\mtx{Z}}
	= \frac{1}{n^2} \fnormsq{\mtx{AA}^* - \Id}
	< \frac{1}{n^2} \fnormsq{\mtx{AA}^*}
	\leq \frac{1}{n} \normsq{\mtx{A}}.
\end{equation}
The strict inequality holds because $\mtx{AA}^*$ has a unit diagonal, which the identity matrix cancels off.  The final bound follows from the interpolation $\fnormsq{\mtx{P}} \leq \trace(\mtx{P}) \norm{\mtx{P}}$, valid when $\mtx{P}$ is positive semidefinite.

Next, let us instate the hypothesis~\eqref{eqn:fit-hyp} from Proposition~\ref{prop:fit}:
$$
\normsq{\mtx{A}} \leq \frac{\cnst{c_{fit}} \cdot n}{\log^3(1+n)}.
$$
Introduce this assumption into the bounds~\eqref{eqn:Y-specnorm} and~\eqref{eqn:Y-frobnorm}, and apply the Hanson--Wright inequality, Proposition~\ref{prop:hw}, to reach
$$
\Prob{ Y \geq t } \leq 2 \exp\left\{ -\frac{\cnst{c_{hw}}}{\cnst{c_{fit}}} \cdot \log^3(1+n)
	\cdot \min\{ t, t^2 \} \right\}.
$$
We may set the constant $\cnst{c_{fit}} := \cnst{c_{hw}} \cnst{c_{inc}}^2 / 3$.  For the choice $t = \cnst{c_{inc}} / \log(1+n)$, we discover that
$$
\Prob{ Y \geq \frac{\cnst{c_{inc}}}{\log(1+n)} } \leq 2 (1 + n)^{-3}.
$$
In other words, it is unlikely that any single pair of rows from $\mtx{W} = \mtx{SA}$ has an inner product much different from its expectation.

To complete the argument, we unfix the pair $(i, \ell)$ of indices.  Forming a union bound over all $n(n+1)/2$ choices where $i \leq \ell$, we conclude that
$$
\Prob{ \max_{i,\ell} \abs{g_{i\ell} - \delta_{i\ell}} \geq \frac{\cnst{c_{inc}}}{\log(1+n)}}
	\leq (1+n)^{-1}.
$$
Therefore, with high probability $\mtx{W}$ is an incoherent matrix that is nearly standardized.


\section*{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 {\sc 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.


%\notate{Deanna: Any support?}



\bibliography{rk-mathscinet}
\bibliographystyle{myalpha}

%\printbibliography

\end{document}

NT-rev.bbl
\newcommand{\etalchar}[1]{$^{#1}$}
\begin{thebibliography}{BHKW88}

\bibitem[AC89]{AC89:Block-Iterative-Projection}
R.~Aharoni and Y.~Censor.
\newblock Block-iterative projection methods for parallel computation of
  solutions to convex feasibility problems.
\newblock {\em Linear Algebra Appl.}, 120:165--175, 1989.

\bibitem[AC09]{AC09:Fast-Johnson-Lindenstrauss}
N.~Ailon and B.~Chazelle.
\newblock The fast {J}ohnson-{L}indenstrauss transform and approximate nearest
  neighbors.
\newblock {\em SIAM J. Comput.}, 39(1):302--322, 2009.

\bibitem[AMT10]{AMT11:Blendenpik-Supercharging}
H.~Avron, P.~Maymounkov, and S.~Toledo.
\newblock Blendenpik: supercharging {L}apack's least-squares solver.
\newblock {\em SIAM J. Sci. Comput.}, 32(3):1217--1236, 2010.

\bibitem[And79]{And79:Extensions-Restrictions}
J.~Anderson.
\newblock Extensions, restrictions, and representations of states on
  {$C\sp{\ast} $}-algebras.
\newblock {\em Trans. Amer. Math. Soc.}, 249(2):303--329, 1979.

\bibitem[Ben09]{Ben09:Key-Moments}
M.~Benzi.
\newblock Key moments in the history of numerical analysis.
\newblock Presented at 2009 SIAM Applied Linear Algebra Conference, Oct. 2009.

\bibitem[BG12]{BG12:Improved-Matrix}
C.~Boutsidis and A.~Gittens.
\newblock Improved matrix algorithms via the {S}ubsampled {R}andomized
  {H}adamard {T}ransform.
\newblock Available at \url{arXiv:1204.0062}, Apr. 2012.

\bibitem[BHKW88]{BHKW88:Matrix-Norm}
K.~Berman, H.~Halpern, V.~Kaftal, and G.~Weiss.
\newblock Matrix norm inequalities and the relative {D}ixmier property.
\newblock {\em Integral Equations Operator Theory}, 11(1):28--48, 1988.

\bibitem[Bj{\"o}96]{Bjo96:Numerical-Methods}
{\AA}.~Bj{\"o}rck.
\newblock {\em Numerical methods for least squares problems}.
\newblock Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
  PA, 1996.

\bibitem[BT87]{BT87:Invertibility-Large}
J.~Bourgain and L.~Tzafriri.
\newblock Invertibility of ``large'' submatrices with applications to the
  geometry of {B}anach spaces and harmonic analysis.
\newblock {\em Israel J. Math.}, 57(2):137--224, 1987.

\bibitem[BT91]{BT91:Problem-Kadison-Singer}
J.~Bourgain and L.~Tzafriri.
\newblock On a problem of {K}adison and {S}inger.
\newblock {\em J. Reine Angew. Math.}, 420:1--43, 1991.

\bibitem[BY93]{BY93:Limit-Smallest}
Z.~D. Bai and Y.~Q. Yin.
\newblock Limit of the smallest eigenvalue of a large-dimensional sample
  covariance matrix.
\newblock {\em Ann. Probab.}, 21(3):1275--1294, 1993.

\bibitem[Byr08]{Byr08:Applied-Iterative}
C.~L. Byrne.
\newblock {\em Applied iterative methods}.
\newblock A K Peters Ltd., Wellesley, MA, 2008.

\bibitem[CD12]{CD12:Invertibility-Submatrices}
S.~Chr{\'e}tien and S.~Darses.
\newblock Invertibility of random submatrices via tail decoupling and a matrix
  {C}hernoff inequality.
\newblock {\em Statist. Probab. Lett.}, 82(7):1479--1487, 2012.

\bibitem[CFM{\etalchar{+}}92]{CFMSS92:New-Variants}
C.~Cenker, H.~G. Feichtinger, M.~Mayer, H.~Steier, and T.~Strohmer.
\newblock New variants of the {POCS} method using affine subspaces of finite
  codimension, with applications to irregular sampling.
\newblock In {\em Proc. SPIE: Visual Communications and Image Processing},
  pages 299--310, 1992.

\bibitem[CHJ09]{CHJ09:Note-Behavior}
Y.~Censor, G.~T. Herman, and M.~Jiang.
\newblock A note on the behavior of the randomized {K}aczmarz algorithm of
  {S}trohmer and {V}ershynin [mr2500924].
\newblock {\em J. Fourier Anal. Appl.}, 15(4):431--436, 2009.

\bibitem[CP12]{CP12:Almost-Sure-Convergence}
X.~Chen and A.~Powell.
\newblock Almost sure convergence of the {K}aczmarz algorithm with random
  measurements.
\newblock {\em J. Fourier Anal. Appl.}, pages 1--20, 2012.
\newblock 10.1007/s00041-012-9237-2.

\bibitem[CT06]{CT06:Kadison-Singer-Problem}
P.~G. Casazza and J.~C. Tremain.
\newblock The {K}adison-{S}inger problem in mathematics and engineering.
\newblock {\em Proc. Natl. Acad. Sci. USA}, 103(7):2032--2039 (electronic),
  2006.

\bibitem[DH01]{DH00:Uncertainty-Principles}
D.~L. Donoho and X.~Huo.
\newblock Uncertainty principles and ideal atomic decomposition.
\newblock {\em IEEE Trans. Inform. Theory}, 47(7):2845--2862, 2001.

\bibitem[EHL81]{EHL81:Iterative-Algorithms}
P.~P.~B. Eggermont, G.~T. Herman, and A.~Lent.
\newblock Iterative algorithms for large partitioned linear systems, with
  applications to image reconstruction.
\newblock {\em Linear Algebra Appl.}, 40:37--67, 1981.

\bibitem[Elf80]{Elf80:Block-Iterative-Methods}
T.~Elfving.
\newblock Block-iterative methods for consistent and inconsistent linear
  equations.
\newblock {\em Numer. Math.}, 35(1):1--12, 1980.

\bibitem[EN11]{EN11:Acceleration-Randomized}
Y.~C. Eldar and D.~Needell.
\newblock Acceleration of randomized {K}aczmarz method via the
  {J}ohnson-{L}indenstrauss lemma.
\newblock {\em Numer. Algorithms}, 58(2):163--177, 2011.

\bibitem[FR12]{FR12:Mathematical-Introduction}
S.~Foucart and H.~Rauhut.
\newblock {\em A mathematical introduction to compressive sensing}.
\newblock Springer, 2012.
\newblock To appear.

\bibitem[FS95]{FS95:Kaczmarz-Based-Approach}
H.~G. Feichtinger and T.~Strohmer.
\newblock A {K}aczmarz-based approach to nonperiodic sampling on unions of
  rectangular lattices.
\newblock In {\em SampTA '95: 1995 Workshop on Sampling Theory and
  Applications}, pages 32--37, Jurmala, Latvia, Sep. 1995.

\bibitem[Gau95]{Gau95:Theoria-Combinationis}
C.~F. Gauss.
\newblock {\em Theoria Combinationis Observationum Erroribus Minimis Obnoxiae,
  Supplementum}.
\newblock Classics in Applied Mathematics. Society for Industrial and Applied
  Mathematics (SIAM), 1995.
\newblock Transl. G.W. Stewart.

\bibitem[GBH70]{GBH70:Algebraic-Reconstruction}
R.~Gordon, R.~Bender, and G.~T. Herman.
\newblock Algebraic reconstruction techniques ({ART}) for three-dimensional
  electron microscopy and {X}-ray photography.
\newblock {\em J. Theoret. Biol.}, 29:471--481, 1970.

\bibitem[GVL96]{GVL96:Matrix-Computations}
G.~H. Golub and C.~F. Van~Loan.
\newblock {\em Matrix computations}.
\newblock Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins
  University Press, Baltimore, MD, third edition, 1996.

\bibitem[HM93]{HM93:Algebraic-Reconstruction}
G.~Herman and L.~Meyer.
\newblock Algebraic reconstruction techniques can be made computationally
  efficient.
\newblock {\em IEEE Trans. Medical Imaging}, 12(3):600--609, 1993.

\bibitem[HMT11]{HMT11:Finding-Structure}
N.~Halko, P.~G. Martinsson, and J.~A. Tropp.
\newblock Finding structure with randomness: probabilistic algorithms for
  constructing approximate matrix decompositions.
\newblock {\em SIAM Rev.}, 53(2):217--288, 2011.

\bibitem[HS78]{HS78:Angles-Null}
C.~Hamaker and D.~C. Solmon.
\newblock The angles between the null spaces of {X}-rays.
\newblock {\em J. Math. Anal. Appl.}, 62(1):1--23, 1978.

\bibitem[HW71]{HW71:Bound-Tail}
D.~L. Hanson and F.~T. Wright.
\newblock A bound on tail probabilities for quadratic forms in independent
  random variables.
\newblock {\em Ann. Math. Statist.}, 42:1079--1083, 1971.

\bibitem[Kac37]{Kac37:Angenaeherte-Aufloesung}
S.~Kaczmarz.
\newblock Angen\"aherte aufl\"osung von systemen linearer gleichungen.
\newblock {\em Bull. Int. Acad. Polon. Sci. Lett. Ser. A}, pages 335--357,
  1937.

\bibitem[KS59]{KS59:Extensions-Pure}
R.~V. Kadison and I.~M. Singer.
\newblock Extensions of pure states.
\newblock {\em Amer. J. Math.}, 81:383--400, 1959.

\bibitem[KT94]{KT94:Some-Remarks}
B.~Kashin and L.~Tzafriri.
\newblock Some remarks on coordinate restriction of operators to coordinate
  subspaces.
\newblock Insitute of Mathematics Preprint~12, Hebrew University, Jerusalem,
  1993--1994.

\bibitem[Lib09]{Lib09:Accelerated-Dense}
E.~Liberty.
\newblock {\em Accelerated Dense Random Projections}.
\newblock Phd dissertation, Yale University, New Haven, CT, 2009.

\bibitem[LL10]{LL10:Randomized-Methods}
D.~Leventhal and A.~S. Lewis.
\newblock Randomized methods for linear constraints: convergence rates and
  conditioning.
\newblock {\em Math. Oper. Res.}, 35(3):641--654, 2010.

\bibitem[Nao11]{Nao11:Sparse-Quadratic}
A.~Naor.
\newblock Sparse quadratic forms and their geometric applications.
\newblock Technical Report No.~1033, S{\'e}minaire Bourbaki, Jan. 2011.

\bibitem[Nat01]{Nat01:Mathematics-Computerized}
F.~Natterer.
\newblock {\em The mathematics of computerized tomography}, volume~32 of {\em
  Classics in Applied Mathematics}.
\newblock Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
  PA, 2001.
\newblock Reprint of the 1986 original.

\bibitem[Nee10]{Nee10:Randomized-Kaczmarz}
D.~Needell.
\newblock Randomized {K}aczmarz solver for noisy linear systems.
\newblock {\em BIT}, 50(2):395--403, 2010.

\bibitem[NW12]{NW12:Two-Subspace-Projection}
D.~Needell and R.~Ward.
\newblock Two-subspace projection method for coherent overdetermined linear
  systems.
\newblock Available at \url{arXiv:1204.0277}, Apr. 2012.

\bibitem[Pop98]{Pop98:Extensions-Block-Projections}
C.~Popa.
\newblock Extensions of block-projections methods with relaxation parameters to
  inconsistent and rank-deficient least-squares problems.
\newblock {\em BIT}, 38(1):151--176, 1998.

\bibitem[Pop99]{Pop99:Block-Projections-Algorithms}
C.~Popa.
\newblock Block-projections algorithms with blocks containing mutually
  orthogonal rows and columns.
\newblock {\em BIT}, 39(2):323--338, 1999.

\bibitem[Pop01]{Pop01:Fast-Kaczmarz-Kovarik}
C.~Popa.
\newblock A fast {K}aczmarz-{K}ovarik algorithm for consistent least-squares
  problems.
\newblock {\em Korean J. Comput. Appl. Math.}, 8(1):9--26, 2001.

\bibitem[Pop04]{Pop04:Kaczmarz-Kovarik-Algorithm}
C.~Popa.
\newblock A {K}aczmarz-{K}ovarik algorithm for symmetric ill-conditioned
  matrices.
\newblock {\em An. \c Stiin\c t. Univ. Ovidius Constan\c ta Ser. Mat.},
  12(2):135--146, 2004.

\bibitem[RR12]{RR12:Beneath-Valley}
B.~Recht and C.~R{\'e}.
\newblock Beneath the valley of the noncommutative arithmetic--geometric mean
  inequality: Conjectures, case studies, and consequences.
\newblock In {\em Proc. 25th Ann. Conf. Learning Theory}, Edinburgh, June 2012.

\bibitem[RT11]{RT11:Iteration-Complexity}
P.~Richt{\'a}rik and M.~Tak{\'a}{\v c}.
\newblock Iteration complexity of randomized block-coordinate descent methods
  for minimizing a composite function.
\newblock Available at \url{arXiv:1107.2848}, Apr. 2011.

\bibitem[Sri10]{Sri10:Spectral-Sparsification}
N.~Srivastava.
\newblock {\em Spectral sparsification and restricted invertibility}.
\newblock Phd dissertation, Yale University, New Haven, CT, 2010.

\bibitem[SS87]{SS87:Incorporation-A-Priori}
M.~I. Sezan and H.~Stark.
\newblock Incorporation of a priori moment information into signal recovery and
  synthesis problems.
\newblock {\em J. Math. Anal. Appl.}, 122(1):172--186, 1987.

\bibitem[SS12]{SS12:Elementary-Proof}
D.~A. Spielman and N.~Srivastava.
\newblock An elementary proof of the restricted invertibility theorem.
\newblock {\em Israel J. Math.}, 190:83--91, 2012.

\bibitem[SV09a]{SV09:Comments-Randomized}
T.~Strohmer and R.~Vershynin.
\newblock Comments on the randomized {K}aczmarz method.
\newblock {\em J. Fourier Anal. Appl.}, 15(4):437--440, 2009.

\bibitem[SV09b]{SV09:Randomized-Kaczmarz}
T.~Strohmer and R.~Vershynin.
\newblock A randomized {K}aczmarz algorithm with exponential convergence.
\newblock {\em J. Fourier Anal. Appl.}, 15(2):262--278, 2009.

\bibitem[Tro08a]{Tro08:Norms-Random}
J.~A. Tropp.
\newblock Norms of random submatrices and sparse approximation.
\newblock {\em C. R. Math. Acad. Sci. Paris}, 346(23-24):1271--1274, 2008.

\bibitem[Tro08b]{Tro08:Random-Paving}
J.~A. Tropp.
\newblock The random paving property for uniformly bounded matrices.
\newblock {\em Studia Math.}, 185(1):67--82, 2008.

\bibitem[Tro09]{Tro09:Column-Subset}
J.~A. Tropp.
\newblock Column subset selection, matrix factorization, and eigenvalue
  optimization.
\newblock In {\em Proceedings of the {T}wentieth {A}nnual {ACM}-{SIAM}
  {S}ymposium on {D}iscrete {A}lgorithms}, pages 978--986, Philadelphia, PA,
  2009. SIAM.

\bibitem[Tro11]{Tro11:Improved-Analysis}
J.~A. Tropp.
\newblock Improved analysis of the subsampled randomized {H}adamard transform.
\newblock {\em Adv. Adapt. Data Anal.}, 3(1-2):115--126, 2011.

\bibitem[Tro12]{Tro12:User-Friendly}
J.~A. Tropp.
\newblock User-{F}riendly {T}ail {B}ounds for {S}ums of {R}andom {M}atrices.
\newblock {\em Found. Comput. Math.}, 12(4):389--434, 2012.

\bibitem[Tse93]{Tse93:Dual-Coordinate}
P.~Tseng.
\newblock Dual coordinate ascent methods for non-strictly convex minimization.
\newblock {\em Math. Programming}, 59(2, Ser. A):231--247, 1993.

\bibitem[Ver01]{Ver01:Johns-Decompositions}
R.~Vershynin.
\newblock John's decompositions: selecting a large part.
\newblock {\em Israel J. Math.}, 122:253--277, 2001.

\bibitem[Ver06]{Ver06:Random-Sets}
R.~Vershynin.
\newblock Random sets of isomorphism of linear operators on {H}ilbert space.
\newblock In {\em High dimensional probability}, volume~51 of {\em IMS Lecture
  Notes Monogr. Ser.}, pages 148--154. Inst. Math. Statist., Beachwood, OH,
  2006.

\bibitem[WLRT08]{WLRT08:Fast-Randomized}
F.~Woolfe, E.~Liberty, V.~Rokhlin, and M.~Tygert.
\newblock A fast randomized algorithm for the approximation of matrices.
\newblock {\em Appl. Comput. Harmon. Anal.}, 25(3):335--366, 2008.

\bibitem[XZ02]{XZ02:Method-Alternating}
J.~Xu and L.~Zikatanov.
\newblock The method of alternating projections and the method of subspace
  corrections in {H}ilbert space.
\newblock {\em J. Amer. Math. Soc.}, 15(3):573--597, 2002.

\bibitem[You12a]{You12:Note-Column}
P.~Youssef.
\newblock A note on column subset selection.
\newblock Available at \url{arXiv:1212.0976}, Dec. 2012.

\bibitem[You12b]{You12:Restricted-Invertibility}
P.~Youssef.
\newblock Restricted invertibility and the {B}anach--{M}azur distance to the
  cube.
\newblock Available at \url{arXiv:1206.0654}, June 2012.

\bibitem[ZF12]{ZF12:Randomized-Extended}
A.~Zouzias and N.~M. Freris.
\newblock Randomized extended {K}aczmarz for solving least-squares.
\newblock Available at \url{arXiv:1205.5770}, May 2012.

\end{thebibliography}
