• 沒有找到結果。

# Gaussian Elimination for Linear Systems

N/A
N/A
Protected

Share "Gaussian Elimination for Linear Systems"

Copied!
56
0
0

(1)

university-logo

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University

October 3, 2011

(2)

university-logo

## Outline

1 Elementary matrices

2 LR-factorization

3 Gaussian elimination

4 Cholesky factorization

5 Error estimation for linear systems

(3)

university-logo

## Elementary matrices

Let A ∈ Cn×n be a nonsingular matrix. We want to solve the linear system Ax = b by

(a) Direct methods (finite steps);

(b) Iterative methods (convergence). (See Chapter 4)

(4)

university-logo

A =



1 1 0 3

2 1 −1 1

3 −1 −1 2

−1 2 3 −1



⇒ A1:= L1A ≡



1 0 0 0

−2 1 0 0

−3 0 1 0

1 0 0 1



 A =



1 1 0 3

0 −1 −1 −5 0 −4 −1 −7

0 3 3 2



⇒ A2:= L2A1



1 0 0 0

0 1 0 0

0 −4 1 0

0 3 0 1



 A1=



1 1 0 3

0 −1 −1 −5

0 0 3 13

0 0 0 −13



= L2L1A

(5)

university-logo

We have

A = L−11 L−12 A2= LR.

where L and R are lower and upper triangular, respectively.

Question

How to compute L−11 and L−12 ?

L1=



1 0 0 0

−2 1 0 0

−3 0 1 0

1 0 0 1



 = I −



 0 2 3

−1



 1 0 0 0 

L2=



1 0 0 0

0 1 0 0

0 −4 1 0

0 3 0 1



 = I −



 0 0 4

−3



 0 1 0 0 

(6)

university-logo

Definition 1

A matrix of the form

I − αxy (α ∈ F, x, y ∈ Fn) is called an elementary matrix.

The eigenvalues of (I − αxy) are {1, 1, . . . , 1, 1 − αyx}. Compute (I − αxy)(I − βxy) = I − (α + β − αβyx)xy.

If αyx − 1 6= 0 and let β =αyαx−1, then α + β − αβyx = 0. We have (I − αxy)−1 = (I − βxy),

where α1+β1 = yx.

(7)

university-logo

Example 1

Let x ∈ Fn, and xx = 1. Let H = {z : zx = 0} and Q = I − 2xx (Q = Q, Q−1= Q).

Then Q reflects each vector with respect to the hyperplane H. Let y = αx + w, w ∈ H. Then, we have

Qy = αQx + Qw = −αx + w − 2(xw)x = −αx + w.

(8)

university-logo

Let y = ei to be the i-th column of the unit matrix and x = li = [0, · · · , 0, li+1,i, · · · , ln,i]T. Then,

I + lieTi =









 1

. ..

1 li+1,i

... . ..

ln,i 1









(1)

Since eTi li= 0, we have

(I + lieTi)−1= (I − lieTi ).

(9)

university-logo

From the equality

(I + l1eT1)(I + l2eT2) = I + l1eT1 + l2eT2 + l1(eT1l2)eT2 = I + l1eT1 + l2eT2 follows that

(I + l1eT1) · · · (I + lieTi) · · · (I + ln−1eTn−1)

= I + l1eT1 + l2eT2 + · · · + ln−1eTn−1

=

 1

l21 . .. 0 ... . .. . .. ln1 · · · ln,n−1 1

. (2)

Theorem 2

A lower triangular with “1” on the diagonal can be written as the product of n − 1 elementary matrices of the form (1).

Remark: (I + l1eT1 + · · · + ln−1eTn−1)−1 = (I − ln−1eTn−1) · · · (I − l1eT1) which can not be simplified as in (2).

(10)

university-logo

## LR -factorization

Definition 3

Given A ∈ Cn×n, a lower triangular matrix L with “1” on the diagonal and an upper triangular matrix R. If A = LR, then the product LR is called a LR-factorization (or LR-decomposition) of A.

(11)

university-logo

Basic problem

Given b 6= 0, b ∈ Fn. Find a vector l1= [0, l21, · · · , ln1]T and c ∈ F such that

(I − l1eT1)b = ce1.

### Solution:

 b1= c,

bi− li1b1= 0, i = 2, . . . , n.

 b1= 0, it has no solution (since b 6= 0), b16= 0, then c = b1, li1= bi/b1, i = 2, . . . , n.

