A

SHORT NOTE
ON
INTRODUCTION TO NUMERICAL LINEAR
ALGEBRA (INLA)

Department of Mathematics

School of CS & AI

SR University, Warangal-506371

Telangana
Topics Covered

➢ Norms of vectors and matrices

➢ singular value decomposition (SVD)
➢ LU factorization

➢ Cholesky factorization
➢ QR factorization

Topics to Recall before studying this Module

➢ Basics of Matrix

➢ Determinant and properties

➢ Transpose of a Matrix
➢ Multiplication of two Matrices
➢ Inverse of a Matrix
➢ Singular and Non-singular Matrix
➢ Solving system of linear equations (Ax=B)
➢ Characteristic Polynomial and Equation
➢ Eigen Values and Eigen vectors
 Module-2: Factorization
2.1. NORMS OF VECTORS AND MATRICES
2.1.1. Introduction
Norms of vectors and matrices are fundamental concepts in linear algebra and functional analysis, playing
a crucial role in various mathematical and computational applications. In linear algebra, a norm is a
mathematical function that assigns a non-negative value to vectors and matrices, eventually measuring their
size or magnitude. It serves as a quantitative measure of distance or size, providing a way to characterize
and compare different mathematical entities. For vectors, norms capture the idea of the length or size of a
vector in a given space. The concept extends naturally to matrices, where matrix norms help in assessing
the behaviour and properties of linear transformations.

2.1.2. Vector Norms

Definition. A vector norm on ℝ𝑛 is a function ‖. ‖: ℝ𝑛 → ℝ satisfies the following conditions:

(N1) ‖𝑥‖ > 0 f o r e v e r y n o n - z e r o 𝑥 a n d ‖𝑥‖= 0 iff x = 0. (positivity)

(N2) ‖𝜆𝑥‖= |𝜆|‖𝑥‖ f o r a l l 𝑥 ∈ ℝ𝑛 and for all scalars 𝜆. (scaling)

(N3) ‖𝑥 + 𝑦‖≤ ‖𝑥‖ + ‖𝑦‖ f o r al l 𝑥, 𝑦 ∈ ℝ𝑛 . (triangle inequality)

Note:
1. ‖−𝑥‖ = ‖𝑥‖
2 . ‖𝑥‖ − ‖𝑦‖ ≤ ‖𝑥 − 𝑦‖ f o r a l l 𝑥, 𝑦 ∈ ℝ𝑛

Different types of Vector Norms.

𝑥1
𝑥2
Let 𝑥 = (𝑥1 , 𝑥2 , … , 𝑥𝑛 ) =
⋮ . Then
[𝑥𝑛 ]
A. ‖𝑥‖1 = |𝑥1 | + |𝑥2 | + ⋯ + |𝑥𝑛 | (Sum norm or one norm).

B. ‖𝑥‖2 = √𝑥12 + 𝑥22 + ⋯ + 𝑥𝑛2 (Euclidean norm or two norm).

C. ‖𝑥‖∞ = max |𝑥𝑖 | (Infinity norm or maximum norm).

𝑖

In general, if p is a real number greater than or equal to 1, then the p-norm or Holder
norm is defined by
 1
‖𝑥‖𝑝 = (|𝑥1 |𝑝 + |𝑥2 |𝑝 + ⋯ + |𝑥𝑛 |𝑝 )𝑝

Example 2.1.1. Let 𝒙 = (𝟏, 𝟏, −𝟐) be a vector. Find the 1-norm, Euclidean norm and infinity
norm.
Solution: The sum norm or 1-norm is given by
‖𝑥‖1 = |1| + |1| + |−2| = 1 + 1 + 2 = 4
The Euclidean norm or two norm is given by
‖𝑥‖2 = √12 + 12 + (−2)2 = √6
The Infinity norm or maximum norm is given by
‖𝑥‖∞ = Max |𝑥𝑖 | = Max{|1|, |1|, |−2|} = Max{1,1, 2} = 2.
𝑖

2.1.4. Matrix Norms

Definition. Let A he an 𝑚 × 𝑛 matrix. Then, analogous to the vector norm, we define a matrix norm
‖𝐴‖ satisfies the following conditions:

(M1) ‖𝐴‖ > 0 f o r e v e r y n o n - z e r o m a t r i x 𝐴 a n d ‖𝐴‖= 0 iff A is a zero matrix.

(M2) ‖𝜆𝐴‖= |𝜆|‖𝐴‖ f or a n y m at r i x A a n d for all scalars 𝜆.

(M3) ‖𝐴 + 𝐵‖≤ ‖𝐴‖ + ‖𝐵‖ f o r al l m at r i ces A a n d B .

Subordinate Matrix Norms

Given a matrix A and a vector norm II·II, a nonnegative number defined by

‖𝐴𝑥‖𝑝
‖𝐴‖𝑝 = Max
𝑥≠0 ‖𝑥‖𝑝

satisfies all the properties of a matrix norm. This norm is called the matrix norm subordinate
to the vector norm or just subordinate matrix norm.
A very useful and frequently used property of a subordinate matrix norm (we shall sometimes call it the
p-norm of a matrix A) is ‖𝐴𝑥‖𝑝 ≤ ‖𝐴‖𝑝 ‖𝑥‖𝑝 .

The two easily computable p-norms are

‖𝐴‖1 = max ∑𝑚
𝑖=1 |𝑎𝑖𝑗 | (maximum column-sum norm),
1≤𝑗≤𝑛

‖𝐴‖∞ = max ∑𝑛𝑗=1 |𝑎𝑖𝑗 | (maximum row-sum norm),

1≤𝑖≤𝑚
 𝟏 −𝟐
Example 2.1.2. Let 𝑨 = ( 𝟑 𝟒 ) be a matrix. Find the maximum column-sum and row-
−𝟓 𝟔
sum norm.
Solution: The column-sum norm is given by
𝑚
‖𝐴‖1 = max ∑|𝑎𝑖𝑗 |
1≤𝑗≤𝑛
𝑖=1

= max(|𝑎1𝑗 | + |𝑎2𝑗 | + |𝑎3𝑗 |)

