MATH4602 Scientific Computing Tutorial 4

Extension to Assignment1

Inverse Matrix and Solution to Linear Equations

When we are calculating the inverse of a matrix, denoted as

An×n , we do the following steps:

A · A−1 = I

Denote A−1 := X, an unknown n-by-n variable. Then

calculating A−1 is equivalent to finding solutions to linear
systems A · X = [e1 , · · · , en ], which can be written as
A · [x1 , · · · , xn ] = [Ax1 , · · · , Axn ] = [e1 , · · · , en ]. As A is
invertible, A is of full column rank, the following linear system
must have n unique solutions.
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Inverse Matrix and Solution to Linear Equations

Ax1 = e1
Ax2 = e2
···
Axn = en

The solutions constitute the inverse matrix.

The relationship between inverse matrix and matrix
factorization.
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Inverse Matrix and Solution to Linear Equations

Human calculate inverse matrix by observation and elementary

row operations. But computers can’t observe independently
and they must follow step-by-step instructions. That’s why they
solve linear equations instead of using elementary row
operations. If A is symmetric, A = L · LT , and
Axk = L · LT xk = L · yk = ek . They firstly solve yk through
substitution, then find LT xk = yk through substitution. If not
symmetric, use LU factorization. (You may also refer to lecture
notes.)
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

The Determinant of Tridiagnal Matrix and

Mathematical Induction

Find the determinant.

α γ 0
β α γ
..
an = |An | = β α .
.. ..
. . γ
0 β α n×n
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

The Determinant of Tridiagnal Matrix and

Mathematical Induction

We solve the above-mentioned determinant by Mathematical

Induction. Clearly, a1 = α, a2 = α2 − βγ.
By Laplace formula

β γ 0
α γ
..
an = α · an−1 − γ β α .
.. ..
. . γ
0 β α (n−1)×(n−1)
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

The Determinant of Tridiagnal Matrix and

Mathematical Induction

α γ 0
..
β α .
= α · an−1 − γβ .. ..
. . γ
0 β α (n−2)×(n−2)
= α · an−1 − γβan−2

Since we already know a1 and a2 , we can solve {an } iteratively.

MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi) Definite Matrice

7. [Positive definite matrix] Determine the values of a such that

the following matrix A is symmetric positive definite where
 
a 1 1 1
 1 a 1 1 
A=  1 1 a 1 .


1 1 1 a

When we encounter a symmetric matrix and need to verify

whether positive (semi)definite, we usually have 3 methods.
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi)Definite Matrice

Solution1:
(Brute force) if n is not large, calculate the eigenvalues and
eigenvectors separately. And let all eigenvalues greater than
(greater or equal to) 0.

Solution2:
Quadratic form = a · (sum of suqares)+ b · (sum)2 , (a, b are
generally the function of unknown parameters we need to deal
with)
Then analyze. (Please refer to Assignment 1 Solution)
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi) Definite Matrice

Solution3:
By observation, since the sum of rows are equal,
v = [1, 1, 1, 1]T is obviously an eigenvector and the
corresponding eigenvalue is λ1 = (a + 3).
 
a 1 1 1
 1 a 1 1 
A·v =  1 1 a 1  · v = (a + 3) · v


1 1 1 a
Similarly by observation, there is an eigenvalue λ2,3,4 = a − 1
and the corresponding eigenvectors are
     
1 1 1
 −1   0   0 
 0  , v3 =  −1  , v4 =  0  .
v2 =      

0 0 −1
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi) Definite Matrice

Hence, let λ1,2,3,4 > 0, we obtain a > 1.

Comments:

1. If rank(A)=r < n, then A has an eigenvalue λ = 0 with

corresponding n − r eigenvectors.
2. If the row sums of A equal, then A has an eigenvector
v =P [1, · · · , 1]T , with corresponding eigenvalue
λ = nj=1 aij
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi)Definite Matrice

Example1:
Given a matrix A with rank(A)=1. Since A can be written as
A = u · v T , A has an eigenvalue λ1,··· ,n−1 = 0, with
corresponding eigenvectors
     
v2 v3 vn
 −v1   0   0 
     
 0   −v1 
v1 =   , v2 =   , · · · , vn−1 =  0  .
 
 ..   ..   .. 
 .   .   . 
0 0 −v1

And an eigenvalue λ = v T · u with corresponding eigenvector u.

(You can easily verify this by
A · u = u · v T · u = u · (v T · u) = (v T · u) · u)
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Positive (Semi) Definite Matrice

