# Presentation: 149793

## Mix and Match----A Multimethod Approach to Solving Linear Systems

- Sanjukta Bhowmick                                    Department of Applied Physics and Applied Mathematics
- (Columbia University)

## In Collaboration With ...

- Adaptive and Composite Solvers
- Lois Curfman McInnes (*Argonne National Laboratory*)
- Boyana Norris (*Argonne National Laboratory*)
- Padma Raghavan (*Pennsylvania State University*)
- Machine Learning Approach
- Victor Eijkhout  (*Texas Austin Computing Centre*)
- Yoav Freund (*University of California San Diego*)
- Erika Fuentes (*University of Tennessee*)
- David Keyes (*Columbia University*)

## Motivation: Extensive Use of Linear Solvers

- The execution time in many PDE-based simulations is dominated by the time to solve linear systems ***Ax=b***
- MHD codes *M3D*, *Nimrod*, *AORSA* in U.S. DOE
- Aerodynamics code *TRANAIR* at Boeing
- Different simulations/applications may have different solution requirements:
* **robustness, accuracy, fast execution time , low memory requirement , parallel scaling,.....*
- The linear solution time (in turn total simulation time) can be reduced if  solvers are chosen to match problem attributes
- Automating solver selection and composition is one of the challenges in scientific computing

## Motivation:  Difficulty of Method Selection

- There exists a vast range of choices for tuning linear solvers
- Solvers:  direct* (**order**), *iterative* (**tolerance, restart  value**), *multilevel/multigrid* (**number of levels, cycles**)*
- Preconditioners : method type, fill parameters, overlap parameters, number of levels, coarsening type, etc.
- Linear solvers typically reflect tradeoffs between metrics,
*robustness, execution time, memory requirement*
*speedup, scalability, efficiency*
- Given methods {A, B, C} it is not hard to find problems arising naturally from PDE discretizations {x, y, z}, such that :
- On x, A > B > C
- On y, B > C > A
- On z, C > A > B

- “>” is an order relation for better performance

- “*The impossibility of  finding the best linear solver for a given problem...is widely appreciated.*”  Towards Polyalgorithmic Linear System Solvers for Nonlinear Elliptic Problems (SIAM Journal of Scientific Computing)

## Opportunity for solver tuning

- Example of performance variations of differently tuned Krylov solvers over a sequence of about 85 systems
- Vertical axis is log of runtime improvements
- Opportunity is potentially orders of magnitude

##  Using Multiple Solvers  

- Multimethod approach uses more than one linear solver to solve a sequence or set of linear systems
- Multimethod Techniques
- Combination: combining  a sequence of multiple solvers   *(**increased robustness**)*_x0013_
- Selection: selecting set of “good” solvers (*faster solution time*)
- Static: previous knowledge  of problem parameters
- Dynamic: observing the problem parameters as they change
- Related Research
- Ern, Giovangili, Keyes and Smooke ----motivates “polyalgorithmic” solvers  through empirical study on nonlinear elliptic PDEs
- Dongarra, Eijkhout and Fuentes ---- Self-Adaptive* *Large-scale Solver Architecture (SALSA), statistical schemes for method classification using performance database
- Barrett, Berry, Dongarra, Eijkhout and Romine ---- multiprocessor implementation, different Krylov methods are applied in parallel to the same system

## Multi-method Approach Combination

- Composite solver - a specifically ordered sequence of distinct methods to improve reliability
- A method fails if it does not converge within “maximum” iterations
- Upon  failure of one method, the next method in the sequence is executed on the *same* system
- The composite fails only if all methods fail
- In the worst case, all methods have to be executed

## Combinatorial Model-I

