# Presentation: 098130

## Subspace Projected Approximate Matrix 			The SPAM MethodM. Minkoff, R. Shepard, A. Wagner

- Outline
- Background on SPAM Approach
- Results
- Cumulative Reaction Probabilities
- Initial Results

- CMCIM02

## Subspace Projected Approximate Matrix: SPAM Method

**• ****Subspace iterative solution to eigenvalue, linear, and nonlinear problems**
**	****e.g. eigenvalue problem  **(**H **–** ******j**** **)**v****j**** **=** 0**
**•**** ****Subspace iterative solutions have form  v****j**** = X****n**** c****j**** **
**	****	****where:****	****X****n**** = x****1****,x****2****,...,x****n**** **
**• ****c****j**** is solved in a subspace:  H****n ****c****j**** = ********n****j ****c****j**** **
**	****	****where:****	**** H****n**** = (W****n****)****T ****X****n**** **
**	****	****	****where****: W****n ****= H X****n****    <---****for N>>n, where all the work is**
**• ****New x****n+1**** vector from residual:  r****n+1**** = (W****n****-********n****j ****X****n****)c****j**
**• ****SPAM gives more accurate or faster converging way to get r****n+1 ****in 3 steps**

- CMCIM02

## Subspace Projected Approximate Matrix: SPAM Method

**	****P****n ****r****n+1**** (=X****n****[X****n****]****T ****r****n+1****)****	****	****Q****n ****r****n+1**** (= r****n+1**** - P****n ****r****n+1****)**
**     ****where**
**       ****n trial vectors     n+1 trial vector  exact         matrix-vector         subspace**
**   ****(already processed)****	****   r****n+1****	****      matrix            product                matrix**
**                ****(X****n****) ****	****	**** ****	**

**S****ubspace ****P****rojected ****A****pproximate ****M****atrix: ****SPAM**** **_**Method**_

**minkoff:**
- The text after “where” is a bit to close to each (horizontally).  Perhaps a smaller font would work.  (e.g. vector and exact are close).

- CMCIM02

## Subspace Projected Approximate Matrix: SPAM Method

**Step 2: Decomposition and Approximation of H r****n+1**
**Hv = (P****n****+Q****n****) H (P****n****+Q****n****) r****n+1**
**	****          ****cheap subspace operations****	****	****     expensive full space**
**	****	****	****	****	****	****            ****operation**

**        ****cheap approximation**
**     ****to full matrix**

**Step 3: Solve approximate subspace problem**
**	****- solution is r****n+1**

- CMCIM02

## Subspace Projected Approximate Matrix: SPAM Method

**SPAM properties (broad view):**
**Relation to other subspace methods (e.g., Davidson)**
**More flexible (sequence of approx. matrices - no sequence => SPAM= Davidson)**
**If iteration tolerances are correct, always no worse than Davidson **
**Relation to multigrid methods**
**SPAM sequences of approx. matrices ~ multigrid sequences of approx. grids**
**Subspace method => solution vector composed of multiple vectors **
**Multigrid methods => single solution vector that is updated**
**Relation to preconditioned conjugate gradient (PCG) methods**
**SPAM has multiple approximations always improving by projection**
**	****PCG has a fixed single preconditioner**
**Deep injection of physical insight into numerics**
**Application experts can design physical approximation sequences**
**SPAM maps approximation sequences onto numerics sequences**
**Projection operators continually improve the approximations**
**	****=> coarse approximations can still be numerically useful**

**minkoff:**
- PCG has a fixed single preconditioner using projection to improve a single vector

- CMCIM02

## Outline of the Davidson Method
Generate an initial vector x1
MAINLOOP: DO n =1
Compute and save  wn = Hxn
Compute the n-th row and column of H:  H1:n,n = wnT X[n]
Compute the subspace eigenvector and value: (H – )c = 0
Compute the residual:  r=W1:nc1:n – X1:nc1:n
Check for convergence using |r|, c, , etc. 
IF (converged) THEN 
EXIT MAINLOOP
ELSE
Generate a new expansion vector  xn+1 from r, , v=Xc, etc. 
ENDIF
ENDDO MAINLOOP