(12)

university-logo

### Construction of LR-factorization:

Let A = A(0)=h

a(0)1 · · · a(0)n

i. Apply basic problem to a(0)1 : If

a(0)11 6= 0, then there exists L1= I − l1eT1 such that (I − l1eT1)a(0)1 = a(0)11e1. Thus

A(1) = L1A(0)

= h

La(0)1 · · · La(0)n

i

=





a(0)11 a(0)12 · · · a(0)1n 0 a(1)22 a(1)2n ... ... ... 0 a(1)n2 · · · a(1)nn





 .

(13)

university-logo

The k-th step:

A(k)= LkA(k−1)= LkLk−1· · · L1A(0) (3)

=















a(0)11 · · · a(0)1n 0 a(1)22 · · · a(1)2n

... 0 . .. ...

... ... . .. a(k−1)kk · · · a(k−1)kn ... ... 0 a(k)k+1,k+1 · · · a(k)k+1,n

... ... ... ... ...

0 0 · · · 0 a(k)n,k+1 · · · a(k)nn















(14)

university-logo

If a(k−1)kk 6= 0, for k = 1, . . . , n − 1, then the method is executable and we have that

A(n−1)= Ln−1· · · L1A(0)= R is an upper triangular matrix. Thus, A = LR.

Explicit representation of L:

Lk= I − lkeTk, L−1k = I + lkeTk

L = L−11 · · · L−1n−1= (I + l1eT1) · · · (I + ln−1eTn−1)

= I + l1eT1 + · · · + ln−1eTn−1 (by (2)).

(15)

university-logo

Theorem 4

Let A be nonsingular. Then A has an LR-factorization (A = LR) if and only if κi:= det(Ai) 6= 0, where Ai is the leading principal matrix of A, i.e.,

Ai=

a11 · · · a1i

. . .

. . . ai1 · · · aii

, for i = 1, . . . , n − 1.

Proof: (Necessity “⇒” ): Since A = LR, we have



a11 · · · a1i

... ... ai1 · · · aii

 =



l11 0

... . ..

li1 · · · lii





r11 · · · r1i

. .. ...

0 rii

 .

From det(A) 6= 0 follows that det(L) 6= 0 and det(R) 6= 0. Thus, ljj 6= 0 and rjj 6= 0, for j = 1, . . . , n. Hence κi= l11· · · liir11· · · rii6= 0.

(16)

university-logo

(Sufficiency “⇐”): From (3) we have

A(0)= (L−11 · · · L−1i )A(i).

Consider the (i + 1)-th leading principle determinant. From (3) we have

a11 · · · ai,i+1

... ...

ai+1 · · · ai+1,i+1

=

1 0

l21 . .. ... . .. . .. ... . .. . .. li+1,1 · · · li+1,i 1

a(0)11 a(0)12 · · · a(0)1,i+1 a(1)22 · · · a(1)2,i+1

. .. ...

a(i−1)ii a(i−1)i,i+1

0 a(i)i+1,i+1

 .

Thus, κi= 1 · a(0)11a(1)22 · · · a(i)i+1,i+16= 0 which implies a(i)i+1,i+16= 0.

Therefore, the LR-factorization of A exists.

(17)

university-logo

Theorem 5

If a nonsingular matrix A has an LR-factorization with A = LR and l11= · · · = lnn= 1, then the factorization is unique.

Proof: Let A = L1R1= L2R2. Then L−12 L1= R2R−11 = I.

Corollary 6

If a nonsingular matrix A has an LR-factorization with A = LDR, where D is diagonal, L and RT are unit lower triangular (with one on the diagonal) if and only if κi6= 0.

(18)

university-logo

Theorem 7

Let A be a nonsingular matrix. Then there exists a permutation P , such that P A has an LR-factorization.

Proof: By construction! Consider (3): There is a permutation Pk, which interchanges the k-th row with a row of index large than k, such that 0 6= a(k−1)kk (∈ PkA(k−1)). This procedure is executable, for k = 1, . . . , n − 1. So we have

Ln−1Pn−1· · · LkPk· · · L1P1A(0) = R. (4) Let P be a permutation which affects only elements k + 1, . . . , n. It holds P (I − lkeTk)P−1 = I − (P lk)eTk = I − ˜lkeTk = ˜Lk, (eTkP−1= eTk) where ˜Lk is lower triangular. Hence we have