- There are *n* methods *M1*, *M2*,....,*Mn.*
- Method *Mi* is associated with metrics:
- normalized execution time *ti  *or *t(i).*
- probability of not converging (failure rate) *fi*  or *f(i);* (reliability: *ri=1-fi* )
- Cumulative failure rate of methods *M1*, *M2*,....,*Mn, *is *f(1**** 2 **** ... **** n) *The cumulative failure rate is lower than the individual failure rates of the methods
- The set ***P*** contains all permutations of *{1,2,..,n}*
*Mi* denotes the i   method in *P******P***
- The composite *C* executes methods in the order specified by *P*

*^*

- th

*^*

*^*

*^*

- _Goal:  Find ordering for the composite with the __least time__ in the __worst case__ (__on average__)_

##  Combinatorial Model-II

-  Reliability is:

* **The reliability is independent of the ordering *
*The reliability is higher than any of the base methods*

-  Execution time (worst case) is :

* **The execution time depends on the ordering*

-  Goal: To determine an optimal permutation

-  Assume that failures are mutually independent

## Composites of 3 Methods

- Method  M2; Time =1.5; Reliability =0.7; Failure=0.3

- Method  M1; Time =1.0; Reliability =0.1; Failure=0.9

- Method  M3; Time =3.0; Reliability =0.8; Failure=0.2

- M3

- M2

- M1

- M3

- M2

- M1

- .9x 1.5

- .9x.3x3.0

- 1

- 3.16

- M2

- M3

- M1

- .2x 1.5

- .2x.3x1.0

- 3

- 3.36

- .3x.2x1.0

- .3x3

- 1.5

- 2.46

*Least Time ***(1,2,3)**

*Least Failure ***(3,2,1)**

*Other    ***(2,3,1)**

- The failure rate of any composite is 0.054

## The Optimal Composite is ...

- Define *u**i**=t**i**/r**i* as the *utility ratio* of *M**i*
- Let *P*** **P** and let *C* be the associated composite
- Theorem: *C* is the optimal composite with
*          **T=min{T:P**** ****P****}* _if and only if_ the sequence  of
-           methods *M**1**,M**2**,...,M**n* is such that

*^*

*~*

*~*

*u**1 **** **u**2** **** ..... **** **u**n-1** **** **u**n*

-  *The optimal composite is one in which the component methods are arranged in increasing order of their utility ratios*

*~*

*~*

*~*

*~*

*~*

*~*

*~*

*~*

*~*

*~*

## Composites as Paths in a Layered Graph

- Method  M2; Time =1.5; Failure=0.3

- Method  M1; Time =1.0; Failure=0.9

- Method  M3; Time =3.0; Failure=0.2

- 1.5

- .3x.2x1.0

- .3x3

- .9x.3x3.0

- .2x.9x1.5

- .2x1

- .9x3

- .2x1.5

- .3x1

- 1

- 3

- 0

- 1,2,3

- 1,3

- 2,3

- 1,2

- 3

- 2

- 1

*The vertices on the shortest path give the optimal composite*

- .9x1.5

## Composites of 3 Methods

- Method  M2; Time =1.5; Reliability =0.7; Failure=0.3

- Method  M1; Time =1.0; Reliability =0.1; Failure=0.9

- Method  M3; Time =3.0; Reliability =0.8; Failure=0.2

- M3

- M2

- M1

- M3

- M2

- M1

- .9x 1.5

- .9x.3x3.0

- 1

- 3.16

- M2

- M3

- M1

- .2x 1.5

- .2x.3x1.0

- 3

- 3.36

- .3x.2x1.0

- .3x3

- 1.5

- 2.46

*Least Time ***(1,2,3)**

*Least Failure ***(3,2,1)**

*Least      Utility Ratio    ***(2,3,1)**

- The failure rate of any composite is 0.054

## Composite Solvers: Experiments

- Driven Cavity with Pseudo-transient Continuation                                              (Petsc SNES example ex27)
- Velocity - vorticity formulation
- Flow driven by lid and/or buoyancy
- Logically regular grid, parallelized with DAs
- Finite difference discretization
- Grashof number and lid velocity determine the degree of nonlinearity
- Problem Parameters
- Mesh Size 128 by 128
- Grashof Number [580,620,660,700,740,780]
- Lid Velocity: [10,13,16,20]
- Total of 24 problem instances

