# Lecture 2- Solution of Sparse Linear Systems

ECE 530 – Analysis Techniques for Large – Scale Electrical Systems
Solution of Sparse Linear Systems

George Gross
Department of Electrical and Computer Engineering

University of Illinois at Urbana-Champaign

ECE 530

1

INTRODUCTION
? We are concerned with the solution of linear systems

Ax ? b
where A is an n ? n matrix with elements aij , and x and b are n - vectors with elements xi and bi respectively ? We are, particularly, interested when n is large and A is sparse
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 2

INTRODUCTION
? A matrix is said to be sparse if a large percentage of its elements aij have zero value ? Sparse linear systems arise naturally in many power systems analysis and power electronics problems ? The best known and most widely used method for solving linear systems of algebraic equations is attributed to Gauss
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 3

GAUSSIAN ELIMINATION
? Gaussian elimination is the elementary procedure in which we use the first equation to eliminate the first variable from the last ( n - 1) equations, then we use the new second equation to eliminate the second variable from the last ( n ? 2) equations,

and so on
? After performing ( n ? 1) such eliminations we end

up with a triangular system which is easily solved
in a backward direction
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 4

EXAMPLE 1
? We need to solve for x in the system

3 ?1 0 ? ? x1 ? ? 20 ? ? 2 ? ? ?? 6 ? 5 ? ? ? 0 2 x2 ? 45 ? ?? ??? ? 6 ? 6 ? ? x3 ? ? ? 3 ? ? 2 ? 5 ?x ? ? ? 4 ? ? 2 2 ? 3 ? ? 4 ? ? 30 ? ?
? The three elimination steps are given below: for

simplicity, we have appended the r.h.s. vector to
the matrix
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 5

EXAMPLE 1
? Eliminate x1 by subtracting row 1 from all the rows below it
multiply row 1 by 1 2

multiply row 1 by 6 and add to row 2

multiply row 1 by ? 2 and add to row 3 multiply row 1 by ? 4 and add to row 4
ECE 530

? ?1 ? ?0 ? ?0 ? ?0 ?

3 2 4 ?8 ?4

?

1 2

0

? 3 ? 7 7

1 2
? 6 ? 3

? 10 ? ? 15? ? ? 23? ? ? 10 ? ?
6

EXAMPLE 1
? Eliminate x2 by subtracting row 2 from all the rows below it
multiply row 2 by 1 4

multiply row 2 by 8 and add to row 3 multiply row 2 by 4 and add to row 4
ECE 530

? ?1 ? ?0 ? ?0 ? ?0 ?

3 2 1 0 0

?

1 2

0

3 ? 4
1 1

1 2
?2 ? 1

? 10 ? ? 15 ? 4 ? 7? ? 5? ?
7

EXAMPLE 1
? Elimination of x3 from row 3 and 4

subtract row 3 from row 4
ECE 530

? 1 ? ? ?0 ? ?0 ? ?0 ?

3 2

?

1 2

0

3 1 ? 4
0 0 1 0

1 2
?2 1

? 10 ? ? 15 ? 4 ? 7? ? ? 2? ?
8

EXAMPLE 1
? Then, we solve for x by “going backwards”, i.e., using backsubstitution:
x4 x3 ? 2x4 ? ?2 ? 7 ? x3 ? 3

3 1 x2 ? x3 ? x 4 4 2 3 1 x1 ? x2 ? x3 2 2
ECE 530

? ?

15 ? 4 10 ?

x2 x1

? ?

7 1
9

EXAMPLE 1
? In this example, we have: ? triangularized the original matrix by Gaussian elimination using column elimination

? then, we used backsubstitution to solve the
triangularized system
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 10

TRIANGULAR DECOMPOSITION
? We give a general scheme for the triangular

decomposition of a matrix by Gaussian
elimination

? We base the development of the system,
A x=b

(1)

where x is the unknown vector and b ? 0
? We assume that A is a nonsingular matrix, i.e., A

is invertible and A exists
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 11

?1

