Understanding Linear Systems and LU Decomposition
Explore the fundamental concepts of linear algebra, including matrix notation, existence of solutions, vector spaces, computation tasks, and LU decomposition techniques. Learn about Gauss elimination, Crout's algorithm, and how to solve linear systems efficiently using LU decomposition.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
E N D
Presentation Transcript
Chapter 2, Linear Systems, Mainly LU Decomposition
Linear Systems In matrix notation: A x = b See, e.g., Matrix Analysis and Applied Linear Algebra , Carl D. Meyer, SIAM, 2000. Or Matrix Computations , Golub and van Loan.
Existence and Solutions M = N, for nonsingular square matrix, a unique solution exists M < N, more variables than equations infinitely many solutions, need to find the null space of A, Ax = 0 (e.g. SVD) M > N, one can ask for a least-square solution, e.g., from the normal equation (ATA) x = (AT b)
Some Concepts in Linear Algebra Vector space, linear independence, and dimension Null space of a matrix A: the set of all x such that Ax = 0 Range of A: the set of all Ax Rank of A: max number of linearly independent rows or columns Rank of A + Null space Dim of A = number of columns of A
Computation Tasks and Pitfalls Solve A x = b Find A-1 Compute det(A) Cramer s rule ??=det(??) det(?) Round off error makes the system singular, or numerical instability makes the answer wrong.
Gauss Elimination Basic facts about linear equations Interchanging any two rows of A and b does not change the solution x Replace a row by a linear combination of itself and any other row does not change x Interchange column permutes the solution Example of Gauss elimination Pivoting LU equivalent to Gauss elimination
LU Decomposition 1 = ii L . U = A Thus A x = (LU) x = L (U x) = b Let y = U x, then solve y in L y = b by forward substitution Solve x in U x = y by backward substitution
Crouts Algorithm = 1 Set For each j= 1, 2, 3, , N (a) for i= 1, 2, , j a ij ij = for all i aij for A ii for L for U 1 i A = LU ik kj L y = b U x = y Ax = LUx = b = 1 k (b) for i = j +1, j +2, , N 1 jj 1 j k = a ij ij ik kj = 1
Order of Update in Crouts Algorithm ???= 1
LU Decomposition (without pivoting) 2 1 3 0 5 1 0 1 0 0 1 11 0 0 12 13 = = = 2 1 2 A LU 21 22 0 23 1 31 32 33 Determining the entries in an order scanning the column from top to bottom and left to right. L and U are stored in the same place as the original A. 1 0 1 1/8 1 0 0 2 0 0 1 4 0 5 4 = 1 1/ 2 4
Pivoting Pivoting is essential for stability Interchange rows of A to get largest ii, this results LU = P A. Implicit pivoting (when comparing for the biggest elements in a column, use the normalized one so that the largest coefficient in an equation is 1)
ludcmp() import math def ludcmp(a, n, indx, d): d[0] = 1.0 vv=indx.copy() for i in range(0,n): big = 0.0 for j in range(0,n): temp = math.fabs(a[i][j]) if (temp > big): big = temp vv[i]=1.0/big Indentations are important in Python
Run Crouts algorithm for j in range(0,n): big = 0.0 for i in range(0,n): if(i<j): l = i else: l = j sum = a[i][j] for k in range(0,l): sum -= a[i][k]*a[k][j] a[i][j] = sum
if(i>=j): dum = vv[i]*math.fabs(sum) if (dum >= big): big = dum imax = i if (j != imax): dum = a[imax] a[imax] = a[j] a[j] = dum d[0] = - d[0] vv[imax] = vv[j] indx[j] = imax dum = 1.0/a[j][j] for i in range(j+1,n): a[i][j] *= dum Need very carefully pay attention to indentation!
Computational Complexity How many basic steps does the LU decomposition program take? Big O( ) notation & asymptotic performance O(N3)
Compute A-1 Let B = A-1 then AB = I (I is identity matrix) or A [b1, b2, ,bN] = [e1, e2, ,eN] or Abj = ej for j= 1, 2, , N where bj is the j-th column of B. I.e., to compute A-1, we solve a linear system N times, each with a unit vector ej.
Compute det(A) Definition of determinant = ( 1) P det( ) A a a a 1, 2, , 1 i i n i 2 n P Properties of determinant (anti-symmetry, Laplace expansion, etc) Since det(LU) = det(L) det(U), thus det(A) = det(U) = 1 j N = jj
Use numpy and scipy numpy: define N-dimensional arrays of various types scipy: efficient numerical routines such as integration, interpolation, linear algebra, and statistics covered in this course. scipy need numpy in most cases.
Solve Ax =b in scipy import numpy as np from scipy import linalg # A is 2x2 matrix A = np.array([[1.03,2.0],[3.5,4.1]]) # b is the right-side vector b = np.array([[5.0],[6.1]]) # x is solution x = np.linalg.solve(A,b)
Use LAPACK Lapack is a free, high quality linear algebra solver package (downloadable at www.netlib.org/lapack/). Much more sophisticated than NR routines. Written in Fortran 90 Calling inside Python made possible through scipy
What Lapack can do? Solution of linear systems, Ax = b Least-square problem, min ||Ax-b||2 Singular value decomposition, A = UsVT Eigenvalue problems, Ax = x
An Example for Lapack Fortran Routine SUBROUTINE DSYEVR(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO) DSYEVR computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues. JOBZ (input) CHARACTER*1 = 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors. RANGE (input) CHARACTER*1 = 'A': all eigenvalues will be found. = 'V': all eigenvalues in the half-open interval (VL,VU] will be found. = 'I': the IL-th through IU-th eigenvalues will be found. For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and DSTEIN are called UPLO (input) CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored. N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if UPLO='L') or the upper triangle (if UPLO='U') of A, including the diagonal, is destroyed. skip
Other matrix decompositions Cholesky: A = LLT, A is real symmetric, positive definite. QR factorization: A=QR, Q is orthogonal QTQ=I, R is upper triangular. Schur: U AU= upper triangular with eigenvalues on the diagonals. Singular value decomposition: A=U V , U and V are unitary U U=I, V V=I , is nonnegative diagonal.
Reading, Reference Read NR, Chap. 2 M. T. Heath, Scientific Computing: an introductory survey For a more thorough treatment on numerical linear algebra computation problems, see G H Golub & C F van Loan, Matrix Computations , as well as Horn & Johnson, Matrix Analysis.
Problems for week 3/chapter 2 (due Thu 12 Sep 2024) 1. Consider the following linear equation 1 2 5 0 1 2 x x x 1 = 1 2 3 6 0 1 2 3 (a) Solve the system (by hand) with LU decomposition, following Crout s algorithm (slide 8, NR p.44) exactly, without pivoting. (b) Find the inverse of the 3 3 matrix, using LU decomposition. (c) Find the determinant of the 3 3 matrix, using LU decomposition. 2. Can we do LU decomposition for all square, real matrices? What is the condition for the existence of an LU decomposition? The answers depend if pivoting is applied or not.
