Machine learning has gained renewed interest in recent years due to advances in computer hardware (processing power and high-capacity storage) and the availability of large amounts of data which can be used to develop accurate, robust models. While hardware improvements have facilitated the development of machine learning models in a single machine, the analysis of large amounts of data still requires parallel computing to obtain shorter running times or where the dataset cannot be stored on a single machine.

In addition to hardware improvements, algorithm redesign is also an important direction to further reduce running times. On modern computer architectures, the cost of moving data (communication) from main memory to caches in a single machine is orders of magnitude more expensive than the cost of performing floating-point operations (computation). On parallel machines the cost of moving data from one processor to another over an interconnection network is the most expensive operation. The large gap between computation and communication suggests that algorithm redesign should be driven by the goal of avoiding communication and, if necessary, decreasing communication at the expense of additional computation.

Many problems in machine learning solve mathematical optimization problems which, in most non-linear and non-convex cases, requires iterative methods. This thesis is focused on deriving communication-avoiding variants of the block coordinate descent method, which is a first-order method that has strong convergence rates for many optimization problems. Block coordinate descent is an iterative algorithm which at each iteration samples a small subset of rows or columns of the input matrix, solves a subproblem using just the chosen rows or columns, and obtains a partial solution. This solution is then iteratively refined until the optimal solution is reached or until convergence criteria are met. In the parallel case, each iteration of block coordinate descent requires communication. Therefore, avoiding communication is key to attaining high performance.

This thesis adapts well-known techniques from existing work on communication-avoiding (CA) Krylov and s-step Krylov methods. CA-Krylov methods unroll vector recurrences and rearrange the sequence of computation in way that defers communication for s iterations, where s is a tunable parameter. For CA-Krylov methods the reduction in communication cost comes at the expense of numerical instability for large values of s.

We apply a similar recurrence unrolling technique to block coordinate descent in order to obtain communication-avoiding variants which solve the L2-regularized least-squares, L1-regularized least-squares, Support Vector Machines, and Kernel problems. Our communication-avoiding variants reduce the latency cost by a tunable factor of s at the expense of a factor of s increase in computational and bandwidth costs for the L2 and L1 least-squares and SVM problems. The CA-variants for these problems require additional computation and bandwidth in order to update the residual vector. For CA-kernel methods the computational and bandwidth costs do not increase. This is because the CA-variants of kernel methods can reuse elements of the kernel matrix already computed and therefore do not need to compute and communicate additional elements of the kernel matrix.

Our experimental results illustrate that our new, communication-avoiding methods can obtain speedups of up to 6.1x on a Cray XC30 supercomputer using MPI for parallel processing. For CA-kernel methods we show modeled speedups of 26x, 120x, and 197x for MPI on a predicted Exascale system, Spark on a predicted Exascale system, and Spark on a cloud system, respectively. Furthermore we also experimentally confirm that our algorithms are numerically stable for large values of s.

Finally, we also present an adaptive batch size technique which reduces the latency cost of training convolutional neural networks (CNN). With this technique we have achieved speedups of up to 6.25x when training CNNs on up to 4 NVIDIA P100 GPUs. Furthermore, we were able to train the ImageNet dataset using the ResNet-50 network with a batch size of up to 524,228 which would allow neural network training to attain a higher fraction of peak GPU performance than training with smaller batch sizes.

The cost of an algorithm includes both arithmetic and communication.

We use "communication" in a general sense to mean data movement,

either between levels of a memory hierarchy ("sequential") or

between processors ("parallel"). Communication costs include both

bandwidth terms, which are proportional to the number of words sent,

and latency terms, which are proportional to the number of messages in

which the data is sent. Communication costs are much higher than

arithmetic costs, and the gap is increasing rapidly for technological

reasons. This suggests that for best performance, algorithms should

minimize communication, even if that may require some redundant

arithmetic operations. We call such algorithms

"communication-avoiding."

Krylov subspace methods (KSMs) are iterative algorithms for solving

large, sparse linear systems and eigenvalue problems. Current KSMs