Example 2:(Householder Matrix)
(a) Let A be a real n × n matrix and v ∈ Rn \ {0}. Show that

2
P=I− vvT = I − 2 (v/∥v∥) (v/∥v∥)T
vT v
is symmetric (i.e., P T = P) and orthogonal (i.e., P T = P −1 ).
Proof.
Direct verification.

(b) Let u be a unit vector and P = I − 2uuT . Find the geometric

meaning of Px for any x ∈ Rn \ {0}. Hence, or otherwise,
determine the eigenvaluse and the corresponding eigenectors
of P.
MATH4602 Scientific Computing Tutorial 4
Extension to Assignment1

Proof.
Let U be the span of u. Note that
h i h i
Px = x − (uT x)u − (uT x)u = x − (uT x)u + −(uT x)u

The first summand is the component of x in U ⊥ (the orthogonal

complement of U) while −(uT x)u is the negative of the
projection of x along the (unit) vector u. Thus, Px is the mirror
reflection of x along the hyperplane U ⊥ .
We then proceed to find the eigensystem of P. Clearly,
Pu = −u so that −1 is an eigenvalue of P and every nonzero
multiple of u is the corresponding eigenvectors. On the other
hand, the restriction of P to U ⊥ is the identity map whereby 1 is
an eigenvalue of P and every nonzero vector in U ⊥ is an
eigenvector corresponding to 1.
MATH4602 Scientific Computing Tutorial 4
Tensor

A review on tensor.

1. A ⊗ (B + C) = A ⊗ B + A ⊗ C
2. (B + C) ⊗ A = B ⊗ A + C ⊗ A
3. (kA) ⊗ B = A ⊗ (kB) = k (A ⊗ B)
4. (A ⊗ B) ⊗ C = A ⊗ (B ⊗ C)
5. A ⊗ 0 = 0 ⊗ A = 0
6. In general
A ⊗ B ̸= B ⊗ A.
7. If A, B, C, D are matrices of such size that one can form the
matrix products AC and BD, then we have

(A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD).

MATH4602 Scientific Computing Tutorial 4
Tensor

N N N
(A B)(C D) = AC BD
Proof.

O O
(A B)(C D)
   
a11 B a12 B · · · a1n B c11 D c12 D · · · c1n D
 .. .
.. ..  ·  ...
. .. .. 
= .
 
. . 
am1 B am2 B · · · amn B cm1 D cm2 D · · · cmn D
 Pn Pn 
k =1 a1k ck 1 BD · · · k =1 a1k ckn BD
=
 .. .. 
Pn . Pn . 
k =1 amk ck 1 BD · · · k =1 amk ckn BD
O
= (AC) (BD)(by definition)
MATH4602 Scientific Computing Tutorial 4
Tensor

Let A be an n × n matrix with eigenvalues λ1 , . . . , λn and B be a

m × m matrix with eigenvalues µ1 , . . . , µm . From Assignment 1,
we know that the eigenvalues of A ⊗ B are λi µj
(i = 1, . . . , n; j = 1, . . . , m). We have
n X
X m Xn Xm
trace(A⊗B) = λi µj = ( λi )( µj ) = trace(A)·trace(B).
i=1 j=1 i=1 j=1

det(A ⊗ B) = λm m n n m n
1 · · · λn · µ1 · · · µm = (det(A)) (det(B)) .
MATH4602 Scientific Computing Tutorial 4
Tensor

Proof.
Perform SVD on A and B. Then A = U1T Σ1 V1 and
B = U2T Σ2 V2 . Obviously,
A⊗B = (U1 Σ1 V1T )⊗(U2 Σ2 V2T ) = (U1 ⊗U2 )(Σ1 ⊗Σ2 )(V1T ⊗V2T ).
Thus rank (A ⊗ B) = rank (Σ1 ⊗ Σ2 ) = rank (A) · rank (B).

Second proof
Proof.
Suppose A has rank r1 and B has rank r2 . Then there exists
invertible matrices P1 , P2 , Q1 and Q2 such that

P1 AQ1 = Ir1 and P2 BQ2 = Ir2

Thus rank (A ⊗ B) = rank ((P1 ⊗ P2 )(A ⊗ B)(Q1 ⊗ Q2 )) =

rank ((P1 AQ1 ) ⊗ (P2 BQ2 )) = rank (Ir1 ⊗ Ir2 ) = r1 r2 =
rank (A) · rank (B).
MATH4602 Scientific Computing Tutorial 4
Exam 2002

(a) Let A be an non-singular n × n matrix, u and v be two n × k