P Lk= ˜LkP.

(19)

university-logo

Now write all Pk in (4) to the right as

Ln−1n−2· · · ˜L1Pn−1· · · P1A(0)= R.

Then we have P A = LR with L−1= Ln−1n−2· · · ˜L1and P = Pn−1· · · P1.

(20)

university-logo

## Gaussian elimination

Given a linear system

Ax = b

with A nonsingular. We first assume that A has an LR-factorization, i.e., A = LR. Thus

LRx = b.

We then (i) solve Ly = b; (ii) solve Rx = y. These imply that LRx = Ly = b. From (4), we have

Ln−1· · · L2L1(A | b) = (R | L−1b).

(21)

university-logo

Algorithm: Gaussian elimination without permutation

1: for k = 1, . . . , n − 1 do

2: if akk= 0 then

3: Stop.

4: else

5: ωj := akj (j = k + 1, . . . , n);

6: end if

7: for i = k + 1, . . . , n do

8: η := aik/akk, aik := η;

9: for j = k + 1, . . . , n do

10: aij:= aij− ηωj, bj := bj− ηbk.

11: end for

12: end for

13: end for

14: xn = bn/ann;

15: for i = n − 1, n − 2, . . . , 1 do

16: xi = (bi−Pn

j=i+1aijxj)/aii.

17: end for

(22)

university-logo

Cost of computation (A flop is a floating point operation):

(i) LR-factorization: 2n3/3 flops;

(ii) Computation of y: n2 flops;

(iii) Computation of x: n2flops.

For A−1: 8/3n3≈ 2n3/3 + 2kn2(k = n linear systems).

(23)

university-logo

Pivoting: (a) Partial pivoting; (b) Complete pivoting.

From (3), we have

A(k−1)=













a(0)11 · · · a(0)1n

0 . .. ...

... a(k−2)k−1,k−1 · · · a(k−2)k−1,n ... 0 a(k−1)kk · · · a(k−1)kn

... ... ... ...

0 . . . 0 a(k−1)nk · · · a(k−1)nn











 .

(24)

university-logo

## Partial pivoting



Find a p ∈ {k, . . . , n} such that

|apk| = maxk≤i≤n|aik| (rk = p)

swap akj with apj for j = k, . . . , n, and bk with bp.

(5)

Replacing stopping step in Line 3 of Gaussian elimination Algorithm by (5), we have a new factorization of A with partial pivoting, i.e., P A = LR (by Theorem 7 and |lij| ≤ 1 for i, j = 1, . . . , n).

For solving linear system Ax = b, we use

P Ax = P b ⇒ L(Rx) = PTb ≡ ˜b.

It needs extra n(n − 1)/2 comparisons.

(25)

university-logo

## Complete pivoting









Find p, q ∈ {k, . . . , n} such that

|apq| ≤ max

k≤i,j≤n|aij|, (rk:= p, ck := q)

swap akj and bk with apj and bp, resp., (j = k, . . . , n), swap aik with aiq(i = 1, . . . , n).

(6)

Replacing stopping step in Line 3 of Gaussian elimination Algorithm by (6), we also have a new factorization of A with complete pivoting, i.e., P AΠ = LR (by Theorem 7 and |lij| ≤ 1, for i, j = 1, . . . , n).

For solving linear system Ax = b, we use

P AΠ(ΠTx) = P b ⇒ LR˜x = ˜b ⇒ x = Π˜x.

It needs n3/3 comparisons.

(26)

university-logo

Let

A =

 10−4 1

1 1



be in three decimal-digit floating point arithmetic.

κ(A) = kAkkA−1k≈ 4. A is well-conditioned.

Without pivoting:

L =

 1 0

f l(1/10−4) 1



, f l(1/10−4) = 104, R =

 10−4 1

0 f l(1 − 104· 1)



, f l(1 − 104· 1) = −104. LR =

 1 0

104 1

  10−4 1 0 −104



=

 10−4 1

1 0

 6= A.

(27)

university-logo

Here a22 entirely “lost” from computation. It is numerically unstable.

Let Ax =

 1 2



. Then x ≈

 1 1

 . But Ly =

 1 2



solves y1= 1 and y2= f l(2 − 104· 1) = −104, Rˆx = y solves ˆx2= f l((−104)/(−104)) = 1,

