Transcript conjugate gradient method - Department of Computer Science

Monica Garika
Chandana Guduru
METHODS TO SOLVE LINEAR SYSTEMS
 Direct methods
Gaussian elimination method
LU method for factorization
Simplex method of linear programming
 Iterative method
Jacobi method
Gauss-Seidel method
Multi-grid method
Conjugate gradient method
Conjugate Gradient method
 The CG is an algorithm for the numerical
solution of particular system of linear equations
Ax=b.
Where A is symmetric i.e., A = AT and positive
definite i.e.,
xT * A * x > 0 for all nonzero vectors
If A is symmetric and positive definite then the
function
Q(x) = ½ x`Ax – x`b + c
Conjugate gradient method
 Conjugate gradient method builds up a solution
x*€ Rn in at most n steps in the absence of round off
errors.
 Considering round off errors more than n steps may be
needed to get a good approximation of exact solution
x*
 For sparse matrices a good approximation of exact
solution can be achieved in less than n steps in also
with round off errors.
Practical Example
In oil reservoir simulation,
The number of linear equations corresponds to the
number of grids of a reservoir
 The unknown vector x is the oil pressure of reservoir
 Each element of the vector x is the oil pressure of a
specific grid of the reservoir
Linear System
Square matrix
Unknown vector
(what we want to find)
Known vector
Matrix Multiplication
Positive Definite Matrix
[ x 1 x 2 … xn ]
>0
Procedure
 Finding the initial guess for solution 𝑥0
 Generates successive approximations to 𝑥0
 Generates residuals
 Searching directions
Conjugate gradient iteration
 x0 = 0, r0 = b, p0 = r0
 for k = 1, 2, 3, . . .

αk = (rTk-1rk-1) / (pTk-1Apk-1) step length

xk = xk-1 + αk pk-1
approximate solution

rk = rk-1 – αk Apk-1
residual

βk = (rTk rk) / (rTk-1rk-1)
improvement

pk = rk + βk pk-1
search direction
 Iteration of conjugate gradient method is of the form
x(t) = x(t-1) + s(t)d(t)
where,
x(t) is function of old value of vector x
s(t) is scalar step size x
d(t) is direction vector
Before first iteration, values of x(0), d(0) and
g(0) must be set
Steps to find conjugate gradient
method
 Every iteration t calculates x(t) in four steps :
