Numerical Linear Algebra

1
Numerical Linear Algebra

• Brian Hickey
• Jason Shields

2
What is Numerical Linear Algebra?

• Ax b
• The same algebra learned in high school.

3
Common terms used in NLA

• Scalars values described by magnitude alone.
• Vectors values described by both magnitude and
  direction

4
Common terms used in NLA (cont.)

• Matrix n x m 2-dimensional array of values

5
Types of Matrices

• Dense Matrices most values are filled

6
Types of Matrices (cont.)

• Band matrices nonzero only for certain
  diagonals
• Sparse matrices many zero entries
• a practical requirement for a family of
  matrices to be usefully' sparse is that they
  have only O(n) nonzero entries - CSEP

7
Common terms used in NLA (cont.)

• Eigenvalues a scalar c such that, given an n x
  n matrix A and a vector x, Ax cx, the vector x
  is known as the Eigenvector
• Least Squares Solution a vector x in
• Ax b that minimizes the length of the vector
  Ax b. A is an m x n matrix and b an m x 1
  vector.

8
Numerical Linear Algebra

• History

9
History of NLA

• Ancient Beginnings in Babylonian times (1900
  B.C.)
• Advanced base 60 mathematical system
• Used tables to calculate simple ax b formulas

10
Hermann Grassman

• 1809 - 1877

11
H. Grassman (cont.)

• Father of Linear Algebra
• Wrote Die lineale Ausdehnungslehre, ein neuer
  Zweig der Mathematik
• Introduced exterior algebra
• Symbols representing points, lines, etc. can be
  manipulated using certain rules.

12
H. Grassman (cont.)

• Never formally studied mathematics
• Taught at gymnasiums (German equivalent of high
  school) in Stettin and Berlin
• Also wrote papers on color vision, botany, tide
  theory, and acoustics
• Wrote a Sanskrit dictionary that is still widely
  used

13
Johann Carl Friedrich Gauss

• 1777 - 1855

14
JCF Gauss (cont.)

• Best known for Gaussian Elimination
• Used to simplify Ax b if A is a dense matrix
• Math prodigy
• Claimed to know the existence of non-Euclidean
  geometry at 15
• Contributed to astronomy

15
Arthur Caley

• 1821 - 1895

16
Arthur Caley (cont.)

• Theory of Algebraic Invariances
• Work with matrices
• Developed foundation for quantum mechanics
• Divided geometry into types
• Euclidean, non-Euclidean, n-dimensional

17
Cholesky Decomposition

• Used to simplify matrices that are symmetric and
  positive definite
• Used in parallel computing

18
Computational Science Matrices

• Pattern of non-zero elements within the matrix
  important for sparse matrix solve
• Reorder the matrix to optimize solve time
• Real-world problems usually exhibit a systematic
  pattern that can be exploited
• Sparse storage overhead vs. dense storage and
  work
• Reduce storage and work by manipulating zero
  elements

19
Sparse Matrices Fill

• New non-zero elements introduced by column
  modification
• Preserving sparsity requires that we limit fill
• Finding optimal row, column ordering to minimize
  fill is NP-complete combinatorial problem

20
(No Transcript)
21
(No Transcript)
22
Iterative Methods

• Use successive approximation to obtain more
  accurate solutions to a linear system at each
  step
• Two basic types, stationary and non-stationary
• Preconditioner transformation matrix used to
  improve convergence of the iteration method

23
Stationary Iterative Methods

• Perform the same operations on the current
  iteration vectors at each step
• Stationary simpler, easier to understand and
  implement, but usually not as effective
• Xk Bx(k-1) c
• where neither B nor c rely on the iteration
  count k

24
The Jacobi Method

• Examine each of the N equations in the system Ax
  b in isolation.
• The updates to the elements can be done
  simultaneously, sometimes called method of
  simultaneous displacements
• Defined by equation

25
The Gauss-Seidel Method

• Similar to Jacobi Method, except the equations
  are examined one at a time sequentially
• Previously computed results are used
• Also called method of successive displacements
  due to the iterates dependence on the ordering of
  equations

26
The Gauss-Seidel Method (Cont)

• Computations are serial, each new iterate depends
  on all previously computed components
• The iterates depend upon the order in which
  equations are examined
• Defined by the equation

27
SOR Method

• Successive Overrelaxation Method
• Applies extrapolation to Gauss-Seidel Method,
  uses weighted average w
• Choosing w to accelerate rate of convergence,
  this method is defined by the equation

28
Non-Stationary Iterative Methods

• Uses iteration-dependent coefficients
• Computation involves variable information that
  changes at each iteration
• Typically computed by taking inner products of
  residuals or other vectors arising from iterative
  method
• Harder to implement, or follow, than stationary
  counterparts

29
Conjugate Gradient Method

• Effective on symmetric positive definite systems
• Generates vector sequences of iterates, residuals
  corresponding to the iterates, and search
  directions used for updating the vector iterates
  and residuals
• Produces large sequences, but only needs to store
  a small number of vectors

30
Conjugate Gradient Method (Cont)