TRIANGULAR DECOMPOSITION
? We form the matrix A a using A and b with
? a11 ? ? ? a21 ? ? ? a31 ? ? ? ? ? a n1 ? a12 a22 a32 a1n a2 n a3 n b1 ? ? ? b2 ? ? ? (2) b3 ? ? ? ? ? bn ? ?
12

Aa

=

[ A b]

=

an 2

ann

and show the steps of the triangularization scheme

TRIANGULAR DECOMPOSITION
? Step 1: normalize the first equation

a

(1) 1j

=

a1 j a 11

j ? 2 , ... ,n

b
ECE 530

(1) 1

=

b1 a 11

? ? ? ? (3) ? ? ?
13

TRIANGULAR DECOMPOSITION
? Step 2: a) eliminate x1 from row 2:

(1) a (1) = a ? a a 2j 2j 21 1 j,

j = 2 , ... ,n

(1) b (1) = b ? a b 2 2 21 1

? ? ?(4) ? ? ? ?
14

ECE 530

TRIANGULAR DECOMPOSITION
b) normalize the second equation

a

(2) 2j

=

a a

(1) 2j (1) 22

,

j = 3 , ... ,n

b (2) 2 =
ECE 530

b (1) 2 a
(1) 22

? ? ? ?(5) ? ? ?
15

TRIANGULAR DECOMPOSITION
and we end up at the end of step 2 with

? 1 ? ? ? 0 ? ? a31 ? ? ? ? a n1 ? ?
ECE 530

a

(1) 12

a a

(1) 13 (2) 23

a

(1) 1n (2) 2n

1 a32

a

a33

a3 n

an 2

an 3

ann

? ? (2) ? b2 ? ? b3 ? (6) ? ? ? ? bn ? ? b
(1) 1
16

TRIANGULAR DECOMPOSITION
? Step 3: a) eliminate x1 and x2 from row 3:
(1) a (1) = a ? a a 3j 3j 31 1j

j = 2 , ... ,n

(1) b (1) = b ? a b 3 3 31 1

(1) (2) a (1) ? a a 3j 32 2j

j = 3 , ... ,n b
(2) 2

b
ECE 530

(2) 3

= b

(1) 3

? a

(1) 32

? ? ? (7) ? ? ? ?
17

TRIANGULAR DECOMPOSITION
b) normalize the third equation:

a

(3) 3j

=

a (2) 3j a
(2) 33

j = 4 , ... ,n

b
ECE 530

(3) 3

=

b a

(2) 3 (2) 33

? ? ? (8) ? ? ? ?
18

TRIANGULAR DECOMPOSITION
and we have the system at the end of step 3

? 1 ? ? 0 ? ? 0 ? ? a 41 ? ? ? ? a n1
ECE 530

a (1) 12 1 0 a 42

a (1) 13 a (2) 21 1 a 43

a (1) 14 a (2) 24 a (3) 34 a 44

a (1) 1n a (2) 2n a (3) 3n a 4n

a n2

a n3

a n4

a nn

? b (2) ? 2 ? b (3) 3 ? ? b4 ? ? ? ? bn ?

? b (1) 1

(9)

19

TRIANGULAR DECOMPOSITION
? In general, we have for

Step

k : a) eliminate x1 , x2 , ... , xk ? 1 from row k :
j = m + 1 , ... ,n ?

(m ?1 ) (m ?1 ) (m) a (m) = a ? a kj kj km a mj

(m ?1 ) (m ?1 ) (m) b (m) = b ? a k k km b m
ECE 530

? ? (10) ? ? m = 1, 2 , ... , k - 1 ? ?
20

