Parallelising the dual revised simplex method

1 Parallelising the dual revised simplex methodJulian Hall1Qi Huangfu2 Miles Lubin31 School of Mathematics, University of Edinburgh2 FICO3 MITC onvex Optimization and BeyondEdinburgh27 June 2014 Parallelising the dual revised simplex method : OverviewBackgroundThree approachesMultiple iteration parallelism for general LPSingle iteration parallelism for general LPData parallelism for stochastic LPConclusionsJulian HallParallelising the dual revised simplex method2 / 42 Linear programming (LP)minimizecTxsubject toAx=b x 0 BackgroundFundamental model in optimaldecision-makingSolution techniques simplex method (1947) Interior point methods (1984)Large problems have 103 107/8variables 103 107/8constraintsMatrixAis (usually) sparseExampleSTAIR: 356 rows, 467 columns and 3856 nonzerosJulian HallParallelising the dual revised simplex method3 / 42 Solving LP problemsminimizefP=cTxmaximizefD=bTysubj ect toAx=b x 0(P)subject toATy+s=c s 0(D)Optimality conditionsFor a partitionB Nof the variable set with nonsingularbasis matrixBinBxB+NxN=bfor (P)and[BTNT]y+[sBsN]=[cBcN]for (D)withxN=0andsB=0 Primal basic variablesxBgiven by b=B 1bDual non-basic variablessNgiven by cTN=cTN cTBB 1 NPartition is optimal if there isPrimal feasibility b 0 dual feasibility cN 0 Julian HallParallelising the dual revised simplex method4 / 42 simplex algorithm: Each iterationRHS aq aTp cTN apq cq bp bNBDual algorithm: Assume cN 0 Seek b 0 scan bi,i B, for a good candidatepto leaveBCHUZRScan cj/ apj,j N, for a good candidateqto leaveNCHUZCU pdate.

2 ExchangepandqbetweenBandNUpdate b:= b p aq p= bp/ apqUPDATE-PRIMALU pdate cTN:= cTN d aTp d= cq/ apqUPDATE-DUALJ ulian HallParallelising the dual revised simplex method5 / 42 Standard simplex method (SSM): ComputationMajor computational componentUpdate of tableau: N:= N 1 apq aq aTpwhere N=B 1 NHopelessly inefficient for sparse LP problemsProhibitively expensive for large LP problemsJulian HallParallelising the dual revised simplex method6 / 42 revised simplex method (RSM): ComputationMajor computational components Tp=eTpB 1 BTRAN aTp= TpNPRICE aq=B 1aqFTRANI nvertBINVERTH yper-sparsityVectorsepandaqare always sparseBmay be highly reducible;B 1may be sparseVectors p, aTpand aqmay be sparseEfficient implementations must exploit these featuresH and McKinnon (1998 2005), Bixby (1999)Clp, Koberstein and Suhl (2005 2008)Julian HallParallelising the dual revised simplex method7 / 42 revised simplex method (RSM): Algorithmic enhancementsRow selection: dual steepest edge (DSE)Weight bibywi: measure of B 1ei 2 Requires additionalFTRANbut can reduce iteration count significantlyColumn selection: Bound-flipping ratio test (BFRT)Minimizes the dual objective whilst remaining dual feasibleDual variables may change sign if corresponding primal variables can flip boundsRequires additionalFTRANbut can reduce iteration count significantlyJulian HallParallelising the dual revised simplex method8 / 42 Exploiting parallelism.

3 BackgroundData parallel standard simplex methodGood parallel efficiencywasachievedOnly relevant for dense LP problemsData parallel revised simplex methodOnly immediate parallelism is in forming TpNWhenn msignificant speed-upwasachievedBixby and Martin (2000)Task parallel revised simplex methodOverlap computational components for different iterationsWunderling (1996), H and McKinnon (1995-2005)Modest speed-upwasachieved on general sparse LP problemsJulian HallParallelising the dual revised simplex method9 / 42 Parallelising the dual revised simplex method : OverviewParallelising the dual revised simplex method : OverviewSingle iteration parallelism for general LPPure dual revised simplexData parallelism:Form TpNTask parallelism:Identify serial computation which can be overlappedMultiple iteration parallelism for general LPDual revised simplex with minor iterations of dual standard simplexData parallelism:Form TpNand update (slice of) dual standard simplex tableauTask parallelism:Identify serial computation which can be overlappedData parallelism for stochastic LPPure dual revised simplex for column-linked block angular LP problemsData parallelism:SolveBT =ep,B aq=aqand form TpNJulian HallParallelising the dual revised simplex method11 / 42 Single iteration parallelismSingle iteration parallelism.