𝑗

= max{(|1| + |2| + |−5|), (|−2| + |4| + |6|)}

= max{ 8, 12} = 12
The row-sum norm is given by
𝑛
‖𝐴‖∞ = max ∑|𝑎𝑖𝑗 |
1≤𝑖≤𝑚
𝑗=1

= max(|𝑎𝑖1 | + |𝑎𝑖2 | + |𝑎𝑖3 |)

𝑖

= max{(|1| + |−2|), (|3| + |4|), (|−5| + |6|)}

= max{ 3, 7,11} = 11
Example 2.1.3.

Solution:
Spectral Norm of a matrix
Another useful p-norm is the spectral norm, denoted by ‖𝐴‖2 a n d i s d ef i ne d b y
‖𝐴𝑥‖2
‖𝐴‖2 = Max .
𝑥≠0 ‖𝑥‖2

It can be shown that

‖𝐴‖2 = √maximum eigenvalue of 𝐴𝑇 𝐴
Note that the eigenvalues of 𝐴𝑇 𝐴 are real and nonnegative.

𝟐 𝟓
Example 2.1.4. Let 𝑨 = ( ) be a matrix. Find the spectral norm of A.
𝟏 𝟑
S o l u t i o n : The eigenvalues of 𝐴𝑇 𝐴 are 𝜆1 = 0.0257 and 𝜆2 = 38.9743. (Verify yourself!!!)
Then the spectral norm of A is given by
‖𝐴‖2 = √maximum eigenvalue of 𝐴𝑇 𝐴
𝑜𝑟, ‖𝐴‖2 = √38.9743
‖𝐴‖2 = 6.2429
The Frobenius Norm
An important matrix norm compatible with the vector norm ‖𝑥‖2 is the Frobenius norm:
1
𝑛 𝑚 2
2
‖𝐴‖𝐹 = [∑ ∑|𝑎𝑖𝑗 | ]
𝑗 𝑖

A matrix norm ‖. ‖𝑀 and vector norm ‖. ‖𝑣 are compatible if ‖𝐴𝑥‖𝑣 ≤ ‖𝐴‖𝑀 ‖𝑥‖𝑣 .
𝟐 𝟓
Example 2.1.5. Let 𝑨 = ( ). Then find the Frobenius Norm of the matrix.
𝟏 𝟑
Solution: The Frobenius Norm of the matrix A is given by
1
𝑛 𝑚 2
2
‖𝐴‖𝐹 = [∑ ∑|𝑎𝑖𝑗 | ] = √30.
𝑗 𝑖

Two important properties of the Frobenius norm

1. For the identity matrix I, ‖𝐼‖𝐹 = √𝑛, w h e r ea s whereas ‖𝐼‖1 = ‖𝐼‖2 = ‖𝐼‖∞ = 1.
2. ‖𝐴‖𝐹 = √𝑡𝑟𝑎𝑐𝑒(𝐴𝑇 𝐴)
Consistent norm
A norm ‖. ‖ is consistent if it satisfies ‖𝐴𝐵‖ ≤ ‖𝐴‖‖𝐵‖ whenever the product 𝐴𝐵 is defined.
The Frobenius norm and all subordinate matrix norms are consistent. The max norm ‖𝐴‖ is not
consistent
Exercise 2.1
1. Define 1-norm, Euclidean norm and infinity norm of a vector.
2. Find the 1-norm, Euclidean norm and infinity norm of the following vectors:
1
(a) [0]
2
1
0
(b) [ ]
−2
2
(c) (5,2,0, −1)
2. Find the matrix norms (maximum column-sum norm, maximum row-sum norm, Spectral Norm and
Frobenius norm) of the following matrices:
1 2
(a) [ ]
0 4
1 0 1
(b) [2 −1 0 ]
3 0 −2
2 −1 0
(c) [ 3 1 −2]
−2 0 2
1 2 −1 0
0 3 1 −2
(d) [ ]
0 0 3 6
0 0 0 2
2.2. LU Factorization
2.2.1. Definition
The 𝐿𝑈 factorization or 𝐿𝑈 decomposition of a matrix is an important matrix factorization useful for solving
a linear system and computing the determinant and the inverse of a matrix. The 𝐿𝑈 factorization of the
matrix 𝐴 can be written as 𝐴 = 𝐿𝑈, where 𝐿 is a lower triangular matrix and 𝑈 is upper triangular matrix.

Example 1 (General procedure of finding the LU factorization of a matrix)

1 2 4 1 0 0 𝑈11 𝑈12 𝑈13
Let 𝐴 = [3 8 14] = 𝐿𝑈, where 𝐿 = [𝐿21 1 0] and 𝑈 = [ 0 𝑈22 𝑈23 ]
2 6 13 𝐿31 𝐿32 1 0 0 𝑈33

1 2 4 1 0 0 𝑈11 𝑈12 𝑈13

or, 𝐴 = [3 8 14] = [𝐿21 1 0] [ 0 𝑈22 𝑈23 ]
2 6 13 𝐿31 𝐿32 1 0 0 𝑈33

1 2 4 𝑈11 𝑈12 𝑈13

or, [3 8 14] [ 21 𝑈11
= 𝐿 𝐿21 𝑈12 + 𝑈22 𝐿21 𝑈13 + 𝑈23 ]
2 6 13 𝐿31 𝑈11 𝐿31 𝑈12 + 𝑈22 𝐿31 𝑈13 + 𝐿32 𝑈23 + 𝑈33

We begin by running along the top row to see that

𝑈11 = 1
𝑈12 = 2
𝑈13 = 4

Now comparing the second row, we get

𝐿21 𝑈11 = 3 𝑜𝑟, 𝐿21 × 1 = 3 or, 𝐿21 = 3

𝐿21 𝑈12 + 𝑈22 = 8 or, 3 × 2 + 𝑈22 = 8, or, 𝑈22 = 2

𝐿21 𝑈13 + 𝑈23 = 14, or, 3 × 4 + 𝑈23 = 14 or, 𝑈23 = 2

Again, comparing the last row, we get