rely on sparse matrix-vector multiply (SpMV) and vector-vector

operations (like dot products and vector sums). All of these

operations are communication-bound. Furthermore, data dependencies

between them mean that only a small amount of that communication can

be hidden. Many important scientific and engineering computations

spend much of their time in Krylov methods, so the performance of many

codes could be improved by introducing KSMs that communicate less.

Our goal is to take s steps of a KSM for the same communication cost

as 1 step, which would be optimal. We call the resulting KSMs

"communication-avoiding Krylov methods." We can do this under certain

assumptions on the matrix, and for certain KSMs, both in theory (for

many KSMs) and in practice (for some KSMs). Our algorithms are based

on the so-called "s-step" Krylov methods, which break up the data

dependency between the sparse matrix-vector multiply and the dot

products in standard Krylov methods. This idea has been around for a

while, but in contrast to prior work (discussed in detail in Section

1.5), this thesis makes the following contributions:

We have fast kernels replacing SpMV, that can compute the

results of s calls to SpMV for the same communication cost as one

call (Section 2.1).

We have fast dense kernels as well, such as Tall Skinny QR (TSQR

-- Section 2.3) and Block Gram-Schmidt (BGS -- Section 2.4), which

can do the work of Modified Gram-Schmidt applied to s vectors for a

factor of Theta(s^2) fewer messages in parallel, and a factor of

Theta(s/W) fewer words transferred between levels of the memory

hierarchy (where W is the fast memory capacity in words).

We have new communication-avoiding Block Gram-Schmidt algorithms

for orthogonalization in more general inner products (Section 2.5).

We have new communication-avoiding versions of the following

Krylov subspace methods for solving linear systems: the Generalized

Minimum Residual method (GMRES -- Section 3.4), both

unpreconditioned and preconditioned, and the Method of Conjugate

Gradients (CG), both unpreconditioned (Section 5.4) and

left-preconditioned (Section 5.5).

We have new communication-avoiding versions of the following Krylov

subspace methods for solving eigenvalue problems, both standard (Ax =

&lambda x, for a nonsingular matrix A) and "generalized" (Ax =

&lambda Mx, for nonsingular matrices A and M): Arnoldi iteration

(Section 3.3), and Lanczos iteration, both for Ax = &lambda x (Section

4.2) and Ax = &lambda M x (Section 4.3).

We propose techniques for developing communication-avoiding versions

of nonsymmetric Lanczos iteration (for solving nonsymmetric eigenvalue

problems Ax = &lambda x) and the Method of Biconjugate Gradients

(BiCG) for solving linear systems. See Chapter 6 for details.

We can combine more stable numerical formulations that use different

bases of Krylov subspaces with our techniques for avoiding

communication. For a discussion of different bases, see Chapter 7.

To see an example of how the choice of basis affects the formulation

of the Krylov method, see Section 3.2.2.

We have faster numerical formulations. For example, in our

communication-avoiding version of GMRES, CA-GMRES (see Section 3.4),

we can pick the restart length r independently of the s-step basis

length s. Experiments in Section 3.5.5 show that this ability

improves numerical stability. We show in Section 3.6.3 that it also

improves performance in practice, resulting in a 2.23 × speedup

in the CA-GMRES implementation described below.

We combine all of these numerical and performance techniques in a

shared-memory parallel implementation of our communication-avoiding

version of GMRES, CA-GMRES. Compared to a similarly highly optimized

version of standard GMRES, when both are running in parallel on 8

cores of an Intel Clovertown (see Appendix A), CA-GMRES achieves 4.3

× speedups over standard GMRES on standard sparse test matrices

(described in Appendix B.5). When both are running in parallel on 8

cores of an Intel Nehalem (see Appendix A), CA-GMRES achieves 4.1

× speedups. See Section 3.6 for performance results and Section

3.5 for corresponding numerical experiments. We first reported

performance results for this implementation on the Intel Clovertown

platform in Demmel et al. [78].

We have incorporated preconditioning into our methods. Note

that we have not yet developed practical communication-avoiding