- Generate an initial vector **x**1
- MAINLOOP: DO *n *=1
- Compute and save  **w***n *= **Hx***n*
- Compute the *n*-th row and column of **H**:  **H**1:*n,n* = **w***n**T ***X**[*n*]
- Compute the subspace eigenvector and value: (**H** – **)**c **= **0**
- Compute the residual:  **r**=**W**1:*n***c**1:*n* – ****X**1:*n***c**1:*n*
- Check for convergence using |**r**|, **c**, **, etc.
- IF (*converged*) THEN
- EXIT MAINLOOP
- ELSE
- Generate a new expansion vector ** x***n+1* from **r**, **, **v**=**Xc**, etc.
- ENDIF
- ENDDO MAINLOOP

- Outline of the SPAM Method
- Generate an initial vector **x**1
- Set *wtype*1=1  *! Start the iterations with approximate products*
- Set *n*0=0; *n*=1
- MAINLOOP: DO
- Compute and save **w***n *= **H**(*wtype**n*,* n*0) **x***n*
- Compute the *n*-th row and column of **H**:  **H**1:*n,n* = **w***n**T ***X**[*n*]
- Compute the subspace eigenvector and value: (**H** – **)**c **= **0**
- Compute the residual:  **r**=**W**1:*n***c**1:*n* – ****X**1:*n***c**1:*n*
- Check for convergence using |**r**|, **c**, **, etc.
- IF (*converged* .AND. *wtype**n*.EQ.0) then
- EXIT MAINLOOP *! Final convergence is achieved*
- ELSEIF (*converged* .AND. *wtype**n*≠0) then
- Contract
- Set *n*=*n*0+1; *n*0=* n*1
- Set *wtype**n*=0  *!  The next product will be exact*
- ELSE
- Set *n**n*+1
- Generate a new expansion vector ** x***n* from **r**, **, **v**=**Xc**, etc.
- Set *wtype**n* =1  *!  The next product will be approximate*
- ENDIF
- ENDDO MAINLOOP

- CMCIM02

## Subspace iterative solution to eigenvalue problem:

- (**H **–** ******j**** **)**v****j**** **=** 0**

** ****subspace iterative solutions have form  v****j**** = X****n**** c****j**** **

** ****where:****	****X****n**** = x****1****,x****2****,...,x****n**** **
**c****j**** is solved in a subspace:  H****n ****c****j**** = ********nj ****c****j****  **

- CMCIM02

## That subspace Hn is of the form:

**That subspace H****n**** is of the form:**

**where:**

** **

- CMCIM02

## For N>>n, work/iteration = H xn 

New x vector from residual:  xn+1 = (Wn-nj Xn)cj
SPAM gives better way to get xn+1 via 4 steps:

• Define Projection Operators P and Q:
  Pn	= Xn ((Xn)T Xn)-1 (Xn)T
= Xn (Xn)T  with orthonormal basis vectors,
   Qn	= (1–Pn)

  P v = the Xn components of v
  Q v = the part of v that is not spanned by Xn 

• Decomposition of H with P and Q:
  H	= (Pn+Qn) H (Pn+Qn)
= PnHPn+ PnHQn+QnHPn+QnHQn
= Xn Hn(Xn)T + Xn(Wn)TQn + QnWn(Xn)T + QnHQn
           blue	=> cheap subspace operations  
           red	=> expensive matrix vector multiply

• Approximate QnHQn with cheap approximate H(1):
  H1 = XnHn(Xn)T+Xn(Wn)TQn+QnWn(Xn)T + QnH(1)Qn

• Iteratively solve for xn+1:
  xn+1 = wj            where (H1 – j )wj = 0 

** **
**New x vector from residual:  xn+1 = (Wn-********nj Xn)cj**
**SPAM gives better way to get xn+1 via 4 steps:**
** **
**• **_**Define Projection Operators P and Q**_**:**
**  ****Pn****	****= Xn ((Xn)T Xn)-1 (Xn)T**
**= Xn (Xn)T  with orthonormal basis vectors,**
**   ****Qn****	****= (1–Pn)**
** **
**  ****P v = the Xn components of v**
**  ****Q v = the part of v that is not spanned by Xn **
** **
**• **_**Decomposition of H with P and Q**_**:**
**  ****H****	****= (Pn+Qn) H (Pn+Qn)**
**= PnHPn+ PnHQn+QnHPn+QnHQn**
**= ****Xn Hn(Xn)T + Xn(Wn)TQn + QnWn(Xn)T**** + ****QnHQn**
**           ****blue****	****=>**** cheap**** ****subspace operations****  **
**           ****red****	****=> ****expensive**** matrix vector multiply**
** **
**• **_**Approximate **__**QnHQn**__** with **__***cheap***__** approximate H(1)**_**:**
**  *****H*****1 = ****XnHn(Xn)T+Xn(Wn)TQn+QnWn(Xn)****T**** + ****Qn****H****(1)****Qn**
** **
**• **_**Iteratively solve for xn+1**_**:**
**  ****xn+1 = wj            where **(***H*****1 **–** ******j **)**wj **=** 0**