𝐿31 𝑈11 = 2, or, 𝐿31 × 1 = 2, or, 𝐿31 = 2

𝐿31 𝑈12 + 𝑈22 = 6, or, 2 × 2 + 𝐿32 × 2 = 6, or, 𝐿32 = 1

𝐿31 𝑈13 + 𝐿32 𝑈23 + 𝑈33 = 13, or, 2 × 4 + 1 × 2 + 𝑈33 = 13, or, 𝑈33 = 3.

Hence the matrix A has the LU factorization as follows:

1 2 4 1 0 0 1 2 4
𝐴 = [3 8 14] = [3 1 0] [0 2 2].
2 6 13 2 1 1 0 0 3
2.2.2. LU Factorization using Gaussian Elimination
Now we will discuss the procedure of obtaining the LU factorization of the matrix A by Gaussian
Elimination method with an example as follows:
Example 2.
2 2 3
Let us consider a square matrix 𝐴 = [4 5 6]
1 2 4
Step 1 (eliminate the entries on the first column of A below the diagonal)
1
Multiply the first row of A by -2 and − and add, respectively, to the second and third rows.
2

2 2 3
𝐴 (1)
= [0 1 0] ; ℎ𝑒𝑟𝑒 𝑚 = −2, 𝑚 = − 1 𝑎𝑟𝑒 𝑡ℎ𝑒 𝑚𝑢𝑡𝑖𝑝𝑙𝑖𝑒𝑟𝑠.
5 21 31 2
0 1 2

Step 2 (eliminate the entries of the second column of 𝑨(𝟏) below the diagonal)
Multiply the second column of 𝑨(𝟏) by −1 and add to the third row.
2 2 3
𝐴 (2)
= [0 1 0] ; ℎ𝑒𝑟𝑒 𝑚 = −1.
5 32
0 0 2
2 2 3
So, 𝑈 = 𝐴 (2)
= [0 1 0],
5
0 0
2

1 0 0 1 0 0
𝐿= 𝐿−1
1 = [−𝑚21 1 0] = [2 1 0].
1
−𝑚31 −𝑚32 1 2
1 1

2.2.3. LU factorization using Gaussian Elimination for a Rectangular Matrix

The above-described Gaussian elimination process for an 𝑛 × 𝑛 matrix A can he easily extended to an
𝑚 × 𝑛 matrix to compute its LU factorization, when it exists. The process is identical. However, the number
of required steps in this case is 𝑘 = min{𝑚 − 1, 𝑛}. We illustrate this with an example.
Example 3
𝟏 𝟐
Consider 𝑳 = [𝟑 𝟒] , 𝒎 = 𝟑, 𝒏 = 𝟐
𝟓 𝟔 𝟑×𝟐
Here 𝑘 = 𝑚𝑖𝑛{𝑚 − 1, 𝑛} = 𝑚𝑖𝑛(3 − 1,2) = 𝑚𝑖𝑛{2,2} = 2.
k =1 (Eliminate the entries in the first column of A below the diagonal)
The multipliers are 𝑚21 = −3, 𝑚31 = −5
 (1)
Update: 𝑎22 = 𝑎22 = 𝑎22 + 𝑚21 𝑎12 = −2
(1)
𝑎32 = 𝑎32 = 𝑎32 + 𝑚31 𝑎12 = −4.
𝟏 𝟐
𝐴 = 𝐴(1) = [𝟎 −𝟐].
𝟎 −𝟒
k = 2. (Eliminate the entries in the second column of 𝐴(1) below the diagonal)
The multipliers are 𝑚32 = −2
(1)
Update: 𝑎32 = 𝑎32 = 0
𝟏 𝟐
𝐴 = 𝐴(2) = [𝟎 −𝟐].
𝟎 𝟎
Note that U in this case is a 3 x 2 upper triangular matrix.
1 0 0 1 0 0
Form 𝐿 = [−𝑚 21 1 0] = [3 1 0].
−𝑚31 −𝑚32 1 5 2 1
1 0 0 𝟏 𝟐 𝟏 𝟐
Verification: 𝐿𝑈 = [3 1 0] [𝟎 −𝟐] = [𝟑 𝟒] (Verify yourself!!!)
5 2 1 𝟎 𝟎 𝟓 𝟔
2.2.3. Solution of systems of linear equations using LU decomposition
Once a matrix A has been decomposed into lower and upper triangular parts it is possible to obtain
the solution to AX = B in a direct way. The procedure can be summarised as follows
Step-1. Given A, find L and U so that A = LU. Hence LUX = B.
Step-2. Let Y = UX so that LY = B. Solve this triangular system for Y.
Step-3. Finally solve the triangular system UX = Y for X.
The benefit of this approach is that we only ever need to solve triangular systems. The cost is that
we have to solve two of them.

Example 4
𝒙𝟏 𝟏 𝟐 𝟒 𝒙𝟏 𝟑
Find the solution of 𝑿 = [𝒙𝟐 ] of the system [𝟑 𝟖 𝟏𝟒 ] [ 𝒙 𝟐 ] = [𝟏𝟑]
𝒙𝟑 𝟐 𝟔 𝟏𝟑 𝒙𝟑 𝟒
Solution
The first step is to calculate the LU decomposition of the coefficient matrix. Here, the factorization has
already been done (see Example 2.2.1).
We found that
 1 0 0 1 2 4