## Experiments: Base and Composite Methods 

- ILU(l)=Incomplete LU with l level of fill
- ILUT=Incomplete LU with drop threshold .01

## Failures and Reliability

- Composites show 100% reliability

## Execution Times for Base and Composite Methods

- CU requires only *56%(58%)* of the *nonlinear(linear)* solution time required by the base method

## Multimethod Approach: Selection

- The graph plots time spent in solving linear systems during each nonlinear iteration using a fixed solver
- Different amounts of work required for each linear system

- _Goal: To __reduce the total execution__ time by selecting “suitable” linear solvers_

## Dynamic Selection:Adaptive Solver

- Adaptive Solvers: sparse linear solver is dynamically changed to match the linear system characteristics
- General method,
- design switching criterion
- select the first method
- select window size (fixed number of consecutive nonlinear iterations)
- switch at each window according to criterion
- Two classes of adaptive heuristics:
- Sequence based
- methods arranged in increasing order of t some heuristic (time per iteration *ti)*
* *assumption:  weaker methods appear earlier in the series
- selection among current method and methods preceding or succeeding it
- suitable for applications where linear systems get progressively more easier (or difficult) to solve
- Non-sequence based
- all methods considered during selection
- suitable for  applications where changes in linear systems are more complex

##  Adaptive Solver I(Sequence Based Selection)

- Switching is based on,*                                           *
- convergence rate
- rate of increase in the number of  linear iterations*                        *
- User defined boundary [-λ λ]  for the convergence rate
- The switching criteria:
- convergence rates are within the boundary: do not switch
- convergence rates exceed the upper bound: switch to preceding solver
-  convergence rates are lower than the lower bound or rate of increase of linear iterations is greater than a threshold β: switch to succeeding solver

## Adaptive Solver II (Sequence Based Selection)

- Switching is based on
- execution time
- Adaptive parameters:
- linear iterations using solver using *Mj* in windows *i *is* I(i,j)*
- Time per iteration for method *Mj *is *tj*
- linear solution time using solver using *Mj* in windows *i *is* T(i,j)=I(i,j)*tj*
- assume *I(i,j+1)=I(i,j)α, *and* I(i,j-1)=I(i,j)β*;
-  *α , β  *are user defined parameters
- The switching criteria:
- switch up (to method *M**j+1*) if* *
- switch down (to method *M**j-1*) if* *
* *

## Adaptive Solver III(Non-sequence Based Selection)

- Switching based on
- predicted time required to solve linear systems from current to final nonlinear iteration
- Predictions based on polynomial interpolation
- calculate number of steps, *m*, required to converge through polynomial interpolation of  rate of convergence
- time for linear solution taken by method *Mj* at nonlinear iteration *t* is given by an n-ordered polynomial
- The switching criteria:
- At iteration x switch to method Mk such that
- i.e., the method with the least cumulative time is selected

## Adaptive Solvers :Experiments

- Application: Driven Cavity with Pseudo-transient Continuation
- Problem Parameters
- Mesh Size 128 by 128
- Grashof Number [400,450,500,580,600,650]
- Lid Velocity: [10,13,15,20,25]
- Total of 30 problem instances

## Adaptive Solvers:Driven Cavity Flow(600:20)

- Linear solution time per nonlinear iteration for the base and adaptive methods

## Execution Time for Base and Adaptive Solvers

- Cumulative time for simulations over 30 instances
- Adaptive solvers are 39%-32% faster than the average execution time of the base methods

## Static Selection:Using Machine Learning 