(k ≤ n) matrices such that the matrix (Ik + v t A−1 u) is
non-singular. Prove that

−1  −1
A + uv t = A−1 − A−1 u Ik + v t A−1 u v t A−1 .

(b) Use the result in 1(a) or otherwise to obtain the solution of

the following n × n linear system:
  
1
 
n+1 1 ··· 1 1 x1
 1 n+2 1 ··· 1  0 
  x2  
  


 .. . . . . .  .
. 
Bx ≡  . .. .. .. ..   ..  = 
  
. 
  
1 2n − 1 1   xn−1   ... 
 
 1 ···

1 ··· ··· 1 2n xn 0

where Bii = n + i and Bij = 1 when i ̸= j.

MATH4602 Scientific Computing Tutorial 4
Exam 2002

(a)
  −1 
−1 −1 t −1 t −1
A + uv t

A −A u I+v A u v A

= In + A−1 uv t −A−1 u(I + v t A−1 u)−1 v t − A−1 u(I + v t A−1 u)−1 v t

= In + A−1 (uv t ) − A−1 u(I + v t A−1 u)−1 (I + v t A−1 u)v t
= In + A−1 (uv t ) − A−1 uIk v t
= In .
MATH4602 Scientific Computing Tutorial 4
Exam 2002

(b) We note that

B = Diag(n, n + 1, . . . , 2n − 1) + u t · u

where u = (1, 1, . . . , 1). Let

1 1 1
z = Diag(n, n + 1, . . . , 2n − 1)−1 u = ( , ,..., )t
n n+1 2n − 1
and
1
y = Diag(n, n + 1, . . . , 2n − 1)−1 (1, 0, . . . , 0)t = ( , 0, . . . , 0)t .
n
MATH4602 Scientific Computing Tutorial 4
Exam 2002

Since
n
X 1
(1 + (1, 1, . . . , 1) · z) = 1 + ̸= 0
n+i −1
i=1

we can apply the result in 1(a). Let

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

then we have
1 − α −α −α t
x = B −1 (1, 0, . . . , 0)t = y − αz = ( , ,..., ).
n n+1 2n − 1
MATH4602 Scientific Computing Tutorial 4
Exam 2003

(a) Consider the following n × n (n ≥ 3) symmetric matrix:

 
2 1 0 ··· 0 1
.. .
. .. 0 
 
 1 2 1
. 
 
 0 1 . . . . . . 0 .. 

A= .
 .. . . . 

.. ..
 . . . . 1 .. 
 
 0 ··· 0 1 2 0 
1 0 · · · · · · 0 2 n×n

Prove that A is positive definite by showing that

xAxt > 0 if and only if x ̸= 0 where x = (x1 , x2 , . . . , xn ) ∈ R n .

MATH4602 Scientific Computing Tutorial 4
Exam 2003

(b) Suppose that the Cholesky factorization of A takes the following

form:
a1 0 · · · ··· ···
  
0 0 a1 l1 0 0 d1
 .. ..  .. .. 
 l1 a2 . . 0  .
  0 a2 l2 . d2
 
.. ..   .. . . ..
 
 . . . . . . . . 
 0
 l2 . . . .  .
 . . . 0 . .

 . . . .  . .
 .. .. .. a ..   .. .. a

 n−2 0  n−2 ln−2 dn−2 

 0 ··· 0 ln−2 an−1 0   0 ··· ··· 0 an−1 dn−1 
d1 d2 · · · dn−2 dn−1 an 0 0 ··· ··· 0 an
show that
 2
 a1 = 2;
2
+ li2 = 2

ai+1 for i = 1, 2, . . . , n − 2;




 ai li = 1 for i = 1, 2, . . . , n − 2;



a1 d1 = 1;


 ai di = −li−1 di−1 for i = 2, . . . , n − 1;


 n−1
X
di2 + an2 = 2.




MATH4602 Scientific Computing Tutorial 4
Exam 2003

(c) Show that

2 1
ai+1 =2− for i = 1, 2, . . . , n − 2.
ai2

Hence prove that

i +1
ai2 = for i = 1, 2, . . . , n − 1
i
and
i
li2 = for i = 1, 2, . . . , n − 2.
i +1
MATH4602 Scientific Computing Tutorial 4
Exam 2003

(d) Prove that

1
di2 = for i = 1, 2, . . . , n − 1
i(i + 1)

and
1
. an2 = 1 +
n
(e) What is the operational cost of the Cholesky factorization in
(b)?
MATH4602 Scientific Computing Tutorial 4
Exam 2003