4 dual revised simplex methodComputational components appear sequentialEach has highly-tuned sparsity-exploiting serial implementationExploit slack in data dependenciesJulian HallParallelising the dual revised simplex method13 / 42 Single iteration parallelism: Computational schemeParallelPRICEto form aTp= TpNOther computational componentsserialOverlap any independent calculationsOnly four worthwhile threads unlessn msoPRICE dominatesMore than Bixby and Martin (2000)Better than Forrest (2012)Julian HallParallelising the dual revised simplex method14 / 42 Single iteration parallelism: 4-coresipvs mean speedup is is generally poor for problems with high hyper-sparsityPerformance is generally good for problems with low hyper-sparsityJulian HallParallelising the dual revised simplex method15 / 42 Multiple iteration parallelismMultiple iteration parallelismsiphas too little work to be performed in parallel to get good speedupPerform standard dual simplex minor iterations for rows in setP(|P| m)Suggested by Rosander (1975) but never implemented efficientlyin serialRHS cTN aTP b bPNBTask-parallel multipleBTRANto form P=B 1ePData-parallelPRICEto form aTp(as required)Data-parallel tableau updateTask-parallel multipleFTRANfor primal, dual and weight updatesJulian HallParallelising the dual revised simplex method17 / 42 Multiple iteration parallelism.

5 8-corepamivs for all problemsGeometric mean speedup is HallParallelising the dual revised simplex method18 / 42 Multiple iteration parallelism: 8-corepamivs mean speedup is than speedup relative to 1-corepamiGeometric mean speed of 1-corepamirelative to 1-corehsolis HallParallelising the dual revised simplex method19 / 42 Multiple iteration parallelism: Performance profile benchmarking12345678910020406080100pamic lpcplexpamiis plainly better thanclppamiis comparable withcplexpamiideas have been incorporated in FICO Xpress (Huangfu 2014)Julian HallParallelising the dual revised simplex method20 / 42 Data parallelism for stochastic LPsStochastic MIP problems: GeneralTwo-stage stochastic LPs have column-linked block angular structureminimizecT0x0+cT1x1+cT2x2+..+cT NxNsubject toAx0=b0T1x0+W1x1=b1T2x0+W2x2= +WNxN=bNx0 0x1 0x2 0 Variablesx0 Rn0arefirst stagedecisionsVariablesxi Rnifori= 1.

6 ,Naresecond stagedecisionsEach corresponds to ascenariowhich occurs with modelled probabilityThe objective is the expected cost of the decisionsIn stochastic MIP problems, some/all decisions are discreteJulian HallParallelising the dual revised simplex method22 / 42 Stochastic MIP problems: For ArgonnePower systems optimization project at ArgonneInteger second-stage decisionsStochasticity comes from availability of wind-generated electricityInitial experiments carried out using model problemNumber of scenarios increases with refinement of probability distribution samplingSolution via branch-and-boundSolve root node using parallel IPM solverPIPSL ubin, Petraet al.(2011)Solve subsequent nodes using parallel dual simplex solverPIPS-SLubin, Het al.(2013)Julian HallParallelising the dual revised simplex method23 / 42 Stochastic MIP problems: GeneralConvenient to permute the LP thus:minimizecT1x1+cT2x2+.

7 +cTNxN+cT0x0subject toW1x1+T1x0=b1W2x2+T2x0= +TNx0=bNAx0=b0x1 0x2 0x0 0 Julian HallParallelising the dual revised simplex method24 / 42 Exploiting problem structure: Basis matrix inversionInversion of the basis matrixBis key to revised simplex efficiencyFor column-linked BALP problemsB= WBiare columns corresponding tonBibasic variables in scenarioi are columns corresponding tonB0basic first stage HallParallelising the dual revised simplex method25 / 42 Exploiting problem structure: Basis matrix inversionInversion of the basis matrixBis key to revised simplex efficiencyFor column-linked BALP problemsB= Bis nonsingular soWBiare tall : full column rank[WBiTBi]are wide : full row rankABis wide : full row rankScope for parallel inversion is immediate and well knownDuff and Scott (2004).

