Reproducible Tall-Skinny QR J. Demmel and H.D. Nguyen - PowerPoint PPT Presentation

Reproducible ¡Tall-­‑Skinny ¡QR ¡ J. ¡Demmel ¡and ¡H.D. ¡Nguyen ¡ ARITH22 ¡ 1 ¡

Plan ¡ 1. MoBvaBon ¡ 2. Algorithms ¡ 3. Experimental ¡results ¡ 4. Conclusions ¡ 2 ¡

MoBvaBon: ¡Reproducibility ¡ • FloaBng-­‑point ¡computaBons ¡are ¡usually ¡“inexact” ¡and ¡non-­‑ determinisBc ¡due ¡to ¡the ¡non-­‑associaBvity ¡of ¡floaBng-­‑point ¡addiBon. ¡ Reproducibility: ¡obtains ¡ bit-­‑wise ¡idenBcal ¡results ¡from ¡different ¡ runs ¡of ¡the ¡program. ¡ • Sources ¡of ¡non-­‑reproducibility: ¡ – Data ¡alignment ¡ – Data ¡parBBoning/ordering ¡ – FMA ¡support ¡ – Intermediate ¡precision ¡(64 ¡bits, ¡80 ¡bits, ¡128 ¡bits, ¡etc) ¡ – Data ¡path ¡(SSE, ¡AVX, ¡GPU ¡warp, ¡etc) ¡ – Number ¡of ¡processors ¡ – Network ¡topology ¡ – ??? ¡ ¡ 3 ¡

Reproducibility ¡at ¡Large ¡Scale ¡ • Performance ¡is ¡increased ¡by ¡increasing ¡the ¡number ¡of ¡ processors ¡ – Highly ¡dynamic ¡scheduling ¡ – Network ¡heterogeneity: ¡reducBon ¡tree ¡shape ¡can ¡vary ¡ – CommunicaBon ¡cost ¡tends ¡to ¡dominate ¡arithmeBc ¡cost. ¡A ¡liale ¡ extra ¡arithmeBc ¡cost ¡is ¡allowed ¡so ¡long ¡as ¡the ¡communicaBon ¡ cost ¡is ¡controlled. ¡ • DOE ¡ASCAC ¡Report ¡Feb ¡10, ¡2014: ¡Reproducibility ¡is ¡listed ¡in ¡ top ¡10 ¡ExaScale ¡Research ¡Challenges. ¡“Reproducibility ¡will ¡ be ¡expensive ¡if ¡not ¡impossible ¡to ¡achieve ¡on ¡exascale ¡ machines. ¡… ¡requirement ¡will ¡need ¡to ¡be ¡relaxed” ¡ 4 ¡

ReproBLAS ¡ • FloaBng-­‑point ¡summaBon ¡can ¡be ¡made ¡reproducible ¡using ¡ either ¡the ¡long ¡accumulator ¡technique ¡or ¡our ¡“Indexed ¡ FloaBng-­‑Point ¡format” ¡ ¡ • ReproBLAS ¡ hap://bebop.eecs.berkeley.edu/reproblas ¡ – Reproducible ¡regardless ¡of ¡#processors, ¡reducBon ¡tree ¡shape, ¡… ¡ – BLAS ¡level ¡1: ¡ ¡{s|d|c|z}sum, ¡{s|d|c|z}asum, ¡{s|d|c|z}dot{c|u}, ¡ {s|d|c|z}norm2 ¡ – gemv, ¡gemm: ¡under ¡development ¡ – Some ¡might ¡be ¡reproducible: ¡trsv, ¡trsm, ¡Cholesky, ¡… ¡ • Not ¡all ¡linear ¡algebra ¡rouBnes ¡are ¡“automaBcally” ¡ reproducible ¡by ¡using ¡a ¡reproducible ¡BLAS ¡library ¡ ¡ 5 ¡