• Two inner products performed to compute update
  scalars that make the sequences satisfy certain
  orthogonality conditions
• These conditions minimize the distance to the
  solution in some norm (for SPD LS)

31
Conjugate Gradient Method (Cont)

• The iterates are updated in each iteration by a
  multiple of the search direction vector
• The residuals are updated also
• Where

32
Conjugate Gradient Method (Cont)

• The search direction is updated
• Where

33
MINRES Method

• Applied to symmetric indefinite systems
• Does not suffer breakdown upon encountering a
  zero pivot (like the Conjugate Gradient Method
  does)
• This method minimizes the residual in the 2-norm

34
Chebyshev Iteration Method

• Non-Stationery iteration method for solving
  non-symmetric problems
• Does not involve computing inner products
• Requires enough knowledge about the spectrum of
  the coefficient matrix A to identify an ellipse
  that envelopes the spectrum and excludes the
  origin

35
Time-Consuming Computational Kernels of Iterative
Schemes

• Inner Products
• Vector Updates
• Matrix Vector Products
• Preconditioning (solve for w in Kw r)
• Preconditioners are transformation matrices that
  improve spectral properties, and hence
  convergence rates

36
IML

• Iterative Methods Library
• C template library of modern iterative methods
• Solves symmetric and non-symmetric systems of
  linear equations
• http//math.nist.gov/iml/

37
Parallelizing the Kernels

• Each processor computes inner product of two
  segments of each vector
• Each processor updates its vector segment
• Split matrix into vector segment strips for
  matrix-vector products (shared memory), or into
  connected nearby node blocks (distributed memory)

38
Parallelism in Preconditioning

• Preconditioning most difficult
• Reordering the computations
• Reordering the unknowns
• Forced parallelism
• Other preconditioners simple diagonal scaling,
  polynomial preconditioning, and truncated Neumann
  series

39
How is NLA used in Computational Sciences?

• Why computers brought upon an interest in NLA
• Parallel computers and there use of matrices

40
The resurgence of NLA from computers

• Before computers, NLA did not get a whole lot of
  attention
• The raw computing power of computers made NLA a
  much more practical tool

41
The use of NLA in parallel computing

• Because of the nature of parallel computing, it
  is very beneficial to use matrices in NLA
• Matrices are easily split up and distributed
  throughout a network of processors

42
Several Numerical Software Libraries for using
NLA in Parallel Computing

• BLAS
• LINPACK
• LAPACK
• SCALAPACK
• TNT (?)

43
BLAS

• Library of basic Linear Algebra Subroutines
• Such as
• matrix vector
• matrix scalar

44
LINPACK

• System Library of widely used Linear Algebra
  Routines
• Such as Gaussian Elimination and Cholesky Decomp.
• Initially written in Fortran 66.
• Now written in Fortran 77
• Written before the advent of vector and parallel
  computers

45
LAPACK

• Linear Algebra Package
• Similar to LINPACK
• Written with parallel architectures such as Cray
  in mind.
• Written in several languages including
• Fortran 77
• with calls to BLAS

46
ATLAS

• Automatically Tuned Linear Algebra Software
• Both a research project and a software package
• Purpose is to provide portably optimal linear
  algebra software
• Optimizes BLAS and a small subset of LAPACK for
  optimal performance for each individual machine
• ATLAS homepage http//math-atlas.sourceforge.net/
  

47
SCALAPACK

• Scalable LAPACK
• Continuation of the LAPACK project
• a subset of LAPACK routines redesigned for
  distributed memory MIMD parallel computers
• currently written in a Single-Program-Multiple-Dat
  a style using explicit message passing for
  interprocessor communication
• Taken from ScalaPACK homepage http//www.netlib.o
  rg/scalapack/scalapack_home.html

48
TNT

• Template Numerical Toolkit
• According to several websites, this is replacing
  SCALAPACK as the standard software library for
  NLA in parallel computing
• Not necessarily reliable
• TNT homepage http//math.nist.gov/tnt/overview.ht
  ml

49
TNT (cont.)

• TNT Data Structures
• C-style arrays
• Fortran-style arrays
• Sparse Matrices
• Vector/Matrix
• TNT Utilities
• array I/O
• math routines (hypot(), sign(), etc.)
• Stopwatch class for timing measurements
• Libraries that utilize TNT
• JAMA a linear algebra library with QR, SVD,
  Choleksy and Eigenvector solvers.
• old (pre 1.0) TNT routines for LU, QR, and
  Eigenvalue problems

50
Sources and Links

• http//csep1.phy.ornl.gov/la/node14.html
• http//www.netlib.org/linalg/html_templates/node9.
  html
• http//www.cise.ufl.edu/davis/sparse/
• http//ipdps.eece.unm.edu/2000/papers/wylin.pdf
• http//math.nist.gov/iml/

51
Sources and Links

• http//www-groups.dcs.st-and.ac.uk/history/Mathem
  aticians/
• http//math.nist.gov/tnt/overview.html
• http//www.netlib.org/scalapack/scalapack_home.htm
  l
• http//www.mgnet.org/douglas/Classes/cs521-s01/in
  dex.html