𝐿 = [3 1 0], 𝑈 = [0 2 2].
2 1 1 0 0 3
𝑦1
The next step is to solve 𝐿𝑌 = 𝐵 for the vector 𝑌 = [𝑦2 ]. That is, we consider
𝑦3
1 0 0 𝑦1 3
𝐿𝑌 = [3 1 0] [𝑦2 ] = [13] = 𝐵
2 1 1 𝑦3 4
which can be solved by forward substitution.
From the top equation we see that 𝑦1 = 3.
The middle equation states that 3𝑦1 + 𝑦2 = 13 and hence 𝑦2 = 4.
Finally, the last one that 2𝑦1 + 𝑦2 + 𝑦3 = 4 from which we see that 𝑦3 = −6.
Now that we have found 𝑌 we finish the procedure by solving 𝑈𝑋 = 𝑌 for 𝑋. That is, we solve
1 2 4 𝑥1 3
𝑈𝑋 = [0 2 2] [𝑥2 ] = [ 4 ] = 𝑌
0 0 3 𝑥3 −6
by using back substitution.
Starting with the bottom equation we see that 3𝑥3 = −6 so clearly 𝑥3 = −2.
The middle equation implies that 2𝑥2 + 2𝑥3 = 4 and it follows that 𝑥2 = 4.
The top equation states that 𝑥1 + 2𝑥2 + 4𝑥3 = 3 and consequently 𝑥1 = 3.
Therefore, the solution to then given system of simultaneous equations is
3
𝑋=[ 4 ]
−2
2.2.4. Existence (Do matrices always have an LU Factorization?)
No. Sometimes it is impossible to write a matrix in the form “lower triangular” × “upper triangular”.
Why not?
An invertible matrix A has an 𝐿𝑈 decomposition provided that all its leading submatrices have non-zero
determinants. The kth leading submatrix of A is denoted 𝐴𝑘 and is the k × k matrix found by looking only
at the top k rows and leftmost k columns. For example, if
1 2 4
𝐴 = [3 8 14]
2 6 13
then the leading submatrices are
1 2 4
1 2
𝐴1 = [1], 𝐴2 = [ ], 𝐴3 = [3 8 14]
3 8
2 6 13
The fact that this matrix A has an LU decomposition can be guaranteed in advance because none of these
determinants is zero:

𝟏 𝟐 𝟑
Example 5. Show that Show that [𝟐 𝟒 𝟓] does not have an LU decomposition
𝟏 𝟑 𝟒
Solution
The second leading submatrix has determinant equal to

which means that an LU decomposition is not possible in this case.

Do the following problem yourself

2.2.5. Difficulties of Gaussian Elimination without Pivoting

As we have seen before, Gaussian elimination without pivoting faits if any of the pivots is zero. However,
it is worse yet if any pivot becomes close to zero: in this case the method can be carried to completion, but
the obtained results may be totally wrong. However, we can avoid this small pivot just by row interchanges.
It is true that with the interchange above, we now obtained an LU factorization of the matrix A1, a permuted
version of the matrix A, and not of the original matrix.
2.2.6. Permutation Matrices
A permutation matrix P is a square matrix obtained by permuting the rows of an identity matrix I. Each row
and column of a permutation matrix contains only one non-zero element, which is equal to 1. The non-zero
elements are arranged such that each row and column have a single 1, and the rest of the elements are zeros.
0 1 0 1 0 0 1 0 0
Example: 𝑃1 = [0 0 1], 𝑃2 = [0 1 0], 𝑃3 = [0 0 1] are all permutation matrices.
1 0 0 0 0 1 0 1 0

2.2.7. LU Factorization using Gaussian elimination with partial pivoting

(GEPP)
The Gaussian elimination with partial pivoting yields a factorization of A in the form PA = LU, where P is
a permutation matrix, L is a unit lower triangular matrix, and U is an upper triangular matrix.
Partial pivoting involves choosing the pivot element (the element used for elimination in each step)
as the maximum absolute value among the elements in the current column. This helps prevent division by
very small numbers, enhancing numerical stability.
➢ Brief description of PLU Factorization of A as 𝑷𝑨 = 𝑳𝑼:
• Permutation Matrix (P): The partial pivoting introduces row exchanges during the Gaussian
elimination process. The permutation matrix P records these row exchanges. It is a square matrix
that, when multiplied with the original matrix A, results in a permuted matrix PA.
• Lower Triangular Matrix (L): The lower triangular matrix L contains the multipliers used during
the elimination process. These multipliers are chosen to eliminate the entries below the main
diagonal of the matrix. The matrix L is unit lower triangular, meaning that its diagonal entries are
all 1.
• Upper Triangular Matrix (U): The upper triangular matrix U is the result of the Gaussian
elimination process. It is a triangular matrix with zeros below the main diagonal, and its diagonal
entries are the pivots obtained during the elimination.
Further the method can be illustrated by the following example

𝟏 𝟐 𝟒
Example 6. Find the LU Factorization of the matrix 𝑨 = (𝟒 𝟓 𝟔) using Gauss elimination
𝟕 𝟖 𝟗
with partial pivoting method.
Solution

Verify that 𝑃𝐴 = 𝐿𝑈. (Verify yourself!!!)

2.2.8. Applications and disadvantages of LU decomposition
Applications of LU Factorization:
1. Solving Linear Systems: LU factorization is commonly used for solving systems of linear
equations. Once a matrix is factorized into LU form, it becomes computationally more efficient to
solve linear systems.
2. Matrix Inversion: LU factorization is employed to find the inverse of a matrix. If a matrix A has
been factorized as A=LU, the inverse of A can be found by solving 𝐴−1 = 𝑈 −1 𝐿−1.
3. Eigenvalue Problems: LU decomposition is used in eigenvalue computations and various
algorithms related to matrix eigenvalues.

Disadvantages of LU Factorization:
1. Pivoting Issues: Without proper pivoting strategies, LU factorization may encounter problems
with matrices that have small pivots or are close to singular. Pivoting is often required to improve
stability.
2. Storage Requirements: The LU factorization process requires additional storage space for the
factors L and U, which can be a concern for large matrices.
3. Computationally Expensive: LU factorization can be computationally expensive for large
matrices.

Exercises 2.2
2.3. Cholesky Factorization
2.3.1. Definition (positive definite matrix)
A symmetric matrix A is positive definite if for every nonzero vector x, 𝒙𝑻 𝑨𝒙 > 𝟎. Let 𝒙 =
𝑥1
𝑥2
(𝒙𝟏 , 𝒙𝟐 , . . . , 𝒙𝒏 ) = [ … ]. Then 𝒙𝑻 𝑨𝒙 = ∑𝒏𝒊,𝒋=𝟏 𝒂𝒊𝒋 𝒙𝒊 𝒙𝒋 is called the quadratic farm associated with 𝐴.
𝑥𝑛
A positive semidefinite matrix is similarly definite. A symmetric matrix 𝐴 is positive semidefinite
if 𝒙𝑻 𝑨𝒙 ≥ 0 for all vectors 𝑥.
A commonly used notation for a symmetric positive definite (positive semidefinite matrix) is 𝐴 >
0 (≥ 0).