Time (normalized by dasum time) 0 1 2 3 4 5 Performance ¡results ¡on ¡1024 ¡proc ¡Cray ¡XC30 ¡ dasum 1 3.1 Algo 3 4.7 Algo 2 4.9 Algo 6 1.2x ¡to ¡3.2x ¡slowdown ¡vs ¡fastest ¡code, ¡for ¡n=1M ¡ dasum 2 3.2 Algo 3 4.7 Algo 2 5.0 Algo 6 dasum 4 3.1 Algo 3 4.4 (IEEE ¡Trans. ¡Comp., ¡ ¡July ¡2015) ¡ ¡ ¡ ¡ ¡ Algo 2 4.6 Algo 6 dasum 8 2.9 Algo 3 4.1 Algo 2 4.3 Algo 6 dasum 16 2.8 Algo 3 3.9 Algo 2 3.9 Algo 6 # Processors dasum 32 3.1 Algo 3 3.7 Algo 2 3.4 Algo 6 dasum 64 3.0 Algo 3 3.6 Algo 2 3.0 Algo 6 dasum 128 2.9 Algo 3 3.1 Algo 2 2.2 Algo 6 dasum 256 2.3 Algo 3 2.3 Algo 2 1.4 Algo 6 local sum Absolute Max communication All-Reduce dasum 512 2.4 Algo 3 2.3 Algo 2 1.2 Algo 6 1024 dasum 2.2 Algo 3 2.2 Algo 2 1.2 Algo 6 6 ¡

MoBvaBon: ¡Reproducible ¡TSQR ¡ • Householder ¡transformaBon ¡can ¡be ¡made ¡reproducible ¡ using ¡the ¡norm2 ¡funcBon ¡from ¡ReproBLAS, ¡however ¡it ¡ is ¡not ¡efficient ¡in ¡terms ¡of ¡communicaBon. ¡ • Tall-­‑Skinny ¡QR ¡(TSQR) ¡is ¡opBmal ¡in ¡terms ¡of ¡ communicaBon ¡but ¡computed ¡results ¡strongly ¡depend ¡ on ¡the ¡reducBon ¡tree ¡shape ¡ – Ex: ¡13x ¡speedup ¡on ¡GPU, ¡… ¡ • CholQR ¡is ¡communicaBon ¡opBmal ¡and ¡can ¡be ¡made ¡ reproducible ¡but ¡not ¡numerically ¡stable. ¡ • Main ¡goal : ¡Implement ¡a ¡reproducible ¡and ¡ ¡stable ¡ ¡ ¡ ¡ Tall-­‑Skinny ¡QR ¡factorizaBon ¡ 7 ¡

cholQR ¡ function [Q,R] = cholQR(A) � � % Require: A is nxb matrix. cond(A) < ε -1/2 � � Z = A T A � � R = chol(Z) � � Q = A / R � end � ProperBes: ¡ • Can ¡be ¡reproducible ¡using ¡a ¡reproducible ¡BLAS ¡library ¡ • Efficient ¡in ¡both ¡arithmeBc ¡and ¡communicaBon ¡cost ¡ • cholQR ¡can ¡have ¡a ¡small ¡residual ¡but: ¡ – Orthogonality ¡of ¡Q ¡is ¡not ¡guaranteed ¡ – Might ¡not ¡run ¡to ¡compleBon ¡when ¡ cond(A)> ε -1/2 ¡ 8 ¡

cholQR ¡ function [Q,R] = cholQR(A) � � % Require: A is nxb matrix. cond(A) < ε -1/2 � � Z = A T A � � R = chol(Z) � � Q = A / R � end � ProperBes: ¡ • Can ¡be ¡reproducible ¡using ¡a ¡reproducible ¡BLAS ¡library ¡ • Efficient ¡in ¡both ¡arithmeBc ¡and ¡communicaBon ¡cost ¡ • cholQR ¡can ¡have ¡a ¡small ¡residual ¡but: ¡ Reconstruct ¡Householder ¡vectors ¡ – Orthogonality ¡of ¡Q ¡is ¡not ¡guaranteed ¡ – Might ¡not ¡run ¡to ¡compleBon ¡if ¡ cond(A)> ε -1/2 ¡ 9 ¡ IteraBve ¡refinement ¡

