Computational Methods for Engineering

Korea Institute of Science and Technology

Imaging Media Research Center
Jinwook Kim
2
Contents

1 Review on Linear Algebra 7

1.1 Vector Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.1 Basic notations . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.2 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.3 Linear independence . . . . . . . . . . . . . . . . . . . . . 8
1.1.4 Subspace . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.5 Rank of a matrix . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.6 Determinant of a square matrix . . . . . . . . . . . . . . . 9
1.1.7 Linear equations . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.8 Inner products and norms . . . . . . . . . . . . . . . . . . 11
1.2 Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.1 Linear operators . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.2 Similarity transform . . . . . . . . . . . . . . . . . . . . . 12
1.2.3 Eigenvalues and eigenvectors . . . . . . . . . . . . . . . . 13
1.2.4 Orthogonal projections . . . . . . . . . . . . . . . . . . . . 13
1.2.5 Quadratic forms . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.6 Singular value decomposition . . . . . . . . . . . . . . . . 15

2 Solving Linear System 19

2.1 Linear System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.1 Glimpse of a numerical method to linear system . . . . . 19
2.1.2 Existence of the solution . . . . . . . . . . . . . . . . . . . 20

3
4 CONTENTS

2.1.3 Least squares problem . . . . . . . . . . . . . . . . . . . . 20

2.1.4 Condition and Stability . . . . . . . . . . . . . . . . . . . 22
2.2 Direct Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2.1 LU factorization . . . . . . . . . . . . . . . . . . . . . . . 24
2.2.2 Solving modified problem . . . . . . . . . . . . . . . . . . 26
2.2.3 Cholesky factorization . . . . . . . . . . . . . . . . . . . . 27
2.2.4 QR factorization . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.5 Computational complexity . . . . . . . . . . . . . . . . . . 33
2.3 Iterative Approach . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3.2 Basics of iterative approach . . . . . . . . . . . . . . . . . 35
2.3.3 Jacobi iteration . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.4 Gauss-Seidel iteration . . . . . . . . . . . . . . . . . . . . 36
2.3.5 Successive over-relaxation . . . . . . . . . . . . . . . . . . 37

3 Representing Rotation 45
3.1 Rotation Group . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Properties of Rotation Matrices . . . . . . . . . . . . . . . . . . . 45
3.3 Exponential Map . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 Euler Angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.5 Quaternions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.5.1 Structure of Quaternions . . . . . . . . . . . . . . . . . . 51
3.5.2 Unit Quaternions . . . . . . . . . . . . . . . . . . . . . . . 52
3.5.3 Interpolation of Quaternioins . . . . . . . . . . . . . . . . 54

4 Introduction to Optimization 59
4.1 Di↵erential Calculus . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2 Unconstrained Optimization . . . . . . . . . . . . . . . . . . . . . 60
4.2.1 Conditions for local minima . . . . . . . . . . . . . . . . . 61
4.2.2 Gradient method . . . . . . . . . . . . . . . . . . . . . . . 62
4.2.3 Steepest descent method . . . . . . . . . . . . . . . . . . . 63
CONTENTS 5

4.2.4 Newton’s method . . . . . . . . . . . . . . . . . . . . . . . 63

4.2.5 Conjugate direction methods . . . . . . . . . . . . . . . . 64
6 CONTENTS
Chapter 1

Review on Linear Algebra

1.1 Vector Space

1.1.1 Basic notations

A column n vector is an array of n numbers, denoted by

2 3
a1
6 7
6 7
6 a2 7
a=6 7
6 .. 7 .
6.7
4 5
an

Denote the set of real numbers by R, and the set of n vectors with real compo-
nents by Rn .

1.1.2 Properties

The vector space with the operations of a vector addition and a scalar multipli-
cation holds the following properties:

1. commutative : a + b = b + a

2. associative : (a + b) + c = a + (b + c)

7
8 CHAPTER 1. REVIEW ON LINEAR ALGEBRA

3. zero vector : 90, such that a + 0 = 0 + a = a

4. distributive: ↵(a + b) = ↵a + ↵b

, where a, b, c 2 Rn , ↵ 2 R.

1.1.3 Linear independence

Definition 1. A set of vectors {a1 , · · · , an } is said to be linearly independent

if the equality
↵1 a1 + ↵2 a2 · · · + ↵n an = 0

implies that all coefficients ↵i , i = 1, · · · , n are equal to zero. If it is not linearly

independent, the set is linearly dependent.

1.1.4 Subspace

Definition 2. A subset V of Rn is called a subspace of Rn if V is closed under

the operations of vector addition and scalar multiplication.

Definition 3. Let a1 , a2 , · · · , ak be arbitrary vectors in Rn . The set of all their

linear combinations is called the span of a1 , a2 , · · · , ak and is denoted
k
X
span(a1 , a2 , · · · , ak ) = ↵i ai | ↵i 2 R .
i=1

Definition 4. Given a subspace V , any set of linearly independent vectors

{ai , a2 , · · · , ak } ⇢ V such that V = span(a1 , a2 , · · · , ak ) is referred to as a basis
of the subspace V . All bases of a subspace V contain the same number of
vectors. This number is called the dimension of V , denoted dim V .

1.1.5 Rank of a matrix

We will represent an m ⇥ n matrix A by

2 3
a11 a12 · · · a1n
6 7
6 7
6 a21 a22 · · · a2n 7
A=6 6 .. .. ..
7
.. 7 .
6 . . . . 7
4 5
am1 am2 ··· amn
1.1. VECTOR SPACE 9

Let us denote the kth column of A by ak , that is,

2 3
a1k
6 7
6 7
6 a2k 7
ak = 6 7
6 .. 7 .
6 . 7
4 5
amk

Definition 5. The maximal number of linearly independent columns of A is

called the rank of the matrix A.

Lemma 1. The rank of a matrix A is invariant under the following operations:

1. Multiplication of the columns of A by non-zero scalars

2. Interchange of the columns

3. Addition to a given column a linear combination of other columns

1.1.6 Determinant of a square matrix

Definition 6. The determinant of a square matrix is a function of its columns

and has the following properties:
2 3
p p p
6 7
6 7
1. The determinant of the matrix A = 6a1 a2 ··· an 7 is a linear func-
4 5
p p p
tion of each column; that is
h i
det a1 , · · · , ak 1 , ↵ak + a0k , ak+1 , · · · , an
h i h i
=↵ det a1 , · · · , ak 1 , ak , ak+1 , · · · , an + det a1 , · · · , ak 0
1 , ak , ak+1 , · · · , an ,

for each ↵, 2 R, a0k 2 Rn .

2. If for some k we have ak = ak+1 , then

h i h i
det A = det a1 , · · · , ak , ak+1 , · · · , an = det a1 , · · · , ak , ak , · · · , an = 0.
10 CHAPTER 1. REVIEW ON LINEAR ALGEBRA

3. Let 2 3
1 0 ··· 0
6 7
h i 6 7
60 1 ··· 07
I n = e 1 , · · · , en =6
6 .. .. ..
7
.. 7 ,
6. . . .7
4 5
0 0 ··· 1

where {e1 , · · · , en } is the natural basis for Rn . Then

det In = 1.

Lemma 2. The determinant of A satisfies the following properties:

h i h i
1. det a1 , · · · , ak 1 , ak + ↵aj , ak+1 , · · · , aj , · · · , an = det a1 , · · · , an .
h i h i
2. det a1 , · · · , ak , ak+1 , · · · , an = det a1 , · · · , ak+1 , ak , · · · , an .

1.1.7 Linear equations

Suppose we are given m equations with n unknowns as

a11 x1 + a12 x2 + · · · + a1n xn = b1

a21 x1 + a22 x2 + · · · + a2n xn = b2

..
.

am1 x1 + am2 x2 + · · · + amn xn = bm .

We can represent the above set of equations as a vector equation

x1 a1 + x2 a2 + · · · xn an = b,

where aj = [a1j , a2j , · · · , amj ]T , b = [b1 , b2 , · · · , bm ]T .

Associated with the above system of equations are the following matrices
2 3
p p p
6 7
6 7
A = 6a1 a2 · · · an 7 ,
4 5
p p p
1.1. VECTOR SPACE 11

and an augmented matrix

2 3
p p p p
6 7
6 7
[ A, b ] = 6a1 a2 ··· an b7 .
4 5
p p p p

We can also represent the above system of equations as

Ax = b,

where x = [x1 , x2 , · · · , xn ]T .

Theorem 1. The system of equations Ax = b has a solution if, and only if,

rank(A) = rank([ A, b ]).

1.1.8 Inner products and norms

Definition 7. The Euclidean inner product is defined by

X
hx, yi = xi yi = xT y.

Definition 8. A general vector norm is a real valued function of a vector, which

satisfies the following properties:

1. Positivity: kxk 0, kxk = 0 if, and only if, x = 0

2. Homogeneity: krxk = |r|kxk, r 2 R

3. Triangle inequality: kx + yk  kxk + kyk.

Definition 9. The Euclidean norm of a vector x is defined as

p
kxk = hx, xi.

Theorem 2. For any two vectors x and y in Rn , the Cauch-Schwarz in-

equality
|hx, yi|  kxkkyk

holds. Furthermore, equality holds if, and only if, x = ↵y for some ↵ 2 R.
12 CHAPTER 1. REVIEW ON LINEAR ALGEBRA

1.2 Transformations

1.2.1 Linear operators

Definition 10. A function L : Rn ! Rm is called a linear operator if

1. L(ax) = aL(x) 8x 2 Rn and a 2 R

2. L(x + y) = L(x) + L(y) 8x, y 2 Rn .

If one chooses natural bases for Rn and Rm , then the linear operator L can
be represented by a matrix A as

L(x) = Ax.

1.2.2 Similarity transform

Let {e1 , e2 , · · · , en } and {e01 , e02 , · · · , e0n } be two bases for Rn and the matrix
T relate these bases as follows:
2 3 1 2 3
p p p p p p
6 7 6 7
6 0 7 6 7
T = 6e1 e02 · · · e0n 7 6e1 e2 ··· en 7 .
4 5 4 5
p p p p p p

We call T the transformation matrix from {e1 , e2 , · · · , en } to {e01 , e02 , · · · , e0n }.

Note that the ith column of T is the vector of coordinates of ei with respect to
the basis {e01 , e02 , · · · , e0n }.
Consider a linear operator L : Rn ! Rn , and let A be its matrix repre-
sentation with respect to {e1 , e2 , · · · , en }, and B its matrix representation
with respect to {e01 , e02 , · · · , e0n }. Then y = Ax or y 0 = Bx0 , where x and
x0 are di↵erent representations of an identical vector with respect to bases
{e1 , e2 , · · · , en } and {e01 , e02 , · · · , e0n }. Similarly, y and y 0 are di↵erent
representations of an identical vector with respect to those two bases. Now we
have a similarity transform of the linear operator L between two bases as

1
A=T BT.
1.2. TRANSFORMATIONS 13

1.2.3 Eigenvalues and eigenvectors

Definition 11. The characteristic polynomial of A is defined as

( ) , det( I A) = 0.

Theorem 3. Suppose the characteristic polynomial of A has n distinct roots

1, 2, · · · , n. Then, there exist n linearly independent vectors v1 , v2 , · · · , vn
such that
Avi = i vi , i = 1, 2, · · · , n.

The roots of the characteristic polynomial, i and their associated vectors vi are
called eigenvalues and eigenvectors, respectively.

Lemma 3. Consider a basis formed by a linearly independent set of eigenvectors

v1 , v2 , · · · , vn . Then the similarity transform of the matrix A with respect to this
basis has a diagonal form

