Computing Matrix Functions on theK Computer
William Dawson
Takahito Nakajima
RIKEN Center for Computational Science

2018年度第2回計算科学フォーラム!1

Computational Molecular Science ResearchTeam• Also in collaboration with Luigi Genovese (French Alternative Energies and Atomic Energy

Commission), Marco Zaccaria (Boston College), Massimo Reverberi (Sapienza University ofRome).

!2

Introduction to Computational Chemistry• In our team, we work to develop software and methods
for understanding molecules and materials fromaccording to the laws of quantum mechanics.
AtomicPositions

QuantumChemistry
Software

SystemProperties

ForceIntegrator• Quantum chemistry software gives us atomic level insight,

allowing us to go beyond the technical limits of experiment.

• Calculation quantities: band gaps, chemical reactions, rate
constants, durability, etc.

!3Introduction - Large Scale Calculations• Large uniform environment withsmall perturbation (e.g. dilutesolutions).

• Large Nanostructures (e.g. CarbonNanotubes).

• Molecules in a realistic environment(e.g. proteins).

• Not only perform calculations on thesystem, but also to gain insight into
the actual chemistry.

• These calculations require cleveralgorithms, and large computationalresources.

Ratcliff, Laura E., Stephan Mohr, Georg Huhs, Thierry Deutsch, Michel Masella, and Luigi Genovese. "Challenges in large scale quantum mechanicalcalculations." Wiley Interdisciplinary Reviews: Computational Molecular Science 7, no. 1 (2017): e1290.

!4

Outline• Introduction to Matrix

Functions• Formal Definition• Methods of Computing
Matrix Functions• Motivating Matrices• NTPoly introduction• Parallelization Techniques• Distributed Memory•• Usability Considerations• Data Distribution• Programming LanguageSupport• Example Applications• Quantum Chemistry• Social Network Analysis• Search Engine
Optimization

Parallelization
On Node Parallelization!5

Introduction to Matrix Functions!6

Introduction to Matrix Functions
Cauchy integral definition:

where f is analytic on and inside a closed contour 𝚪 that encloses thespectrum of A.
Simple definition: We are all familiar with functions of a single variable f(x).
In the matrix function case, just replace the variable x with a matrix A.
Standard Function

Matrix Function

Interpretationf(x) = x2

f(A) = A2

Matrix Productf(x) = 1/x

f(A) = A-1

Matrix Inverse

f(x) = exf(A) = eAMatrix Diﬀerential Equationf(x) = sign(x)f(A) = sign(A)Projection on to Subspace

Higham, Nicholas J. Functions of matrices: theory and computation. Vol. 104. Siam, 2008.

!7

Motivating Applications• Solving the generalized eigenvalue equation.
Ax = 𝝀Bx => B-1/2AB-1/2x = 𝝀x (if B is SPD).

• Constructing good preconditioners.
Ax = c => BAx = Bc, where B ≅ A-1.

• Computing centrality measures of a network.
A: the adjacency matrix of a graph.
Kantz centrality: (I - 𝝰A)-1
Estrada Centrality: eβA• Solution to Sylvester equation, algebraic Ricatti, etc.
!8Motivating Applications - Chemistry• Diagonalization free methods for quantum chemistry.
• Given the hamiltonian matrix H, we wish to construct thedensity matrix D.

• Usually we do this by solving the eigenvalue equation:
HV = 𝝀V (where V is a nbasis x nelectrons matrix).

• And compute the density matrix as D = VVT.
• Instead, we can compute the density matrix directly using
the fermi function:

D = 1/(eβ(H-μ) + I).

!9Methods for Computing Matrix Functions1. Diagonalization (if possible):
A = ZDZ-1.
f(A) = Zf(D)Z-1.
2. Schur Decomposition (explicit formulas exist for upper
triangular matrices exist).
3. Taylor series expansion:cos(A) = I - A2/2! + A4/4! + A6/6! …
4. Polynomial Approximation (and Rational Approximation).
5. From each function’s definition (A-1 : A*X = I).

!10Motivating Matrices• In many domains, the problem of interest can be represented using a sparse,
hermitian matrix.

• Under certain conditions, not only is the matrix A sparse, but also some matrix
functions f(A) are sparse.

• Estrada matrix exponential eβA contains a scaling factor which might be interpreted
as a unit of edge weight.

• For certain values of β, the matrix exponential of small world matrices is alsosparse.

102Sparsity(%)

101=1= 12100= 14= 18
1= 16
1= 32

29

210

211

212
Nodes

213

214

!11

Motivating Matrices - Chemistry• For insulating systems (and metals at high temperature), it
is known that the density matrix is sparse whenrepresented in a localized basis.

• Example: the Hamiltonian and Density Matrix of 1920
water molecules computed using the BigDFT code.

Hamiltonian

Density Matrix

http://bigdft.org/ !12

Sparsity Aware Matrix Function Calculation• From the list of methods for computing matrix functions, wewill select calculations based on matrix polynomials.
e.g. Chebysehv polynomials: f(A) ≅ ∑ciTi(x).
T0(A) = I

T2(A) = 2A2 - I