ˆ

x1= f l((1 − 1)/10−4) = 0.

We have an erroneous solution with cond(L), cond(R) ≈ 108.

(28)

university-logo

### Partial pivoting:

L =

 1 0

f l(10−4/1) 1



=

 1 0

10−4 1

 ,

R =

 1 1

0 f l(1 − 10−4)



=

 1 1 0 1

 . L and R are both well-conditioned.

(29)

university-logo

T

## -factorizations

### Algorithm 2

[Crout’s factorization or compact method]

For k = 1, . . . , n,

for p = 1, 2, . . . , k − 1, rp:= dpapk, ωp:= akpdp, dk:= akk−Pk−1

p=1akprp, if dk = 0, then stop; else

for i = k + 1, . . . , n, aik := (aik−Pk−1

p=1aiprp)/dk, aki:= (aki−Pk−1

p=1ωpapi)/dk. Cost: n3/3 flops.

With partial pivoting: see Wilkinson EVP pp.225-.

Advantage: One can use double precision for inner product.

(30)

university-logo

Theorem 8

If A is nonsingular, real and symmetric, then A has a unique LDLT factorization, where D is diagonal and L is a unit lower triangular matrix (with one on the diagonal).

Proof: A = LDR = AT = RTDLT. It implies L = RT.

(31)

university-logo

Theorem 9

If A is symmetric and positive definite, then there exists a lower triangular G ∈ Rn×n with positive diagonal elements such that A = GGT.

Proof:

A is symmetric positive definite

⇔ xTAx ≥ 0 for all nonzero vector x ∈ Rn

⇔ κi≥ 0 for i = 1, . . . , n

⇔ all eigenvalues of A are positive

From Corollary 6 and Theorem 8 we have A = LDLT. From L−1AL−T = D follows that

dk= (eTkL−1)A(L−Tek) > 0.

Thus, G = Ldiag{d1/21 , · · · , d1/2n } is real, and then A = GGT.

(32)

university-logo

Derive an algorithm for computing the Cholesky factorization A = GGT: Let

A≡ [aij] and G =

g11 0 · · · 0 g21 g22 . .. ... ... ... . .. 0 gn1 gn2 · · · gnn

 .

Assume the first k − 1 columns of G have been determined after k − 1 steps.

Bycomponentwise comparisonwith

[aij] =

g11 0 · · · 0 g21 g22 . .. ... ... ... . .. 0 gn1 gn2 · · · gnn

g11 g21 · · · gn1

0 g22 · · · gn2

... . .. . .. ... 0 · · · 0 gnn

 ,

one has

akk=

k

X

j=1

gkj2 ,

(33)

university-logo

which gives

gkk2 = akk

k−1X

j=1

gkj2 .

Moreover,

aik= Xk j=1

gijgkj, i = k + 1, . . . , n,

hence the k-th column of G can be computed by

gik =

aik

k−1X

j=1

gijgkj

gkk, i = k + 1, . . . , n.

(34)

university-logo

Cholesky Factorization

Input: n × n symmetric positive definite matrix A.

Output: Cholesky factorization A = GGT.

1: Initialize G = 0.

2: for k = 1, . . . , n do

3: G(k, k) =q

A(k, k) −Pk−1

j=1G(k, j)G(k, j)

4: for i = k + 1, . . . , n do

5: G(i, k) =

A(i, k) −Pk−1

j=1G(i, j)G(k, j)  G(k, k)

6: end for

7: end for

In addition to n square root operations, there are approximately Xn

k=1

[2k − 2 + (2k − 1)(n − k)] = 1 3n3+1

2n2−5 6n

floating-point arithmetic required by the algorithm.

(35)

university-logo

For solving symmetric, indefinite systems: See Golub/ Van Loan Matrix Computation pp. 159-168. P APT = LDLT, D is 1 × 1 or 2 × 2 block-diagonal matrix, P is a permutation and L is lower triangular with one on the diagonal.

(36)

university-logo

## Error estimation for linear systems

Consider the linear system

Ax = b, (7)

and the perturbed linear system

(A + δA)(x + δx) = b + δb, (8)

where δA and δb are errors of measure or round-off in factorization.

Definition 10

Let k · k be an operator norm and A be nonsingular. Then