2.3.2. Some characterizations and properties of positive definite matrices

1. A symmetric matrix 𝐴 is positive definite if and only if all its eigenvalues are positive.

2. A symmetric matrix A is positive definite if and only if all its leading principal minors are positive. i.e.,
the determinant of 𝐴𝑘 (the upper left 𝑘 × 𝑘 submatrix) should be positive.

3. If 𝐴 = (𝑎𝑖𝑗 ) is symmetric positive definite, then 𝑎𝑖𝑖 > 0 for all 𝑖.

4. If 𝐴 = (𝑎𝑖𝑗 ) is symmetric positive definite, then the largest element (in magnitude) of the whole matrix
must lie on the diagonal.

5. The sum of two symmetric positive definite matrices is symmetric positive definite.

𝟗 𝟔
Example 1: Let 𝑨 = [ ] be a 𝟐 × 𝟐 matrix.
𝟔 𝒂
A is symmetric
|A1|=9>0
9 6
|A2|= | | = 9𝑎 − 36
6 𝑎
If 𝑎 > 4, then 𝐴 is positive definite.
If 𝑎 = 4, then 𝐴 is positive semi definite.
If 𝑎 < 4, then 𝐴 is not positive definite.
 2 −𝟏 𝟎
Example 2: Is the matrix 𝑨 = [−𝟏 2 −1] positive definite?
𝟎 −1 𝟐
Solution
We will check whether all of its leading principal minors |𝐴𝑘 | are positive or not.
Now, |𝐴1 | = 2 > 0
𝟐 −𝟏
|𝐴2 | = | |=4−1=3>0
−𝟏 𝟐
𝟐 −𝟏 0
|𝐴3 | = |−𝟏 𝟐 −𝟏| = 2(4 − 1) + 1(−2 − 0) + 0(1 − 0) = 6 − 2 = 4 > 0
𝟎 −𝟏 𝟐
Therefore, all of its leading principal minors |𝐴𝑘 | are positive. So, the given matrix is positive definite.

2.3.3. Definition (Cholesky factorization)

The Cholesky factorization is a method to decompose a positive definite symmetric matrix into the product
of a lower triangular matrix and its transpose.
Given a symmetric positive definite matrix A, there exists a lower triangular matrix 𝐻 (called
Cholesky factor) with positive diagonal entries such that 𝐴 = 𝐻𝐻 𝑇
Note: LU decomposition can be used more generally. However, in the special case where 𝐴 is a
symmetric positive definite matrix, one can use Cholesky Factorization method that is more
efficient in time and space.

2.3.4. The Cholesky Factorization procedure for a 3×3 matrix

We now show how the Cholesky factorization can be computed directly from 𝐴 = 𝐻𝐻 𝑇 , when n = 3.
The general case is analogous.
𝑎11 𝑎12 𝑎13 ℎ11 0 0 ℎ11 ℎ21 ℎ31
[𝑎21 𝑎22 𝑎23 ] = [ℎ21 ℎ22 0 ][ 0 ℎ22 ℎ32 ]
𝑎31 𝑎32 𝑎33 ℎ31 ℎ32 ℎ33 0 0 ℎ33
A H HT
Step-1. Compute the first column of H. Compare the corresponding entries of the first column of both sides:
2
𝑎11 = ℎ11 ⇒ ℎ11 = √𝑎11
𝑎21
𝑎21 = ℎ11 ℎ21 ⇒ ℎ21 =
ℎ11
𝑎31
𝑎31 = ℎ11 ℎ31 ⇒ ℎ31 =
ℎ11
Step-2. Compute the second column of H. Compare the second and third entries of the second column of
both sides:

2 2 2
𝑎22 = ℎ21 + ℎ22 ⇒ ℎ22 = √𝑎22 − ℎ21

𝑎32 − ℎ21 ℎ31

𝑎32 = ℎ21 ℎ31 + ℎ22 ℎ32 ⇒ ℎ32 =
ℎ22
Step-3. Compute the third column of H. Compare the third entry of the third column of both sides:

2 2 2 2 2
𝑎33 = ℎ31 + ℎ32 + ℎ33 ⇒ ℎ33 = √𝑎33 − ℎ31 − ℎ32

In general, compute recursively the first through nth columns of H by comparing the entries of the respective
column of both sides of 𝐴 = 𝐻𝐻 𝑇 .
Remarks.
1. The matrix H above is computed column by column.
2. The positive definiteness of A will make the quantities under the square root signs positive.
𝟏 𝟏 𝟏
Example 2: Find the Cholesky factorization of the matrix 𝑨 = [𝟏 𝟓 𝟓]
𝟏 𝟓 𝟏𝟒
Solution
𝟏 𝟏 𝟏 𝑎11 𝑎12 𝑎13
𝑎
We have 𝑨 = [𝟏 𝟓 𝟓 ] = [ 21 𝑎22 𝑎23 ] (𝑠𝑎𝑦)
𝟏 𝟓 𝟏𝟒 𝑎31 𝑎32 𝑎33
ℎ11 0 0
Then the Cholesky factorization of A is given by 𝐴 = 𝐻𝐻 𝑇 , where 𝐻 = [ℎ21 ℎ22 0 ]
ℎ31 ℎ32 ℎ33
Now we will find the Cholesky Factor H in the following the steps:
Step-1 (Computation of the 1st column of H)
ℎ11 = √𝑎11 = √1 = 1
𝑎21 1
ℎ21 = = =1
ℎ11 1
𝑎31 1
ℎ31 = = =1
ℎ11 1

Step-2 (Computation of the 2nd column of H)

2
ℎ22 = √𝑎22 − ℎ21 = √5 − 12 = 2
 𝑎32 − ℎ21 ℎ31 5 − 1.1