T4(A) = 8A4 - 8A2 + I

T1(A) = A

T3(A) = 4A3 - 3A…

• Computing a matrix polynomial requires only two coreroutines: matrix addition, matrix multiplication.

• Easy to parallelize.
• Many functions can be tuned through just two routines.
• In the case of sparse matrices, we replace these kernels withsparse matrix addition, and sparse matrix multiplication.

!13

NTPoly - A Library for Computing Matrix
Functions• General Polynomials• Standard Polynomials• Chebyshev Polynomials• Hermite Polynomials• Transcendental Functions• Trigonometric Functions• Exponential and Logarithm• Matrix Roots• Square Root and Inverse• Quantum Chemistry• Density Matrix Minimization• Density Matrix Purification• Chemical Potential Calculation• Density Matrix Extrapolation• Other• Matrix Inverse (and Moore-• Matrix pth Root and Inversepth root

• Parallel File I/O• MIT License (available onSquare Root

Penrose Inverse)• Sign Function/Polar
DecompositionGithub)Dawson, William, and Takahito Nakajima. "Massively parallel sparse matrix function calculations with NTPoly."
Computer Physics Communications 225 (2018): 154-165.

!14

Parallelization Techniques!15Matrix Multiplication Parallelization - 2.5D
1D Algorithm: Each processor has a matrix sliceStart With A SliceReceive Columns

Compute Element

2D Algorithm: Each processor has a matrix block
Start With A Block

Receive Blocks

Compute Block

Schatz, Martin D., Robert A. Van de Geijn, and Jack Poulson. "Parallel matrix multiplication: A systematic journey." SIAM Journal on
Scientific Computing 38, no. 6 (2016): C748-C781.

!16

Matrix Multiplication Parallelization - 2.5D
2.5D Algorithm: Duplicate In Z Direction
Start With A Block

Gather Blocks

Partial Contribution

Block Copy

Gather Partial

Partial 2Reduce Sum

Solomonik, Edgar, and James Demmel. "Communication-optimal parallel 2.5 D matrix multiplication and LU factorization algorithms." In European
Conference on Parallel Processing, pp. 90-109. Springer, Berlin, Heidelberg, 2011.

!17

Matrix Multiplication Parallelization - OpenMP• Important to have a hybrid OpenMP/MPI implementation to target
future architectures.

• Main idea: thread parallel over matrix blocks.
• Local blocked matrix multiply works like dense multiply.
Local Matrix A

Local Matrix B

Local Matrix C• Can also block the communication, allowing for overlapping ofcommunication and computation.

• Little overhead for blocking.
!18

Matrix Multiplication Parallelization - OpenMP• OpenMP loop parallelism doesn’t work well withoverlapping communication. Instead we use OpenMPtask framework.

• Creating a task manager, and dependency graph.
Local Gather Row Ai

Transpose Column Bj

Transpose Column Bk

Global Gather Row Ai

Local Gather Column Bj

Local Gather Column Bk

Transpose Row Ai

Global Gather Column Bj

Global Gather Column Bk

Compute Block Cij

Compute Block Cik

Sum Z Direction Cij

Sum Z Direction Cik!19

Usability Considerations!20

Usability Considerations• Challenge: Integrate With Parallel Programs Using aVariety of Diﬀerent Data Layouts.

• Solution 1: Parallel File I/O through the standard matrix
market format for rapid prototyping.

• Solution 2: Arbitrary Data Remapping Routines.
120

2D All-To-All GatherInvert
Read From File
Redistribute

Sum Blocks
10080Time(s)Start With Arbitrary Elements

60

40

20

00

10000

20000

30000
Cores

40000

50000

60000

!21

Usability Considerations - 2
Challenge: Integration with codes written in a variety of
programming language.
Solution: Programming Language Wrapper Hierarchy.
Fortran Implementation Base
Handle Based Fortran Wrapper (Only Primitive Types)

C++ Object Oriented Wrapper

C InterfaceSWIG Wrapper GeneratorPython

RubyETC!22

Programming Language Support - Details
Using complex data types makes life easier in Fortran, but
makes it harder to call from other languages.
! Complex Fortran Data TypeTYPE :: DistributedSparseMatrix_t! Simple data
INTEGER :: matrix_dimension
INTEGER :: start_column, end_column…
! Variety of members
TYPE(ProcessGrid_t) :: grid
! Also contains allocatable subtypes
TYPE(LocalMatrix_t), DIMENSION(:,:), ALLOCATABLE :: local_data…
END TYPESUBROUTINE ComputeExponential(InputMat, OutputMat)
TYPE(DistributedSparseMatrix_t), INTENT(in):: InputMat

TYPE(DistributedSparseMatrix_t), INTENT(inout) :: OutputMat! Solver Logic
END SUBROUTINE

Pletzer, Alexander, Douglas McCune, Stefan Maszala, Srinath Vadlamani, and Scott Kruger. "Exposing Fortran derived types to C
and other languages." Computing in Science & Engineering 10, no. 4 (2008): 86-92.

!23