κ ≡ κ(A) = kAkkA−1k is a condition number of A corresponding to k k.

(37)

university-logo

Theorem 11 (Forward error bound)

Let x be the solution of (7) and x + δx be the solution of the perturbed linear system (8). If kδAkkA−1k < 1, then

kδxk

kxk ≤ κ

1 − κkδAkkAk

kδAk

kAk +kδbk kbk

 .

Proof: From (8) we have

(A + δA)δx + Ax + δAx = b + δb.

Thus,

δx = −(A + δA)−1[(δA)x − δb]. (9)

(38)

university-logo

Here, Corollary 2.7 implies that (A + δA)−1 exists. Now, k(A + δA)−1k = k(I + A−1δA)−1A−1k ≤ kA−1k 1

1 − kA−1kkδAk. On the other hand, b = Ax implies kbk ≤ kAkkxk. So,

1

kxk ≤kAk

kbk. (10)

From (9) follows that kδxk ≤ 1−kAkA−1−1kkδAkk (kδAkkxk + kδbk). By using (10), the inequality (11) is proved.

Remark 1

If κ(A) is large, then A (for the linear system Ax = b) is called ill-conditioned, else well-conditioned.

(39)

university-logo

## Error analysis for Gaussian algorithm

A computer in characterized by four integers:

(a) the machine base β;

(b) the precision t;

(c) the underflow limit L;

(d) the overflow limit U .

Define the set of floating point numbers.

F = {f = ±0.d1d2· · · dt× βe| 0 ≤ di< β, d16= 0, L ≤ e ≤ U } ∪ {0}.

Let G = {x ∈ R | m ≤ |x| ≤ M } ∪ {0}, where m = βL−1and

M = βU(1 − β−t) are the minimal and maximal numbers of F \ {0} in absolute value, respectively.

(40)

university-logo

We define an operator f l : G → F by

f l(x) = the nearest c ∈ F to x by rounding arithmetic.

One can show that f l satisfies

f l(x) = x(1 + ε), |ε| ≤ eps,

where eps = 12β1−t. (If β = 2, then eps = 2−t). It follows that f l(a ◦ b) = (a ◦ b)(1 + ε)

or

f l(a ◦ b) = (a ◦ b)/(1 + ε), where |ε| ≤ eps and ◦ = +, −, ×, /.

(41)

university-logo

Given x, y ∈ Rn. The following algorithm computes xTy and stores the result in s.

s = 0,

for k = 1, . . . , n, s = s + xkyk.

Theorem 12

If n2−t≤ 0.01, then

f l(

Xn k=1

xkyk) = Xn k=1

xkyk[1 + 1.01(n + 2 − k)θk2−t], |θk| ≤ 1

Proof of Theorem 12

(42)

university-logo

Let the exact LR-factorization of A be L and R (A = LR) and let ˜L, ˜R be the LR-factorization of A by using Gaussian Algorithm (without pivoting). There are two possibilities:

(i) Forward error analysis: Estimate |L − ˜L| and |R − ˜R|.

(ii) Backward error analysis: Let ˜L ˜R be the exact LR-factorization of a perturbed matrix ˜A = A + E. Then E will be estimated, i.e.,

|E| ≤ ?.

Theorem 13

The LR-factorization ˜L and ˜R of A using Gaussian Elimination with partial pivoting satisfies

L ˜˜R = A + E, (2.6)

where

kEk≤ n2ρkAk2−t (2.7)

with

ρ = max

i,j,k

a(k)ij

.

kAk. (2.8)

(43)

university-logo

Applying Theorem 12 to the linear system ˜Ly = b and ˜Rx = y, respectively, the solution x satisfies

( ˜L + δ ˜L)( ˜R + δ ˜R)x = b or

( ˜L ˜R + (δ ˜L) ˜R + ˜L(δ ˜R) + (δ ˜L)(δ ˜R))x = b. (2.9) Since ˜L ˜R = A + E, substituting this equation into (2.9) we get

[A + E + (δ ˜L) ˜R + ˜L(δ ˜R) + (δ ˜L)(δ ˜R)]x = b.

The entries of ˜L and ˜R satisfy

|elij| ≤ 1, and |erij| ≤ ρkAk.

(44)

university-logo

Therefore, we get





















k ˜Lk≤ n, k ˜Rk≤ nρkAk,