preconditioners; this is future work. We have

accomplished the following:

* We show (in Sections 2.2 and 4.3) what the s-step basis should

compute in the preconditioned case for many different types of Krylov

methods and s-step bases. We explain why this is hard in Section 4.3.

* We have identified two different structures that a preconditioner

may have, in order to achieve the desired optimal reduction of

communication by a factor of s. See Section 2.2 for details.

We present a detailed survey of related work, including s-step KSMs

(Section 1.5, especially Table 1.1) and other techniques for reducing

the amount of communication in iterative methods (Section 1.6).

Dense linear algebra computations are essential to nearly every problem in scientific computing and to countless other fields. Most matrix computations enjoy a high computational intensity (i.e., ratio of computation to data), and therefore the algorithms for the computations have a potential for high efficiency. However, performance for many linear algebra algorithms is limited by the cost of moving data between processors on a parallel computer or throughout the memory hierarchy of a single processor, which we will refer to generally as communication. Technological trends indicate that algorithmic performance will become even more limited by communication in the future. In this thesis, we consider the fundamental computations within dense linear algebra and address the following question: can we significantly improve the current algorithms for these computations, in terms of the communication they require and their performance in practice?

To answer the question, we analyze algorithms on sequential and parallel architectural models that are simple enough to determine coarse communication costs but accurate enough to predict performance of implementations on real hardware. For most of the computations, we prove lower bounds on the communication that any algorithm must perform. If an algorithm exists with communication costs that match the lower bounds (at least in an asymptotic sense), we call the algorithm communication optimal. In many cases, the most commonly used algorithms are not communication optimal, and we can develop new algorithms that require less data movement and attain the communication lower bounds.