Reconstruct ¡Householder ¡vectors ¡ YTY ¡format: ¡ Q = (I – YTY T ) � R ¡ T ¡ Y T ¡ A ¡ = ¡ × ¡ I- � Y ¡ ¡ A = QR = (I – YTY T ) R � 10 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ YTY ¡format: ¡ Q = (I – YTY T ) � R ¡ T ¡ R ¡ Y T ¡ A ¡ = ¡ - � Y ¡ ¡ A = QR = R – YTY T R � 11 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ YTY ¡format: ¡ Q = (I – YTY T ) � -­‑TY 1 T R ¡ R ¡ A ¡ = ¡ - � Y ¡ ¡ T R) � A-R = Y(-TY 1 12 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ YTY ¡format: ¡ Q = (I – YTY T ) � -­‑TY 1 T R ¡ R ¡ R ¡ A ¡ A ¡ = ¡ = ¡ LU ¡ - � Y ¡ - � ¡ T R) = LU(A-R) � A-R = Y(-TY 1 13 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ YTY ¡format: ¡ Q = (I – YTY T ) � -­‑TY 1 T R ¡ R ¡ R ¡ A ¡ A ¡ = ¡ = ¡ modLU ¡ - � Y ¡ - � ¡ Flipping ¡signs ¡to ¡ avoid ¡cancellaBon ¡ T R,S] = modLU(A,R) = LU(A,S*R) � [Y,-TY 1 14 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ function Y = r2y(A,R) ; % A = [A1;A2] � � (Y1,W) = modLU(A1,R) % on root only � � Y2 = A2/W % Broadcast + local trsm � end � function Y = repQR(A) � � R = chol(A T A) %require reduction operation � � Y = r2y(A,R) � end � • Q ¡computed ¡from ¡Y: ¡Q ¡= ¡(I ¡-­‑ ¡YTY T ) ¡is ¡always ¡orthogonal ¡ • Cholesky ¡might ¡not ¡run ¡to ¡compleBon ¡if ¡cond(A) ¡> ¡ε -­‑1/2 ¡ • LU ¡is ¡also ¡not ¡stable ¡when ¡A-­‑R ¡is ¡ill ¡condiBoned ¡ • SoluBon: ¡Using ¡iteraBve ¡refinement ¡to ¡improve ¡stability ¡ 15 ¡

Reconstruct ¡ ¡Householder ¡vectors ¡ function Y = r2y(A,R) ; % A = [A1;A2] � � (Y1,W) = modLU(A1,R) % on root only � � Y2 = A2/W % Broadcast + local trsm � end � function Y = repQR(A) � � R = chol(A T A) %require reduction operation � � Y = r2y(A,R) � end � • Q ¡computed ¡from ¡Y: ¡Q ¡= ¡(I ¡-­‑ ¡YTY T ) ¡is ¡always ¡orthogonal ¡ • Cholesky ¡might ¡not ¡run ¡to ¡compleBon ¡if ¡cond(A) ¡> ¡ε -­‑1/2 ¡ • LU ¡is ¡also ¡not ¡stable ¡when ¡A-­‑R ¡is ¡ill ¡condiBoned ¡ • SoluBon: ¡Using ¡iteraBve ¡refinement ¡to ¡improve ¡stability ¡ 16 ¡

Recursive ¡Cholesky ¡QR ¡ • cholQR ¡might ¡not ¡run ¡to ¡compleBon ¡if ¡A ¡is ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ill-­‑condiBoned: ¡cond(A)>ε -­‑1/2 , ¡or ¡cond(A T A)>ε -­‑1 . ¡ • [R,p] ¡= ¡chol(A T A): ¡ ¡ – p=0: ¡Cholesky ¡runs ¡to ¡compleBon ¡ – 0<p<b: ¡R ¡is ¡just ¡the ¡Cholesky ¡factor ¡of ¡the ¡ p × p ¡top-­‑lex ¡ submatrix ¡of ¡A T A � R ¡is ¡the ¡R ¡factor ¡of ¡A1=A(:,1:p). ¡Similarly ¡to ¡Householder ¡ factorizaBon ¡we ¡can ¡use ¡Y1=r2y(A(:,1:p), ¡R) ¡to ¡update ¡ A2=A(:,p+1,:) ¡and ¡conBnue ¡the ¡factorizaBon ¡of ¡trailing ¡ matrix. ¡ 17 ¡