kδ ˜Lkn(n+1)2 1.01 · 2−t, kδ ˜Rkn(n+1)2 1.01ρ2−t.

(2.10)

In practical implementation we usually have n22−t<< 1. So it holds kδ ˜Lkkδ ˜Rk≤ n2ρkAk2−t.

Let

δA = E + (δ ˜L) ˜R + ˜L(δ ˜R) + (δ ˜L)(δ ˜R). (2.11) Then, from (2.7) and (2.10) we get

kδAk≤ kEk+ kδ ˜Lkk ˜Rk+ k ˜Lkkδ ˜Rk+ kδ ˜Lkkδ ˜Rk

≤ 1.01(n3+ 3n2)ρkAk2−t (2.12)

(45)

university-logo

Theorem 14

For a linear system Ax = b the solution x computed by Gaussian Elimination with partial pivoting is the exact solution of the equation (A + δA)x = b and δA satisfies (2.12).

Remark: The quantity ρ defined by (2.9) is called a growth factor. The growth factor measures how large the numbers become during the process of elimination. In practice, ρ is usually of order 10 for partial pivot selection. But it can be as large as ρ = 2n−1, when

A =

1 0 · · · 0 1

−1 1 0 · · · 0 1

... −1 . .. . .. ... 1 ... ... . .. . .. 0 1

−1 −1 · · · −1 1 1

−1 −1 · · · −1 1

 .

(46)

university-logo

Better estimates hold for special types of matrices. For example in the case of upper Hessenberg matrices, that is, matrices of the form

A =





× · · · ×

× . .. ... ...

. .. ... ...

0 × ×





the bound ρ ≤ (n − 1) can be shown. (Hessenberg matrices arise in eigenvalus problems.)

(47)

university-logo

For tridiagonal matrices

A =









α1 β2 0

γ2 . .. ...

. .. ... ...

. .. ... βn

0 γn αn









it can even be shown that ρ ≤ 2 holds for partial pivot selection. Hence, Gaussian elimination is quite numerically stable in this case.

(48)

university-logo

For complete pivot selection, Wilkinson (1965) has shown that

|akij| ≤ f (k) max

i,j |aij| with the function

f (k) := k21h

21312413 · · · k(k−1)1 i12 . This function grows relatively slowly with k:

k 10 20 50 100

f (k) 19 67 530 3300

(49)

university-logo

Even this estimate is too pessimistic in practice. Up until now, no matrix has been found which fails to satisfy

|a(k)ij | ≤ (k + 1) max

i,j |aij| k = 1, 2, ..., n − 1, when complete pivot selection is used. This indicates that Gaussian elimination with complete pivot selection is usually a stable process.

Despite this, partial pivot selection is preferred in practice, for the most part, because:

(i) Complete pivot selection is more costly than partial pivot selection.

(To compute A(i), the maximum from among (n − i + 1)2 elements must be determined instead of n − i + 1 elements as in partial pivot selection.)

(ii) Special structures in a matrix, i.e. the band structure of a tridiagonal matrix, are destroyed in complete pivot selection.

(50)

university-logo

### Iterative Improvement:

Suppose that the linear system Ax = b has been solved via the LR-factorization P A = LR. Now we want to improve the accuracy of the computed solution x. We compute



r = b − Ax, Ly = P r, Rz = y, xnew = x + z.

(2.13)

Then in exact arithmatic we have

Axnew= A(x + z) = (b − r) + Az = b.

Az = r by using P A = LR.

(51)

university-logo

Unfortunately, r = f l(b − Ax) renders an xnew that is no more accurate than x. It is necessary to compute the residual b − Ax with extended precision floating arithmetic.

### Algorithm 4

Compute P A = LR (t-digit) Repeat: r := b − Ax (2t-digit)

Solve Ly = P r for y (t-digit) Solve Rz = y for z (t-digit) Update x = x + z (t-digit)

From Theorem 14 we have (A + δA)z = r, i.e.,

A(I + F )z = r with F = A−1δA. (2.14)

(52)

university-logo

Theorem 15

Let {xk} be the sequence constructed by Algorithm 4 for solving Ax = b and x= A−1b be the exact solution. Assume Fk in (2.14) satisfying kFkk ≤ σ < 1/2 for all k. Then {xk} → x.

Proof of Theorem 15

Corollary 16 If

1.01(n3+ 3n2)ρ2−tkAk kA−1k < 1/2, then Algorithm 4 converges.