TRIANGULAR DECOMPOSITION
where, the superscript m denotes the order of the
) derived system: for row k , a ( m is the value of the kj

element in column j after the elements in columns
1,2, ... , m have become 0, i.e., after eliminating
th k the x1 , x2 , ... , xm variables from the equation
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 21

TRIANGULAR DECOMPOSITION
Step

k : b) normalize the k equation:
?1) a (k kj

th

a (k) kj =

a

(k ? 1 ) kk

j = k + 1 , ... ,n

b (k) k
ECE 530

?1) b (k k = ?1) a (k kk

? ? ?(11) ? ? ? ?
22

TRIANGULAR DECOMPOSITION
and proceed in this manner until we obtain the
upper triangular matrix (the n th derived system):

?1 ? ?0 ? ?0 ? ?0 ? ? ? ?0
ECE 530

(1) a 12

(1) a 13

(1) a 14

(1) a 1n

1 0 0

a (2) 21 1 0

a (2) 24 a (3) 34 1

a (2) 2n a (3) 3n a (4) 4n

0

0

0

1

? b (2) 2 ? ? (3) b3 ? ? (4) b 4 ? ? ? (n) ? b n?

? b (1) 1

(12)

23

TRIANGULAR DECOMPOSITION
? Note that in the scheme presented, unlike in the example, we triangularly decompose the system

by eliminating row – wise rather than column –
wise: in successive rows we eliminate (reduce to 0) each element to the left of the diagonal rather

than those below the diagonal
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 24

TRIANGULAR DECOMPOSITION
? To compute x in (1), we perform backsubstitution

xn ? bn
xn-1 =

? n?

bn-1 ? an-1 ,n xn

? n ? 1?

? n ? 1?

xk

=

bk ?

?k ?

j=k+1

?a

n

?k ?
k, j

xj

? ? ? ? ? ? (13) ? ? ? ? k = n ? 1, n ? 2 , ... , 1 ? ? ?
25

ECE 530

UPPER TRIANGULAR MATRIX
? The triangular decomposition scheme applied to the matrix A results in the upper triangular

matrix U in (12) with the elements

? 1 ? ? (i) u ij = ? a ij ? ? 0 ?
ECE 530

i = j j > i j < i
26

(14)

UPPER TRIANGULAR MATRIX
? The matrix U plays a key role in the so – called factorization of A ? The following theorem is an important result that is the basis for the development of the computational scheme
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 27

THEOREM 1 STATEMENT
? Any nonsingular matrix A has the following

factorization:

A = L U
where U is given by (14) and L is defined by

(15)

ij

? a (jij ? 1 ) ? = ? ? 0 ?

j ? i

(16)

j ? i
28

ECE 530

THEOREM 1 : PROOF
? To verify (15), let

C = LU
? Then,

ck, j =

m=1

?

n

min? k, j ? k,m

um, j =

m=1

?

a k,m u m, j

? m ? 1?

(17)

? From (14) and (10), we have that
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 29

THEOREM 1 : PROOF
?1) (m ? 1 ) (m) (m ? 1 ) (m) a (m u = a a = a ? a k,m m, j k,m m, j k, j k, j

for j ? m and m ? k ? 1 ? If j > k , then min(k, j) = k so that (17) becomes
k ?1

ck, j =

m=1

(m ? 1 ) (k ? 1 ) a u + a ? k,m m, j k, j =

(m ? 1 ) (m) ? (k ? 1 ) ? a ? a + a = ak, j ? k, j k, j k, j ? ? m=1 (18)
30

k ?1

ECE 530

THEOREM 1 PROOF
? However, if j ? k , then min(k, j) = j so that (17) becomes

?1) (m ? 1 ) (j ? 1 ) (j) ck, j = ? a (m u = a u + a a k, j = ak, j (19) ? k,m m, j k,m m, j k,k m=1 m=1 ?1) a (j k, j

j

j -1

? Thus, (15) is verified
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 31

THEOREM 1
? This theorem shows that every nonsingular matrix A may be expressed as the product of a

lower triangular matrix L and an upper triangular
matrix U

? We rewrite (1) as

A x = LU x = L y= b y

(20)
32

THEOREM 1
? Once L and U are known for a given A , then we can solve for x by first solving

L y = b

(21)

for y using forward substitution and then solve for
.

x using backsubstitution

U x = y
? We may think of L as a record of the forward operations performed on b

(22)

33

THEOREM 2 STATEMENT
? The vector,
(2) (n) ? b = ?b (1) ,b , ... ,b n ? ? 1 2 ^ T

(23)

is related to the original b by the transformation

b = L ?1 b

^

(24)
34

THEOREM 2 : PROOF
? Let

c = L b
? Then for

^

k = 1, 2 , ... , n
ck ?

m=1

?l

k

k,m

bm?

^

m=1

?a

k

(m ? 1 ) k,m

bm ?
(m)

k ?1

m=1

?a

(m ? 1 ) k,m

b m ? a k,k
(m)

(k ? 1 )

bk

(k)

ECE 530

35

THEOREM 2 : PROOF
? Now, from (10) and (11) we have:
?1) (m) (m ? 1 ) (m) a (m b = b ? b , k,m m k k

m = 1,2 , ... ,k ? 1

and
? 1 ) (k) (k ? 1 ) a (k b = b k,k k k

? Then, for k = 1, 2 , ... , n
ck =
ECE 530

(m ? 1 ) (m) (k ? 1 ) (0 ) ? ? b ? b ? b = b = bk ? k k ? k k ? m =1

k ?1

THEOREM 2 : PROOF
? Thus,

c= b

and the theorem is proved since L is nonsingular
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 37

THEOREM 2
? We note that by having a record of the decomposition of A we can solve (1) for any b without

repeating the triangularization
? We store the L and U triangular factors into a composite array which we term the table of

factors denoted by
ECE 530

? L \U ? :
38

THEOREM 2 : PROOF
? ? ? ? ? ? L \U ? = ? ? ? ? ? ? ? ?
11

u12

u13 u23

21

22

n1

n2

n3

u1n ? ? ? u2 n ? ? ? ? ? ? ? ? nn ? ?

with
? L \U ?i, j
ECE 530

? i, j ? ? = ? ?u ? ? i, j

i? j

(25)
i< j
39

THEOREM 2 IMPLICATIONS
? The diagonal terms of U are not stored since they are known to be 1’s by the steps of the

Gaussian elimination
? The record in (25) makes it unnecessary to include the vector b since the table of factors has all the

required information
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 40

THEOREM 2 IMPLICATIONS
? We can further decompose the array L as
L= LD,

(26)

where D is a diagonal matrix with elements
? ? d ij = ? ? ?
ECE 530

ii

i? j

(27)
i? j
41

0

THEOREM 2 IMPLICATIONS
? Then, L is a unit lower rectangular matrix with
elements
~

ij

=

? ? ? ? ? ? ? ? ?

0 1

i < j i = j
(28)
ij

d jj

i > j
42

ECE 530

THEOREM 2 IMPLICATIONS
? Then we may write

A = L DU
which provides the L DU decomposition of A

~

(29)

? A special case of (29) is obtained for the case of a
symmetric matrix A :
A = AT

(30)
43

THEOREM 2 IMPLICATIONS
? Since
~ ~

A = L DU = U

T

D LT = A

T

? Since the transpose of a lower triangular matrix is an upper triangular matrix and vice versa, and the .L D U factorization is unique it follows that

U = LT
and, therefore,
A= U
ECE 530

~

T

DU

(31)
44

SYMMETRIC MATRIX FACTORIZATION
? In the case of a symmetric A, only the diagonal
and upper triangular elements of the table need to be stored, thereby reducing storage requirements down by nearly a factor of 2 ? For a symmetric matrix, in effect, it is unnecessary

to perform any operations to the left of the
diagonal
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 45

SYMMETRIC MATRIX FACTORIZATION
? We can obtain the elements a (j ? 1 ) , k ? j k, j

using the relation

a k, j

(j ? 1 )

=

k, j

= d

~ j, j k, j

=

j, j

?1) (j) u j,k = a (j a j, j j,k

(31a)
j = 1,2 , ... ,k ? 1
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 46

SYMMETRIC MATRIX FACTORIZATION
?1 ) ? In fact, we can obtain a (j for k ? j from the k, j

unnormalized row j since (31a) and (11) imply that

?1) (j ? 1 ) (j) (j ? 1 ) a (j = a a = a k, j j, j j,k j,k
ECE 530

j = 1,2 , ... ,k ? 1 (31b)
47

SYMMETRIC MATRIX FACTORIZATION
? We use the result in (31a) to modify our scheme
so as to bring about an efficiency improvement for symmetric matrices, thereby saving on the number of multiplication/division operations in (31a)

? Thus, after processing row j each term is left in
? 1) its unnormalized form a ( jj,k , k = j , ... , n
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 48

SYMMETRIC MATRIX FACTORIZATION
? 1) , j = 1, 2, ... , k ? 1 is needed ? Now, the term a ( jj ,k

as a multiplying factor for row j to process row .k and in the modified scheme is obtained directly from the unnormalized row k , i.e., from the term

a

( j ? 1) j,k
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 49

SYMMETRIC MATRIX FACTORIZATION
? Thus, the operation given in (10) for row k ,i.e.,
= j ? 1 , ... ,n
(j ? 1 ) (j ? 1 ) (j) a (j) = a ? a a k k kj j

j = 1 , 2 , ... ,k ? 1

is replaced in the modified scheme using (11) by
a a
(j ? 1 ) jl (j ? 1 ) jj

= j ? 1 , ... ,n j = 1 , 2 , ... ,k ? 1
50

a kl = a
ECE 530

(j)

(j ? 1 ) kl

?a

(j ? 1 ) kj

SYMMETRIC MATRIX FACTORIZATION
(j ? 1 ) ? ? (j ? 1 ) a kj (j ? 1 ) = a kl ? ? (j ? 1 ) ? a j ? ? a jj ? ?

(31c)
(j ? 1 ) ? ? (j ? 1 ) a jk (j) (j ? 1 ) a k = a k ? ? (j ? 1 ) ? a j j =1,2 , ... ,k ? 1 = k,k ? 1, ... ,n ? ? a jj ? ? ? The bracketed term is simply a (j) j k – the

?1 ) normalized value – and replaces a (j as soon as jk

it has been used in (31c)
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 51

SYMMETRIC MATRIX FACTORIZATION
? It is important to note that in this modified scheme

we perform no operations on any elements to the

left of the diagonal
ECE 530 ? 2002 - 2010 George Gross, University of Illinois at Urbana-Champaign, All Rights Reserved. 52

SUMMARY
? We reproduce the key results and relations developed in this section

? To solve

Ax = b
? We decompose A into LU and L DU factors, so
that we solve A x ? b with
~

LU x = b
ECE 530

or

L DU x = b
53

~

SUMMARY
? 1 ? ? u ij = ? a (i) ij ? ? 0 ?
? ? ii ? ? 0 ?
ECE 530

i = j j > i j < i
i = j i ? j
ij

~ ij

? ? =? ? ?

0 1
ij

i< j i= j i> j

d jj

d ij =

?1 ) ? a (j ? ij = ? ? ? 0

j ? i j > i
54

SUMMARY
? The key basic relations are
(m ? 1 ) (m ? 1 ) (m) a (m) = a ? a a kj kj km mj

(m ? 1 ) (m ? 1 ) (m) b(m) = b ? a bm k k km

j = m ? 1 , ... ,n ? ? ? (10) ? m = 1,2 , ... ,k ? 1 ? ? ?
b (k) k =
?1) b(k k

a (k) kj =

?1) a (k kj

a

(k ? 1 ) kk
ECE 530

j = k ? 1, ... ,n

a

(k ? 1 ) kk

(11)
55

SUMMARY
? An important special case is the symmetric

matrix A ,i.e., A ? A and the decomposition

T

reduces to

A = LU = L DU = U D L ? L = U
T
ECE 530

~

T

~

~

T

56

SUMMARY
? 1 ? ? (i) u ij = ? a ij ? ? 0 ? i = j j > i j < i
ij ?1) ? a (j ? ij = ? ? ? 0

j ? i j > i

? ? ii d ij = ? ?0 ?
ECE 530

i = j i ? j

~
i, j

= u j,i