ℎ32 = = =2
ℎ22 2

Step-3 (Computation of the 3rd column of H)

2 2
ℎ33 = √𝑎33 − ℎ31 − ℎ32 = √14 − 12 − 22 = ℎ33 = √9 = 3

1 0 0
Thus 𝐻 = [1 2 0]
1 2 3
1 1 1
Hence the Cholesky Factorization of the matrix 𝐴 = [1 5 5 ] is given by
1 5 14
1 1 1 1 0 0 1 1 1
𝐴 = [1 5 5 ] = 𝐻𝐻 𝑇 = [1 2 0] [0 2 2]
1 5 14 1 2 3 0 0 3

2.3.5. Application: Solving the system of equations 𝑨𝒙 = 𝒃 with positive definite

matrix 𝑨
Having the Cholesky factorization 𝐴 = 𝐻𝐻 𝑇 at hand, the positive definite linear system 𝐴𝑥 = 𝑏 can now
be solved by solving the lower trlangular system 𝐻𝑦 = 𝑏 first, followed by the upper triangular system
𝐻 𝑇 𝑥 = 𝑦.
𝟏 𝟏 𝟏 𝟏
Example 3: Solve the system 𝑨𝑿 = 𝑩 𝐰𝐡𝐞𝐫𝐞 𝑨 = [𝟏 𝟓 𝟓 ] and 𝑩 = [𝟏𝟏]
𝟏 𝟓 𝟏𝟒 𝟐𝟎
Solution
Step-1 (Find the Cholesky factorization 𝑨 = 𝑯𝑯𝑻 )
1 1 1
We have already found the Cholesky Factorization of the matrix 𝐴 = [1 5 5 ] in the previous example
1 5 14
which is given by
1 1 1 1 0 0 1 1 1
𝐴 = [1 5 5 ] = 𝐻𝐻 𝑇 = [1 2 0] [0 2 2]
1 5 14 1 2 3 0 0 3
So that the LHS of the given equation becomes 𝑨𝑿 = 𝐻𝐻 𝑇 𝑋 = 𝑯𝒀 where 𝑌 = 𝐻 𝑇 𝑋
Step-2 (To solve the system 𝑯𝒀 = 𝑩)
1 0 0 𝑦1 1
𝐻𝑌 = 𝐵 ⇒ [1 2 0] [𝑦2 ] = [11]
1 2 3 𝑦3 20
𝑦1 1
𝑦
Solving [ 2 ] = [5]
𝑦3 3
Step-3 (To solve the system 𝑯𝑻 𝑿 = 𝒀 )
1 1 1 𝑥1 1
𝑇 𝑥
𝐻 𝑋 = 𝑌 ⇒ [0 2 2] [ 2 ] = [ 5]
0 0 3 𝑥3 3
𝑥1 1
3
Solving [𝑥2 ] = [ 2 ]
𝑥3 −
3
2

Exercise 2.3
1. Find the Cholesky Factorization of the following symmetric positive definite matrices:
𝟐𝟓 𝟏𝟓 −𝟓
(a) [ 𝟏𝟓 𝟏𝟖 𝟎 ]
−𝟓 𝟎 𝟏𝟏
𝟗 −𝟔 𝟔
(b) [−𝟔 𝟓 −𝟏]
𝟔 −𝟏 𝟏𝟓
2. Solve the following systems of linear equations using Cholesky Factorization:
(a) 𝟐𝟓𝒙 + 𝟏𝟓𝒚 − 𝟓𝒛 = 𝟓
𝟏𝟓𝒙 + 𝟏𝟖𝒚 + 𝟎𝒛 = 𝟑
−𝟓𝒙 + 𝟎𝒚 + 𝟏𝟏𝒛 = 𝟏𝟐

(b) 𝟗𝒙 − 𝟔𝒚 + 𝟔𝒛 = 𝟐𝟏
−𝟔𝒙 + 𝟓𝒚 − 𝒛 = 𝟏
𝟔𝒙 − 𝒚 + 𝟏𝟓𝒛 = 𝟎
2.4. QR Factorization
2.4.1. Definition (Orthogonal matrix)
Recall that a square matrix 𝑄 is said to be an orthogonal matrix if 𝑄𝑄 𝑇 = 𝑄 𝑇 𝑄 = 𝐼,

2.4.2. Definition (QR factorization)

Given an 𝑚 × 𝑛 matrix 𝐴, there exist an 𝑚 × 𝑚 orthogonal matrix 𝑄 and an 𝑚 × 𝑛 upper triangular matrix
R such that 𝐴 = 𝑄𝑅. Such a factorization of A is called the 𝑄𝑅 factorization.

The process of obtaining the QR factorization involves the Gram-Schmidt orthogonalization or other
methods like Householder reflections or Givens rotations. In this course we will focus only on QR
factorization using Gram-Schmidt orthogonalization.

2.4.3. QR factorization using Classical Gram-Schmidt (CGS) process

The CGS process gives a reduced QR factorization of a matrix A whose columns are 𝑎1 , 𝑎2 , … , 𝑎𝑛
; that is, 𝐴 = [𝑎1 𝑎2 … 𝑎𝑛 ] ∈ ℝ𝑚×𝑛 by generating an orthogonal set of vectors 𝑞1 , 𝑞2 , … , 𝑞𝑛 .

The CGS process is described briefly as follows:

𝑣
Step-1: Set 𝑣1 = 𝑎1. Then compute 𝑞1 = ‖𝑣1‖.
1

Step-2: For 𝑖 = 2, 3, … , 𝑛:

Compute 𝑣𝑖 = 𝑎𝑖 − ∑𝑖−1
𝑗=1(𝑞𝑗 . 𝑎𝑖 ). 𝑞𝑗

i.e., 𝑣𝑖 = 𝑎𝑖 − (𝑞1 . 𝑎𝑖 ). 𝑞1 − (𝑞2 . 𝑎𝑖 ). 𝑞2 − ⋯ − (𝑞𝑖−1 . 𝑎𝑖 ). 𝑞𝑖−1 .