In this thesis, we develop both new communication lower bounds and new algorithms, tightening (and in many cases closing) the gap between best known lower bound and best known algorithm (or upper bound). We consider both sequential and parallel algorithms, and we asses both classical and fast algorithms (e.g., Strassen's matrix multiplication algorithm).

In particular, the central contributions of this thesis are: proving new communication lower bounds for nearly all classical direct linear algebra computations (dense or sparse), including factorizations for solving linear systems, least squares problems, and eigenvalue and singular value problems; proving new communication lower bounds for Strassen's and other fast matrix multiplication algorithms; proving new parallel communication lower bounds for classical and fast computations that set limits on an algorithm's ability to perfectly strong scale; summarizing the state-of-the-art in communication efficiency for both sequential and parallel algorithms for the computations to which the lower bounds apply; developing a new communication-optimal algorithm for computing a symmetric-indefinite factorization (observing speedups of up to 2.8x compared to alternative shared-memory parallel algorithms); developing new, more communication-efficient algorithms for reducing a symmetric band matrix to tridiagonal form via orthogonal similar transformations (observing speedups of 2--6x compared to alternative sequential and parallel algorithms); and developing a new communication-optimal parallelization of Strassen's matrix multiplication algorithm (observing speedups of up to 2.84x compared to alternative distributed-memory parallel algorithms).

Modern commodity processors such as GPUs may execute up to about a thousand of physical threads per chip to better utilize their numerous execution units and hide execution latencies. Understanding this novel capability, however, is hindered by the overall complexity of the hardware and complexity of typical workloads. In this dissertation, we suggest a better way to understand modern multithreaded performance by considering a family of synthetic workloads, which use the same key hardware capabilities – memory access, arithmetic operations, and multithreading – but are otherwise as simple as possible.

One of our surprising findings is that prior performance models for GPUs fail on these workloads: they mispredict observed throughputs by factors of up to 1.7. We analyze these prior approaches, identify a number of common pitfalls, and discuss the related subtleties in understanding concurrency and Little’s Law. Also, we help to further our understanding by considering a few basic questions, such as on how different latencies compare with each other in terms of latency hiding, and how the number of threads needed to hide latency depends on basic parameters of executed code such as arithmetic intensity. Finally, we outline a performance modeling framework that is free from the found limitations.

As a tangential development, we present a number of novel experimental studies, such as on how mean memory latency depends on memory throughput, how latencies of individual memory accesses are distributed around the mean, and how occupancy varies during execution.

Obtaining good performance for tensor problems requires computing an architecture-specific mapping from the algorithm to the target hardware. Choosing a mapping requires optimizing an objective function over a large, nonconvex space; this objective function represents a performance metric which may be modeled, simulated, or (if possible) measured. Each such objective function incurs different tradeoffs in terms of speed, accuracy, and the strength of results that can be formally proven, and as a result requires its own optimization methods. In this dissertation, we describe optimization approaches for several performance objectives.

• In a simple abstract memory hierarchy model, we derive an unconditional communication lower bound for "projective" tensor operations (which includes most dense linear algebra) and convolutions. We show that this can always be attained (up to a constant factor) by solving a mathematical optimization problem, and that these methods provide significant benefits in practice, halving the communication volume in one case.

• We extend these lower bound techniques - and corresponding algorithms - to handle combinations of communication and convolution for randomized matrix multiplication.

• We then describe how to incorporate support for additional mapping decisions and to account for more architectural information by extending our optimization-based approach to support analytical performance models. While these performance models are too complex to prove strong lower bounds for, empirical results show that optimization-based methods can find more performant mappings using significantly less runtime and sample complexity than pure search-based methods.

As computing devices evolve with successive technology generations, many machines target either the mobile or high-performance computing/datacenter environments. In both of these form factors, energy consumption often represents the limiting factor on hardware and software efficiency. On mobile devices, limitations in battery technology may reduce possible hardware capability due to a tight energy budget. On the other hand, large machines such as datacenters and supercomputers have budgets directly related to energy consumption and small improvements in energy efficiency can significantly reduce operating costs. Such challenges have influenced research upon the impact of applications, operating and runtime systems upon energy consumption. Until recently, little consideration was given to the potential energy efficiency of algorithms themselves.

A dominant idea within the high-performance computing (HPC) community is that applications can be decomposed into a set of key computational problems, called kernels. Via automatic performance tuning and new algorithms for many kernels, researchers have successfully demonstrated performance improvements on a wide variety of machines. Motivated by the large and increasingly growing dominant cost (in time and energy) of moving data, algorithmic improvements have been attained by proving lower bounds on the data movement required to solve a computational problem, and then developing communication-optimal algorithms that attain these bounds.

This thesis extends previous research on communication bounds and computational kernels by presenting bounds on the energy consumption of a large class of algorithms. These bounds apply to sequential, distributed parallel and heterogeneous machine models and we detail methods to further extend these models to larger classes of machines. We argue that the energy consumption of computational kernels is usually predictable and can be modeled via linear models with a handful of terms. Thus, these energy models (and the accompanying bounds) may apply to many HPC applications when used in composition.

Given energy bounds, we analyze the implications of such results under additional constraints, such as an upper bound on runtime, and also suggest directions for future research that may aid future development of a hardware/software co-tuning process. Further, we present a new model of energy efficiency, Cityscape, that allows hardware designers to quickly target areas for improvement in hardware attributes. We believe that combining our bounds with other models of energy consumption may provide a useful method for such co-tuning; i.e. to enable algorithm and hardware architects to develop provably energy-optimal algorithms on customized hardware platforms.

In this paper, we present the main algorithmic features in the software package SuperLU_DIST, a distributed-memory sparse direct solver for large sets of linear equations. We give in detail our parallelization strategies, with focus on scalability issues, and demonstrate the parallel performance and scalability on current machines. The solver is based on sparse Gaussian elimination, with an innovative static pivoting strategy proposed earlier by the authors. The main advantage of static pivoting over classical partial pivoting is that it permits a priori determination of data structures and communication pattern for sparse Gaussian elimination, which makes it more scalable on distributed memory machines. Based on this a priori knowledge, we designed highly parallel and scalable algorithms for both LU decomposition and triangular solve and we show that they are suitable for large-scale distributed memory machines.

Advancements in the field of high-performance scientific computing are necessary to address the most important challenges we face in the 21st century. From physical modeling to large-scale data analysis, engineering efficient code at the extreme scale requires a critical focus on reducing communication -- the movement of data between levels of memory hierarchy or between processors over a network -- which is the most expensive operation in terms of both time and energy at all scales of computing. Achieving scalable performance thus requires a dramatic shift in the field of algorithm design, with a key area of innovation being the development of communication-avoiding algorithms.

Solvers for sparse linear algebra problems, ubiquitous throughout scientific and mathematical applications, often limit application performance due to a low computation/communication ratio. Among iterative methods, Krylov subspace methods are the most general and widely-used. To alleviate performance bottlenecks, much prior work has focused on the development of communication-avoiding Krylov subspace methods, which can offer asymptotic performance improvements over a set number of iterations.

In finite precision, the convergence and stability properties of classical Krylov methods are not necessarily maintained by communication-avoiding Krylov methods. Depending on the parameters used and the numerical properties of the problem, these communication-avoiding variants can exhibit slower convergence and decreased accuracy compared to their classical counterparts, making it unclear when communication-avoiding Krylov subspace methods are suitable for use in practice.

Until now, the literature on communication-avoiding Krylov methods lacked a detailed numerical stability analysis, as well as both theoretical and practical comparisons with the stability and convergence properties of standard implementations. In this thesis, we address this major challenge to the practical use of communication-avoiding Krylov subspace methods. We extend a number of theoretical results and algorithmic techniques developed for classical Krylov subspace methods to communication-avoiding Krylov subspace methods and identify constraints under which these methods are competitive in terms of both achieving asymptotic speedups and meeting application-specific numerical requirements.

Solving square linear systems of equations Ax=b is one of the primary workhorses in scientific computing. With asymptotically and practically small amounts of extra calculation and higher precision, we can render solution techniques dependable. We produce a solution with tiny error for almost all systems where we should expect a tiny error, and we correctly flag potential failures.

Our method uses a proven technique: iterative refinement.

We extend prior work by applying extra precision not only in

calculating the residual b-A y_i of an intermediate solution

y_i but also in carrying that intermediate solution

y_i. Analysis shows that extra precision in the

intermediate solutions lowers the limiting backward error

(measuring perturbations in the initial problem) to levels that

produce a forward error (measuring

perturbations in the solution) not much larger than the

precision used to store the result. We also demonstrate that

condition estimation is not necessary for determining success,

reducing the computation in refinement substantially.

This basic, dependable solver applies to typical dense LU

factorization methods using partial pivoting as well as methods

that risk greater failure by choosing pivots for non-numerical

reasons. Sparse factorization methods may choose pivots to

promote structural sparsity or even choose pivots before

factorization to decouple the phases. We show through

experiments that solutions using these restrictive pivoting

methods still have small error so long as an estimate of

factorization quality, the growth factor, does not grow too

large. Our refinement algorithm dependably flags such

failures. Additionally, we find a better choice of heuristic

for sparse static pivoting than the defaults in Li and Demmel's

SuperLU package.

Static pivoting in a distributed-memory setting needs an

algorithm for choosing pivots that does not rely on fitting the

entire matrix into one memory space. We investigate a set of

algorithms, Bertsekas's auction algorithms, for choosing a

static pivoting via maximum weight perfect bipartite matching.

Auction algorithms have a natural mapping to distributed memory

computation through their bidding mechanism.

We provide an analysis of the auction algorithm fitting it

comfortably in linear optimization theory and characterizing

approximately maximum weight perfect bipartite matches.

These approximately maximum weight perfect matches work well as

static pivot choices and can be computed much more quickly than

the exact maximum weight matching.

Finally, we consider the performance of auction algorithm implementations on a suite of real-world sparse problems.

Sequential performance is roughly equivalent to existing

implementations like Duff and Koster's MC64, but varies widely

with different parameter and input settings. The parallel

performance is even more wildly unpredictable. Computing

approximately maximum weight matchings helps performance

somewhat, but we still conclude that the performance is too

variable for a black-box solution method.