Programming Language Support - 2To simplify things, we will only expose handles to data objects.
TYPE :: DistributedSparseMatrix_wrp ! Handle Datatype
TYPE(DistributedSparseMatrix_t), POINTER :: DATAEND TYPESUBROUTINE ConstructMatrix_wrp(ih_this)INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) ! SIZE_wrp is size of a pointer struct.
TYPE(DistributedSparseMatrix_wrp) :: this
ALLOCATE(this%data)
ih_this = TRANSFER(this,ih_this) ! Convert between handle and integer.
END SUBROUTINESUBROUTINE ComputeExponential_wrp(ih_InputMat, ih_OutputMat) bind(c,name=“ComputeExponential_wrp")
INTEGER(kind=c_int), INTENT(in):: ih_InputMat(SIZE_wrp)

INTEGER(kind=c_int), INTENT(inout) :: ih_OutputMat(SIZE_wrp)
TYPE(DistributedSparseMatrix_wrp) :: InputMat
TYPE(DistributedSparseMatrix_wrp) :: OutputMat
InputMat = TRANSFER(ih_InputMat,InputMat)
OutputMat = TRANSFER(ih_OutputMat,OutputMat)
CALL ComputeExponential(InputMat%data, OutputMat%data)
END SUBROUTINE!24

Programming Language Support - 3
C Interface is now simple to expose.
// C Routine To Callvoid ConstructMatrix_wrp(int *ih_this);void ComputeExponential_wrp(const int *ih_Input, int *ih_Output);

C++ Uses the same interface
class DistributedSparseMatrix {public:DistributedSparseMatrix() {ConstructMatrix_wrp(this->handle);}
int handle[SIZE_wrp];};
void ComputeExponential(const DistributedSparseMatrix &InputMat,DistributedSparseMatrix &OutputMat) {ComputeExponential_wrp(InputMat.handle, OutputMat.handle);}!25

Example Applications!26

Quantum Chemistry• Standard eigensolvers can make limited use of the sparsity of a matrix,
but will be outperform by matrix function based approaches.

• Calculation of water clusters of various sizes, 631G basis set, using theTRS2 density matrix method to approximate the fermi function.

!27

Quantum Chemistry - 2• Use of communication avoiding algorithms and task

based openmp parallelization allows for calculations
using tens of thousands of cores.

!28

Social Network Analysis• Estrada’s Scaled Matrix Exponential Metric: eβA• Example, social networks:

Node B

Node F

Node A

Node H

Node D

Node E

Node C

Node G

Node I

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks." Physics reports 514, no. 3 (2012): 89-119.

!29

Social Network Analysis• Estrada’s Scaled Matrix Exponential Metric: eβA• Example, social networks:

Node B

Node F

Node A

Node H

Node D

Node E

Node C

Node G

Node I!30

Social Network Analysis• Estrada’s Scaled Matrix Exponential Metric: eβA• Example, social networks:

Node B

Node F

Node A

Node H

Node D

Node E

Node C

Node G

Node I!31

Social Network Analysis• Estrada’s Scaled Matrix Exponential Metric: eβA• Example, social networks:

Node B

Node F

Node A

Node H

Node D

Node E

Node C

Node G

Node I!32

Social Network Analysis - In Practice• Network Resiliency calculations:• Data Set: Israeli Social Network “TheMarker Cafe”• Nodes: 69413. Sparsity: 0.04%.
• Procedure:• Remove a node from the graph• Random Node
• Node with largest degree• Compute the matrix exponential• Compute the sparsity• Repeat50
45

Average Number of Edges

40
3530
25
20
15Random
Largest Degree

100500

1000

1500
2000
Nodes Removed

2500

3000!33

Search Engine Optimization• Finding important nodes in a directed network.
• Going beyond page rank by ranking nodes as Authorities andHubs. Authorities are important nodes, hubs point toauthorities.

• Matrix functions to compute:• Hub similarity matrix: H = cosh(√(AAT)).
• Authority similarity matrix: A = cosh(√(ATA)).
• Calculation Method:• Scaling and squaring method with Chebyshev polynomials:
cosh(√4x) = 2cosh2√x - 1.

Benzi, Michele, Ernesto Estrada, and Christine Klymko. "Ranking hubs and authorities using matrix functions." Linear Algebra and its Applications 438,
no. 5 (2013): 2447-2474.

!34

Search Engine Optimization - 2• Authority similarity for the High Energy Physics
Phenomenology citation graph.

• 34,546 papers, 421,578 citations.
• Largest eigenvalue is > 3000, so some scaling factor isnecessary.

• Sparsity:• A: ~0.04%.
• AAT: ~1%• cosh(√(βAAT)): 31%.

!35

Conclusion• Matrix functions have a number of diﬀerent applications,motivating the creation of new libraries.

• Using sparse matrix algebra techniques, low order scaling

techniques exist to compute the functions of sparse matrices.

• A combination of communication avoiding algorithms andtask based parallelization enable strong scaling, massively
parallel calculations.

• Automated data redistribution techniques, parallel file i/ousing standard formats, and support for many programming
languages leads to an easy to use library.

• Applications to quantum chemistry, graph analysis, and more.
• https://william-dawson.github.io/NTPoly/!36