Direct Methods for Solving Linear Systems

(1)

Chapter 6 Direct Methods for Solving Linear Systems

Hung-Yuan Fan (范洪源)

Department of Mathematics, National Taiwan Normal University, Taiwan

Spring 2016

(2)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Section 6.1

Linear Systems of Equations

(線性系統; 線性聯立⽅程組)

(3)

Objective

To solve a system of n linear equations with n unknowns:

E

₁: a₁₁x₁+ a₁₂x₂+· · · + a1nx_n= b₁

E

₂: a₂₁x₁+ a22x₂+· · · + a2nx_n= b2

...

E

_n: a_n1x₁+ a_n2x₂+· · · + annx_n= b_n.

(1)

Definition

The vectorx = [x₁, x₂,· · · , xn]^T∈ Rⁿ is called a solution to the linear system (1).

(4)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Matrix-Vector Forms (矩陣-向量形式)

The linear system (1) can be rewritten as the followg matrix-vector form:

Ax = b,

where coefficient matrix A∈ Rⁿ^×n and right-hand side vector b∈ Rⁿ are defined by

A =







a₁₁ a₁₂ · · · a1n

a₂₁ a₂₂ · · · a2n

... ... ... a_n1 a_n2 · · · ann





, b =





 b₁ b₂ ... b_n





.

Assumption: The matrix A is nonsingular throughout the

context.

(5)

The Linear Solvers (線性系統的數值⽅法)

1 Direct Methods: (直接法)

Used for solving small- or medium-sized linear systems with full and dense coefficient matrices. (適⽤於求解中⼩型線性系統)

Floating-point operation count (簡稱 flop)≈ O(n³).

Gaussian Elimination (⾼斯消去法，簡稱 GE) is an efficient

and stable algorithm for solving this type of linear systems.

2 Iterative Methods: (迭代法)

Used for solving large and sparse linear systems with problem size n≥ 10⁴. (適⽤於求解⼤型稀疏線性系統)

flop≈ O(n) per iteration if A is a sparse matrix.

Jacobi, Gauss-Seidel, SOR and CG-based methods,. . ..

(6)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Elementary Row Operations (基本列運算)

Three Row Operations

In order to simplify linear system (1), we use

1 (λEi)→ (Ei): eq. Ei is replaced by λ· Ei for any λ̸= 0.

2 (E_i+ λE_j)→ (Ei): eq. E_j is multiplied by any λ∈ R and added to eq. E_i.

3 (Ei)↔ (Ej): exchange eqs. Ei and E_j for i̸= j.

(7)

Gaussian Elimination (GE)

The Procedure of GE

Given a linear system Ax = b with A∈ Rⁿ^×n and b∈ Rⁿ.

1 Form the augmented matrix ˜A⁽¹⁾= [A| b] ∈ Rⁿ^×(n+1).

2 Applying row operations continuously, we obtain a finite sequence of augmented matrices, i.e.,

A˜⁽¹⁾ → Ã⁽²⁾→ · · · → Ã⁽ⁿ⁾, where Ã⁽ⁿ⁾ is upper triangular. (上三⾓矩陣)

3 Use backward substitution (向後代入) to obtain x_n, xn−1,· · · , x2, x₁.

(8)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

From ˜ A

⁽¹⁾

to ˜ A

⁽²⁾

From the augmented matrix ˜A⁽¹⁾= [A| b] = [a⁽¹⁾ij ], we have

A˜⁽¹⁾ =







a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾1,n+1

a⁽¹⁾₂₁ a⁽¹⁾₂₂ · · · a⁽¹⁾_2,n+1 ... ... ... a⁽¹⁾_n1 a⁽¹⁾_n2 · · · a⁽¹⁾n,n+1





 .

(9)

From ˜ A

⁽¹⁾

to ˜ A

⁽²⁾

(Conti’d)

If a⁽¹⁾₁₁ ̸= 0, do(

E_i− ^a⁽¹⁾ⁱ¹

a⁽¹⁾₁₁E₁)

→ (Ei) for i = 2, 3, . . . , n⇒

A˜⁽²⁾ =