1
V AV = diag{ 1, 2, · · · , n },

h p p p
i
where V = v1 v2 ··· vn .
p p p

Note that all matrix can not be diagonalized. If the eigenvalues are not all
distinct, we may find only m(< n) distinct eigenvectors.

Corollary 1. All eigenvalues of a symmetric matrix are real.

Corollary 2. Any real symmetric n ⇥ n matrix has a set of n eigenvectors that

are mutually orthogonal.

Corollary 3. Eigenvalues of a matrix are invariant under any similarity trans-

form.

1.2.4 Orthogonal projections

V is a subspace of Rn if x1 , x2 2 V ) ↵x1 + x2 2 V, 8↵, 2 R. Furthermore,

the dimension of a subspace V is equal to the maximum number of linearly
independent vectors in V .
14 CHAPTER 1. REVIEW ON LINEAR ALGEBRA
*6
x
x Px

-
Px

Figure 1.1: Geometric meaning of an orthogonal projector.

Definition 12. The orthogonal complement of V , denoted V ? , defined as fol-

lows:
V ? = x | hv, xi = 0 8v 2 V .

Lemma 4. Together, V and V ? span Rn in the sense that every vector x 2 Rn

can be represented uniquely as

x = x1 + x2 ,

where x1 2 V, x2 2 V ? .

We call the representation the orthogonal decomposition of x onto the sub-

spaces V and V ? . Also we say that x1 and x2 are orthogonal projections of x
onto the subspaces V and V ? , respectively.

Definition 13. A linear transformation P is an orthogonal projector onto V if,

for all x 2 Rn , we have

Px 2 V and x P x 2 V ?.

Note that an orthogonal projector P is not an orthogonal matrix.

Definition 14. The range and nullspace of A are defined as

R(A) = Ax | x 2 Rn

N (A) = x 2 Rn | Ax = 0 .

Theorem 4. (The Fundamental Theorem of Linear Algebra) Suppose A 2

Rm⇥n . Then

1. N (A)? = R(AT ) and R(AT )? = N (A);

1.2. TRANSFORMATIONS 15

2. N (AT )? = R(A) and R(A)? = N (AT ).

Corollary 4. A matrix P is an orthogonal projector onto the subspace V =

R(P ) if, and only if, P 2 = P = P T .

Corollary 5. If a matrix P is an orthogonal projector, then also I P is an

orthogonal projector.

1.2.5 Quadratic forms

A quadratic form f : Rn ! R is a function

f (x) = xT Qx,

where Q is an n ⇥ n real symmetric matrix.

Definition 15. A quadratic form xT Qx is said to be positive definite if xT Qx >

0 8x 6= 0. Is is positive semi-definite if xT Qx 0 8x. Similarly, we define the
quadratic form to be negative definite or negative semi-definite, if xT Qx <
0 8x 6= 0, or xT Qx  0 8x, respectively.

Theorem 5. A symmetric matrix Q is positive definite (or positive semi-

definite) if, and only if, all eigenvalues of Q are positive (or non-negative).

1.2.6 Singular value decomposition

Theorem 6. Every m ⇥ n real matrix A of a rank r can be decomposed into

the form 2 3
⌃ 0
A=U4 5V T,
0 0

where U T U = Im⇥m , V T V = In⇥n , and ⌃ = diag{ 1, 2, · · · , r} with i >

0, i = 1, · · · , r.

Corollary 6. The rank of a matrix A is same to the number of non-zero singular

values.
16 CHAPTER 1. REVIEW ON LINEAR ALGEBRA

1u1
v2

A
v1 2u2

Figure 1.2: A maps a unit circle to an ellipse.

Corollary 7. If rank(A) = r, then the set of column vectors of U associated

with non-zero singular values spans R(A) and the set of column vectors of V
associated with zero singular values spans N (A).

Corollary 8. Matrix norms of A is characterized by its singular values.

X
kAk2F = Tr(AT A) = Tr(AAT ) = 2
i (A)
i
q
kAk2 = T A)
max (A) = max (A

Let us discuss the geometric meaning of SVD. Using the singular value de-
composition, a matrix A can be represented as
X
T
A= i ui v i .
i

Hence A maps vi to i ui . That is, A maps a unit hyper-sphere to a hyper-ellipse

which principal axes are determined by singular values of A.

Problem Set
1. Suppose that {v1 , v2 , · · · , vn } is a basis for Rn . Show that if A is an n ⇥ n
invertible matrix, then {Av1 , Av2 , · · · , Avn } is also a basis for Rn .

2. Consider a linear operator T : R3 ! R2 ,

T (x, y, z) = (3x z, 4y + 2z).

1.2. TRANSFORMATIONS 17

Let A1 , A2 be two bases for R3 and B1 , B2 be two bases for R2 .

A1 = {[1, 1, 1]T , [1, 1, 0]T , [1, 0, 0]T }, B1 = {[0, 1]T , [1, 2]T }

A2 = {[2, 1, 3]T , [1, 0, 2]T , [0, 1, 0]T }, B2 = {[1, 1]T , [2, 1]T }

(a) Find a matrix representation of T1 with respect to A1 and B1 .

(b) Find transformation matrices from A1 to A2 and from B1 to B2 .

(c) Using the similarity transform, find a matrix representation of T2

with respect to A2 and B2 .

3. Let w be a column vector in Rn , and let A = I + wwT , where I denotes

the identity matrix. Find the eigenvalues and eigenvectors of A. Is A
diagonalizable?

4. If all eigenvalues of a symmetric matrix A are distinct, show that

X 1
1
(sI A) = vi viT ,
i
s i

where i and vi represent the ith eigenvalue and eigenvector respectively.

5. The exponential of a symmetric real matrix A is defined as

X1
1 n
eA = A .
n=0
n!

Show that det(eA ) = eTr(A) .

[ Hint: Use matrix properties of Tr(AB) = Tr(BA), det(AB) = det(A) det(B).]

[ Note: Recall that the Taylor series expansion of the exponential function
P1 n
gives us ex = n=0 xn! for x 2 R. Hence we can regard the above defini-
tion as a matrix version of the exponential function. In this problem, we
assume that A is real and symmetric. However the above definition and
property hold for an arbitrary square matrix.]

6. Show that eigenvalues of any orthogonal projector are 0 or 1.

18 CHAPTER 1. REVIEW ON LINEAR ALGEBRA

7. Consider a linear transform T : R2 ! R3 ,

T (x, y) = (x y, 2x 4y, 2x + 2y) .

Find a basis for R(T )? .

8. Show that AT A is positive semi-definite for all matrix A 2 Rm⇥n .

9. Consider subspaces U and V of R4 ,

2 3 2 3 2 3 2 3 2 3
2 5 1 3 0
6 7 6 7 6 7! 6 7 6 7!
6 7 6 7 6 7 6 7 6 7
6 0 7 6 17 617 6 17 6 0 7
U = span 6 7 6 7 6 7
6 7 , 6 7 , 6 7 , V = span
6 7,6 7 .
6 7 6 7
6 17 6 37 607 6 27 6 1 7
4 5 4 5 4 5 4 5 4 5
3 7 2 4 1

Find dim(U ? \ V ? ), where U ? \ V ? ⌘ {x 2 R4 |x 2 U ? , x 2 V ? }.

10. Find directions of principal axes and their magnitude of the ellipse

5x2 4xy + 8y 2 = 1.
Chapter 2

Solving Linear System

2.1 Linear System

2.1.1 Glimpse of a numerical method to linear system

By applying row manipulations successively, we can convert a given linear sys-

tem into the form of
2 32 3
2 3
a011 a012 ··· a01n x1 b0
6 7 6 7 6 17
6 76 7 6 7
6 0 a022 ··· a02n 7 6 x2 7 6 b02 7
6 76 7 6 7
6 .. .. .. 7 6 .. 7 = 6 .. 7 .
6 0 . . . 7 6 7 6 7
4 54 . 5 4 . 5
0 0 ··· a0nn xn b0n

It is called the Gaussian elimination procedure. For the system with a full
rank matrix A, a solution can be computed efficiently by the back-substitution
procedure.

1 ⇣ 0 ⌘
n
X
xi = b a0ij xj , i = n, · · · , 1.
a0ii i j=i+1

Note that the above back-substitution holds only when the diagonal entries a0ii
are not zero.

19
20 CHAPTER 2. SOLVING LINEAR SYSTEM

for i = n:1 do
x(i) = b(i)
for j = i+1:n do
x(i) = x(i) - a(i,j) * x(j)
end
x(i) = x(i) / a(i,i)
end
Algorithm 1: Back-substitution

2.1.2 Existence of the solution

If there exists unique solution to a given linear system, then it is called an

exactly-determined system. If there exist infinitely many solutions or no solution,
then it is called an under-determined or an over-determined system, respectively.
Any linear system Ax = b falls into the above categorization as follows:

A b x

rank(A) = n b 2 span(A) exactly-determined

rank(A) < n b 2 span(A) under-determined
rank(A)  n b2
/ span(A) over-determined

, where A 2 Rm⇥n .

2.1.3 Least squares problem

Definition 16. A least squares solution of Ax = b is to find x̃ that minimizes

a norm of the residual vector r = b Ax, that is,

x̃ = arg min ||b Ax||2 .

x

Assume that A 2 Rm⇥n . In the most usual case m n and rank(A) = n,

the solution to the linear least squares problem is unique, and the problem it
also referred to as finding a least squares solution to an over-determined system
of linear equations.
2.1. LINEAR SYSTEM 21

Theorem 7. A least squares solution to an over-determined system Ax = b is

given by x̃ = A† b, where A† = (AT A) 1
AT .

Proof. Minimizing the objective function f (x) = (b Ax)T (b Ax) = bT b

2xT AT b + xT AT Ax,

@f (x)
= 2(AT Ax AT b) = 0,
@x

x̃ = (AT A) 1
b.

Geometric interpretation of the least square problem is to find the point

Ax̃ 2 R(A) closest to b, which means that x̃ will minimize the magnitude of
residual kr = Ax̃ bk. It happens when hAx̃, ri = 0.

xT AT (Ax b) = 0

xT (AT Ax AT b) = 0 , 8x 6= 0,

x = (AT A) 1
b.

When m < n and rank(A) = m, there are an infinite number of solutions x

which exactly satisfy b Ax = 0. In this case, it is often useful to find the unique
solution x which minimizes ||x||2 , and the problem is referred to as finding a
minimum norm solution to an under-determined system of linear equations.

Theorem 8. A minimum norm solution to an under-determined system Ax = b

is given by x̃ = A† b, where A† = AT (AAT ) 1
.

Proof. Let xm = AT (AAT ) 1

b be the minimum norm solution. Think of x 6= xm
satisfying Ax = b. Then

(x xm )T xm = (x xm )T AT (AAT ) 1
b

= (Ax Axm )T (AAT ) 1

b

= 0.
22 CHAPTER 2. SOLVING LINEAR SYSTEM

Therefore

kxk2 = kx xm + xm k2

= kx xm k2 + kxm k2 + 2(x xm ) T xm

= kx xm k2 + kxm k2 kxm k2 .

In general, when we may have rank(A) < min(m, n) – in other words, A may
be rank-deficient – we seek the minimum norm least squares solution x which
minimizes both ||x||2 and ||b Ax||2 .

Corollary 9. Given ↵ 2 R, a solution to arg min ||Ax b||2 + ↵2 ||x||2 is

x
x̃ = (AT A + ↵2 I) 1
AT b.

It is called the Tikhonov regularization, which is one of the most commonly

used methods to regularize ill-posed problems. Note that it improves the condi-
tioning of the problem, thus enabling a numerical solution even when cond(A)
is quite big considering a machine precision.

2.1.4 Condition and Stability

In a numerical approach, existence and uniqueness of the solutions may be

problematic due to the machine precision. Consider a following problem:
2 32 3 2 3
1 2 x 4
4 54 5 = 4 5 (2.1)
2 3.9 y 7.9

It is obvious that the solution is [ x y ]T = [ 2 1 ]T . Now let us perturb

matrices and see how the solution changes.
2 32 3 2 3 2 3 2 3
1 2 x 4.1 x 3.9
4 54 5 = 4 5 ) 4 5=4 5
2 3.9 y 7.8 y 4
2 32 3 2 3 2 3 2 3
1 2 x 4 x 6
4 54 5 = 4 5 ) 4 5=4 5
2 4.1 y 7.9 y 1
2.1. LINEAR SYSTEM 23

