Finite and Infinite Dimensional

Matrices as Operators

• We have discussed matrices. These are “operators” in finite (n) dimensional linear space.

• You can generalize to infinite dimensional space.

— e.g., space of “suitable” functions f(x)

• \frac{d}{dx} is an operator!

• \frac{d}{dx}f(x) = g(x) is another function.

— e.g., another operator

\frac{d}{dx} as an Operator and Its Scalar Products

• Consider the Integral

I = ∫_{a}^{b} g(x)\frac{d}{dx}f(x)dx = ⟨g|D|f⟩

where D is operator \frac{d}{dx}.

• Integrate by parts.

I = − ∫_{a}^{b} f(x) \frac{d}{dx} g(x) dx
+ [f(x)g(x)]_{a}^{b}
= f(b) g(b) − f(a) g(a)

• The last term [f(x)g(x)]_{a}^{b} = f(b) g(b) − f(a) g(a)

• We need to get rid of this surface (endpoint) term.

i\frac{d}{dx} as a Hermitean Operator

• Suppose we consider periodic functions

f(b) = f(a)
g(b) = g(a)

— i.e., impose particular boundary conditions

• Then surface term vanishes and

∫_{a}^{b} g(x)\frac{d}{dx}f(x) dx = − ∫_{a}^{b} f(x)\frac{d}{dx}g(x) dx

• The operator L = i\frac{d}{dx} with \underline{boundary conditions} f(a) = f(b) is Hermitean. The i in L cancels minus sign.

• Must include boundary conditions when discussing operators in infinite dimensional space

The Laplacian as a (Matrix) Operator

• ∇^{2} = \frac{d^{2}}{dx^{2}} = −L L

• (∇^{2})^{†} = −L^{†}L^{†} = −L L = ∇^{2}

• Note: (A B)^{†} = B^{†}A^{†} in general

• Thus, ∇^{2} (\frac{d^{2}}{dx^{2}} in one dimension) is Hermitean (symmetric) operator

Mapping of Function Spaces

to a Finite Number of Dimensions

• In solving differential equations numerically, we take infinite dimensional operators on space of functions and convert them to finite matrix operators in space of n dimensional vectors.

• The different numerical approaches correspond to different ways of mapping the infinite dimensional space to a finite space—there is no unique way and some are much better than others.

Mathematical and Pseudo Code Form

of Gauss Seidel Iteration Method

Choose an initial guess x^{(0)} to the solution x. for k = 1,2,…

for i = 1,2,…,n σ = 0 for j = 1,2,…,i − 1

σ = σ + a_{i,j} x_{j}^{(k)} end for j = i + 1,…,n

σ = σ + a_{i,j} x_{j}^{(k−1)} end x_{i}^{(k)} = (b_{i} − σ) /a_{i,i}

end check convergence; continue if necessary

end

Pseudo code for Gauss-Seidel Method

Mathematical (Matrix) Form of Gauss Seidel

• Mathematically, Gauss Seidel is

x_{i}^{(k)} = \frac{[b_{i} − _{j<i}^{∑} a_{i,j} x_{j}^{(k)} − _{j>i}^{∑} a_{i,j} x_{j}^{(k−1)}]}{a_{i,i}}

• In matrix language \underline{x}^{(k)} = (D − L)^{−1} (U\underline{x}^{(k−1)} + \underline{b})

• where D=Diagonal, L=lower, and U=Upper parts of matrix A

• A= D − L − U

Parallelism in Gauss-Seidel Iteration

• Conventional Gauss-Seidel canNOT be easily parallelized because it has sequential update order

x_{n}^{(k)} requires x_{n−1}^{(k)} …
x_{n−1}^{(k)} requires x_{n−2}^{(k)} …
⋮
x_{2}^{(k)} requires x_{1}^{(k)}, etc.

• There are two parallelization strategies

1. Keep same algorithm—pipelined wavefront method

2. Change algorithm—red/black two-phase

6 × 6 Matrix Example Stencil

• We can exploit special structure of matrix of A (e.g., x_{16}^{(k)} requires x_{15}^{(k)}, x_{17}^{(k−1)}, x_{10}^{(k)}, x_{22}^{(k−1)}, etc.).

CPS615J, G.C.Fox, Sept. 94  9  gcf@npac.syr.edu

6 × 6 Matrix—Wavefront Parallelism

for Gauss Seidel

• Here number l labels l = 1 … 11 phases

• Each phase l ONLY depends on \underline{x}^{(k)} components from (l − 1) th phase

