Module IV: Solution of linear system

Amlan Dutta∗
Department of Metallurgical and Materials Engineering,
Indian Institute of Technology Kharagpur, West Bengal 721302, India
(Dated: October 7, 2025)

1
 Abstract
A linear system refers to a set of linear equations in multiple variables. While the standard
methods of linear algebra offer techniques to solve such systems, it turns out such methods cannot
be directly applied in a practical scenario. It is due to the unfavorable time–complexity, which
renders such methods infeasible for a computer if the number of variables is even moderately large.
In this module, we shall learn ways to overcome this issue through more efficient numerical means
with reduced time–complexity.

1. LINEAR SYSTEMS

In applied sciences and engineering, we often encounter linear equations of the form,

Ax̄ = b̄, (1)

where A is an (n×n) coefficient matrix and b̄ is an n–dimensional vector (column matrix).

The goal is to find the value of the vector variable, x̄, which satisfies Eq. (1). The above
equation may be expanded as a set of linear equations:

a11 x1 + a12 x2 + ... + a1n xn = b1 ;

a21 x1 + a22 x2 + ... + a2n xn = b2 ;
.. (2)
.
an1 x1 + an2 x2 + ... + ann xn = bn ,
or simply,

n
X
aij xj = bi . (3)
j=1

Equation (1) is solvable iff |A| =

̸ 0, i.e., A is non–singular. In that case the system may
be solved simply by applying the Cramer’s rule,

|Ai |
xi = , i = 1, 2, ..., n, (4)
|A|

∗
[Link]@[Link]

2
where Ai is the matrix obtained by replacing the i-th column of A by b̄. Although the
Cramer’s method appears to be a direct and simple way of solving Eq. (1), in practice, it
loses its utility for n ≳ 15. It is because if the determinants are evaluated using Laplace
expansion, the Cramer’s rule would need about 3(n + 1)! floating point operations. Using
some sophisticated methods of computing the determinants, the time complexity can be
brought down to O(n3.8 ), which is still not enough to be put into practical use.

2. THE LU FACTORIZATION METHOD

Let us assume that the matrix, A ∈ Rn×n , can be factorized as,

A = LU, (5)

with L and U are lower and upper triangular matrices, respectively, both of size n × n.
Since A is non–singular, so are L and U , which further implies that all their diagonal ele-
ments are non–zero.

By substituting Eq. (4) in Eq. (1), we obtain,

LU x̄ = b̄. (6)

The above equation can be split into two linear equations,

Lȳ = b̄; (7)

U x̄ = ȳ. (8)

We begin with Eq. (7), which can be rewritten as,

n
X
lij yj = bi . (9)
j=1

Since L is lower triangular, lij = 0, ∀j > i. Therefore, Eq. (9) becomes,

i
X
lij yj = bi , (10)
j=1

3
 i.e.,

l11 y1 = b1 ;
l21 y1 + l22 y2 = b2 ;
.. (11)
.
ln1 y1 + an2 y2 + ... + ann yn = bn .
Equation (10) (or (11)) can be rearranged as,

Pi−1 n o
1
bi − i−1
P
j=1 lij yj + lii yi = bi , =⇒ yi = lii j=1 lij y j . (12)

Thus, we arrive at the so called forward substitution iterations,


 b1 ;

i = 1,
l
yi = 11 n o (13)
 1 bi − Pi−1 lij yj ; i ≥ 2.

lii j=1

The iterations given above require n2 floating point operations.

Once the vector y has been evaluated using the iterative scheme of Eq. (13), we can
proceed to solve Eq. (8). To this end, we use a procedure similar to that given in Eqs. (9)–
(13) but in the reverse direction starting from i = n and iterating sequentially for i =
n − 1, n − 2, ..., 2, 1. We arrive at the following scheme (detail shown in the class):

yn
; i = n,


unn
xi = n o (14)
1
yi − nj=i+1 uij xj ; i = n − 1, n − 2, ..., 1.

 P
uii

The above sequence is known as backward substitution and also required n2 operations
similar to the forward substitution. Thus, the original problem, Ax̄ = b̄ reduces to finding
the matrices L and U , such that A = LU .

3. FACTORIZATION OF MATRIX A

The technique of factorizing the coefficient matrix can be understood by first considering
a 2 × 2 matrix, A, such that LU = A.
    
l11 0 u11 u12 a11 a12
∴  = ,
l21 l22 0 u22 a21 a22

4
i.e.,

l11 u11 = a11 ; (15) l11 u12 = a12 ; (16)

l21 u11 = a21 ; (17) l21 u12 + l22 u22 = a22 . (18)
This system is essentially underdetermined, for we have four equations against six un-
known variables. Hence, we are free to choose arbitrary values of any two of the variables
and solve for the remaining four. Let us assign the value of unity to two diagonal elements
of L, i.e., l11 = l22 = 1. Then, from Eqs. (15) and (16), we have,

u11 = a11 /l11 = a11 ; u12 = a12 /l11 = a12 .

Now Eqs. (17) and (18) can be solved in sequence as,

l21 = a21 /u11 ; u22 = a22 − l21 u12 .

Thus, we have been able to factorize the matrix A ∈ Rn×n into L and U .
Let us try a similar strategy for a 3 × 3 system. In this case,
    