A small perturbation in A or b leads to a relatively huge change in the solution

x. To investigate what is the reason of instability, let us recall a matrix norm
first.

Definition 17. The 2-norm of an m ⇥ n matrix A is defined as follows:

kAxk
kAk2 , maxn .
x2R kxk
x6=0

Theorem 9. A matrix 2-norm is the maximum of singular values of the matrix.

P T
P T
P
Proof. Since A = i ui vi , Ax = i ui v i x = ↵i ui , where ↵i = i hvi , xi.
P P P
Then kAxk2 = ↵i2 = 2 2 2
i hvi , xi  max hvi , xi2 = max 2
kxk2 . Therefore
kAk2 = max kAxk
kxk = max (A).

Now let us consider how the perturbation of b a↵ects the change of solution
x by comparing Ax = b with A(x + x) = b + b.

1
x=A b , b = Ax
1
k xk  kA kk bk , kbk  kAkkxk
k xk k bk max k bk
)  kA 1
kkAk = .
kxk kbk min kbk

Definition 18. A condition number of a nonsingular matrix A is defined by

max
cond(A) = .
min

Similarly, consider a case of (A + A)(x + x) = b.

1
x=A Ax
1 k Ak
k xk  kA kk Akkxk = cond(A) kxk
kAk
k xk k Ak
 cond(A) .
kxk kAk

Definition 19. If small changes either in b or in A lead to small changes in

x, the linear system Ax = b is called well-conditioned. Otherwise, it is called
ill-conditioned.
24 CHAPTER 2. SOLVING LINEAR SYSTEM
⇥1 2
⇤
Back to the previous example, ( 2 3.9 ) = {4.9203, 0.0203}. So the con-
dition number is higher than 242.38, which means that a small error as in the
example might be magnified a lot. In practice, most errors do not get magnified
so much like the example. However, once you encounter the disaster, it is always
too late. Predict first, and be cautious about it!

2.2 Direct Approach

2.2.1 LU factorization

Gaussian elimination process is equivalent to multiplying A by a sequence of

lower-triangular matrices Lk on the left to get an upper-triangular matrix.

Ln 1 Ln 2 · · · L1 A = U

For a general n ⇥ n matrix, let xk denote the kth column of the matrix at
the beginning of step k. Then Lk can be chosen such that it transforms all the
entries below the (k + 1)th element of xk to be zero.
2 32 3 2 3
1 x1k x1k
6 76 . 7 6 . 7
6 .. 76 . 7 6 . 7
6 . 76 . 7 6 . 7
6 76 7 6 7
6 76 7 6 7
6 1 7 6 xkk 7 6xkk 7
Lk xk = 6
6
76
76
7=6
7 6
7,
7
6 lk+1,k 1 7 6xk+1,k 7 6 0 7
6 76 7 6 7
6 .. .. 76 . 7 6 . 7
6 . . 7 6 .. 7 6 .. 7
4 54 5 4 5
lnk 1 xnk 0

where
xjk
ljk = .
xkk
Lemma 5. Let lk = [0, · · · , 0, lk+1,k , · · · , lnk ]T . Then

Lk = I lk eTk

Lk 1 = I + lk eTk

Lk 1 Lk+1
1
= I + lk eTk + lk+1 eTk+1 .
2.2. DIRECT APPROACH 25

Hence we can construct L as follows:

2 3
1
6 7
6 7
6 l21 1 7
6 7
1 1 1 6 7
L = L1 L2 · · · Ln 1 = 6 l31 l32 1 7.
6 7
6 .. .. .. .. 7
6 . . . . 7
4 5
ln1 ln2 ··· ln,n 1 1

initialization: U = A, L = I
for k = 1 : n 1 do
for j = k + 1 : n do
l(j, k) = u(j, k)/ u(k, k)
u(j, k : n) = u(j, k : n) l(j, k) ⇤ u(k, k : n)
end
end
Algorithm 2: LU factorization

A solution of Ax = b using LU factorization can be computed by solving two

triangular systems.

1. Solve Ly = b. (forward-substitution)

2. Solve U x = y. (back-substitution)

Consider a case when akk is near to zero. Obviously it may induce a large
round-o↵ error because the Gaussian elimination requires a division of elements
by akk . But there is no inherent reason to eliminate against the kth row. We
introduce a permutation matrix P which is simply an identity matrix with its
rows or columns permuted. A permutation matrix is useful to permute rows or
columns of a matrix.
h1 0 0i
For example, if P = 001 , P A is A with the second and third rows per-
010
muted and AP is A with second and third columns permuted. Note that a
permutation matrix is not symmetric in general.
26 CHAPTER 2. SOLVING LINEAR SYSTEM

Now we modify Gaussian elimination using the permutation matrix to avoid

zero-division. So the kth step of Gaussian elimination with pivoting can be
summed up as
2 3 2 3 2 3
⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥
6 7 6 7 6 7
6 7 6 7 6 7
6 ⇥ ⇥ ⇥ ⇥7 6 xik ⇤ ⇤ ⇤7 6 xik ⇥ ⇥ ⇥7
6 7 6 7 6 7
6 7 Pk 6 7 Lk 6 7
6 ⇥ ⇥ ⇥ ⇥7 ! 6 ⇥ ⇥ ⇥ ⇥7 ! 6 0 ⇤ ⇤ ⇤7.
6 7 6 7 6 7
6 7 6 7 6 7
6 xik ⇥ ⇥ ⇥7 6 ⇤ ⇤ ⇤ ⇤7 6 0 ⇤ ⇤ ⇤7
4 5 4 5 4 5
⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ 0 ⇤ ⇤ ⇤

After n 1 steps, A is transformed into an upper-triangular matrix U :

Ln 1 Pn 1 · · · L 2 P2 L 1 P1 A = U

Now we reorder the elementary operations as follows:

L0n 1 · · · L02 L01 Pn 1 · · · P2 P1 A,

where L0k = Pn 1 · · · Pk+1 Lk Pk+1 · · · Pn 1. Note that L0k has the same structure
as Lk , that is, L0k is simply Lk with non trivial entries permuted. This property
is important because it makes L0n 1 · · · L01 lower triangular.
Finally we have
P A = LU,
1 1 1
where L = L0 1 L0 2 · · · L0 n 1 and P = Pn 1 · · · P2 P1 . This process is called
the partial pivoting and one of the most popular criteria to choose the pivot is
to find the largest entries from xk .

2.2.2 Solving modified problem

Assume that we have solved Ax = b already using LU factorization. Then

solving slightly modified system of equations can be performed efficiently using
the previous results.
1 1
Theorem 10. (Woodbury formula) If there exist A and C , and A + U CV
has full rank,

1 1 1 1 1 1 1
(A + U CV ) =A A U (C +VA U) VA ,
2.2. DIRECT APPROACH 27

initialization: U = A, L = I, P = I
for k = 1 : n 1 do
select i k maximizing |u(i, k)|
u(k, k : n) u(i, k : n)
l(k, 1 : k 1) l(i, 1 : k 1)
p(k, :) p(i, :)
for j = k + 1 : n do
l(j, k) = u(j, k)/ u(k, k)
u(j, k : n) = u(j, k : n) l(j, k) ⇤ u(k, k : n)
end
end
Algorithm 3: LU factorization with partial pivoting

where A 2 Rn⇥n , U 2 Rn⇥k , C 2 Rk⇥k , V 2 Rk⇥n . In the special case where

C is the 1 ⇥ 1 unit matrix, this identity reduces to the Sherman-Morrison
formula:
1
(A + uv T ) 1
=A 1
1 ui
A 1
u vT A 1
,
1 + hv, A
where u, v 2 Rn .

Using the Sherman-Morrison formula, the solution of (A + uv T )x = b is

x = (A + uv T ) 1
b
1 1 1
=A b A uv T A 1
b
1 + hv, A 1 ui
1 hv, A 1 bi 1
=A b A u.
1 + hv, A 1 ui
Note that we don’t need an explicit inversion of (A + uv T ) to solve the modified
problem.

2.2.3 Cholesky factorization

Theorem 11. If A is real and positive definite, then there exists a unique
factorization of A = RT R, where R is an upper triangular matrix with strictly
positive diagonal entries. It is called the Cholesky factorization.
28 CHAPTER 2. SOLVING LINEAR SYSTEM

Consider a symmetric matrix

2 3
1 wT
A=4 5.
w K

Applying the Gaussian elimination gives us

2 32 3
1 0 1 wT
A=4 54 5.
w I 0 K wwT

To maintain symmetry, we apply another elementary transformation to the RHS

of the above as follows:
2 32 32 3
1 0 1 0 1 wT
A=4 54 54 5.
w I 0 K wwT 0 I

The idea behind Cholesky factorization is to continue this process. We can

generalize this for any a11 > 0 as
2 3
a11 wT
A=4 5
w K
2 32 32 3
1
↵ 0 1 0 ↵ ↵ w5
=4 54 54
1 1 T
↵w I 0 K a11 ww 0 I

= R1T A1 R1 ,

p 1 T
where ↵ = a11 . If the first diagonal entry of K a11 ww is positive, we can
factor it as A1 = R2T A2 R2 and hence A = R1T R2T A2 R2 R1 . By repeating this
process until An = I, we finally have

A = RT R,

where R = Rn · · · R2 R1 .

Definition 20. Let A be an n ⇥ n matrix. A is diagonally dominant if for all

i = 1, 2, · · · , n
X
|aii | |aij |
j6=i

and strictly diagonally dominant if the inequality is strict for all i.

2.2. DIRECT APPROACH 29

Proposition 1. A matrix A is symmetric positive definite i↵ all the diagonal

entries of A are positive and A is strictly diagonally dominant.

Lemma 6. If a matrix A is symmetric positive definite, then any principal

submatrix of A is also symmetric positive definite.

P
From aij = k rki rkj , we can derive rij as follows:
8 q Pi 1 2
< aii k=1 rki if i = j
rij = ⇣ ⌘
: 1 aij Pi 1 rki rkj if i < j .
rii k=1

for i = 1 : n do
r(i, i : n) = a(i, i : n) r(1 : i 1, i)T ⇤ r(1 : i 1, i : n)
p
r(i, i : n) = r(i, i : n)/ r(i, i)
end
Algorithm 4: Cholesky factorization

In the above algorithm, elements with zero indices are assumed to be zero.
A solution of Ax = b using the Cholesky factorization can be computed by
solving two triangular systems.

1. Solve RT y = b. (forward-substitution)

2. Solve Rx = y. (back-substitution)

2.2.4 QR factorization

Gram-Schmidt orthogonalization

Consider a linearly independent set S = {a1 , a2 , · · · , an }. The Gram-Schmidt