• Can execute each update in a given phase concurrently

• Some Load Imbalance

• Speedup of 2.5 in example of three processors

The Red-Black Two Phase Parallel

Gauss Seidel Iteration

Red ⊙ Black •

• Note calculation of ⊙ only depends on • points and vice versa

• So two phase algorithm. Calculate all ⊙—Phase I. Then all •—Phase II.

Analysis of Parallel Red Black

Gauss Seidel

• This can be parallelized as easily as Jacobi!

• Phase I follows by Phase II.

• Red black is a different ordering from original Gauss Seidel.

• For (some) elliptic problems, can show red black as good as and maybe better than original ordering.

— i.e., red black does not need more iterations than pure Gauss Seidel.

• For related hyperbolic problems—see NAS benchmark SSOR discussion in CPS713—not so clear.

• For higher order differencing, need more phases (more colors) but principle the same.

Eigenvalues of Gauss Seidel

Iteration Matrix

• It was relatively easy to find the eigenvalues and eigenfunctions of Jacobi’s method for a simple geometry.

• Gauss-Seidel it is not so easy as matrix M^{−1}N is not symmetric. There is a nice discussion in

“Numerical Computation of Internal and External Flows,” C. Hirsch, Wiley, 1988.

• Eigenvalues λ_{GS} of Gauss-Seidel Iteration Matrix are related to those of Jacobi Iteration λ_{J} by λ_{GS} = λ_{J}^{2}

• We showed that

1 − ϵ < ∼ λ_{J} < ∼ −1 + ϵ ϵ ∼ \frac{π^{2}}{2N^{2}}

• Therefore,

|λ_{GS}|_{max} = (|λ_{J}|_{max})^{2} ∼ 1 − 2ϵ

Comparison of Convergence of

Gauss-Seidel and Jacobi Iteration

• Convergence can be characterized by a given fractional error (say .001).

• This correponds to

|λ|_{max}^{k} ∼ .001

where k number of required iterations

• exp [−k (1 − |λ|_{max})] = .001 i.e., value of k ∝ \frac{1}{(1−|λ|_{max})}

• We showed Gauss Seidel doubles value of

1/(1 − |λ|_{max}),

i.e., it converges in half as many iterations as Jacobi

Successive Overrelaxation Iteration

Method (SOR)

• This is relaxation (acceleration) method applied to Gauss-Seidel.

• In Jacobi, the optimal case was ω = 1 (i.e., no gain in using relaxation).

• However, for Gauss-Seidel, we will find ω ∼ 2 is best and does improve the method.

• The spectral radius (|λ|_{max}) of Jacobi is given by

ρ_{J} = 1 − \frac{π^{2}}{2N^{2}}

and for conventional (unaccelerated) Gauss Seidel by

ρ_{GS} = 1 − \frac{π^{2}}{N^{2}}

• For choice

ω_{opt} ∼ 2 (1 − \frac{π}{N} )

one can show that

ρ_{SOR} ∼ 1 − \frac{π}{N}

Convergence of SOR Compared to

Jacobi and Gauss Seidel

• Remember number of iterations is proportional to

\frac{1}{1 − ρ} = \frac{2N^{2}}{π^{2}} Jacobi
= \frac{N^{2}}{π^{2}} Gauss Seidel
= \frac{N}{π}  SOR

— i.e., number of iterations decreases in SOR by a factor N/π, which is large!

• Unfortunately, the above analysis is only valid in simple cases

• simple geometry

• some orderings in Gauss-Seidel

Estimate of Over Relaxation

Parameter ω

• Thus, in general for practice problems, one cannot precalculate optimal overrelaxation parameter.

• Find ω by numerical experimentation.

• However, note that one can expect a choice 1 < ω < 2 to be better than natural ω = 1.

• One can use the estimate

ω_{opt} = 2 − \frac{constant}{N}
= 2 − constant · h

(h is grid spacing)

Pseudo Code for SOR—Successive

Over Relaxation

Choose an initial guess x^{(0)} to the solution x. for k = 1,2,…

for i = 1,2,…,n σ = 0 for j = 1,2,…,i − 1

σ = σ + a_{i,j} x_{j}^{(k)} end for j = i + 1,…,n

σ = σ + a_{i,j} x_{j}^{(k−1)} end σ = (b_{i} − σ) /a_{i,i} x_{i}^{(k)} = x_{i}^{(k−1)} + ω (σ − x_{i}^{(k−1)})

end check convergence; continue if necessary

end

The SOR Method