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

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 X^k = 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-Data style using explicit message passing for interprocessor communication Taken from ScalaPACK homepage: http://www.netlib.org/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.html

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/Mathematicians/ http://math.nist.gov/tnt/overview.html http://www.netlib.org/scalapack/scalapack_home.html http://www.mgnet.org/~douglas/Classes/cs521-s01/index.html