16 r 1 = 1.6 = 1 + 2/16 + 3/162 + … 16 r 2 = 9.6 = 2 + 3/16 + 4/162 + … Hence, a number with m digits in one base may not have terminal representation in another base. It is not just irrationals that are a problem (e.g., consider ). Consider r = .115 if b = 10 and m = 2, then r = .11 chopping r = .12 symmetric rounding (r+.5bc-m-1 and then chop) Most computers chop instead of round off. IEEE compliant CPUs can do both and there may be a system call to switch, which is usually not user accessible. Note: When the rounding changes, almost all nontrivial codes break.
Google operates the world’s largest Cloud/GRID system. The Top 500 list provides an ongoing list of the fastest computers willing to be measured. It is not a comprehensive list and Google, Yahoo!, and many governments and companies do not participate. Data needs to be distributed sensibly among the p processors. Where the data needs to be can change, depending on the operation, and communication is usual. Algorithms that essentially never need to communicate are known as embarrassingly parallel. These algorithms scale wonderfully and are frequently used as examples of how well so and so’s parallel system scales. Most applications are not in this category, unfortunately. To do parallel programming, you need only a few functions to get by: Initialize the environment and find out processor numbers i. Finalize or end parallel processing on one or all processors. Send data to one, a set, or all processors.
Continuing basic definitions… If is an n-vector (i.e., a n1 matrix), then . Let ei be a n-vector with all zeroes except the ith component, which is 1. Then I = [ e1, e2, …, en ] is the nn identity matrix. Further, if A is nn, then IA=AI=A. The nn matrix A is said to be nonsingular if ! x such that Ax=b, b. Tests for nonsingularity:
Let 0n be the zero vector of length n. A is nonsingular if and only if 0n is the only solution of Ax=0n. A is nonsingular if and only if det(A)0. Lemma: ! A-1 such that A-1A=AA-1=I if and only if A is nonsingular. Proof: Suppose C such that CA-1, but CA=AC=I. Then C=IC=(A-1A)C=A-1(AC)=A-1I=A-1. Diagonal matrices: . Triangular matrices: upper , strictly upper . lower , strictly lower . 3b. Gaussian elimination
Theorem 3.1: Let A be such that Gaussian elimination yields nonzero diagonal elements , k=1, 2, …, n. Then A is nonsingular and . Also, is upper triangular and A has the factorization , where is lower triangular with elements The vector . Proof: Note that once is proven, , so follows.
Now we prove . Set . Then (since L and U are triangular and A(k) is satisfied for k=n) . From the definitions of we get and recall that . Thus, if ij, then . When i>j, . Finally, we prove . Let . So, .
Generalized Gaussian elimination Order of elimination arbitrary. Set . Select an arbitrary as the first pivot element. We can eliminate from all but the i1-st equation. The multipliers are . The reduced system is now . Select another pivot and repeat the elimination. If , then the remaining equations are degenerate and we halt. Theorem 3.2: Let A have rank r. Then we can find a sequence of distinct row and column indices (i1,j1), (i2,j2), …, (ir,jr) such that corresponding pivot elements in A(1), A(2), …, A(r) are nonzero and . Define permutation matrices (whose columns are unit vectors) and , where {ik} and {jk} are permutations of {1,2,…,n}. Then
Using the Ly=f and Ux=y approach, computing xi requires (n-i) multiplies plus 1 divide. Hence, only multiplies are required to solve the triangular systems. Lemma: operations are required to solve m systems , j=1, …, m by Gaussian elimination. Note: To compute A-1 requires n3 operations. In general, n2 operations are required to compute A-1f(j). Thus, to solve m systems requires mn2 operations. Hence, n3+mn2 operations are necessary to solve m systems. Thus, it is always more efficient to use Gaussian elimination instead of computing the inverse! We can always compute A-1 by solving Axi=ei, i=1,2,…,n and then the xi’s are the columns of A-1.
Theorem 3.3: If A is nonsingular, P such that PA=LU is possible and P is only a permutation of the rows. In fact, P may be found such that for i>k, Theorem 3.4: Suppose A is symmetric. If A=LU is possible, then the choice of lkk=ukklik=uki. Hence, U=LT. Variants of Gaussian elimination LDU factorization: L and U are strictly lower and upper triangular and D is diagonal. Cholesky: A=AT, so factor A=LLT. Fun example: is symmetric, but cannot be factored into LU form. Definition: A is positive definite if .
. Thus, . (lhs) Note that since , or . Further, . Combining the two inequalities gives us the lhs. QED Theorem 3.7: Suppose x and x satisfy , where x and x are perturbations. Let A be nonsingular and A be so small that . Then for we have . Note: Theorem 3.7 implies that when x is small, small relative changes in f and A cause small changes in x. Iterative Improvement
3c (i) Splitting or Relaxation Methods Let A=ST, where S is nonsingular. Then . Then the iterative procedure is defined by To be useful requires that is easy to compute. in a reasonable amount of time. Example: Let A=D-L-U, where D is diagonal and L and U are strictly lower and upper triangular, respectively. Then S=D and T=L+U: both are easy to compute, but many iterations are required in practice.
3d. Sparse Matrix Methods We want to solve Ax=b, where A is large, sparse, and NN. By sparse, A is nearly all zeroes. Consider the tridiagonal matrix, . If N=10,000, then A is sparse, but if N=4 it is not sparse. Typical sparse matrices are not just banded or diagonal matrices. The nonzero pattern may appear to be random at first glance. There are a small number of common storage schemes so that (almost) no zeroes are stored for A, ideally storing only NZ(A) = number of nonzeroes in A: Diagonal (or band) Profile
The compressed column storage scheme is defined similarly. Modified compressed row storage scheme (new Yale Sparse Matrix Package format) 2 vectors: IJA, AM, each of length NZ(A)+N+1. Assume A = D + L + U, where D is diagonal and L and U are strictly lower and upper triangular, respectively. Let . Then column index of jth element in AM is arbitrary The modified compressed column storage scheme is defined similarly.
Very modified compressed column storage scheme (Bank-Smith format) Assumes that A is either symmetric or nearly symmetric. Assume A = D + L + U, where D is diagonal and L and U are strictly lower and upper triangular, respectively. Let that will be stored. Let . 2 vectors: IJA, AM with both aij and aji stored if either is nonzero. row index of jth element in AM is arbitrary If , then AM contains first D, an arbitrary element, UT, and then possibly L.
Example: Then D and column “pointers” UT Optional L index 1 2 3 4 5 6 7 8 9 IJA 6 6 6 7 8 1 2 AM a11 a22 a33 a44 a13 a24 a31 a42 Compute Ax or ATx Procedure MULT( N, A, x, y ) do i = 1:N y(i)=A(i)x(i) enddo Lshift=0 if L is not stored and IJA(N+1)-IJA(1) otherwise
k=LIST(k) LIST(ksave)=0 Sum+=U(j)y(k) enddo // j y(i)=b(i)Sum LIST(k)=0 The i loop ends after this substitution. Solving Ux=z follows the same pattern, but columns are processed in the reverse order. Adding Lshift and Ushift parameters allows the same code handle both cases A=AT and AAT equally easily. R.E. Bank and R.K. Smith, General sparse elimination requires no permanent integer storage, SIAM J. Sci. Stat. Comp., 8 (1987), pp. 574-584 and the SMMP and Madpack2 packages in the Free Software section of my home web. 4. Solution of Nonlinear Equations
, which is a contradiction. QED Note that the condition on must be continuous. Algorithm: Let be arbitrary and set Note that after n steps, . Theorem 4.3: Let and . For , the sequence converges to the fixed point s and the nth error satisfies .
Note that Theorem 4.3 is a nonlocal convergence theorem because s is fixed, a known interval I is assumed, and convergence is for any . Proof: (convergence) Recall that s is unique. For any n, between xn-1 and s such that . Repeating this gives . Since , . Error bound: Note that Since , . QED
Let , I a closed interval, and g satisfies a local theorem’s requirements on I. The Taylor series of g about x=s is = = , where = . If , and , then kthorder convergence. The important k’s are 1 and 2. If we only have 1st order convergence, we can speed it up using quadratic interpolation: given , fit a 2nd order polynomial p to the data such that . Use p to get the next guess. Let
Properties: 1st order convergence Convergence if xn+1 can be found (always) Can obtain error estimates Newton’s Method Choose . Let s be such that . Then = = If exists in and , then 2nd order convergence. So, . What if exists? Then , where exists. So, =
, or (better) Solve Set Quadratic convergence if exists for nonsingular x0 sufficiently close to s 1D Example: . To reduce , Bisection: 20 steps False position: 7 steps Secant: 6 steps Newton: 5 steps Zeroes of Polynomials
Let . Properties: Computing and are easy. Finding zeroes when n4 is by use of formulas, e.g., when n=2, When n5, there are no formulas. Theorem 4.5 (Fundamental Theorem of Algebra): Given with n1, there exists at least one r (possibly complex) such that . We can uniquely factor p using Theorem 4.5: = , degree polynomial = , degree polynomial =
We can prove by induction that there exist no more than n roots. Suppose that such that . Then , where . Theorem 4.6 (Division Theorem): Let and be polynomials of degree n and m, where 1mn. Then of degree nm and a unique polynomial of degree m1 or less such that . Evaluating Polynomials How do we compute ? We may need to make a change of variables: , which leads to . Using Taylor’s Theorem we know that
Let be any linear interpolatory process that is exact for polynomials of degree m, i.e., if is a polynomial of degree m, . For a given function , Taylor’s Theorem says that . Define so that . Even More Error Analysis Define .
Quick review of probability Event B is a set of possible outcomes that has probability . The set of all events is denoted by and particular outcomes are . Hence, . Suppose . Then represents events in both B and C. Similarly, represents events that are in B or C. Some axioms of probability are The conditional probability that a C outcome is also a B outcome is given by Bayes formula, . Frequently, we already know both and and use Bayes formula to calculate .
The covariance matrix is positive semidefinite. Common random variables The standard uniform random variable U has a probability density of We can create a random variable in [a,b] by . The PDF for Y is The exponential random variable T with rate constant has a PDF