Proof: From (2.14) and (2.12) follows that

kFkk ≤ 1.01(n3+ 3n2)ρ2−tκ(A) < 1/2.

(53)

university-logo

## Appendix

Proof of Theorem 12: Let sp= f l(Pp

k=1xkyk) be the partial sum in Algorithm 41. Then

s1= x1y1(1 + δ1) with |δ1| ≤ eps and for p = 2, . . . , n,

sp= f l[sp−1+ f l(xpyp)] = [sp−1+ xpyp(1 + δp)](1 + εp) with |δp|, |εp| ≤ eps. Therefore

f l(xTy) = sn = Xn k=1

xkyk(1 + γk),

where (1 + γk) = (1 + δk)Qn

j=k(1 + εj), and ε1≡ 0. Thus, f l(

Xn k=1

xkyk) = Xn k=1

xkyk[1 + 1.01(n + 2 − k)θk2−t].

The result follows immediately from the following useful Lemma.

(54)

university-logo

Lemma 7.1 If (1 + α) =Qn

k=1(1 + αk), where |αk| ≤ 2−t and n2−t≤ 0.01, then Yn

k=1

(1 + αk) = 1 + 1.01nθ2−twith |θ| ≤ 1.

Proof: From assumption it is easily seen that

(1 − 2−t)n≤ Yn k=1

(1 + αk) ≤ (1 + 2−t)n.

Expanding the Taylor expression of (1 − x)n as −1 < x < 1, we get (1 − x)n = 1 − nx +n(n − 1)

2 (1 − θx)n−2x2≥ 1 − nx.

Hence

(1 − 2−t)n ≥ 1 − n2−t.

(55)

university-logo

Now, we estimate the upper bound of (1 + 2−t)n: ex= 1 + x +x2

2! +x3

3! + · · · = 1 + x +x

2x(1 +x 3 +2x2

4! + · · · ).

If 0 ≤ x ≤ 0.01, then

1 + x ≤ ex≤ 1 + x + 0.01x1

2ex≤ 1 + 1.01x

(Here, we use the fact e0.01< 2 to the last inequality.) Let x = 2−t. Then the left inequality of (55) implies

(1 + 2−t)n≤ e2−tn

Let x = 2−tn. Then the second inequality of (55) implies e2−tn ≤ 1 + 1.01n2−t

From (55) and (55) we have

(1 + 2−t)n≤ 1 + 1.01n2−t.

Return

(56)

university-logo

Proof of Theorem 15: From (2.14) and rk= b − Axk we have A(I + Fk)zk= b − Axk.

Since A is nonsingular, we have (I + Fk)zk= x− xk.

From xk+1= xk+ zk we have (I + Fk)(xk+1− xk) = x− xk, i.e., (I + Fk)xk+1= Fkxk+ x. (2.15) Subtracting both sides of (2.15) from (I + Fk)x we get

(I + Fk)(xk+1− x) = Fk(xk− x).

Then, xk+1− x= (I + Fk)−1Fk(xk− x). Hence, kxk+1− xk ≤ kFkkkxk− xk

1 − kFkk ≤ σ

1 − σkxk− xk.

Let τ = σ/(1 − σ). Then

kxk− xk ≤ τk−1kx1− xk.

But σ < 1/2 follows τ < 1. This implies convergence of Algorithm 4.

Return

In this Learning Unit, students should be able to use Cramer’s rule, inverse matrices and Gaussian elimination to solve systems of linear equations in two and three variables, and

Abstract Based on a class of smoothing approximations to projection function onto second-order cone, an approximate lower order penalty approach for solving second-order cone

Based on a class of smoothing approximations to projection function onto second-order cone, an approximate lower order penalty approach for solving second-order cone

where L is lower triangular and U is upper triangular, then the operation counts can be reduced to O(2n 2 )!.. The results are shown in the following table... 113) in

Since the subsequent steps of Gaussian elimination mimic the first, except for being applied to submatrices of smaller size, it suffices to conclude that Gaussian elimination

Since the subsequent steps of Gaussian elimination mimic the first, except for being applied to submatrices of smaller size, it suffices to conclude that Gaussian elimination

Department of Mathematics – NTNU Tsung-Min Hwang November 30, 2003... Department of Mathematics – NTNU Tsung-Min Hwang November