Step 1: Compute the gradient
g(t) = Ax(t-1) – b
Step 2: Compute direction vector
d(t) = - g(t) + [g(t)` g(t) / g(t-1)` g(t-1)] d(t-1)
Step 3: Compute step size
s(t) = [- d(t)` g(t)]/d(t)’ A d(t)]
Step 4: Compute new approximation of x
x(t) = x(t-1) + s(t)d(t)
Sequential Algorithm
1) 𝑥0 = 0
2) 𝑟0 := 𝑏 − 𝐴𝑥0
3) 𝑝0 := 𝑟0
4) 𝑘 := 0
5) 𝐾𝑚𝑎𝑥 := maximum number of iterations to be done
6) if 𝑘 < 𝑘𝑚𝑎𝑥 then perform 8 to 16
7) if 𝑘 = 𝑘𝑚𝑎𝑥 then exit
8) calculate 𝑣 = 𝐴𝑝k
9) αk : = rkT rk/pTk v
10) 𝑥k+1:= 𝑥k + ak pk
11) 𝑟k+1 := 𝑟k − ak v
12) if 𝑟k+1 is sufficiently small then go to 16
end if
13) 𝛽k :=( r T k+1 r k+1)/(rTk 𝑟k)
14) 𝑝k+1 := 𝑟k+1 +βk pk
15) 𝑘 := 𝑘 + 1
16) 𝑟𝑒𝑠𝑢𝑙𝑡 = 𝑥k+1
Complexity analysis
 To Identify Data Dependencies
 To identify eventual communications
 Requires large number of operations
 As number of equations increases complexity also
increases .
Why To Parallelize?
 Parallelizing conjugate gradient method is a way to
increase its performance
 Saves memory because processors only store the
portions of the rows of a matrix A that contain nonzero elements.
 It executes faster because of dividing a matrix into
portions
How to parallelize?
 For example , choose a row-wise block striped
decomposition of A and replicate all vectors
 Multiplication of A and vector may be performed
without any communication
 But all-gather communication is needed to replicate
the result vector
 Overall time complexity of parallel algorithm is
Θ(n^2 * w / p + nlogp)
Row wise Block Striped Decomposition of a
Symmetrically Banded Matrix
(a)
(b)
Dependency Graph in CG
Algorithm
of a Parallel CG on each Computing Worker (cw)
1) Receive 𝑐𝑤𝑟0,𝑐𝑤 𝐴,𝑐𝑤 𝑥0
2) 𝑐𝑤𝑝0 =𝑐𝑤 𝑟0
3) 𝑘 := 0
4) 𝐾𝑚𝑎𝑥 := maximum number of iterations to be done
5) if 𝑘 < 𝑘𝑚𝑎𝑥 then perform 8 to 16
6) if 𝑘 = 𝑘𝑚𝑎𝑥 then exit
7) 𝑣 =𝑐𝑤 𝐴𝑐𝑤𝑝𝑘
8) 𝑐𝑤𝑁𝛼𝑘
9) 𝑐𝑤𝐷𝛼𝑘
10) Send 𝑁𝛼𝑘,𝐷𝛼𝑘
11) Receive 𝛼𝑘
12) 𝑥𝑘+1 = 𝑥𝑘 + 𝛼𝑘𝑝𝑘
13) Compute Partial Result of 𝑟𝑘+1: 𝑟𝑘+1 =𝑟𝑘 − 𝛼𝑘𝑣
14) Send 𝑥𝑘+1, 𝑟𝑘+1
15) Receive 𝑠𝑖𝑔𝑛𝑎𝑙
16) if 𝑠𝑖𝑔𝑛𝑎𝑙 = 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛𝑟𝑒𝑎𝑐ℎ𝑒𝑑 go to 23
17) 𝑐𝑤𝑁𝛽𝑘
18) 𝑐𝑤𝐷𝛽𝑘
19) Send 𝑐𝑤𝑁𝛽𝑘 , 𝑐𝑤𝐷𝛽𝑘
20) Receive 𝛽𝑘
21) 𝑐𝑤𝑝𝑘+1 =𝑐𝑤 𝑟𝑘+1 + 𝛽𝑘 𝑐𝑤𝑝𝑘
22) 𝑘 := 𝑘 + 1
23) Result reached
Speedup of Parallel CG on Grid
Versus Sequential CG on Intel
Communication and Waiting Time
of the Parallel CG on Grid
We consider the difference between f at the solution x and any other vector p:
i)
ii )
f p  
1
f x 
1
p Ap  b p  c
t
t
2
x Ax  b x  c
t
t
2
1 t
 1 t

t
t
ii  i ) f ( x )  f  p    x A x  b x  c    p A p  b p  c 
2
 2

sim plify ) f ( x )  f  p  
1
1
2
2
replace b ) f ( x )  f  p  
1
x Ax  x Ax 
x Ax  b x 
t
t
t
t
2

1
1
2
positive-definitess ) f ( x )  f  p   0
t
1
t
p Ap  x Ap
t
t
2
1
x Ax 
t
2

p Ap  b p
p Ap  x Ap
t
2
x  p
t
A x  p
t
Parallel Computation Design
Parallel Computation Design
 – Iterations of the conjugate gradient method can be executed
 only in sequence, so the most advisable approach is to
parallelize the computations, that are carried out at each
iteration
 The most time-consuming computations are the multiplication
of matrix A by the vectors x and d
 – Additional operations, that have the lower computational
complexity order, are different vector processing procedures
(inner product, addition and subtraction, multiplying by a scalar).
 While implementing the parallel conjugate gradient method,
it can be used parallel algorithms of matrix-vector
multiplication,

“Pure” Conjugate Gradient Method (Quadratic Case)
0 - Starting at any x0 define d0 = -g0 = b - Q x0 , where gk is the column
vector of gradients of the objective function at point f(xk)
1 - Using dk , calculate the new point xk+1 = xk + ak dk , where
gkTd k
ak=- T
d k Qd k
2 - Calculate the new conjugate gradient direction dk+1 , according to:
dk+1 = - gk+1 + bk dk
where
gk+1 TQd k
bk=
d kTQd k
ADVANTAGES
Advantages:
 1) Gradient is always nonzero and linearly
independent of all previous direction vectors.
 2) Simple formula to determine the new direction.
Only slightly more complicated than steepest descent.
 3) Process makes good progress because it is based
on gradients.
ADVANTAGES
 Attractive are the simple formulae for updating the
direction vector.
 Method is slightly more complicated than steepest
descent, but converges faster.
 Conjugate gradient method is an Indirect Solver
 Used to solve large systems
 Requires small amount of memory
 It doesn’t require numerical linear algebra
Conclusion
 Conjugate gradient method is a linear solver tool
which is used in a wide range of engineering and
sciences applications.
 However, conjugate gradient has a complexity
drawback due to the high number of arithmetic
operations involved in matrix vector and vector-vector
multiplications.
 Our implementation reveals that despite the
communication cost involved in a parallel CG, a
performance improvement compared to a sequential
algorithm is still possible.
References
 Parallel and distributed computing systems by Dimitri




P. Bertsekas, John N.Tsitsiklis
Parallel programming for multicore and cluster
systems by Thomas Rauber,Gudula Runger
Scientific computing . An introduction with parallel
computing by Gene Golub and James M.Ortega
Parallel computing in C with Openmp and mpi by
Michael J.Quinn
Jonathon Richard Shewchuk, ”An Introduction to the
Conjugate Gradient Method Without the Agonizing
Pain”, School of Computer Science, Carnegie Mellon
University, Edition 1 1/4