  We give a stochastic optimization algorithm that solves a dense $n\times n$
real-valued linear system $Ax=b$, returning $\tilde x$ such that $\|A\tilde
x-b\|\leq \epsilon\|b\|$ in time: $$\tilde
O((n^2+nk^{\omega-1})\log1/\epsilon),$$ where $k$ is the number of singular
values of $A$ larger than $O(1)$ times its smallest positive singular value,
$\omega < 2.372$ is the matrix multiplication exponent, and $\tilde O$ hides a
poly-logarithmic in $n$ factor. When $k=O(n^{1-\theta})$ (namely, $A$ has a
flat-tailed spectrum, e.g., due to noisy data or regularization), this improves
on both the cost of solving the system directly, as well as on the cost of
preconditioning an iterative method such as conjugate gradient. In particular,
our algorithm has an $\tilde O(n^2)$ runtime when $k=O(n^{0.729})$. We further
adapt this result to sparse positive semidefinite matrices and least squares
regression.
  Our main algorithm can be viewed as a randomized block coordinate descent
method, where the key challenge is simultaneously ensuring good convergence and
fast per-iteration time. In our analysis, we use theory of majorization for
elementary symmetric polynomials to establish a sharp convergence guarantee
when coordinate blocks are sampled using a determinantal point process. We then
use a Markov chain coupling argument to show that similar convergence can be
attained with a cheaper sampling scheme, and accelerate the block coordinate
descent update via matrix sketching.