- Machine learning methods  are used for selecting a method based on problem features such as;
- application parameters, matrix condition number or sparsity pattern, departure of Jacobian from normality, symmetry, etc.
- Machine Learning has been used to harvest large datasets to associate features of an event with outcomes
- medical diagnoses, retail purchases, hand-written characters, etc.
- Many applications are solved over code lifetime of decades
- We can obtain a vast amount of information about the systems and corresponding solvers
- This information (mountains of performance data) can give us pointers for selecting solvers for future runs
- _A “trained” classifier can reliably decide how to respond to such features in subsequent encounters_
- Input: Database of matrix characteristics X performance with candidate solvers
- Output: Mechanism to predict a set of suitable solvers for a given matrix

## Applying  Machine Learning

- _x0016__x0016__x0016__x0016__x0016__x0003__x0016_

- _x0016__x0016__x0016__x0016__x0016__x0003__x0016_

- _x0016__x0016__x0016__x0016__x0016__x0003__x0016_

- _x0016__x0016__x0016__x0016__x0016__x0003__x0016_

- _x0016__x0016__x0016__x0016__x0016__x0003__x0016_

- Linear Systems B

- Linear Systems A

- Matrix Properties

- Matrix Properties

- Performance Parameters

- List of Solvers

- Database A     (seen)

- Database B (unseen)

- Prediction Rules  (Classifiers)

- Set of Suitable Solvers

- Petsc

- Anamod

- Anamod

- Combine

- Combine

- (Machine Learning Tool)

- (Perl script)

- (Perl script)

- MLJava

- MLJava

- Training Set 2/3 Testing Set 1/3

## Alternating Decision Trees

- An alternating decision tree is a mapping from instances to real numbers which is defined in terms of a set of base rules

- -0.7

- a<4.5

- b>1

- -0.2

- +0.4

- +0.5

- -0.7

- -0.1

- -0.2

- +0.4

- +0.2

- -0.6

- +0.3

- +0.1

- +0.5

- a<1

- b>1

- b>0

- a<4.5

- +0.2

- Decision Tree Instance: a=b=0.5 Classification= -1

- General Alternating Tree            Instance: a=b=0.5                     Classification= sign(.5+.3-.7-.2-.1)= -1

- Alternating Tree            Instance: a=b=0.5   Classification= sign(.5-.7-.2)= -1

- Y

- Y

- Y

- Y

- Y

- Y

- N

- N

- N

- N

- N

- N

- Decision node:
-  If(precondition) then
-       If(condition)
-         then output p1
-       Else output p2
-  Else output 0

## Boosting Algorithm: AdaBoost

- Boosting refers to a general and provably effective method of producing a very accurate prediction by combining rough and moderately inaccurate rules of thumb---A Short Introduction to Boosting
- AdaBoost implements the boosting algorithm as follows;
- Given database (*x**1**,y**1*),....,(*x**m**,y**m*) where *x**i** **** X* and *y**i** **** Y=*{-1,+1}
- Initialize distribution *D**1**(i)=1/m*
- For iterations *t *= 1,...,T
- Train weak learner using distribution Dt
- Get weak hypothesis ht: X{-1,+1}
- Update distribution Dt+1(i) giving more weight to wrongly classified entries
- Final hypothesis H(x) is formed by combining the all the ht
- Applying Boosting to ADTrees:
- Alternating decision trees can be defined as a sum of base rules
- It is easy to apply boosting to ADTrees setting base rules as weak learners
- MLJava implements AdaBoost on ADTrees

**Notes:**
H(x)=sum at*ht  at=.5ln(1-err)/err err Prob(false results).

##  Accuracy:ROC Curves

- ROC (Receiver Operator Characteristic)  curves give the accuracy of the classifier.

- The closer the curve follows the left-hand border and then the top border
-  the more accurate the test.

##   Driven Cavity Flow with Pseudo-transient Continuation(Petsc SNES example ex27)