𝑣
Then find 𝑞𝑖 = ‖𝑣𝑖 ‖.
𝑖
Therefore, from the above steps, it can easily be seen that the vectors 𝑎𝑖 can be expressed as a
linear combination of 𝑞𝑖 as follows:

𝑎1 = ‖𝑣1 ‖𝑞1 = (𝑞1 . 𝑎1 )𝑞1

𝑎2 = (𝑞1 . 𝑎2 )𝑞1 + ‖𝑣2 ‖𝑞2
𝑎3 = (𝑞1 . 𝑎3 )𝑞1 + (𝑞2 . 𝑎3 )𝑞2 + ‖𝑣3 ‖𝑞3
⋮ ⋮ ⋮
𝑎𝑛 = (𝑞1 . 𝑎𝑛 )𝑞1 + (𝑞2 . 𝑎𝑛 )𝑞2 + (𝑞3 . 𝑎𝑛 )𝑞3 + ⋯ + ‖𝑣𝑛 ‖𝑞𝑛
The resulting QR Factorization is as follows:

‖𝑣1 ‖ (𝑞1 . 𝑎2 ) … (𝑞1 . 𝑎𝑛 )

0 ‖𝑣2 ‖ … (𝑞2 . 𝑎𝑛 )
𝐴 = [𝑎1 𝑎2 … 𝑎𝑛 ] = [𝑞1 𝑞2 … 𝑞𝑛 ] [ ] = 𝑄𝑅
⋮ ⋮ ⋱ ⋮
0 0 … ‖𝑣𝑛 ‖
‖𝑣1 ‖ (𝑞1 . 𝑎2 ) … (𝑞1 . 𝑎𝑛 ) 𝑟11 𝑟12 … 𝑟1𝑛
0 ‖𝑣2 ‖ … (𝑞2 . 𝑎𝑛 ) 0 𝑟22 … 𝑟2𝑛
where 𝑄 = [𝑞1 𝑞2 … 𝑞𝑛 ] and 𝑅 = [ ]=[ ⋮ ⋮ ⋱ ⋮ ].
⋮ ⋮ ⋱ ⋮
0 0 … ‖𝑣𝑛 ‖ 0 0 … 𝑟𝑛𝑛

𝟏 𝟐 𝟎
Example 1. Compute the QR factorization for the matrix 𝑨 = [𝟎 𝟏 𝟏].
𝟏 𝟎 𝟏

Solution

1
Step-1: 𝑣1 = 𝑎1 = [0] and ‖𝑣1 ‖ = √12 + 02 + 12 = √2.
1

𝟏
𝑣1 1 1 √𝟐
𝑞1 = = [ 0] = 𝟎
‖𝑣1 ‖ √2 𝟏
1
[√𝟐]

Step-2: 𝑣2 = 𝑎2 − (𝑞1 . 𝑎2 ) 𝑞1

2 1 1
√2
= [1] − 2 [0] = [ 1 ]
√
0 1 −1
 √2
This calculation required that (𝑞1 . 𝑎2 ) = = 1.
√2

Moreover 𝑟22 = ‖𝑣2 ‖ = √12 + 12 + (−1)2 = √3.

𝟏
1 √𝟑
𝑣 1 𝟏
and 𝑞2 = ‖𝑣2 ‖ = [ 1 ]=
2 √3 √𝟑
−1 −𝟏
[√𝟑 ]

Step-3: 𝑣3 = 𝑎3 − (𝑞1 . 𝑎3 ) 𝑞1 − (𝑞2 . 𝑎3 ) 𝑞2

1 0
1 1 1
First, we compute (𝑞1 . 𝑎3 ) = [0] . [1] = (1.0 + 0.1 + 1.1) =
√2 √2 √2
1 1

1 1 0 1
(𝑞2 . 𝑎3 ) = [ 1 ] [1] =
. (1.0 + 1.1 − 1.1) = 0
√3 −1 1 √3

1
0 1 1 0 1 − −1
1 1 1 1 2 1
Therefore, 𝑣3 = [1] − 2 [ 0 ] − 0. [ 1 ] = [1 ] − [0 ] = [ 1 ] = [ 2]
√ √2 √3 2 1 2
1 1 −1 1 1 1
2

1 √6
Also, ‖𝑣3 ‖ = 2 √(−1)2 + 22 + 12 = 2
.

−𝟏
−1 −1 √𝟔
𝑣3 2 1 1 𝟐
and 𝑞3 = = [ 2 ] = 6[ 2 ] =
‖𝑣3 ‖ √6 2 √ √𝟔
1 1 𝟏
[√𝟔 ]

Collecting all of the information we get

𝟏 𝟏 −𝟏
√𝟐 √𝟑 √𝟔
𝟏 𝟐
𝑄 = [𝑞1 |𝑞2 |𝑞3 ] = 𝟎
√𝟑 √𝟔
𝟏 −𝟏 𝟏
[√𝟐 √𝟑 √𝟔 ]
 𝟏
‖𝑣1 ‖ (𝑞1 . 𝑎2 ) (𝑞1 . 𝑎3 ) √𝟐 √𝟐 √𝟐
and 𝑅 = [ 0 ‖𝑣2 ‖ (𝑞2 . 𝑎3 )] = 𝟎 √𝟑 𝟎
0 0 ‖𝑣3 ‖ √𝟔
[𝟎 𝟎 𝟐]

We can easily verify that 𝑄𝑅 = 𝐴

𝟏 𝟏 −𝟏
𝟏
√𝟐 √𝟑 √𝟔 √𝟐 √𝟐
𝟏 𝟐 √𝟐 𝟏 𝟐 𝟎
𝑖. 𝑒. , 𝟎 𝟎 √𝟑 𝟎 = [𝟎 𝟏 𝟏]
√𝟑 √𝟔 𝟏 𝟎 𝟏
𝟏 −𝟏 𝟏 √𝟔
[ 𝟎 𝟎
[√𝟐 𝟐]
√𝟑 √𝟔 ]

