1 / 32

Avoiding Communication in Sparse Matrix-Vector Multiply ( SpMV )

Avoiding Communication in Sparse Matrix-Vector Multiply ( SpMV ). Sequential and shared-memory performance is dominated by off-chip communication Distributed-memory performance is dominated by network communication The problem: SpMV has low arithmetic intensity.

wynn
Download Presentation

Avoiding Communication in Sparse Matrix-Vector Multiply ( SpMV )

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. Avoiding Communication in Sparse Matrix-Vector Multiply (SpMV) • Sequential and shared-memory performance is dominated by off-chip communication • Distributed-memory performance is dominated by network communication The problem: SpMV has low arithmetic intensity

  2. SpMV Arithmetic Intensity (1) dimension: n = 5 number of nonzeros: nnz = 3n-2 (tridiagonalA) Assumption: A is invertible ⇒ nonzero in every row ⇒ nnz ≥ n overcounts flops by up to n(diagonal A)

  3. SpMV Arithmetic Intensity (2) O( lg(n) ) O( 1 ) O( n ) more flops per byte A r i t h m e t i c I n t e n s i t y SpMV, BLAS1,2 FFTs Stencils (PDEs) Dense Linear Algebra (BLAS3) Lattice Methods Particle Methods • Arithmetic intensity := Total flops / Total DRAM bytes • Upper bound: compulsory traffic • further diminished by conflict or capacity misses

  4. 256.0 128.0 64.0 32.0 16.0 8.0 4.0 2.0 1.0 0.5 1/8 1/4 1/2 1 2 4 8 16 SpMV Arithmetic Intensity (3) • In practice, A requires at least nnzwords: • indexing data, zero padding • depends on nonzero structure, eg, banded or dense blocks • depends on data structure, eg, • CSR/C, COO, SKY, DIA, JDS, ELL, DCS/C, … • blocked generalizations • depends on optimizations, eg, index compression or variable block splitting Opteron 2356 (Barcelona) peak double precision floating-point rate attainable gflop/s stream bandwidth How to do more flops per byte? Reuse data (x, y, A) across multiple SpMVs • 2 flops per word of data • 8 bytes per double • flop:byteratio ≤ ¼ • Can’t beat 1/16 of peak!  actual flop:byteratio

  5. Combining multiple SpMVs (1) kindependent SpMVs • (1) used in: • Block Krylov methods • Krylov methods for multiple systems (AX = B) (2) k dependent SpMVs • (2) used in: • s-step Krylov methods, • Communication-avoiding Krylov methods • …to compute kKrylov basis vectors (3) k dependent SpMVs, in-place variant Def. Krylov space (given A, x, s): • (3) used in: • multigridsmoothers, power method • Related to Streaming Matrix Powers optimization for CA-Krylov methods • What if we can amortize cost of reading A over kSpMVs ? • (k-fold reuse of A)

  6. (1) k independent SpMVs(SpMM) • SpMMoptimization: • Compute row-by-row • Stream A only once =

  7. (2) kdependent SpMVs(Akx) Naïve algorithm (no reuse): • Akx (Akx) optimization: • Must satisfy data dependencies while keeping working set in cache

  8. (2) kdependent SpMVs(Akx) Akx algorithm(reuse nonzeros of A):

  9. (3) kdependent SpMVs, in-place (Akx, last-vector-only) • Last-vector-only Akx optimization: • Reuses matrix and vector k times, instead of once. • Overwrites intermediates without memory traffic • Attains O(k) reuse, even when nnz < n • eg, A is a stencil (implicit values and structure)

  10. Combining multiple SpMVs(summary of sequential results)

  11. 256.0 128.0 64.0 32.0 16.0 8.0 4.0 2.0 1.0 0.5 1/8 1/4 1/2 1 2 4 8 16 Avoiding Serial Communication • Reduce compulsory misses by reusing data: • more efficient use of memory • decreased bandwidth cost (Akx, asymptotic) • Must also consider latency cost • How many cachelines? • depends on contiguous accesses • When k = 16 ⇒ compute-bound? • Fully utilize memory system • Avoid additional memory traffic like capacityand conflict misses • Fully utilize in-core parallelism • (Note: still assumes no indexing data) • In practice, complex performance tradeoffs. • Autotune to find best k Opteron 2356 (Barcelona)  peak DP attainable gflop/s Stream Bandwidth ?  actual flop:byteratio

  12. On being memory bound • Assume that off-chip communication (cache to memory) is bottleneck, • eg, that we express sufficient ILP to hide hits in L3 • When your multicore performance is bound by memory operations, is it because of latency or bandwidth? • Latency-bound: expressed concurrency times the memory access rate does not fully utilize the memory bandwidth • Traversing a linkedlist, pointer-chasing benchmarks • Bandwidth-bound: expressed concurrency times the memory access rate exceeds the memory bandwidth • SpMV, stream benchmarks • Either way, manifests as pipeline stalls on loads/stores (suboptimal throughput) • Caches can improve memory bottlenecks – exploit them whenever possible • Avoid memory traffic when you have temporal or spatial locality • Increase memory traffic when cache line entries are unused (no locality) • Prefetcherscan allow you to express more concurrency • Hide memory traffic when your access pattern has sequential locality (clustered or regularly strided access patterns)

  13. Distributed-memory parallel SpMV • Harder to make general statements about performance: • Many ways to partition x, y, and A to P processors • Communication, computation, and load-balance are partition-dependent • What fits in cache? (What is “cache”?!) • A parallel SpMV involves 1 or 2 rounds of messages • (Sparse) collective communication, costly synchronization • Latency-bound (hard to saturate network bandwidth) • Scatter entries of x and/or gather entries of y across network • kSpMVs cost O(k) rounds of messages • Can we do kSpMVs in one round of messages? • k independent vectors? SpMM generalizes • Distribute all source vectors in one round of messages • Avoid further synchronization • k dependent vectors? Akx generalizes • Distribute source vector plus additional ghost zone entries in one round of messages • Avoid further synchronization • Last-vector-only Akx ≈ standard Akx in parallel • No savings discarding intermediates

  14. Distributed-memory parallel Akx Example: tridiagonal matrix, k = 3, n = 40, p = 4 Naïve algorithm: k messages per neighbor Akx optimization: 1 message per neighbor

  15. Polynomial Basis for Akx • Today we considered the special case of the monomials: • Stability problems - tends to lose linear independence • Converges to principal eigenvector • Given A, x, k > 0, compute where pj(A) is a degree-j polynomial in A. • Choose p for stability.

  16. Tuning space for Akx • DLP optimizations: • vectorization • ILP optimizations: • Software pipelining • Loop unrolling • Eliminate branches, inline functions • TLP optimizations: • Explicit SMT • Memory system optimizations: • NUMA-aware affinity • Software prefetching • TLB blocking • Memory traffic optimizations: • Streaming stores (cache bypass) • Array padding • Cache blocking • Index compression • Blocked sparse formats • Stanza encoding • Distributed memory optimizations: • Topology-aware sparse collectives • Hypergraph partitioning • Dynamic load balancing • Overlapped communication and computation • Algorithmic variants: • Compositions of distributed-memory parallel, shared memory parallel, sequential algorithms • Streaming or explicitly buffered workspace • Explicit or implicit cache blocks • Avoiding redundant computation/storage/traffic • Last-vector-only optimization • Remove low-rank components (blocking covers) • Different polynomial bases pj(A) • Other: • Preprocessing optimizations • Extended precision arithmetic • Scalable data structures (sparse representations) • Dynamic value and/or pattern updates

  17. Krylov subspace methods (1) • Want to solve Ax = b (still assume A is invertible) • How accurately can you hope to compute x? • Depends on condition number of A and the accuracy of your inputs A and b • condition number with respect to matrix inversion • cond(A) – how much A distorts the unit sphere (in some norm) • 1/cond(A) – how close A is to a singular matrix • expect to lose log10(cond(A)) decimal digits relative to (relative) input accuracy • Idea: Make successive approximations, terminate when accuracy is sufficient • How good is an approximation x0 to x? • Error: e0 = x0 - x • If you know e0, then compute x = x0 - e0 (and you’re done.) • Finding e0 is as hard as finding x; assume you never have e0 • Residual: r0 = b – Ax0 • r0 = 0 ⇔ e0 = 0, but they do not necessarily vanish simultaneously • cond(A) small ⇒ (r0small ⇒ e0 small)

  18. Krylov subspace methods (2) • Given approximation xold, refine by adding a correction xnew= xold + v • Pick v as the ‘best possible choice’ from search space V • Krylov subspace methods: V := • 2. Expand V by one dimension • 3. xold =xnew. Repeat. • Once dim(V) = dim(A) = n, xnew should be exact • Why Krylov subspaces? • Cheap to compute (via SpMV) • Search spaces V coincide with the residual spaces - • makes it cheaper to avoid repeating search directions • K(A,z) = K(c1A - c2I, c3z) ⇒ invariant under scaling, translation • Without loss, assume |λ(A)|≤ 1 • As s increases, Ks gets closer to the dominant eigenvectors of A • Intuitively, corrections v should target ‘largest-magnitude’ residual components

  19. Convergence of Krylov methods • Convergence = process by which residual goes to zero • If A isn’t too poorly conditioned, error should be small. • Convergence only governed by the angles θm between spaces Km and AKm • How fast does sin(θm) go to zero? • Not eigenvalues! You can construct a unitary system that results in the same sequence of residuals r0, r1, … • If A is normal, λ(A) provides bounds on convergence. • Preconditioning • Transforming A with hopes of ‘improving’ λ(A) or cond(A)

  20. Conjugate Gradient (CG) Method Given starting approximation x0 to Ax = b, let p0 := r0 := b - Ax0. For m = 0, 1, 2, … until convergence, do: • Vector iterates: • xm = candidate solution • rm = residual • pm = search direction Correct candidate solution along search direction Update residual according to new candidate solution Expand search space • Communication-bound: • 1 SpMV operation per iteration • 2 dot productsper iteration Reformulate to use Akx Do something about the dot products

  21. Applying Akx to CG (1) • Ignore x, α, and β, for now • Unroll the CG loop s times (in your head) • Observe that: • ie, two Akx calls • This means we can represent rm+j and pm+jsymbolically as linear combinations: • And perform SpMV • operations symbolically: • (same holds for Rj-1) vectors of length n vectors of length 2j+1 CG loop: For m = 0,1,…, Do

  22. Applying Akx to CG (2) 6. Now substitute coefficient vectors for vector iterates (eg, for r) SpMV performed symbolically by shifting coordinates: CG loop: For m = 0,1,…, Do

  23. Blocking CG dot products 7. Let’s also compute the 2j+1-by-2j+1 Gram matrices: Now we can perform all dot products symbolically: CG loop: For m = 0,1,…, Do

  24. CA-CG Given approximation x0 to Ax = b, let Take s steps of CG without communication For j = 0 to s - 1, Do { For m = 0, s, 2s, …, until convergence, Do { Expand Krylov basis, using SpMM and Akx optimizations: Represent the 2s+1 inner products of length n with a 2s+1-by-2s+1 Gram matrix Communication (sequential and parallel) } End For Recover vector iterates: Represent SpMV operation as a change of basis (here, a shift): Communication (sequential only) CG loop: For m = 0,1,…, Do } End For Represent vector iterates of length n with vectors of length 2s+1 and 2s+2:

  25. CA-CG complexity (1)

  26. CA-CG complexity (2)

  27. CA-CG complexity (3) • Using Gram matrix and coefficient vectors have additional costs for CA-CG: • Dense work (besides dot products/Gram matrix) does not increase with s: • CG ≈ 6sn • CA-CG ≈ 3(2s+1)(2s+n) ≈ (6s+3)n • Sequential memory traffic decreases (factor of 4.5) • CG ≈ 6sn reads, 3sn writes • CA-CG ≈ (2s+1)nreads, 3n writes b = cacheline size, p = number of processors

  28. CA-CG tuning • Performance optimizations: • 3-term recurrence CA-CG formulation: • Avoid the auxiliary vector p. • 2x decrease in Akx flops, bandwidth cost and serial latency cost (vectors only - A already optimal) • 25% decrease in Gram matrix flops, bandwidth cost, and serial latency cost • Roughly equivalent dense flops • Streaming Akx formulation: • Interleave Akx and Gram matrix construction, then interleave Akx and vector reconstruction • 2x increase in Akx flops • Factor of s decrease in Akx sequential bandwidth and latency costs (vectors only - A already optimal) • 2x increase in Akx parallel bandwidth, latency • Eliminate Gram matrix sequential bandwidth and latency costs • Eliminate dense bandwidth costs other than 3n writes • Decreases overall sequential bandwidth and latency by O(s), regardless of nnz. • Stability optimizations: • Use scaled/shifted 2-term recurrence for Akx: • Increase Akx flops by (2s+1)n • Increase dense flops by an O(s2) term • Use scaled/shifted 3-term recurrence for Akx: • Increase Akx flops by 4sn • Increase dense flops by an O(s2) term • Extended precision: • Constant factor cost increases • Restarting: • Constant factor increase in Akx cost • Preconditioning: • Structure-dependent costs • Rank-revealing, reorthogonalization, etc… • Can also interleave Akx and Gram matrix construction without the extra work - • decreases sequential bandwidth and latency by 33% rather than a factor of s

More Related