(a) We note that

n−2
X
t
xAx = (xi + xi+1 )2 + (x1 + xn )2 + xn−1
2
+ xn2 ≥ 0.
i=1

Moreover if xAxt = 0 then xn = xn−1 = . . . = x1 = 0.

MATH4602 Scientific Computing Tutorial 4
Exam 2003

(b) We note that LLt =

 2
a1 a1 l1 ··· 0 a1 d1
..
a 1 l1 l12 + a22 a2 l2 . l1 d1 + a2 d2


 .. .. ..

 0 a2 l2 . . .
 .. .. .. ..

 . . . an−2 ln−2 .
2 2
 0 ··· an−2 ln−2 ln−2 + an−1 ln−2 dn−1 + an−1 dn−1
Pn−1 2 2
a1 d1 l1 d1 + a2 d2 ··· ln−2 dn−2 + an−1 dn−1 i=1 di + an

By comparing coefficients with A, we have



 a12 = 2;
2 2

 ai+1 + li = 2


 for i = 1, 2, . . . , n − 2;
 a i li = 1


 for i = 1, 2, . . . , n − 2;
a1 d1 = 1;


 ai di = −li−1 di−1 for i = 2, . . . , n − 1;
n−1


 X
di2 + an2 = 2.




MATH4602 Scientific Computing Tutorial 4
Exam 2003

2 + l 2 = 2 and a l = 1 for i = 1, 2, . . . , n − 2, we
(c) Since ai+1 i i i
have
2 1
ai+1 = 2 − 2 for i = 1, 2, . . . , n − 2.
ai
1+1
By using induction on i we have a12 = 2 = 1 . Assume that
ai2 = (i + 1)/i then

2 i i +2
ai+1 =2− = .
i +1 i +1
Hence by the principle of Mathematical Induction we have
i +1
ai2 = for i = 1, 2, . . . , n − 1.
i
Since li2 ai2 = 1, we have
i
li2 = for i = 1, 2, . . . , n − 1.
i +1
MATH4602 Scientific Computing Tutorial 4
Exam 2003

(d) We have a12 d12 = 1 and ai2 di2 = li−1

2 d 2 for i = 2, . . . , n − 1
i−1
and
1 l2 i −1 2
d12 = 2 , and di2 = i−1 2
di−1 = d .
ai ai2 i + 1 i−1
Therefore
i −1 i −2 i −3 i −4 2 1 2
di2 = × × × ×· · ·× × ×d12 = ×d 2 .
i +1 i i −1 i −2 4 3 i(i + 1) 1

Since d12 = 12 , the result follows.

Moreover,
n−1 n−1 n−1
X X 1 X 1 1 1
an2 = 2− di2 = 2− = 2− ( − ) = 1+ .
i(i + 1) i i +1 n
i=1 i=1 i=1
MATH4602 Scientific Computing Tutorial 4
Exam 2003

(e) The operational cost for the Cholesky factorization is O(n).

Because from 1(b)-1(d) the total cost for the obtaining ai2 , li2
and di2 is O(n). Hence the Cholesky factorization takes O(n) in
this case.
MATH4602 Scientific Computing Tutorial 4
Exam 2010

(a) Show that the inverse of a lower triangular Toeplitz matrix is

still a lower triangular Toeplitz matrix.
(b) Let T be a lower triangular Toeplitz matrix given as follows:
 
1 0 ··· 0 0
 p
 1 0 ··· 0  
 2 . . . . . . .
. 
T = p . . . . 
.
 .. . .. p 
 . 1 0 
pn−1 · · · p2 p 1

Find T −1 .
MATH4602 Scientific Computing Tutorial 4
Exam 2010

(a) For a lower triangular Toeplitz matrix

 
t1 0 · · · 0 0
 t2 t1
 0 ··· 0 

 . .. . .. ... .. 
T =  t3 . .

 .. . . 
 . . t2 t1 0 
tn · · · t3 t2 t1

We can obtain the i column of T −1 solving T xi = ei where ei is

the unit column vector with the ith entry being equal to 1. The
linear system can be solved by using forward substitution.

Suppose that x1 = (x1 , x2 , . . . , xn )T , then one can see that

x2 = (0, x1 , x2 , . . . , xn−1 )T , x3 = (0, 0, x1 , x2 , . . . , xn−2 )T · · · ,
xn = (0, 0, . . . , 0, x1 )T . Eventually, T −1 is a lower triangular
matrix.
MATH4602 Scientific Computing Tutorial 4
Exam 2010