Exercise 2.4
1. Find the QR Factorization of the following matrices:
𝟏 −𝟏
(a) [ ]
−𝟏 𝟎
𝟐 𝟏
(b) [ ]
𝟏 𝟏
𝟏 𝟏 𝟎
−𝟏 𝟎 𝟏
(c) [ ]
𝟎 𝟏 𝟏
𝟏 −𝟏 𝟎
𝟏 𝟏 𝟎
−𝟏 𝟎 𝟏
(d) [ ]
𝟎 𝟏 𝟏
𝟎 𝟎 𝟏
𝟏 𝟏 𝟏
𝟏 𝟏 𝟎
(e) [ ]
𝟏 𝟎 𝟎
𝟎 𝟎 𝟎
 2.5. Singular Value Decomposition (SVD)

In this section, we will discuss an extension of the diagonalization theory for symmetric matrices to gen-
eral matrices. The results that we will develop in this section have applications for compression, storage,
and transmission of digitized information and form the basis for many of the best computational algorithms
that are currently available for solving linear systems.

2.5.1. SINGULAR VALUE DECOMPOSITION (SVD):

Suppose A is a m×n matrix. Then the singular value decomposition of A, is of the form A = U ΣV T ,
where
• U is an m × m orthogonal matrix
• Σ is a m × n matrix having singular values of A in the main diagonal and other entries are
zeros
• V is an n × n orthogonal matrix
The columns of U are called the left singular vectors, and those of V are called the right singular
vectors.

Singular Values:
If
√ is any√matrix order √
A m × n and λ1 , λ2 , . . . λn are eigenvalues of AT A then the numbers σ1 =
λ1 , σ2 = λ2 , . . . , σn = λn are called singular values of A.

Note A common convention is to list the singular values in descending order.

 
1 1 0
Example 1. Find the singular value decomposition for the matrix A =
0 0 1
Solution:
The SVD of A is of the form A2x3 = U2×2 2×3 V3×3
T
P

First, compute
   
1 0   1 1 0
1 1 0
AT A =  1 0  =  1 1 0  = B (say)
0 0 1
0 1 0 0 1
The characteristic equation of B is |B − λI| = 0

1−λ 1 0
1 1−λ 0 =0
0 0 1−λ
λ3 − 3λ2 + 2λ = 0
λ(λ − 1)(λ − 2) = 0 ⇒ λ = 0, 1, 2
Hence 0, 1, 2 are the eigenvalues of the matrix B
    
1−λ 1 0 x 0
consider the system  1 1−λ 0   y = 0 
 
0 0 1−λ z 0

Eigenvector corresponding to λ = 2:

1
     
−1 1 0 x 0
put λ = 2 in the system(1) we get  1 −1 0   y  =  0 
0 0 −1 z 0
R2 → R2 + R1
    
−1 1 0 x 0
 0 0 0  y  =  0 
0 0 −1 z 0
−z = 0 ⇒ z = 0
−x + y = 0 ⇒ x = y
choose y = k then x = k
     
x k 1
∴  y  =  k  = k 1 , where k ̸= 0 is a scalar

z 0 0
 
1
Hence  1  is an eigenvector of B corresponding to the eigenvalue λ = 2
0

Eigenvector corresponding to λ = 1    
0 1 0 x 0
put λ = 1 in the system (1) we get  1 0 0   y = 0 
 
0 0 0 z 0
R1 ↔ R2
    
1 0 0 x 0
 0 1 0  y  =  0 
0 0 0 z 0
y = 0, x = 0
Choose z = k
     
x 0 0
∴  y  =  0  = k  0 , where k ̸= 0 is a scalar
z k 1
 
0
Hence  0  is an eigenvector of B corresponding to the eigenvalue λ = 1
1
    
1 1 0 x 0
Eigenvector corresponding to λ = 0 put λ = 0 in the system (1) we get 1 1 0
   y = 0 
 
0 0 1 z 0
R2 → R2 − R1 ↔ R3
         
1 1 0 x 0 1 1 0 x O
 0 0 0  y  =  0   0 0 1  y  =  O 
0 0 1 z 0 0 0 0 z O
z=O
x + y = 0 =⇒ x = −y
Choose y = k then x = −k
     
x −k −1
∴  y  =  k  = k  1  , where k ̸= 0 is a scalar
z 0 0
  
−1
Hence  1  is an eigenvector of B corresponding to the eigenvalue λ = 0
0
The eigenvectors corresponding to λ = 2, 1, 0, respectively, are
     
1 0 −1
X1 =  1  , X2 =  0  and X3 =  1 
0 1 0
Since X1T X2 = X2T X3 = X3T X1 = 0
Therefore, the eigenvectors X1 , X2 , X3 are orthogonal.
The normalized eigenvectors are
 √     √ 
1/√2 0 −1/√ 2
v1 =  1/ 2  , v2 =  0  and v3 =  1/ 2 
0 1 0
√ √ √
and the singular values of B are σ1 = 2, σ2 = 1 and σ3 = 0
 √ √ 
1/√2 0 −1/√ 2  √ 
2 0 0
Thus V = v1 v2 v3 = 1/ 2 0 1/ 2 and Σ =
 
 
0 1 0
0 1 0
To find U , we compute
 √ 
 1/ 2
1 1 0  √ 
  
1 1 1
u1 = Av1 = √ 1/ 2 =
σ1 2 0 0 1 0
0
 
  0  
1 1 1 1 0   0
and u2 = Av2 = 0 =
σ2 1 0 0 1 1
1

Therefore,
 
  1 0
U= u1 u2 =
0 1
We can easily verify that, U ΣV T = A
√ √
 √
 
  1/ 2 1/ 2 0
1 0 2 0 0 
U ΣV T = 0√ 0√ 1 
0 1 0 1 0
−1/ 2 1/ 2 0
    
1 0 1 1 0 1 1 0
= = =A
0 1 0 0 1 0 0 1

Exercise 2.5

1. Find the singular value decomposition for the following matrices

 
1 1 0
(a) A =
0 1 1
 
1 2
(b) A =  2 3 
3 4
  
1 1 0
(c) A =
0 0 1