orthogonalization is a process generating an orthogonal set S 0 = {qi , q2 , · · · , qn }
that spans the same subspace as S.
30 CHAPTER 2. SOLVING LINEAR SYSTEM

The Gram-Schmidt orthogonalization works as follows:

v1
v1 = a1 q1 =
||v1 ||
v2
v2 = a2 ha2 , q1 iq1 q2 =
||v2 ||
.. ..
. .
n
X1 vn
vn = an hqi , an iqi qn =
i=1
||vn ||

Theorem 12. Any matrix A 2 Rm⇥n can be factorized as A = QR, where

Q 2 Rm⇥n and R 2 Rn⇥n are orthogonal and upper triangular, respectively. It
is called the reduced QR factorization.

Corollary 10. Every A 2 Rm⇥n (m n) of full rank has a unique reduced QR

factorization with rii > 0.

for i = 1 : n do
vi = ai
for j = 1 : i 1 do
rji = hqj , ai i
vi = vi rji qj
end
rii = ||vi ||
vi
qi =
rii
end
Algorithm 5: Classical Gram-Schmidt orthogonalization

A full QR factorization of A 2 Rm⇥n appends an additional (m n) or-

thonormal columns to Q to make it an m ⇥ m orthonormal matrix.
To deal with instability of the classical Gram-Schmidt orthogonalization,
Pi 1 T
consider a series of orthogonal projectors Pi = I j=1 qj qj such as

P1 a 1 P2 a 2 Pn a n
q1 = , q2 = , ··· , qn = .
||P1 a1 || ||P2 a2 || ||Pn an ||
2.2. DIRECT APPROACH 31

Note that hqi , Pi ai = 0 for all a 2 Rm , which means Pi 2 Rm⇥m is an orthogonal

projector of rank m (i 1) that projects Rm onto the space orthogonal to
span(q1 , q2 , · · · , qi 1 ). It is known that this process introduces smaller round-o↵
errors than the classical one.

for i = 1 : n do
vi = ai
end
for i = 1 : n do
rii = ||vi ||
vi
qi =
rii
for j = i + 1 : n do
rij = hqi , vj i
vj = vj rij qi
end
end
Algorithm 6: Stable Gram-Schmidt orthogonalization

Householder factorization

Definition 21. Given a unit vector v 2 Rn , the matrix

vv T
H=I 2
vT v

is called the Householder matrix.

Proposition 2. The householder matrix H is orthogonal and reflects every

vector x 2 Rn in the hyper-plane span(v)? .

Note that a Householder matrix looks similar to the orthogonal projector

associated with v. From the geometry, it is straightforward to figure out. That
is why H is often called the Householder reflector.

Corollary 11. A Householder matrix H satisfying Hx = ±||x||e1 is associated

with v = x ⌥ ||x||e1
32 CHAPTER 2. SOLVING LINEAR SYSTEM

Among two possibilities of H, we choose v = x + sgn(x1 )||x||e1 for better

numerical stability, that is, it will give us further Hx from x. A Householder
matrix is basically an orthogonal matrix. Hence we can find a Householder
matrix that rotates a vector so as that it is aligned to the first natural basis. This
is a basic idea of the Householder factorization to perform QR decomposition.

for i = 1 : n do
x = A(i : m, i)
vi = x + sgn(x1 )||x||e1
vi
vi =
||vi ||
R(i : m, i : n) = A(i : m, i : n) 2vi viT A(i : m, i : n)
end
Algorithm 7: Householder factorization

It is known that Householder factorization is more stable than Gram-Schmidt

orthogonalization for QR factorization. Note that when we solve Ax = b using
QR factorization, what we need is not Q itself but QT b, because Ax = QRx =
b , Rx = QT b. Finally the solution x can be computed by back-substitution
using R.

for i = 1 : n do
b(i : m) = b(i : m) 2vi viT b(i : m)
end
Algorithm 8: Calculation of QT b

Computing range and nullspace using QR factorization

Consider a full QR factorization of A 2 Rm⇥n , where m > n. If A has rank r

then Q and R take the following forms:
2 3
h i R1 R2
Q = Q1 Q2 , R = 4 5,
0 0
2.3. ITERATIVE APPROACH 33

where Q1 2 Rm⇥r , Q2 2 Rm⇥(n r)

, R1 is an r ⇥ r upper triangular matrix with
no zeros in its diagonal entries, thus invertible, and R2 is an r ⇥ (n r) general
matrix.

Corollary 12. The column vectors of Q1 spans R(A).

Corollary 13. The column vectors of

2 3
R1 1 R 2
V =4 5
I

form a basis of the null space of A, where I represents an (n r) ⇥ (n r)

identity matrix.

2.2.5 Computational complexity