(b) Using the method described in (a)(i) we can solve

 
1 0 ··· 0 0
 −p 1 0 ··· 0 
 
−1
 .. .. .. .. 
T = 0 . . . . .
 .. . . . −p 1 0 
 . 
0 · · · 0 −p 1
MATH4602 Scientific Computing Tutorial 4
Exam 2014

(a) Show that the following symmetric matrix B is positive

definite where  
4 1 1 1
 1 4 1 1 
B=  1 1 4 1 .


1 1 1 4
(b) Consider the following linear system of equations:
    
4 1 1 1 x1 0
 1 4 1 1   x2   0 
Bx =   1 1 4 1   x3  =  0  .
   

1 1 1 4 x4 21

Apply the Gaussian elimination method (show all your steps) to

solve the above linear system.
(c) Using the results in (b) or otherwise, find the Doolittle’s LU
factorization of B and hence the Cholesky’s factorization of B.
MATH4602 Scientific Computing Tutorial 4
Exam 2014

(a) We note that

4
X
T
[x1 , x2 , x3 , x4 ]B[x1 , x2 , x3 , x4 ] = 3 xi2 +(x1 +x2 +x3 +x4 )2 = 0
i=1

if and only if x1 = x2 = x3 = x4 = 0.
MATH4602 Scientific Computing Tutorial 4
Exam 2014

(b) With multiplier 1/4 we have

 
4 1 1 1
 0 15/4 3/4 3/4 
B1 =  0 3/4 15/4 3/4
.

0 3/4 3/4 15/4
Then with multiplier 1/5, we have
 
4 1 1 1
 0 15/4 3/4 3/4 
B2 =   0
.
0 18/5 3/5 
0 0 3/5 18/5
Finally, with multiplier 1/6, we have
 
4 1 1 1
 0 15/4 3/4 3/4 
B3 =   0
.
0 18/5 3/5 
0 0 0 7/2
MATH4602 Scientific Computing Tutorial 4
Exam 2014

Hence we can solve

    
4 1 1 1 x1 0
 0 15/4 3/4 3/4   x2   0 
B3 = 
 0   x3  =  0
   
0 18/5 3/5 
0 0 0 7/2 x3 21

by backward substitution, [x1 , x2 , x3 , x4 ] = [−1, −1, −1, 6].

MATH4602 Scientific Computing Tutorial 4
Exam 2014

(c) Using the multiplier obtained in the Gaussian elimination

process, we have B = LU where
   
1 0 0 0 4 1 1 1
 1/4 1 0 0   0 15/4 3/4 3/4 
L=  1/4 1/5 1 0  and U =  0
  .
0 18/5 3/5 
1/4 1/5 1/6 1 0 0 0 7/2

Now we have B =
   
1 0 0 0 4 0 0 0 1 1/4 1/4 1/4
15
 1/4 1
 0 0  0
 
4 0 0   0 1 1/5 1/5
  
.
 1/4 1/5 1 0  0 18
0 5 0  0 0 1 1/6 
1/4 1/5 1/6 1 7 0 0 0 1
0 0 0 2
MATH4602 Scientific Computing Tutorial 4
Exam 2014

√ √
Then B ≡ L · D · LT = (L · D)(L · D)T = P · P T where

2 √0 0 0
 
√  1/2
P =L· D = √ 15/2 p 0 0 .
 1/2 15/10 18/5 0 
√ p p
1/2 15/10 1/10 7/2
MATH4602 Scientific Computing Tutorial 4
Exam 2014

Let  
5 0 0 2 0 0

 0 4 0 0 1 0 

 0 0 5 0 0 2 
A= .

 1 0 0 4 0 0 

 0 3 0 0 3 0 
0 0 1 0 0 5
With the help of the function diag, define the matrix A by using
MATLAB. Write a short MATLAB script to permute the columns
of A to get  
5 2 0 0 0 0
1 4 0 0 0 0
 
0 0 4 1 0 0
à = 
 
0 0 3 3 0 0

0 0 0 0 5 2
0 0 0 0 1 5
MATH4602 Scientific Computing Tutorial 4
Exam 2014

A sample script is

d11=diag([5,4,5]);d12=diag([2,1,2]);
d21=diag([1,3,1]);d22=diag([4,3,5]);
A=[d11 d12;d21 d22];

A sample script is

n=size(A);n=n(1);
for i=1:(n/2) p(2*i-1)=i;
p(2*i)=n/2+i;
end
tildeA=A(p,p);