- CMCIM02

## What might H(1) be?

** **
**Anything to make ****matrix vector multiplies cheaper**
** **
**•****	****Sparser**
**•****	****Smaller underlying basis**
**•****	****Lower-order expansion of matrix elements**
**•****	****Coarser underlying grid **
**•****	****Lower-order difference equation**
**•****	****Tensor-product approximation **
**•****	****Operator approximation **
** **
** **
_**SPAM has great properties**_**:**
** **
**•****	****projection operators => convergence from any H****(1)**** **
**•****	****multi-level SPAM with dynamic tolerances**
**	****Q****n****HQ****n**** approximated by ****Q****n****H****(1)****Q****n**** **
**	****   ****Q****n****H****(1)****Q****n**** approximated by ****Q****n****H****(2)****Q****n**
**	****      ****Q****n****H****(2)****Q****n**** approximated by ****Q****n****H****(3)****Q****n****...**
**•****	****tens of lines of additional code to existing iterative subspace eigensolver**
**•****	****applicable to any subspace problem (linearsolves?)**

- CMCIM02

## SPAM Example 1:

** **
**Banded Matrix:**
**H characterized by:**
**	****N = matrix dimension**
**	****∆ ****= off diagonal "increment"**
**	****W****0**** = width of band**
** **
**H****(i)**** has W****i**** < W****i-1**

** **

- CMCIM02

## Slide 12

- CMCIM02

## Slide 13

- CMCIM02

## Tensor-Product Convergence Results for Multiple Eigenvectors

- Convergence summaries of the lowest 10 roots of the *m*=8 and *m*=10 perturbed-tensor-product matrices described in the text.  The initial vectors in all cases are the eigenvectors of the tensor-product matrices, which were computed as tensor-products of the eigenvectors of the 4x4 component matrices.  The matrix-vector product counts are the totals for all 10 roots.  For the *m*=8 calculations, *N*=65,536, **=10, and |**r***j*|<10–1.  For the *m*=10 calculations, *N*=1,048,576, **=100, and |**r***j*|<100.

- CMCIM02

##  How Fast Do Molecules React? MPP Simulation of Exact Cumulative Reaction Probabilities