8 Julian HallParallelising the dual revised simplex method26 / 42 Exploiting problem structure: Basis matrix inversionEliminate sub-diagonal entries in eachWBi(independently)Apply elimination operations to eachTBi(independently)Accumulate non-pivoted rows from theWBiwithABandcomplete eliminationJulian HallParallelising the dual revised simplex method27 / 42 Exploiting problem structure: Basis matrix inversionAfter Gaussian elimination, have invertible representation ofB= = SCRV SpecificallyLiUi=Siof dimensionnBi Ci=L 1iCi Ri=RiU 1iLU factors of the Schur complementM=V RS 1 Cof dimensionnB0 Scope for parallelism since each GE applied to[WBi|TBi]is independentJulian HallParallelising the dual revised simplex method28 / 42 Exploiting problem structure: SolvingBx=bFTRANforBx=bSolve[S CR V][x x0]=[b b0]as1 Liyi=bi,i= 1.

9 ,N2zi= Riyi,i= 1,..,N3z=b0 N i=1zi4Mx0=z5 Uixi=yi Cix0,i= 1,..,NAppears to be dominated by parallelizableSolvesLiyi=biandUixi=yi Cix0 Products Riyiand Cix0 Curse of exploiting hyper-sparsityIn simplex ,b is from constraint columnEither or, more likely, 0wiq0 In latter case, theyiinherit structureOnly oneLiyi=wiqOnly one RiyiLess scope for parallelism than HallParallelising the dual revised simplex method29 / 42 Exploiting problem structure: SolvingBTx=bBTRANforBTx=bSolve[STRTCTVT] [x x0]=[b b0]as1 UTiyi=bi,i= 1,..,N2zi= CTiyi,i= 1,..,N3z=b0 N i=1zi4 MTx0=z5 LTixi=yi RTix0,i= 1,..,NAppears to be dominated by parallelizableSolvesUTiyi=biandLTixi=yi RTix0 Products CTiyiand RTix0 Curse of exploiting hyper-sparsityIn simplex ,b=epAt most one solveUTiyi=biAt most one CTiyiLess scope for parallelism than HallParallelising the dual revised simplex method30 / 42 Exploiting problem structure: Forming TpNPRICE forms[ T1 TN T0] =[ T1WN1 TNWNN T0AN+N i=1 TiTNi]Dominated by parallelizable products TiWNiand TiTNiJulian HallParallelising the dual revised simplex method31 / 42 Exploiting problem structure.

10 UpdateUpdate of the invertible representation ofBis second major factor in revisedsimplex efficiencyEach iteration columnaqof the constraint matrix replaces columnBepofBB =B[I+ ( aq ep)eTp]Unfortunately, the structure ofBis not generally maintainedPIPS-Suses standardproduct formupdateB 1= [I+ ( aq ep)eTp] 1B 1=E 1B 1whereE 1=I 1 apq( aq ep)eTpJulian HallParallelising the dual revised simplex method32 / 42 Exploiting problem structure: Applying updateFormx=E 1basxp= bp apqthenxp =bp + aqxpUpdates{Ek}Kk=1ofB0toBKrequire{ aqk}Kk=1andP={pk}Kk=1,|P| mExploit parallelism when formingx=E 11bthusComputexPseriallyComputexP as a parallel matrix-vector productxP =bP +[ aqK]xPSimilar trick for parallelisingxT=bTE 11 Lubin, Het al.(2013)Julian HallParallelising the dual revised simplex method33 / 42 ResultsResults: Stochastic LP test problemsTest1st Stage2nd-Stage ScenarioNonzero ElementsProblemn0m0nimiAWiTiStorm1211851 ,2595286963,220121 SSN891706175892,28489UC123,132056,53259, 4360163,8393,132UC246,2640113,064118,872 0327,9396,264 StormandSSNare publicly availableUC12andUC24are stochastic unit commitment problems developed at ArgonneAim to choose optimal on/off schedules for generators on the power grid of the stateof Illinois over a 12-hour and 24-hour horizonIn practice each scenario corresponds to a weather simulationModel problem generates scenarios by normal perturbationsZavala (2011)Julian HallParallelising the dual revised simplex method35 / 42 Results.