- Problem Parameters:
- Lid Velocity:  5 10 15 20 25
- Grashof Numbers: 100 500 1000
- Solvers
- KSP: bcgs, tfqmr, fgmres(5,30,60), gmres(5,30,60)
- Preconditioner: asm0
- Subdomain Preconditioner: ilu0,ilu1,ilu2,lu,jacobi,none
- Machine Learning Parameters
- Default Solver: gmres30 with ilu0
- Performance Metric: linear solution time
- Threshold1: 1.5 times better than default
- Threshold2: 2.0 times better than default
- Unseen Dataset: LidVelocity 10, 20 (8736 entries)
- Seen DataSet :Lid Velocity 5,15,25
- Training Set ( 10458 entries)
- Testing Set  (5255 entries)

**Notes:**
Train set 10458
Test set 5255
Trial set 8736

## Results: ROC Curves 

## Results: Performance

##  Linear System from M3D (MHD code)(w3.pppl.gov/~jchen)

- Matrix Properties:
- Number of Non-zeros (of the order 10  ):  7.6 12.6 13.1 18.5  24.2
- Matrices were grouped by their diagonal dominance
- Weakly Diagonally Dominant
- Seen Dataset: nnz 7.6,13.1,18.5
- Training Set ( 524  entries)
- Testing Set ( 298 entries)
- Unseen Dataset: nnz  24.6 (matrix 1),12.6 (matrix 2) (700 entries)
- - Moderately Diagonally Dominant
- Seen Dataset: nnz 7.6,13.1,18.5
- Training Set ( 1644  entries)
- Testing Set ( 870 entries)
- Unseen Dataset: nnz 12.6, 24.6 (1788 entries)
- Strongly Diagonally Dominant
- Seen Dataset: nnz 7.6,18.5
- Training Set ( 1144 entries)
- Testing Set ( 596 entries)
- Unseen Dataset: nnz 12.6, 24.6 (3342 entries)

- 6

**Notes:**
Train set 10458
Test set 5255
Trial set 8736

##   Linear Systems from M3D code

- Solvers
- Direct Methods: SuperLU, Spooles
- Iterative Methods : bcgs, tfqmr, fgmres(5,30,60), gmres(5,30,60)
- Preconditioners:
- Domain Decomposition:  asm0,asm1,asm2
- Subdomain Preconditioner: ilu0,ilu1,ilu2, ilu3, icc0, icc1, icc2, icc3, lu,jacobi,sor,none
- Algebraic Multigrid (Hypre): parasails,boomeramg,euclid,pilut
- Machine Learning Parameters
- Classifier 1
- Performance Metric : Convergence
- Classifier 2
- Performance Metric: linear solution time
- Default Solver: gmres30 with asm0 and ilu0
- Threshold: 1.5 times better than default

**Notes:**
Train set 10458
Test set 5255
Trial set 8736

## Results: ROC CurvesWeakly Diagonally Dominant Matrices 

- test-convergence

## Results: PerformanceWeakly Diagonally Dominant Matrices

- False Positives: Unseen Matrix1  gmres30-boomeramg--diverged

## Results: ROC CurvesModerately Diagonally Dominant Matrices 

- train-time

## Results: PerformanceModerately Diagonally Dominant Matrices

- False Positives: Unseen Matrix2  gmres60-asm0-icc2--memory overflow

## Results: ROC CurvesStrongly Diagonally Dominant Matrices 

## Results: PerformanceStrongly Diagonally Dominant Matrices

## Research Plans

- Variations on Composite Solvers
- Subset of composites
- Composites where failure  rates are dependent
- More Adaptive Heuristics
- Machine Learning
- Applying different machine learning techniques and their comparison
- Classifiers highlight the most important matrix characteristics that determine solver selection. Machine learning can be used to study relationship between linear system properties and corresponding solvers
- In the spirit of the multimethod  approach the selection and combination techniques can be further combined with each other
- Use machine learning to obtain initial set of solvers to be used in adaptive or composite methods
- Create adaptive solvers by switching among different composite solvers
- Use different mutlimethod techniques in different stages of the simulation
**WANTED---more applications!!!**

## Adaptive Solvers: FUN3D

- Linear solution time at each nonlinear iteration for the base and adaptive solvers