l11 0 0  u11 u12 u13  a11 a12 a13 
∴ l21 l22 0   0 u22 u23  = a21 a22 a23  ,
    
    
l31 l32 l33 0 0 u33 a31 a32 a33
i.e.,

l11 u11 = a11 ; (19) l11 u12 = a12 ; (20) l11 u13 = a13 ; (21)

l21 u11 = a21 ; (22) l21 u12 + l22 u22 = a22 ; (23) l21 u13 + l22 u23 = a23 ; (24)

l31 u11 = a31 ; (25) l31 u12 + l32 u22 = a32 ; (26) l31 u13 + l32 u23 + l33 u33 = a33 . (27)

As this system is also underdetermined, we set all the diagonal elements of L to 1, thereby
obtaining the equations,

u11 = a11 ; (28) u12 = a12 ; (29) u13 = a13 ; (30)

l21 = a21 /a11 ; (31) u22 = a22 − l21 u12 ; (32) u23 = a23 − l21 u13 ; (33)

l31 = a31 /a11 ; (34) l32 = (a32 − l31 u12 ) /u22 ; (35) u33 = a33 − l31 a13 − l32 u23 . (36)

Hence, on a matrix, A ∈ Rn×n , we proceed as follows:

5
 1. The elements of L and U satisfy the equations,
min(i,j)
X
lir urj = aij ; i, j = 1, 2, ..., n. (37)
r=1

2. The above system has n2 + n unknowns for n2 equations, thereby making it underde-
termined.

3. By setting the diagonal elements of L to 1, the system becomes determined, which is

then solved using the Guass LU factorization algorithm. We start by setting, A(1) = A,
(1)
i.e., aij = aij ∀i, j = 1, 2, ..., n. Then we proceed as follows:

for k = 1, 2, ..., n − 1
for i = k + 1, ..., n
(k)
aik
lik = (k)
akk
for j = k + 1, ..., n
(k+1) (k) (k)
aij = aij − lik akj

(k)
Clearly all the akk values must be non-zero and are called pivot elements. At the end
n o
(k+1)
of each iteration, k, the matrix A(k+1) ≡ aij has n − k rows and columns updated.
A special feature of the Gauss’ algorithm is that it eliminates the need of storing the lower
and upper triangular matrices separately. Instead the values of A are overwritten in a way
(i)
that at the end of the procedure, the upper triangular matrix, U , is given as, uij = aij , for
i = 1, 2, ..., n and j = i, ..., n. Similarly, matrix L is given by, lij , and is overwritten on the
lower triangular part of A. Observe that there is no computation of the diagonal elements
(i = j) of L, for we know them to be unity. The iterations for factorizing the matrix, A,
requires about 2n3 /3 operations.

4. CONDITIONS OF FACTORIZABILITY

The LU factorization does not exist for every non–singular matrix, A. In this context,
the following result holds:
For a given matrix, A ∈ Rn×n , its LU factorization exists and is unique iff all the principal
submatrices, Ai of A, of order i = 1, 2, ..., n − 1 are non–singular. This result also holds for
A ∈ Cn×n .

6
 For example, consider the matrices,
  
1 1−ϵ 3 5−ϵ
   
A = 2 2 2 , and b̄ =  6  .
   
   
3 6 4 13
 
1
One can verify that, x̄ = 1, is a solution of the system Ax̄ = b̄ ∀ϵ. If ϵ = 1, A
 
 
1
factorizes into,
   
1 0 0 1 0 3 
L = 2 1 0 , and U = 0 2 −4 .
   
   
3 3 1 0 0 7
However, if we set ϵ = 0, the Gauss LU decomposition algorithm will fail to produce the
(k)
aik
result even though A is non–singular. This is because in the step, lik = (k) , the pivoting
akk
element at the denominator becomes 0.

4.1. Special cases

There are special cases, where it is guaranteed that all the principal submatrices of A
would be non–singular, thereby confirming unique LU factorization. These cases are as
follows:

1. Strictly diagonally dominant matrix

A is said to be diagonally dominant if,

n
X
|aii | ≥ |aij | , i = 1, 2, ..., n (dominant by row)
j=1
j̸=i

and/or if,
n
X
|aii | ≥ |aji | , i = 1, 2, ..., n (dominant by column)
j=1
j̸=i

If in the above conditions, ≥ is replaced by >, A becomes strictly diagonally dominant.

It is also valid for A ∈ Cn×n .

7
 2. Real symmetric and positive definite matrices

A ∈ Rn×n is positive definite if,

x̄T Ax̄ > 0; ∀x̄ ∈ Rn with x̄ ̸= 0̄.

3. Complex positive definite matrices

A ∈ Cn×n is positive definite if,

x̄H Ax̄ > 0; ∀x̄ ∈ Cn with x̄ ̸= 0̄.

5. CHOLESKY FACTORIZATION

If A ∈ Rn×n is symmetric and positive definite, it is possible to obtain a special case of

factorization,

A = RT R,

which is known as the Cholesky factorization with R as an upper triangular matrix. It

requires ∼ n3 /3 operations, which is about half of the regular LU factorization. Due to
symmetry, only the upper part of A and R are required to be stored. The following iterations
implement the Cholesky factorization:
√
Set, r11 ← a11
for j = 2, 3, ..., n
for i = 1, 2, ..., j − 1
 Pi−1 
1
rij = rii aij − k=1 rki rkj
q Pj−1 2
rjj = ajj − k=1 rkj