a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1,n+1 0 a⁽²⁾₂₂ · · · a⁽²⁾_2,n+1

... ... ... 0 a⁽²⁾_n2 · · · a⁽²⁾n,n+1





 .

(10)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

From ˜ A

⁽²⁾

to ˜ A

⁽³⁾

Next, if a⁽²⁾₂₂ ̸= 0, do(

E_i− ^a⁽²⁾ⁱ²

a⁽²⁾₂₂ E₂)

→ (Ei) for i = 3, 4, . . . , n

⇒

A˜⁽³⁾ =







a⁽¹⁾₁₁ a⁽¹⁾₁₂ a⁽¹⁾₁₃ · · · a⁽¹⁾_1,n+1 0 a⁽²⁾₂₂ a⁽²⁾₂₃ · · · a⁽²⁾_2,n+1 ... 0 a⁽³⁾₃₃ · · · a⁽³⁾_3,n+1 ... ... ... ... 0 0 a⁽³⁾_n3 · · · a⁽³⁾n,n+1





 .

(11)

From ˜ A

^(k⁻¹⁾

to ˜ A

^(k)

, k ≥ 2

For each k≥ 2, suppose augmented matrix ˜A^(k⁻¹⁾ has the form

A˜^(k⁻¹⁾=







a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1,k₋₂ a⁽¹⁾_1,k₋₁ a⁽¹⁾_1k · · · a⁽¹⁾_1,n+1 0 a⁽²⁾₂₂ · · · a⁽²⁾_2,k₋₂ a⁽²⁾_2,k₋₁ a⁽²⁾_2k · · · a⁽²⁾_2,n+1

... 0 ... ... ... ...

... ... . .. a^(k−2)_k_−2,k−2 a^(k−2)_k_−2,k−1 a^(k−2)_k_−2,k · · · a^(k−2)_k_−2,n+1 ... .

.. 0 a^(k_k_−1,k−1⁻¹⁾ a^(k_k_−1,k⁻¹⁾ · · · a^(k_k_−1,n+1⁻¹⁾ ... .

.. .

.. a^(k_k,k−1⁻¹⁾ a^(k_k,k⁻¹⁾ · · · a^(k_k,n+1⁻¹⁾

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

0 0 · · · 0 a^(k_n,k⁻¹⁾₋₁ a^(k_n,k⁻¹⁾ · · · a^(k_n,n+1⁻¹⁾







(12)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

From ˜ A

^(k⁻¹⁾

to ˜ A

^(k)

, k ≥ 2 (Conti’d)

If a^(k_k_−1,k−1⁻¹⁾ ̸= 0, do (

E_i− ^a^(k^i,k−1⁻¹⁾

a^(k_k−1,k−1⁻¹⁾ E_k₋₁)

→ (Ei) for k≤ i ≤ n ⇒

A˜^(k) =







a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1,k−2 a⁽¹⁾_1,k−1 a⁽¹⁾_1k · · · a⁽¹⁾_1,n+1 0 a⁽²⁾₂₂ · · · a⁽²⁾_2,k₋₂ a⁽²⁾_2,k₋₁ a⁽²⁾_2k · · · a⁽²⁾_2,n+1

... 0 ... ... ... ...

... .

.. . .. a^(k_k−2,k−2⁻²⁾ a^(k_k−2,k−1⁻²⁾ a^(k_k−2,k⁻²⁾ · · · a^(k_k−2,n+1⁻²⁾ ... .

.. 0 a^(k−1)_k_−1,k−1 a^(k−1)_k_−1,k · · · a^(k−1)_k_−1,n+1 ... .

.. .

.. 0 a^(k)_k,k · · · a^(k)_k,n+1

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

0 0 · · · 0 0 a^(k)_n,k · · · a^(k)_n,n+1







(13)

At Final Stage of GE

When k = n, the GE will produce an augmented matrix in upper

triangular form, i.e., (省略右上標記號)

A˜⁽ⁿ⁾ ≡







a

₁₁ a₁₂ · · · a_1,n₋₁ a_1n | a_1,n+1 0

a

₂₂ · · · a_2,n₋₁ a_2n | a_2,n+1

... 0 . .. ... ... | ... ... ... . .. an−1,n−1 a_n_−1,n | an−1,n+1

0 0 · · · 0

a

_nn | a_n,n+1





 .

Note: Entries a

_ii̸= 0 are called pivot elements (軸元) for 1≤ i ≤ n.

(14)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Backward Substitution (向後代入)

From the special form of ˜A⁽ⁿ⁾, we obtain the solution x to linear system (1) as follows:

Firstly, compute x_n= an,n+1/ann.

Then compute x_n₋₁, xn−2,· · · , x1 successively by

x_i =

a_i,n+1− ∑ⁿ

j=i+1

a_ijx_j

a_ii , i = n− 1, n − 2, . . . , 1.

Note: This process is called the backward substitution of GE.

(15)

Backward Substitution with n = 3

Applying Backward Substitution for the 3× 3 linear system a₁₁x₁+ a12x₂+ a13x₃= b1

a₂₂x₂+ a₂₃x₃ = b₂ a₃₃x₃ = b₃,

the unique solution to above linear system is computed via x₃= b3/a33,

x₂= b₂− a23x₃ a₂₂ , x₁= b₃− a12x₂− a13x₃

a₁₁ .

(16)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Will GE break down? (1/2)

Example 2, p. 363 (GE 無法執⾏的例⼦) Use GE to find the solution of the linear system

E₁ : x₁− x2+ 2x3− x4=−8, E₂ : 2x1− 2x2+ 3x3− 3x4 =−20, E₃ : x₁+ x₂+ x₃ =−2,

E₄ : x₁− x2+ 4x₃+ 3x₄ = 4.

Sol: Form the augmented matrix ˜

A⁽¹⁾ as

A˜⁽¹⁾ =







1 −1 2 −1 | −8

2 −2 3 −3 | −20

1 1 1 0 | −2

1 −1 4 3 | 4





 .

(17)

Will GE break down? (2/2)

Since a⁽¹⁾₁₁ = 1, do (Ei− a⁽¹⁾i1 E₁)→ Ei for i = 2, 3, 4, we have

A˜⁽²⁾ =







1 −1 2 −1 | −8

0

−1 −1 | −4

0 2 −1 1 | 6

0 0 2 4 | 12





 .

Because a⁽²⁾₂₂ = 0, GE will break down here and STOP!

How to fix it? Partial Pivoting Strategy!

(18)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Partial Pivoting (1/2)

If we do (E₂)↔ (E3), then ˜A⁽²⁾ becomes

A˜⁽³⁾ =







1 −1 2 −1 | −8

0

2

−1 1 | 6

0 0 −1 −1 | −4

0 0 2 4 | 12





 .

Applying (E₄+ 2E3)→ (E4) =⇒

A˜⁽⁴⁾ =







1 −1 2 −1 | −8

0

2

−1 1 | 6

0 0 −1 −1 | −4

0 0 0 2 | 4





 .

(19)

Partial Pivoting (2/2)

Use backward substitution ⇒ x₄= 4

2 = 2

x₃= [−4 − (−1)x4]

−1 = 2

x₂= [6− (−1)x3− x4]

2 = 3

x₁= [−8 − (−1)x2− 2x3− (−1)x4]

1 =−7.

It seems that GE with partial pivoting works well in this case!

(20)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Pseudocode of Gaussian Elimination

To solve the n× n linear system (1).

Algorithm 6.1: GE with Backward Substitution INPUT dimension n; augmented matrix A = [a_ij]∈ Rⁿ^×(n+1). OUTPUT solution x₁, x₂,· · · , xn.

Step 1 For i = 1, . . . , n− 1 do Steps 2–4 Step 2 Find smallest i≤ p ≤ n s.t. api̸= 0.

If not, OUTPUT(‘No unique solution exists.’); STOP.

Step 3 If p̸= i, perform (Ep)↔ (Ei).

Step 4 For j = i + 1, . . . , n do Steps 5–6 Step 5 Set mji= aji/aii.

Step 6 Perform (Ej− mjiEi)→ (Ej).

Step 7 If a_nn = 0, OUTPUT(‘No unique solution exists.’); STOP.

Step 8 Set x_n= a_n,n+1/a_nn. (Start backward substitution.) Step 9 For i = n− 1, . . . , 1 set xi= [ai,n+1−∑n

j=i+1a_ijx_j]/aii. Step 10 OUTPUT(x₁, x₂,· · · , xn); STOP.

(21)

Operation Counts in Steps 5 and 6

Multiplications/Divisions:

n−1

∑

i=1

[(n− i) + (n − i)(n − i + 1)] =

n−1

∑

i=1

(n− i)(n − i + 2)

=

n−1

∑

i=1

(n²− 2ni + i²+ 2n− 2i) =

2n

³+ 3n²− 5n

6

.

Additions/Subtractions:

n−1

∑

i=1

(n− i)(n − i + 1) =

n−1

∑

i=1

(n²− 2ni + i²+ n− i) =

n

³− n

3

.

(22)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Operation Counts in Steps 8 and 9

Multiplications/Divisions:

1 +

n−1

∑

i=1

[(n− i) + 1] = 1 +

n−1

∑

i=1

(n− i) + (n − 1) =

n

²+ n

2

.

Additions/Subtractions:

n−1

∑

i=1

[(n− i − 1) + 1] =

n−1

∑

i=1

(n− i) =

n

²− n

2

.

(23)

Tota Number of Arithmetic Operations

Multiplications/Divisions:

2n³+ 3n²− 5n

6 +n²+ n 2 =

n

³

3

+ n²−

n 3

.

Additions/Subtractions:

n³− n

3 +n²− n 2 =

n

³

3

+

n

²

2

−

5n

6

. Total flop of GE ≈O(²₃n³) (as n→ ∞).

(24)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Section 6.2

Pivoting Strategies

(25)

Occurrence of Small Pivot Elements (1/2)

Example 1, p. 372

Apply GE to solve the linear system

E₁ :

0.003000x

₁+ 59.14x2 = 59.17 E₂ : 5.291x₁− 6.130x2 = 46.78,

using 4-digit rounding arithmetic and the exact solution is

x

₁ = 10.00 and x2 = 1.000.

Sol: Note that m

₂₁= fl(_0.003000^5.291 ) = fl(1763.6¯6) = 1764. Then do (E2− m21E₁)→ (E2)⇒

fl(−6.130 + fl(−1764 · 59.14))x2= fl(46.78 + fl(−1764 · 59.17)) or−104300x2=−104400 and hence x2 ≈ 1.001.

(26)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Occurrence of Small Pivot Elements (2/2)

Substituting x₂ = 1.001 into E1, we obtain x₁ ≈ fl(59.17− (59.14)(1.001)

0.003000

)

=−10.00, which gives a totally wrong answer for x₁! Why?

(27)

Partial Pivoting (部分選軸元)

For each i = 1, 2, . . . , n− 1, find smallest integer p ≥ i s.t.

|a⁽ⁱ⁾_pi| = max

i≤j≤n|a⁽ⁱ⁾_ji |.

Then perform the row operation

(Ei)↔ (Ep) if p̸= i.

This strategy is also called maximal column pivoting.

(28)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Improvement of Accuracy for Example 1

Example 2, p. 373 (改進例題 1 的計算精度)

Apply GE with partial pivoting to solve the linear system E₁ :

0.003000x

₁+ 59.14x₂ = 59.17 E₂ : 5.291x1− 6.130x2 = 46.78, using 4-digit rounding arithmetic.

Sol: Since

|a11| < |a21|, we first perform (E1)↔ (E2):

E₁ : 5.291x₁− 6.130x2 = 46.78 E₂ : 0.003000x₁+ 59.14x₂= 59.17.

Then m₂₁= fl(^0.003000_5.291 ) = 0.0005670. Do

(E₂− m21E₁)→ (E2)⇒ 59.14x2≈ 59.14 and hence x2= 1.000.

Moreover, x₁ = 10.00 is correct now!

(29)

Algorithm 6.2: GE with Partial Pivoting

INPUT dimension n; augmented matrix A = [a_ij]∈ Rⁿ^×(n+1). OUTPUT solution x₁, x₂,· · · , xn.

Step 1 For i = 1, . . . , n− 1 do Steps 2–5

Step 2 Find smallest i≤ p ≤ n s.t.|api| = maxi≤j≤n|aji|.

Step 3 Ifapi= 0, OUTPUT(‘No unique solution exists.’); STOP.

Step 9 Set x_n= an,n+1/ann. (Start backward substitution.) Step 10 For i = n− 1, . . . , 1 set xi= [ai,n+1−∑_n

j=i+1a_ijx_j]/aii. Step 11 OUTPUT(x₁, x2,· · · , xn); STOP.

(30)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Occurrence of Coefficients of Large Magnitude

Example 1’

Apply GE with partial pivoting to solve the linear system E₁ :

30.00x

₁+ 591400x₂ = 591700 E₂ : 5.291x₁− 6.130x2 = 46.78,

using 4-digit rounding arithmetic. Exact solution is x₁ = 10.00 and x₂= 1.000.

Sol: Note that

m₂₁= fl(5.291

30.00) = 0.1764.

Perform (E₂− m21E₁)→ (E2)⇒ −104300x2 ≈ −104400. Hence,

x

₂ ≈ 1.001 and x1 ≈ −10.00! Why?

(31)

Scaled Partial Pivoting

At the start of GE, compute n scale factors once as follows.

s_i = max

1≤j≤n|aij|, i = 1, 2, . . . , n.

For each i = 1, 2, . . . , n− 1, find smallest integer p ≥ i s.t.

|a⁽ⁱ⁾pi|

s_p = max

i≤j≤n

|a⁽ⁱ⁾ji | s_j . If i̸= p, perform (Ei)↔ (Ep).

This strategy is also called scaled-column pivoting.

(32)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Example 2’

Apply GE with scaled partial pivoting to solve the linear system E₁ :

30.00x

₁+ 591400x₂ = 591700

E₂ : 5.291x1− 6.130x2 = 46.78, using 4-digit rounding arithmetic.

Sol: First compute scale factors s

₁= 591400 and ss2 = 6.130. For i = 1, we see that

|a11|

s₁ = fl( 30.00

591400) = 0.5073× 10⁻⁴< |a21|

s₂ = fl(5.291

6.130) = 0.8631.

So, perform (E₂)↔ (E1)⇒ we obtain the correct solution x₁ = 10.00 and x₂ = 1.000!

(33)

Algorithm 6.3: GE with Scaled Partial Pivoting INPUT dimension n; augmented matrix A = [a_ij]∈ Rⁿ^×(n+1). OUTPUT solution x₁, x2,· · · , xn.

Step 1 For i = 1, . . . , n set s_i = max1≤j≤n|aij|;

Ifs_i= 0, OUTPUT(‘No unique solution exists.’); STOP.

Step 2 For i = 1, . . . , n− 1 do Steps 3–6

Step 3 Find smallest i≤ p ≤ n s.t. ^|a_s^pi_p^|= max_i_≤j≤n^|a_s^ji^|

j .

Step 4 Ifapi= 0, OUTPUT(‘No unique solution exists.’); STOP.

Step 10 Set x_n= an,n+1/ann. (Start backward substitution.) Step 11 For i = n− 1, . . . , 1 set xi= [ai,n+1−∑_n

j=i+1a_ijx_j]/aii.

(34)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Complete Pivoting

For each k = 1, 2, . . . , n− 1, find integers k ≤ p, q ≤ n s.t.

|a^(k)pq| = max

k≤i,j≤n|a^(k)_ij |.

If p̸= i or q ̸= i, row and/or column interchanges are performed to bring a^(k)_pq to the pivot position a^(k)_kk.

This strategy is also called the maximal pivoting at the kth step.

(35)

Section 6.5

Matrix Factorization

(矩陣分解)

(36)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

Motivation

For solving a linear system Ax = b, it requires O(¹₃n³) arithmetic operations to determine x∈ Rⁿ.

If the right-hand vector b∈ Rⁿ is changed to another vector ˜b (and coeff. matrix A is unchanged), how can we solve this linear system efficiently using some matrix factorization of A generated from GE?

In fact, if A has been factored into the triangular form A = LU,

where L is lower triangular and U is upper triangular, then the operation counts can be reduced toO(2n²)!

(37)

Comparison of Arithmetic Calculations

The relative rate of reduction of the operation counts O(2n²) compared with O(¹₃n³) becomes larger and larger for n = 10, 10² and 10³, respectively. The results are shown in the following table.

(38)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

The 1st Step of GE

Let A⁽¹⁾= A∈ Rⁿ^×n and b⁽¹⁾ = b∈ Rⁿ for a linear system.

If a⁽¹⁾_1,1 ̸= 0, do(

E_i− ^a⁽¹⁾^i,1

a⁽¹⁾_1,1E₁)

→ (Ei) for i = 2, 3, . . . , n⇒

A⁽²⁾ =







a⁽¹⁾_1,1 a⁽¹⁾_1,2 · · · a⁽¹⁾_1,n 0 a⁽²⁾_2,2 · · · a⁽²⁾_2,n ... ... ... 0 a⁽²⁾_n,2 · · · a⁽²⁾n,n







, b⁽²⁾ =





 b⁽²⁾₁ b⁽²⁾₂ ... b⁽²⁾_n





 .

The corresponding multipliers are given by

m_i,1= a⁽¹⁾_i,1

a⁽¹⁾₁₁ , i = 2, 3, . . . , n.

(39)

The 1st Step of GE (Conti’d)

This is equivalent to

A⁽²⁾= M⁽¹⁾A⁽¹⁾ and b⁽²⁾= M⁽¹⁾b⁽¹⁾, where the first Gaussian transformation matrix M⁽¹⁾∈ Rⁿ^×n is defined by

M⁽¹⁾=







1 0 · · · 0

−m_2,1 1 0 · · · 0

−m_3,1 0 . .. ... ...

... ... . .. ... 0

−m_n,1 0 · · · 0 1





 .

(40)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

At the kth Step of GE, 2 ≤ k ≤ n − 1

If a^(k)_k,k ̸= 0, do( E_i− ^a

(k) i,k

a^(k)_k,kE_k)

→ (Ei) for k + 1≤ i ≤ n ⇒

A^(k+1) =







a⁽¹⁾_1,1 · · · a⁽¹⁾_1,k a⁽¹⁾_1,k · · · a⁽¹⁾_1,n 0 . ..

.. .

.. . ..

. a^(k)_k,k a^(k)_k,k+1 · · · a^(k)_k,n ..

. 0 a^(k+1)_k+1,k+1 · · · a^(k+1)_k+1,n ..

.

.. .

.. . 0 · · · 0 a^(k+1)_n,k+1 · · · a^(k+1)n,n







, b^(k+1)=





b^(k+1)₁ b^(k+1)₂

.. . b^(k+1)n



 .

The corresponding multipliers are given by

m_i,k= a^(k)_i,k

a^(k)_k,k, i = k + 1, k + 2, . . . , n.

(41)

At the kth Step of GE, 2 ≤ k ≤ n − 1 (Conti’d)

This is equivalent to

A^(k+1) = M^(k)A^(k) and b^(k+1)= M^(k)b^(k), where the kth Gaussian transformation matrix M^(k)∈ Rⁿ^×n is defined by

M^(k)=







1

0 · · · · · · 0

0 . .. . .. ...

...

1

0 ...

... −mk+1,k

1

. .. ...

... ... . .. 0

0 −mn,k

1





 .

(42)

. . . . . .

. . . . . . . .

. . .

. .

. . . . .

The Inverse Matrix of M

^(k)

It is easily seen that the inverse matrix of M^(k) is given by

L^(k)= [M^(k)]⁻¹=







1

0 · · · · · · 0

0 . .. . .. ...

...

1

0 ...

... mk+1,k

1

. .. ...

... ... . .. 0

0 mn,k

1







(2)

for each k = 1, 2, . . . , n− 1.

Check . . .

M^(k)L^(k)= I or L^(k)M^(k)= I, where I denotes the n× n identity matrix.