- 	• Combustion modeling requires accurate determination of rate constants which are approximated by statistical simulation techniques.  Exact calculations of rate constants can be obtained from Cumulative Reaction Probabilities (CRP), however this can be done for only a few degrees of freedom (DOF) problems.  Initial computational results suggest that the statistical methods may increasing differ with exact calculations as DOF increases.  Larger problems (up to 10 DOF) can be addressed only via MPP calculations.  (This work is a collaborative project in chemical kinetics modeling (ANL/CHM) and massively parallel simulation (ANL/MCS)._**	**_
- 	   • CRP simulation involves calculation of a few eigenvalues of the in an outer iteration with an inner iteration evaluation of a Greens function dealing the solution of two Hamiltonian (linear) systems.
- 	• PETSc provides the capability to develop rapid-prototyping for the MPP software development of CRP simulations and an software environment for developing novel preconditioners for these and other linear systems and iterative eigenvalue solvers.

- CMCIM02

## Advanced Software for the Calculation of Thermochemistry, Kinetics and DynamicsParallelization of Cumulative Reaction Probabilities 

**Cumulative Reaction Probabilities (CRP)**
**	****- computational core of reaction rate constants**
**	****- exact computation computational intensive**
**	****- approximate computation underlies **
**	****  ****all major reaction rate theories in use**
**	****=> efficient exact CRP code will**
**	****	****- give exact rates (if the computed forces are accurate)**
**	****	****- calibrate ubiquitous approximate rate methods**
_**Time Independent**_** (Miller and Manthe, 1994, and others)**
**Highly parallel approaches to both methods being pursued**

- CMCIM02

## Parallelization of Cumulative Reaction Probabilities Time Independent Approach

**• ****N(E) = ********k ****p****k****(E,J)**
**  ****where p****k****(E,J) = eigenvalues of Probability Operator:**
**   ****	****P(E) = 4 ********r****1/2**** (H+i********-E)****-1********p**** (H-i********-E)****-1********r****1/2**
**	****where i********x**** = absorbing diagonal potentials**
**                      ****= imaginary potentials**
**	****	****   ****H = hamiltonian (differential operator)**
**	****	****      ****=> for realistic problems**
**	****	****           ****size ~10****5****x10****5**** or much larger**
**	****	****          ****number of eigenvalues < 100**
**• ****iterative approach**
**	****macrocycle of iteration for eigenvalue**
**	****microcycle of iteration for**
**	****    ****action of Green’s function (H-i********-E)****-1**** **
**	****    ****=> linear solve**
**	****	****	**** **

**minkoff:**
- Can the text in the figure be bolder or in a different font?

- CMCIM02

## Parallelization of Cumulative Reaction Probabilities Time Independent Approach

**• ****PETSc: data structure, GMRES linear solve, preconditioners**
**  ****USER: Lanczos method for eigensolve**
**• ****Future: user supplied preconditioners**

**Parallelization of ****C****umulative ****R****eaction ****P****robabilities **_**Time Independent Approach**_

- CMCIM02

## Chemical Dynamics Theory3 angles, 3 stretches6 degrees of freedom

- CMCIM02

## Chemical Dynamics Theory

- Probability Operator and It’s Inverse
- Using probability method calculates a few large eigenvalues via iterative methods.  The iterative evaluation involves the action of two Green’s function.
- Using inverse probability method involves a direct calculation each iteration to obtain a few smallest eigenvalues.  At each iteration the action of a vector by the Green’s function is required.  This leads to solving linear systems involving the Hamiltonian.

- CMCIM02

## Parallelization of Cumulative Reaction Probabilities Time Independent Approach

**Performance:**
**model problem**
** ****Optional number of dimensions**
**  • ****Eckhart potential along rxn coord.**
**  • ****parabolic potential perpendicular**
**    ****to reaction coord.**
** ****DVR representation of H**
**Algorithm options**
** ****Diagonal preconditioner**
** ****other PETSc preconditioners slower**** **
**   ****for this problem**
**Computers**
** ****NERSC SP**
** ****others include SGI Power Chanllenge,**
**  ****Cray T3E**

**Performance vs. processors**
**for model with increasing dimensions**

**minkoff:**
- Other PETSc preconditioners slower for this problem
- Correct spelling of Power Challenge
- On the figure should the text be bolder and also should the lines be bolder

- CMCIM02

## Parallelization of Cumulative Reaction Probabilities Time Independent Approach

**Future Preconditioners:**
**• ****SPAM (see next poster)****	****	****             ****great ****	****	****       ****	****         scales**
**	****     ****Exact H****	****	****	****            performance ****	****	****	****         poorly**
**	****	****	****	****	****	****	****	****	****        ****- storage**
**	****	****	****	****	****	****	****	****	****        ****- banded**
**	****	****	****	****	****	**** ****	****	****	****          ****matrix**
**	****	****	****	****	****	**** ****	****	****	****          ****inversion**
**• ****Sparse optimal similarity transforms **
**	****IF Q =                         THEN find optimal Q such that QHQT =**
**	****	****	****	****	****	*****optimal block diagonal***** ->**

**fat band**
**precond.**
**SPAM**
**precond.**

**can SPAM get scalable performance?**

**minkoff:**
- Can the text in the figure be larger?

- CMCIM02

## Combustion Kinetics

**Parallel computing dramatically lowers time-to-solution for rate constants **
**ANL software libraries (PETSc) important**
**Parallel computing is only way to solve large problem:  more processors => larger problems **
**Results obtained on NERSC IBM/SP.**

- CMCIM02

## Extrapolated Storage Requirements

- CMCIM02