We will count operations of multiplication and division only because these op-
erations demand much more computational resources than the other operations
such as addition, subtraction, and substitution.
Pn 1
• LU factorization: k=1 (n k)(n k + 2) = 13 O(n3 )
Pn
• back-substitution: k=1 (n k + 1) = 12 O(n2 )
Pn
• Cholesky factorization: k=1 k(n k + 1) = 16 O(n3 )
Pn
• Gram-Schmidt orthogonalization: k=1 2m(k 1) + n + 2 = O(mn2 )
Pn
• Householder factorization: k=1 1 + (m k + 1)(3 + 2(n k + 1) =
2 1 3
O(mn 3n )
Pn
• Calculating QT b : k=1 (m k + 1)(m k) = O( 13 n3 + m2 n mn2 )

2.3 Iterative Approach

2.3.1 Motivation

A second-order partial di↵erential equation, i.e., one of the form

Auxx + 2Buxy + Cuyy + Dux + Euy + F = 0

34 CHAPTER 2. SOLVING LINEAR SYSTEM
⇥A B⇤
is called elliptic if the matrix B C is positive definite. Elliptic partial equations
have numerous applications in physics. The basic example of an elliptic PDE is
Laplace’s equation
r2 u = 0

in n-dimensional Euclidean space, where the Laplacian r2 is defined by

X @2
r2 = .
@x2i

Other examples of elliptic equations include the non-homogeneous Poisson’s

equation
r2 u = f (x).

In computational fluid dynamics, the solution of an incompressible flow problem

satisfies the Poisson’s equation, r2 u = f (µ, v), where u and v represent a
pressure field and a velocity field, respectively. We handle this problem by
applying the finite di↵erence of the Poisson’s equation. Consider a 2 dimensional
discrete Poisson’s equation on a 5 ⇥ 5 grid as an example;

1
(r2 u)ij = ui+1,j + ui 1,j + ui,j+1 + ui,j 1 4uij = fij ,
h2

where ui0 , ui4 , u0j , and u4j are boundary conditions.

It will result in a linear system Ax = b, where
2 3 2 3 2 3
4 1 0 1 0 0 0 0 0 u11 h2 f11 + u01 + u10
6 7 6 7 6 7
6 7 6 7 6 7
6 1 4 1 0 1 0 0 0 07 6u21 7 6 h2 f21 + u20 7
6 7 6 7 6 7
6 7 6 7 6 2 7
60 1 4 0 0 1 0 0 0 7 6 u 7 6 h f31 + u41 + u30 7
6 7 6 31 7 6 7
6 7 6 7 6 7
6 1 0 0 4 1 0 1 0 07 6u12 7 6 h2 f12 + u02 7
6 7 6 7 6 7
6 7 6 7 6 2 7
A=6 0 1 0 1 4 1 0 1 0 7 , x = 6 u 7 , b = 6 h f22 7.
6 7 6 22 7 6 7
6 7 6 7 6 7
60 0 1 0 1 4 0 0 17 6u32 7 6 h2 f32 + u42 7
6 7 6 7 6 7
6 7 6 7 6 2 7
60 0 0 1 0 0 4 1 0 7 6 u 7 6 h f13 + u03 + u14 7
6 7 6 13 7 6 7
6 7 6 7 6 7
60 0 0 0 1 0 1 4 17 6u23 7 6 h2 f23 + u24 7
4 5 4 5 4 5
0 0 0 0 0 1 0 1 4 u33 h2 f33 + u43 + u34

Note that A is sparse and strictly diagonally dominant.

2.3. ITERATIVE APPROACH 35

In this approach, we will have a giant matrix A, as our problem domain

increases. As a result, we often need efficient, yet stable methods to solve a big
linear system and direct methods discussed before may not be suitable for these
classes of problems.

2.3.2 Basics of iterative approach

Direct methods for solving linear systems theoretically give the exact solution in
a finite number of steps. Unfortunately, this is rarely true in many applications
because of rounding errors: an error made in one step spreads further in all
following steps. Contrary to direct methods, iterative methods construct a
series of solution approximations such that it converges to the exact solution of
a system. In this section, we seek iterative approaches to solve Ax = b.

Definition 22. A fixed point of a function is a point that is mapped to itself

by the function, that is, x⇤ satisfying f (x⇤ ) = x⇤ .

Definition 23. The spectral radius, ⇢(A), of a matrix A is defined by

⇢(A) = max | i |,

where i is the ith eigenvalue of A.

Lemma 7. lim An converges to zero, i↵ ⇢(A) < 1.

n!1

A general technique for linear iterative methods is to split A as

A=P N,

where P is a pre-conditioning matrix such that det(P ) 6= 0. We will denote

1
B=P N as an iteration matrix.

Proposition 3. The sequence of vectors { x(k+1) = Bx(k) + P 1

b } converges
to the solution of Ax = b for any choice of x(0) , i↵ ⇢(B) < 1, where B = P 1
N.
36 CHAPTER 2. SOLVING LINEAR SYSTEM

Proof. Let us define an error at the kth step e(k) = x(k) x⇤ , where x⇤ is the
solution of Ax = b. Then

e(k+1) = x(k+1) x⇤

= Bx(k) + P 1
b x⇤

= B(x(k) x⇤ + x⇤ ) + P 1
b x⇤

= Be(k) + (B I)x⇤ + P 1
b

= Be(k) + P 1
(N P )x⇤ + b

= Be(k) .

Hence e(k) converges to zero, in another word, x(k) converges to the solution x⇤
as k ! 1, i↵ ⇢(B) < 1.

2.3.3 Jacobi iteration

Let L, D, and U be a strictly lower triangular, a diagonal and an strictly upper

triangular sub-matrices of A, that is A = L + D + U . If all the diagonal entries
of A are non-zero, the Jacobi iteration chooses D as a pre-conditioning matrix,
that is, P = D, N = D A. Then we have

x(k+1) = D 1
(D A)x(k) + D 1
b
1
=D b (A D)x(k)

or
1 ⇣ ⌘
n
X
(k+1) (k)
xi = bi aij xj .
aii j=1
j6=i

Proposition 4. The Jacobi iteration converges if A is strictly diagonally dom-

inant:
X
|aii | > |aij |.
j6=i

2.3.4 Gauss-Seidel iteration

In the Gauss-Seidel iteration, the splitting of A is

P = D + L, N= U.
2.3. ITERATIVE APPROACH 37

It results in

x(k+1) = (D + L) 1
b U x(k)
1 1 1
= (I + D L) D b U x(k) .

It is equivalent to

1
(I + D L)x(k+1) = D 1
b U x(k)

or
x(k+1) = D 1
b Lx(k+1) U x(k) .

Finally we have

1 ⇣ ⌘
i 1
X n
X
(k+1) (k+1) (k)
xi = bi aij xj aij xj , i = 1, 2, · · · , n.
aii j=1 j=i+1

Unlike the Jacobi iteration, Gauss-Seidel iteration uses previously computed

(k+1)
results as soon as they are available. That is, the computation of xi requires
(k+1) (k) (k+1)
xj , where 1  j < i, as well as xi . Thus all the elements of x cannot
be computed simultaneously as in the Jacobi iteration.

Proposition 5.
⇢(BG ) = ⇢(BJ )2 .

2.3.5 Successive over-relaxation

Corollary 14. The sequence of vectors { x(k+1) = x(k) + P 1 (k)

r } is identical
to { x(k+1) = P 1
N x(k) + P 1
b }, where the residual vector r(k) at the kth step
is defined as r(k) = b Ax(k) .

Proof.

x(k+1) = x(k) + P 1 (k)

r

= x(k) + P 1
b Ax(k)
1
= (I P A)x(k) + P 1
b
1
= I P (P N ) x(k) + P 1
b
1
=P N x(k) + P 1
b.
38 CHAPTER 2. SOLVING LINEAR SYSTEM

1 (k)
You can consider P r as a relaxation term to correct the approxima-
tion. So it could be possible to introduce a relaxation factor w to control the
convergence rate. A generalized relaxation method has the form

x(k+1)
w = x(k) + wP 1 (k)
r .

If 0 < w < 1 the method is called under-relaxation, while if w > 1 it is called

over-relaxation.

(k+1) (k+1)
Proposition 6. xw = (1 w)x(k) + wx(k+1) , where xw and x(k+1) de-
note update of x at the (k + 1)th step with relaxation and without relaxation,
respectively.

Proof.

x(k+1)
w = x(k) + wP 1 (k)
r

= x(k) + wP 1
(b Ax(k) )
1
= (I wP A)x(k) + wP 1
b

= (1 w)x(k) + w(I P 1
A)x(k) + wP 1
b

= (1 w)x(k) + w P 1
N x(k) + P 1
b

= (1 w)x(k) + wx(k+1) .

Successive over-relaxation iteration or so called SOR iteration is a general-

ization of the Gauss-Seidel iteration.

w⇣ ⌘
i 1
X n
X
(k+1) (k) (k+1) (k)
xi = (1 w)xi + ai aij xj aij xj .
aii j=1 j=i+1

Corollary 15. SOR iteration converges i↵ 0 < w < 2.

Proof. SOR iteration corresponds to

x(k+1) = (1 w)x(k) + wD 1
b Lx(k+1) U x(k) .
2.3. ITERATIVE APPROACH 39

Arranging it,

1
(I + wD L)x(k+1) = (1 w)I wD 1
U x(k) + wD 1
b.

1 1 1
Hence B = (I + wD L) (1 w)I wD U and we check its convergence
condition.

1 1
B = (I + wD (1 w)I wD 1 U
L)
⇣ ⌘
1
= D 1 (D + wL) D 1 (1 w)D wU
⇣ ⌘
= (D + wL) 1 D D 1 (1 w)D wU
1
= (D + wL) (1 w)D wU .

det (1 w)D wU (1 w)n det D

Thus det B = = = (1 w)n , where n de-
det(D + wL) det D
notes the dimension of x. Here we use the property of det(D+L) = det(D+U ) =
det(D), where D, L and U are a diagonal, a strictly lower triangular, and a
strictly upper triangular matrices, respectively. Since
Y
i = det B = (1 w)n  max | i |n = ⇢(P 1
N ) < 1,
i

the condition for SOR iteration to converge is |1 w| < 1 or equivalently

0 < w < 2.

Proposition 7. The optimal choice of w is

2
wopt = p .
1+ 1 ⇢(BJ )2

Note that in order to compute wopt , we should know the spectral radius of the
Jacobi iteration, which requires an impractical amount of computation. However
inexpensive rough estimates of wopt (for example, from the power method) can
improve the convergence rate.

Corollary 16.
!2
⇢(B )
⇢(Bopt ) = p J .
1 + 1 ⇢(BJ )2
40 CHAPTER 2. SOLVING LINEAR SYSTEM

Problem Set
1. Suppose we are given a set of vectors {x1 , · · · , xp }, xi 2 Rn , i = 1, · · · , p.
Find the vector x̄ 2 Rn such that the average squared distance between x̄
and x1 , · · · , xp ,
p
1X
||x̄ xi ||2 ,
p i=1
is minimized.

2. Consider a least square problem of min ||Ax b||2 , where a rank of A 2

x
Rm⇥n is full. Find the least square solution x̃ using QR decomposition
and SVD.

3. Find the minimum norm least squares solution of Ax = b using QR fac-

torization when A 2 Rm⇥n has not full rank.
h i
R1 R2
[ Hint: Assume rank(A) = r < n. Then A = [ Q1 Q2 ] 0 0 , where
R1 2 R r⇥r
]

4. Given ↵ 2 R, find x̃ such that x̃ = arg minx ||Ax b||2P + ↵2 ||x x0 ||2Q ,
where ||x||P stands for the weighted norm, that is, ||x||P = xT P x.

[ Note: It is called the generalized Tikhonov regularization.]

5. Given a matrix 2 3
1 1 1 0
6 7
6 7
6 2 1 2 17
A=6
6
7,
7
6 2 0 1 37
4 5
2 1 1 3

(a) Find L1 and A1 satisfying A = L1 A1 , where L1 is a lower triangular

matrix resulting the first row of A1 is same to A and the first column
of A1 2is zero except3the first2element, that 3
is,
1 0 0 0 · · · ·
6 7 6 7
6 7 6 7
6⇥ 1 0 07 60 ⇥ ⇥ ⇥7
L1 = 66
7 , A1 = 6
7 6
7.
7
6⇥ 0 1 07 60 ⇥ ⇥ ⇥7
4 5 4 5
⇥ 0 0 1 0 ⇥ ⇥ ⇥
2.3. ITERATIVE APPROACH 41

(b) Find L2 and A2 satisfying A1 = L2 A2 , where

2 3 2 3
1 0 0 0 · · · ·
6 7 6 7
6 7 6 7
60 1 0 07 60 · · ·7
L2 = 6
6
7 , A2 = 6
7 6
7.
7
60 ⇥ 1 07 60 0 ⇥ ⇥7
4 5 4 5
0 ⇥ 0 1 0 0 ⇥ ⇥
(c) Find L3 and A3 satisfying A2 = L3 A3 , where
2 3 2 3
1 0 0 0 · · · ·
6 7 6 7
6 7 6 7
60 1 0 07 60 · · ·7
L3 = 6
6
7 , A3 = 6
7 6
7.
7
60 0 1 07 60 0 · ·7
4 5 4 5
0 0 ⇥ 1 0 0 0 ⇥
(d) Find a LU decomposition of A.

6. For a matrix 2 3
25 5 10
6 7
6 7
A=6 5 17 107 ,
4 5
10 10 62

(a) Verify that A is positive definite.

(b) Find the Cholesky factorization of A.

7. Consider a set of vectors, {(1, 1, 1), (1, 0, 1), (1, 1, 2)} in R3 .

(a) Show that given vectors span R3 .

(b) Apply Gram-Schmidt orthogonalization to get orthogonal basis.

(c) Apply Householder transform to get orthogonal basis.

8. A Householder matrix is given as

vv T
H=I 2 .
vT v
1
(a) Show that H = H T = H.

(b) When evaluating Hx, compare the operation counts of

⇣ vv T ⌘
I 2 x
vT v
42 CHAPTER 2. SOLVING LINEAR SYSTEM

and
vT x
x 2 v,
vT v
where x 2 Rn .

9. Consider the simple example of Richardson’s iteration,

x(k+1) = x(k) + ↵(b Ax(k) ),

where ↵ is a non-negative scalar.

(a) Express the iteration in a form x(k+1) = Bx(k) + b0 .

(b) Assume that i, the ith eigenvalue of A are all real and such that,
min  i  max . Show that µi , the ith eigenvalues of B satisfies

1 ↵ max  µi  1 ↵ min .

(c) Show that if min < 0 and max > 0, then ⇢(B) > 1 for any ↵, which
means that the method will always diverge.

(d) Assume that all eigenvalues of A are positive, i.e., min > 0. Show
that the following conditions must be satisfied for the method to
converge:

1 ↵ min <1

1 ↵ max > 1,

or equivalently
2
0<↵ .
max

(e) Show that the spectral radius of B is

⇢(B) = max{|1 ↵ min |, |1 ↵ max |}.

(f) Show that ⇢(B) is minimized at the intersection ⇢ = |1 ↵ min | and

⇢ = |1 ↵ max |.
2.3. ITERATIVE APPROACH 43

(g) Show that the optimal value ↵opt for the parameter ↵, i.e., the value
of ↵ which minimizes ⇢(B) is

2
↵opt = .
max + min

(h) Show that the spectral radius corresponding to ↵opt is

max min
⇢opt (G) = .
max + min

[ Note: It depicts the difficulty with the presence of small and large
eigenvalues. The convergence rate can be extremely small for realistic
problems. In addition, to achieve good convergence, eigenvalue esti-
mates are require to obtain the optimal or near-optimal ↵, and this
may be difficult. These observations are common to many iterative
methods that depends on an acceleration parameter. ]

10. Consider an iterative solver of (I + A)x = b

x(k+1) = b Ax(k) .

What is the condition of A for the above interative solver to converge?

44 CHAPTER 2. SOLVING LINEAR SYSTEM
Chapter 3

Representing Rotation

3.1 Rotation Group

We describe the orientation of the body by giving the relative orientation be-
tween a coordinate frame attached to the body and a global coordinate frame
which is uniquely fixed. We assume that all coordinate frames are right-handed.
Let A be the global frame, B the body frame, and r1 , r2 , r3 2 R3 the coordinates
h p p p i
of the principal axes of B relative to A. Here the 3 ⇥ 3 matrix R = r1 r2 r3 is
p p p
called a rotation matrix.
Generalizing a rotation matrix, we define the special orthogonal group in n
dimension as
SO(n) = {R 2 Rn⇥n |RRT = I, det R = 1}.

We are interested in SO(3) specifically because it represents an orientation in 3

dimensional world where we live.

3.2 Properties of Rotation Matrices

We define the vector space of n ⇥ n skew-symmetric matrices, so(n) as

so(n) = {S 2 Rn⇥n |S T + S = 0}

45
46 CHAPTER 3. REPRESENTING ROTATION

and a skew-symmetric operator [ · ] : R3 ! so(3) as

2 3
0 wz wy
6 7
6 7
[ w ] = 6 wz 0 wx 7 ,
4 5
wy wx 0

where w = [ wx , wy , wz ] 2 R3 . Note that skew-symmetric operator corresponds

to the cross product, that is, [ w ] v = w ⇥ v, where w, v 2 R3 .

Lemma 8. Given R 2 SO(3) and v, w 2 R3 , the following property holds:

R[ w ]RT = [ Rw ].

Proof. Without
 loss of generality, we prove that RT [ w ]R = [ RT w ]. Assume
p p p
that R = r1 r2 r3 .
p p p

" #  " #
r1T p p p r1T p p p
T
R [ w ]R = r2T [w] r 1 r 2 r3 = r2T w⇥r1 w⇥r2 w⇥r3
p p p p p p
r3T r3T
2 3
hr1 , w ⇥ r1 i hr1 , w ⇥ r2 i hr1 , w ⇥ r3 i
6 7
6 7
= 6hr2 , w ⇥ r1 i hr2 , w ⇥ r2 i hr2 , w ⇥ r3 i7
4 5
hr3 , w ⇥ r1 i hr3 , w ⇥ r2 i hr3 , w ⇥ r3 i
2 3
hw, r1 ⇥ r1 i hw, r2 ⇥ r1 i hw, r3 ⇥ r1 i
6 7
6 7
= 6hw, r1 ⇥ r2 i hw, r2 ⇥ r2 i hw, r3 ⇥ r2 i7
4 5
hw, r1 ⇥ r3 i hw, r2 ⇥ r3 i hw, r3 ⇥ r3 i
2 3
0 hw, r3 i hw, r2 i
6 7
6 7 T
= 6 hw, r3 i 0 hw, r1 i7 = [ R w ].
4 5
hw, r2 i hw, r1 i 0

Proposition 8. A rotation R 2 SO(3) is a rigid body transformation; that is,

1. R preserves distance: kRq Rpk = kq pk for all q, p 2 R3 .

2. R preserves orientation: R(v ⇥ w) = (Rv) ⇥ (Rw) for all v, w 2 R3 .

3.3. EXPONENTIAL MAP 47

Proof.

T
kRq Rpk2 = R(q p) R(q p) = (q p)T RT R(q p) = kq pk2 .

R(v ⇥ w) = R[ v ]w = R[ v ]RT Rw = [ Rv ](Rw) = (Rv) ⇥ (Rw).

3.3 Exponential Map

One of the most mathematically elegant parameterization of rotation matrices
will be the exponential map. We will see that the exponential of so(3) results
in SO(3). Recall that the exponential of a general square matrix is exp(A) =
1 2 1 3
I + A + 2! A + 3! A + · · · . In our case, the matrix A has a very special structure
of skew-symmetry. It could provide us a simpler form of the exponential map.

Lemma 9. Given [ w ] 2 so(3), the following relations hold:

2
[ w ] = wwT kwk2 I
3
[w] = kwk2 [ w ]

and higher powers of [ w ] can be calculated recursively.

Theorem 13. Rodrigues’ formula

2
exp(✓[ w ]) = I + sin ✓[ w ] + (1 cos ✓)[ w ] ,

where ✓ 2 R, kwk = 1.

Proof.

1 2 2 1 3 1 4
exp(✓[ w ]) = I + ✓[ w ] + ✓ [ w ] + ✓3 [ w ] + ✓4 [ w ] + · · ·
2! 3! 4!
1 3 1 1 4 2
= I + (✓ ✓ + · · · )[ w ] + ( ✓2 ✓ + · · · )[ w ]
3! 2! 4!
2
= I + sin ✓[ w ] + (1 cos ✓)[ w ]
48 CHAPTER 3. REPRESENTING ROTATION

Proposition 9. Given a skew-symmetric matrix [ w ] 2 so(3),

exp([ w ]) 2 SO(3).

1
Proof. Let R = exp([ w ]). It is sufficient to show that R = RT and det R = 1.

1 T
exp([ w ]) = exp( [ w ]) = exp([ w ] ) = exp([ w ])T .

It is obvious that det exp([ 0 ]) = det I = 1. Using the continuity of the deter-
minant as a function of the entries of a matrix, we conclude that det R = 1.

Since we have the exponential mapping which maps so(3) to SO(3), we can
use a skew-symmetric matrix to parameterize a rotation matrix, if we have the
inverse of the mapping.

Proposition 10. Given R 2 SO(3), there exists w 2 R3 , kwk = 1 and ✓ 2 R

such that R = exp(✓[ w ]).

Proof. Defining v✓ = 1 cos ✓, c✓ = cos ✓, s✓ = sin ✓, expand the Rodrigues’

formula.

2
exp(✓[ w ]) = I + sin ✓[ w ] + (1 cos ✓)[ w ]
2 3
wx2 v✓ + c✓ w x w y v✓ wz s✓ wx wz v✓ + wy s✓
6 7
6 7
= 6wx wy v✓ + w z s✓ wy2 v✓ + c✓ wy wz v✓ wx s✓ 7 .
4 5
wx w z v✓ wy s✓ w y w z v✓ + w x s✓ wz2 v✓ + c✓

Since Tr(R) = wx2 v✓ + c✓ + wy2 v✓ + c✓ + wz2 v✓ + c✓ = 1 + 2 cos ✓, we have

✓ ◆
1 Tr(R) 1
✓ = cos .
2
Note that there is an ambiguity in the value of ✓, in the sense that ✓ ± 2⇡n or
✓ ± 2⇡n could be chosen as well. Now, equating the o↵-diagonal terms of R
and exp(✓[ w ]), we get

r32 r23 = 2wx s✓

r13 r31 = 2wy s✓

r21 r12 = 2wz s✓ .

3.4. EULER ANGLES 49

If ✓ 6= 0, we choose 2 3
r32 r23
1 6 7
6 7
w= 6r13 r31 7 .
2 sin ✓ 4 5
r21 r12

Note that if 2⇡ ✓ had been chosen earlier, the axis of rotation would have
been w. Indeed, the exponential map is a surjective(many-to-one) map from
R3 onto SO(3). If R 6= I, the above construction shows that there are two
distinct w and ✓ 2 [ 0, 2⇡) such that R = exp([ w ]).
Figure 3.1 illustrates rotating a vector p by an angle ✓ about an axis w,
where kwk = 1. Then the rotated vector q is

q = hw, piw + cos ✓ p hw, piw + sin ✓(w ⇥ p)

= (1 cos ✓)hw, piw + cos ✓p + sin ✓[ w ]p

= (1 cos ✓)wwT p + cos ✓p + sin ✓[ w ]p

2
= (1 cos ✓) I + [ w ] p + cos ✓p + sin ✓[ w ]p
2
= p + sin ✓[ w ]p + (1 cos ✓)[ w ] p
2
= I + sin ✓[ w ] + (1 cos ✓)[ w ] p

= exp(✓[ w ])p .

Now we can figure out the geometric meaning of exp(✓[ w ]). That is why the
parameterization using the exponential map has aliases of an angle-axis repre-
sentation or an equivalent axis representation.

3.4 Euler Angles

One method of describing the orientation of a coordinate frame {B} relative to
another coordinate frame {A} is as follows: start with a frame {B} coincident
with a frame {A}. First, rotate the {B} frame about z-axis of frame {B} by an
angle ↵, then rotate about the (new) y-axis of frame {B} by an angle , and
50 CHAPTER 3. REPRESENTING ROTATION

Figure 3.1: Rotation of p.

then rotate about the (one again, new) x-axis of frame {B} by an angle . This
yields a net orientation R(↵, , ) and the the triple of angles (↵, , ) is called
the ZYX Euler angles.
The rotation matrix R(↵, , ) is given as

R(↵, , ) = Rz (↵)Ry ( )Rx ( )

2 32 32 3
cos ↵ sin ↵ 0 cos 0 sin 1 0 0
6 76 76 7
6 76 76 7
= 6 sin ↵ cos ↵ 07 6 0 1 0 7 60 cos sin 7
4 54 54 5
0 0 1 sin 0 cos 0 sin cos
2 3
cos ↵ cos cos ↵ sin sin sin ↵ cos cos ↵ sin cos + sin ↵ sin
6 7
6 7
= 6 sin ↵ cos sin ↵ sin sin + cos ↵ cos sin ↵ sin cos cos ↵ sin 7
4 5
sin cos sin cos cos
h r11r12 r13
i
Given R = r21 r22 r23
r31 r32 r33
, then the ZYX Euler angles are

↵ = atan2(r21 , r11 )
q
= atan2( r31 , 2 + r2 )
r11 21

= atan2(r32 , r33 ) .
h 0 cos ↵ sin sin ↵ cos cos ↵ cos +sin ↵ sin i
Note that when = ⇡/2, R = 0 sin ↵ sin +cos ↵ cos sin ↵ cos cos ↵ sin =
1 0 0
h 0 sin(↵ ) cos(↵ ) i
0 cos(↵ ) sin(↵ ) , which means that there exist infinitely many choices
1 0 0
of ↵ and for the given R. Similar things happen when = ⇡/2. This
singularity problem is often called a gimbal lock. Any kinds of Euler angles
3.5. QUATERNIONS 51

has the singularity problems. So we have to be careful about it. It has been
proven mathematically that there exists no single parameterization of the rota-
tion group which is free from singularity.

3.5 Quaternions

3.5.1 Structure of Quaternions

Quaternions, denoted by H, generalize complex numbers and can be used to

represent rotations in much the same way as complex numbers on the unit circle
can be used to represent planar rotations. Unlike the other parameterizations,
quaternions give a global parameterization of SO(3), at the cost of using four
numbers instead of three to represent a rotation.
Formally, a quaternion is a vector quantity of the form

q = s + v1 i + v2 j + v3 k, s, vi 2 R, i = 1, · · · , 3,

where s is the scalar component of q and v = (v1 , v2 , v3 ) is the vector component.

A shorthand notation is q = (s, v) with s 2 R, v 2 R3 . Multiplication is
distributive and associative, but not commutative; it satisfies the relations

ai = ia aj = ja ak = ka a2R

i·i=j·j=k·k=i·j·k= 1

i·j= j·i=k j·k= k·j=i k·i= i · k = j.

Proposition 11. Given two quaternions q = (s, v) and p = (t, w), the product
between two quaternions has a simple form in terms of the inner and cross
products between vectors in R3 as follows:

qp = (st hv, wi, sw + tv + v ⇥ w).

The conjugate of a quaternion q = (s, v) is given by q ⇤ = (s, v) and the

magnitude of a quaternion is defined as follows;

kqk2 = q · q ⇤ = s2 + kvk2 .
52 CHAPTER 3. REPRESENTING ROTATION

Lemma 10. For q, p 2 H, kqpk = kqkkpk.

Proof.
kqpk2 = (qp)(qp)⇤ = (qp)(p⇤ q ⇤ ) = qkpk2 q ⇤ = kqk2 kpk2 .

Since we have a definition of a multiplication operator for quaternions, we

can find the identity element for quaternion multiplication, 1 2 H, which is a
1
quaternion satisfying q1 = 1q = q, 8q. Also we can find q , the inverse of q
1 1
which satisfies qq =q q = 1.

Lemma 11.

1 = (1, 0)
1 q⇤
q = .
kqk2

3.5.2 Unit Quaternions

Among general quaternions, we are interested especially in the unit quaternions

H1 which is the subset of all q 2 H such that kqk = 1. Note that for unit
quaternion, q 1
= q⇤ .

Proposition 12. Let q = (s, v) 2 H1 . Then there exists v 0 2 R3 and ✓ 2 ( ⇡, ⇡]

such that q = (cos ✓, sin ✓v 0 ).

Proof. kqk2 = s2 +kvk2 = s2 +k 2 = 1, where k = kvk. Hence we find ✓ 2 ( ⇡, ⇡]

satisfying s = cos ✓, k = sin ✓. If we define v 0 = v/kvk, then q = (cos ✓, sin ✓v 0 ),
where kv 0 k = 1.

Proposition 13. Let q 2 H1 , p = (s, v) 2 H. Then qpq ⇤ = p0 , where p0 = (s, v 0 )

with kvk = kv 0 k.

Proof. Below we write S(q) = 12 (q +q ⇤ ) for the real part of q. The proof consists
of three steps. We first show S(p0 ) = S(p) for p 2 {(s, 0)|s 2 R} and then for
p 2 {(0, v)|v 2 R3 }. Finally these results are combined to prove the proposition
for p 2 H.
3.5. QUATERNIONS 53

1. Let p = (s, 0). Then p0 = q(s, 0)q ⇤ = sqq ⇤ = s. Hence S(p0 ) = S(p) = s.

2. Let p = (0, v). Then 2S(p0 ) = p0 + p0⇤ = qpq ⇤ + (qpq ⇤ )⇤ = q(p + p⇤ )q ⇤ =

2S(p) = 0.

3. p0 = qpq ⇤ = q (s, 0) + (0, v) q ⇤ = q(s, 0)q ⇤ + q(0, v)q ⇤ = (s, 0) + (0, v 0 ) =

(s, v 0 ). Finally kp0 k = kqpq ⇤ k = kpk results that kv 0 k = kvk.

Theorem 14. Let q = (cos ✓, sin ✓w) 2 H1 , p = (0, r) 2 H. Then p0 = qpq ⇤

corresponds to rotating r by an angle 2✓ about the axis w.

Proof. Let q = (s, v) = (cos ✓, sin ✓w).

qpq ⇤ = (s, v)(0, r)(s, v)

= 0, (s2 kvk2 )r + 2hv, riv + 2s(v ⇥ r)

= 0, (cos2 ✓ sin2 ✓)r + 2 sin2 ✓hw, riw + 2 cos ✓ sin ✓(w ⇥ r)

= 0, cos 2✓r + (1 cos 2✓)hw, riw + sin 2✓(w ⇥ r)

= 0, exp(2✓[ w ])r .

Given a rotation matrix R = exp(✓[ w ]), its associated unit quaternion is

Q = cos(✓/2), sin(✓/2) w ,

where w 2 R3 represents the axis of rotation and ✓ 2 R represents the angle of

rotation. A detailed calculation shows that if qab represents a rotation between
frame {A} and frame {B}, and qbc represents a rotation between frame {B} and
frame {C}, then the rotation between {A} and {C} is given by the quaternion

qac = qab qbc .

Thus, the group operation on unit quaternions directly corresponds to the group
operation for rotations. Given a unit quaternion q = (s, v), we can extract the
54 CHAPTER 3. REPRESENTING ROTATION

corresponding rotation by setting

8
>
< v
1 sin(✓/2) if ✓ 6= 0,
✓ = 2 cos s , w= ,
>
:0 otherwise,

and R = exp(✓[ w ]).

Since the group structure for quaternions directly corresponds to that of
rotations, quaternions provide an efficient representation for rotations which do
not su↵er from singularities.

3.5.3 Interpolation of Quaternioins

First we fine an exponential and log functions of quaternions.

Definition 24. Let q = (cos ✓, sin ✓w) 2 H1 . Then log(q) is defined by

log(q) = (0, ✓w).

Note that log(q) is not in general a unit quaternion.

Definition 25. For a quaternion of the form q = (0, ✓w), ✓ 2 R, w 2 R3 , kwk =

1, the exp(q) is defined by

exp(q) = (cos ✓, sin ✓w).

Note that the exponential and logarithmic functions are mutually inverse,
and the exponential function maps into H1 .
From the above definitions, we can define exponentiation for q 2 H1 , t 2 R:

Definition 26. Let q 2 H1 , t 2 R. Exponentiation q t is defined by

q t = exp t log(q) .

Proposition 14. Let q = (cos ✓, sin ✓w) 2 H1 , a, b 2 R. Then

q a q b = q a+b .
3.5. QUATERNIONS 55

Proof.

q a q b = exp(a log q) exp(b log q)

= exp a(0, ✓w) exp b(0, ✓w)

= (cos a✓, sin a✓w)(cos b✓, sin b✓w)

= cos (a + b)✓ , sin (a + b)✓ w

⇣ ⌘
= exp 0, (a + b)✓w

= exp (a + b) log(q)

= q a+b .

Proposition 15. Let q, p = (sin ✓, cos ✓w) 2 H1 , t 2 R. Then qpt q ⇤ = (qpq ⇤ )t .

Proof.

qpt q ⇤ = q(cos t✓, sin t✓w)q ⇤

= (cos t✓, sin t✓w0 ) (by Proposition 6)

= exp t(0, ✓w0 )

⇣ ⌘
= exp t log (cos ✓, sin ✓w0 )

= exp t log(qpq ⇤ )

= (qpq ⇤ )t .

Lemma 12. If q 2 H1 , then q t 2 H1 where t 2 R.

Proof. Let q = (cos ✓, sin ✓w). Then kq t k = k exp t log(q) k = k exp t(0, ✓w) k =
k exp (0, t✓w) k = k(cos(t✓), sin(t✓)w)k = 1.

Definition 27. For q, p 2 H, t 2 [0, 1], the linear quaternion interpolation is

given as
Lerp(q, p, t) = (1 t)q + tp.
56 CHAPTER 3. REPRESENTING ROTATION

Figure 3.2: Linear interpolation and spherical interpolation on H1 .

/ H1 for q, p 2 H1 . It is because H1 is not a vector

Note that Lerp(q, p, t) 2
space. Therefore you may need to normalize r = Lerp(q, p, t) by r/krk.

Definition 28. For q, p 2 H1 , t 2 [0, 1], the spherical quaternion interpolation

is given as

Slerp(q, p, t) = q(q ⇤ p)t .

Lemma 13. For q, p 2 H1 , t 2 [0, 1], Slerp(q, p, t) 2 H1 .

Proof. It is obvious that q ⇤ p 2 H1 since kq ⇤ pk = kq ⇤ kkpk = 1, and (q ⇤ p)t

is therefore in H1 from the previous lemma. As a result, kSlerp(q, p, t)k =
kq(q ⇤ p)t k = kqkk(q ⇤ p)t k = 1.

Unlike the linear interpolation, the spherical quaternion interpolation maps

into H1 . Moreover the interpolation curve corresponds to a geodesic on a unit
quaternion, which is similar to the straight line on a vector space or the great arc
on a sphere. In Figure 3.2, green arrows show Lerp(q, p, 0.5) and Slerp(q, p, 0.5).

Lemma 14. For q, p 2 H1 , t 2 [0, 1],

Slerp(q, p, t) = Slerp(p, q, 1 t).

3.5. QUATERNIONS 57

Proof.

Slerp(q, p, t) = q(q ⇤ p)t

= q(p⇤ q) t

= pp⇤ q(p⇤ q) t

= p(p⇤ q)1 t

= Slerp(p, q, 1 t).

Lemma 15. Slerp(q, p, t) corresponds to Rq exp t log(RqT Rp ) , where Rp and

Rq represents rotation matrices associated with q and p respectively.

Problem Set
n Rn + Rn T o
1. Show that a sequence Rn+1 = converges to a rotation ma-
2
trix, where R0 2 R3⇥3 is an arbitrary matrix. Implement the algorithm
and check the result numerically. You can use any computer languages(C,
C++, JAVA, matlab, ...)

2. Show that R = (I [ a ]) 1
(I + [ a ]) 2 SO(3), where a 2 R3 .

3. Let us rotate a vector p by an angle ✓ about an axis w.

(a) Represent p as a projection onto w and its complement.

(b) Rotate the complement vector about w by ✓.

(c) Show that the rotated vector q is represented as follows:

2
q = I + sin ✓[ w ] + (1 cos ✓)[ w ] p.

4. Implement a spherical linear interpolation Slerp(Ra , Rb , t) = Ra exp(t log(RaT Rb ))

for SO(3).

(a) Generate two random rotation matrices Ra and Rb .

58 CHAPTER 3. REPRESENTING ROTATION

(b) draw a 3D curve of log(Slerp(Ra , Rb , t)) for the generated random

rotation matrices, where t 2 [0, 1].

(c) compare the result with the following interpolatoin method

Llerp(Ra , Rb , t) = exp((1 t) log Ra + t log Rb )

5. Derive the maps between a rotation matrix and a corresponding quater-

nion.

6. Show that unit quaternions q and q represent the same rotation.

Chapter 4

Introduction to
Optimization

4.1 Di↵erential Calculus

Definition 29. Let f : Rn ! Rm be a di↵erentiable function.
2 3
f1 (x)
6 7
6 . 7
f (x) = 6 .. 7 ,
4 5
fm (x)
where x 2 Rn . Then the Jacobian matrix of f is defined as
2 3
@f1 @f1
· · ·
6 @x1 @xn 7
6 . .. .. 7
J = 6 .. . . 7 .
4 5
@fm
@x1 · · · @f
@xn
m

Definition 30. If f : Rn ! R is di↵erentiable, then the function rf defined

by 2 3
@f
(x)
6 @x1 7
6 .. 7
rf (x) = 6 . 7 = Df (x)T
4 5
@f
@xn (x)
is called the gradient of f .

59
60 CHAPTER 4. INTRODUCTION TO OPTIMIZATION

Definition 31. Given f : Rn ! R, if Df is di↵erentiable, we write the deriva-

tive of Df as
2 3
@2f @2f @2f
2 @x2 @x1 ··· @xn @x1
6 @x2 1 7
6 @ f @2f @2f 7
6 @x1 @x2 @x22
··· @xn @x2 7
D2 f = 6 6 7
.. .. .. .. 7.
6 . . . . 7
4 5
@2f @2f @2f
@x1 @xn @x2 @xn ··· @x2n

The matrix D2 f (x) is called the Hessian matrix of f at x.

Definition 32. The level set of a function f : Rn ! R at level c is the set of

points

S = {x|f (x) = c}.

Theorem 15. The vector rf (x0 ) is orthogonal to the tangent vector to an

arbitrary smooth curve passing through x0 on the level set determined by f (x) =
f (x0 ).

4.2 Unconstrained Optimization

We consider the optimization problem

minimize f (x)

subject to x 2 ⌦.

Definition 33. Suppose f : Rn ! R is a real-valued function defined on a set

⌦ ⇢ Rn . A point x⇤ 2 ⌦ is called a local minimizer of f over ⌦ if there exists
✏ > 0 such that f (x) f (x⇤ ) for all x 2 ⌦ \ {x⇤ } and ||x x⇤ || < ✏.

Definition 34. A point x⇤ 2 ⌦ is called a global minimizer of f over ⌦ if

f (x) f (x⇤ ) for all x 2 ⌦ \ {x⇤ }.

In general, global minimizers are difficult to find. Therefore, in practice, we

often have to be satisfied with finding local minimizers.
4.2. UNCONSTRAINED OPTIMIZATION 61

4.2.1 Conditions for local minima

Definition 35. A directional derivative of f in the direction d, denoted by

@f /@d, is the real-valued function defined by

@f f (x + ↵d) f (x)
= lim .
@d ↵!1 ↵

If ||d|| = 1, then @f /@d is the rate of increase of f at x in the direction d.

Lemma 16.
@f
= dT rf (x)
@d

Proof.

@f d
= f (x + ↵d)
@d d↵ ↵=0
T
= rf (x) d

= dT rf (x)

Theorem 16. First-Order Necessary Condition(FONC) Let ⌦ be a subset

of Rn and f 2 C 1 a real-valued function on ⌦. If x⇤ is a local minimizer of f
over ⌦, then
rf (x⇤ ) = 0.

Proof. In the neighborhood of x⇤ , f (x⇤ )  f (x⇤ + d) for all d 2 Rn , 2 R.

Then choosing > 0, we have

f (x⇤ + d) f (x⇤ )
0 .

By taking limits as ! 0,

0  dT rf (x⇤ ), 8d 2 Rn .

Since d is arbitrary, we can replace it with d, which results in 0 dT rf (x⇤ ).

Therefore rf (x⇤ ) = 0.
62 CHAPTER 4. INTRODUCTION TO OPTIMIZATION

Theorem 17. Second-Order Necessary Condition(SONC) Let ⌦ ⇢ Rn , f 2

C 2 a function of ⌦, x⇤ a local minimizer of f over ⌦, and d a feasible direction
at x⇤ . If dT rf (x⇤ ) = 0, then

dT H(x⇤ )d 0,

where H is the Hessian of f .

From FONC, all the local minimizer satisfy rf = 0, which is non-linear

systems of equations in general. Hence it is impractical to solve rf = 0 directly
to get the local minimizers. Instead we consider iterative algorithms which
produce points

x(k+1) = x(k) + ↵k g (k) , with f (x(k+1) ) < f (x(k) ),

where g (k) 2 Rn is called the step direction or search direction, ↵k is called the
step size. Given a starting point x 2 Rn , general descent methods repeat the
following procedures until stopping criterion is satisfied.

1. Determine a search direction g (k) .

2. Line search: choose a step size ↵k > 0.

3. Update: x(k+1) = x(k) + ↵k g (k) .

Then our questions will be how to choose the search direction and the step size.

4.2.2 Gradient method

A natural choice for the search direction is the negative gradient

g= rf (x).

The resulting algorithm is called the gradient algorithm. The gradient method
often exhibits linear convergence and the convergence rate depends greatly on
the condition number of the Hessian. Consequently, it is very simple, but rarely
used in practice due to its slow convergence.
4.2. UNCONSTRAINED OPTIMIZATION 63

4.2.3 Steepest descent method

Consider a quadratic function

1 T
f (x) = x Qx bT x,
2
where Q 2 Rn⇥n is a symmetric positive definite matrix, b, x 2 Rn . Then the
search direction is g (k) = rf (x(k) ) = b Qx(k) and the search step is

↵k = arg min f x(k) + ↵g (k)

↵ 0
✓ ◆
1 (k) T
= arg min x + ↵g (k) Q x(k) + ↵g (k) bT x(k) + ↵g (k)
↵ 0 2
✓ ◆
1 (k) T (k) 2 T T
= arg min g Qg ↵ + g (k) Qx(k) g (k) b ↵ + · · · .
↵ 0 2
Hence T
g (k) g (k)
↵k = T
.
g (k) Qg (k)
Consider a case of f (x1 , x2 ) = x21 + x22 . Then, starting from an arbitrary ini-
tial point x(0) we arrive at the solution x⇤ = 0 in only one step. However if
x21
f (x1 , x2 ) = 5 + x22 , then the method of steepest descent shu✏es ine↵ectively
back and forth when searching for the minimizer in a narrow valley. this example
illustrates a major drawback in the steepest descent method.

4.2.4 Newton’s method

Consider a quadratic approximation of f 2 C 2 around x(k) :

1
f (x) ⇡ f (x(k) ) + (x x(k) )T g (k) + (x x(k) )T H(x(k) )(x x(k) ) = q(x),
2
where, for simplicity, we denote g (k) = rf (x(k) ). Applying the FONC yields

0 = rq(x) = g (k) + H(x(k) )(x x(k) ).

If H(x(k) ) is positive definite, then q(x) achieves a minimum at

x(k+1) = x(k) H(x(k) ) 1 (k)

g .

This recursive formula represents Newton’s method. A drawback of Newton’s

method is that evaluation of H(x(k) ) for large n can be computationally expen-
sive. Furthermore, we have to solve the set of n linear equations H(x(k) ) 1 (k)
g .
64 CHAPTER 4. INTRODUCTION TO OPTIMIZATION

4.2.5 Conjugate direction methods

Definition 36. Let Q be a real symmetric n ⇥ n matrix. The directions

d(1) , d(2) , · · · , d(m) are Q-conjugate if, for all i 6= j,

T
d(i) Qd(j) = 0.

Lemma 17. Let Q be a symmetric positive definite n ⇥ n matrix. If the direc-

tions d(1) , d(2) , · · · , d(k) 2 Rn , k  n are non-zero and Q-conjugate, then they
are linearly independent.

Proof. Let ↵1 , ↵2 , · · · , ↵k be scalars such that

↵1 d(1) + ↵2 d(2) + · · · + ↵k d(k) = 0.

T
Premultiplying the above equality by d(j) Q, 1  j  k, yields

T
↵j d(j) Qd(j) = 0

T
because all other terms d(j) Qd(i) = 0, i 6= j, by Q-conjugacy. Here Q is positive
definite, so ↵j = 0, j = 1, · · · , k. Therefore d(1) , d(2) , · · · , d(k) , k  n are linearly
independent.

Now we present the conjugate direction algorithm for minimizing the quadratic
function of n variables
1 T
f (x) = x Qx xT b,
2
where Q is positive definite, x 2 Rn . Note that since Q > 0, the function f has
a global minimizer that can be found by solving Qx = b.
[Basic Conjugate Direction Algorithm] Given a starting point x(1) , and
Q-conjugate directions d(1) , d(2) , · · · , d(n) ; for k 0,

g (k) = rf (x(k) ) = b Qx(k) ,

T
g (k) d(k)
↵k = T
,
d(k) Qd(k)
x(k+1) = x(k) + ↵k d(k) .
4.2. UNCONSTRAINED OPTIMIZATION 65

Theorem 18. For any starting point x(1) , the basic conjugate direction al-
gorithm converges to the unique x⇤ which solves Qx = b in n steps; that is,
x(n+1) = x⇤ .

Proof. Consider x⇤ x(0) 2 Rn . Because the d(i) are linearly independent, there
exist constant i, i = 0, · · · , n, such that

x⇤ x(1) = 1d
(1)
+ ··· nd
(n)
.

T
Now pre-multiply both sides of the above equation by d(k) Q, 1  k  n, to
obtain
T (k) T
d(k) Q(x⇤ x(1) ) = kd Qd(k) ,
T
where the terms d(k) Qd(i) = 0, k 6= i, by the Q-conjugate property. Hence
T
d(k) Q(x⇤ x(1) )
k = T
.
d(k) Qd(k)

Now, we can write

x(k) = x(1) + ↵1 d(1) + · · · + ↵k d(k) .

Therefore,
x(k) x(1) = ↵1 d(1) + · · · + ↵k d(k) .

So writing
x⇤ x(1) = x⇤ x(k) + x(k) x(1)
T
and pre-multipying the above by d(k) Q, we obtain

T T T
d(k) Q(x⇤ x(1) ) = d(k) Q(x⇤ x(k) ) = d(k) g (k)

since g (k) = Qx(k) b and Qx⇤ = b. Thus

T
d(k) g (k)
k = T
= ↵k
d(k) Qd(k)

and x⇤ = x(n+1) , which completes the proof.

66 CHAPTER 4. INTRODUCTION TO OPTIMIZATION

Note that basic conjugate direction algorithm requires priori information of

Q-conjugate directions, which is not preferred. Hence we consider the conjugate
gradient algorithm which does not use pre-specified conjugate directions, but
instead calculates the directions as the algorithm progresses. At each stage of
the algorithm, the direction is calculated as a linear combination of the previous
direction and the current gradient, in such a way that all the directions re
mutually Q-conjugate. It is why it is called the conjugate gradient algorithm.

Lemma 18. In the conjugate direction algorithm,

T
g (k+1) d(i) = 0

for all k, 1  k  n, and 1  i  k.

Proof. Note that

Q(x(k+1) x(k) ) = Q(x(k+1) b (Q(x(k) b) = g (k+1) g (k)

since g (k) = Qx(k) b. Thus

g (k+1) = g (k) + ↵k Qd(k) .

We prove the lemma by induction. The result is true for k = 0 because

(1) T (i)
g d = 0, as shown before. We now show that if the result is true for
T T
k 1[that is g (k) d(i) = 0.i  k 1] then it is true for k[that is g (k+1) d(i) =
T
0.i  k]. By the induction hypothesis, g (k) d(i) = 0. Since

g (k+1) = g (k) + ↵k Qd(k) ,

T
then by Q-conjugacy, d(k) Qd(i) = 0, and hence

T T T
g (k+1) d(i) = g (k) d(i) + ↵k d(k) Qd(k) = 0.

It remains to be shown that

T
g (k+1) d(k) = 0.
4.2. UNCONSTRAINED OPTIMIZATION 67

Indeed,

T
g (k+1) d(k) = (Qx(k+1) b)T d(k)
✓ T ◆ !T
(k) g (k) d(k) (k)
= x T
d Qd(k) bT d(k)
d(k) Qd(k)
T
= (x(k) Q b)T d(k) g (k) d(k)

=0

Since Qx(k) b = g (k) . Therefore, by induction, for all 0  k  n 1 and

0  i  k,
T
g (k+1) d(i) = 0.

In the conjugate gradient algorithm, we choose d(k+1) as a linear combination

of g (k+1) and d(k) , that is,

d(k+1) = g (k+1) kd
(k)
, k = 1, 2, · · · .

The coefficient k are chosen in such a way that d(k+1) is Q-conjugate to

d(1) , · · · , d(k) , that is,
T
g (k+1) Qd(k)
k = T
.
d(k) Qd(k)
Summarizing the conjugate gradient algorithm,

1. Set k = 1; select the initial point x(1) .

2. g (1) = rf (x(1) ). If g (1) = 0, stop, else set d(1) = g (1) .

T
g (k) d(k)
3. ↵k = T
.
d(k) Qd(k)

4. x(k+1) = x(k) + ↵k d(k) .

5. g (k+1) = rf (x(k+1) ). If g (k+1) = 0, stop.

T
g (k+1) Qd(k)
6. k = T
.
d(k) Qd(k)
68 CHAPTER 4. INTRODUCTION TO OPTIMIZATION

7. d(k+1) = g (k+1) kd
(k)
.

8. Set k := k + 1; goto 3.

Proposition 16. In the conjugate gradient algorithm, the direction d(1) , · · · , d(n)
are Q-conjugate.
T
Proof. We use induction. We first show d(1) Qd(2) = 0.
T T
d(1) Qd(2) = d(1) Q g (2) + 1 d(1)
✓ T ◆
T g (2) Qd(1) (1)
= d(1) Q g (2) T
d
d(1) Qd(1)
T
(1) T (2) g (2) Qd(1) T
=d Qg T
d(1) Qd(1)
d(1) Qd(1)
= 0.

Now assume that d(1) , · · · , d(k) , k < n are Q-conjugate directions. From Lemmm 18,
T
we have g (k+1) d(i) = 0, j = 1, · · · , k. Thus g (k+1) is orthogonal to each of the
directions d(1) , · · · , d(k) . We now show that
T
g (k+1) g (j) = 0, j = 1, · · · , k.

Fix j 2 1, · · · , k. We have d(j) = g (j) j 1d

(j 1)
.
Substituting this equation into the previous one yields
T T (k+1) T (j 1)
g (k+1) d(j) = 0 = g (k+1) g (j) j 1g d .
T T
Because g (k+1) d(j 1)
= 0, it follows that g (k+1) g (j) = 0.
T
We are now ready to show that d(k+1) Qd(j) = 0, j = 1, · · · , k. We have
T (k) T
d(k+1) Qd(j) = g (k+1) kd Qd(j) .
T
If j < k, then d(k+1) Qd(j) = 0 by virtue of the inductive hypothesis. Hence,
we have
T T
d(k+1) Qd(j) = g (k+1) Qd(j) .
T
But g (j+1) = g (j) ↵j Qd(j) . Since g (k+1) g (i) = 0, i = 1, · · · , k, then
T T (g (j+1) g (j) )
d(k+1) Qd(j) = g (k+1) = 0.
↵j
4.2. UNCONSTRAINED OPTIMIZATION 69

T
It remains to be shown that d(k+1) Qd(k) = 0. We have

T
d(k+1) Qd(k) = ( g (k+1) + kd
(k) T
) Qd(k) .

T
Using the expression for k, we get d(k+1) Qd(k) = 0, which completes the
proof.