Description
The goal of this homework is to illustrate applications of the singular value decomposition (SVD) and
Tikhonov regularization to image deblurring.
Mathematical Background:
Consider a linear system of equations
Ax = d (1)
where d ∈ R
n
is a given vector (observed/received data), x ∈ R
n
is an unknown vector (transmitted data),
and A ∈ R
n×n
is a given nonsingular matrix (transformation of data).
In many applications the observed/received data is corrupted by unknown noise or measurement errors,
such that we only have access to
dˆ = d + ξ
where ξ ∈ R
n
is an unknown noise vector. The problem is then formulated as follows:
Given the matrix A and a vector dˆ, provide an approximation of the unknown vector x.
Practical difficulties arise when A is ill-conditioned. In such case, an attempt to provide an approximation
xˆ ≈ x by solving
Axˆ = dˆ (2)
will not work since xˆ will be drastically corrupted by noise. Regularization techniques are required to
address this issue.
Truncated SVD. Consider the singular value decomposition of the matrix A
A = USVT ∈ R
n×n
(3)
where U = [u1 u2 . . . un], S = diag(σ1, σ2, . . . σn), V = [v1 v2 . . . vn]. The solution to the system (2) is
expressed as
xˆ = [USVT
]
−1dˆ = VS−1UTdˆ =
Xn
i=1
u
T
i dˆ
σi
vi (4)
For very small σi
, the terms (u
T
i dˆ)/σi will result in noise amplification.
In the truncated SVD regularization, an approximation xˆk to the solution is obtained by truncating the
sum in (4) at a selected index k,
xbk =
X
k
i=1
u
T
i dˆ
σi
vi (5)
The main difficulty is to find an appropriate truncation index k such that xbk is a good approximation to
the true value x.
Tikhonov Regularization. In this approach, an approximation to x is obtained as the solution to the
least squares minimization problem:
min
xb∈Rn
J(xb), where J(xb)
def
= kAxb − dbk
2 + λ
2
kxbk
2
(6)
In (6), k · k denotes the Euclidean vector norm and λ is a scalar regularization parameter that controls
the smoothness of the solution. If λ = 0 then no regularization is applied and the solution is (4). If λ
is large, then λ
2kxbk
2 has a significant contribution to the cost (6) and the solution xb can not be a good
approximation of x. We view the solution to (6) as a function of the parameter λ and denote it xbλ. The
main difficulty is to find an appropriate value of the regularization parameter λ such that xbλ is a good
approximation to the true value x.
The solution xbλ to the minimization problem (6) is obtained by solving the linear system
ATA + λ
2
I
xbλ = ATdb (7)
and may be expressed using the SVD (3) of A as
xbλ =
Xn
i=1
fi(λ)
u
T
i db
σi
vi (8)
where the filter factors fi(λ) are defined as
fi(λ) = σ
2
i
σ
2
i + λ
2
, i = 1, 2, . . . n (9)
The main difficulty in this approach is to select an appropriate value of the regularization parameter λ.
Homework content
Consider an image represented as a matrix X ∈ R
n×m and a blurring process represented by the matrix
A ∈ R
n×n
. Whereas the true image is the solution to the matrix equation
AX = D
we want to reconstruct an approximation Xb to the image X given the blurring matrix A ∈ R
n×n and a
matrix of received noisy data
Db = AX + ξ
The regularization procedure (5) or (8-9) is used to find an approximate solution Xb , one column at a time:
Algorithm
[U, S, V ] = svd(A)
for i = 1 : m
ˆd = Dˆ(:, i)
evaluate ˆxi
from (5) – TSVD solution or from (8-9) – Tikhonov solution
Xˆ = [ˆx1 xˆ2 . . . xˆm
For this assignment n = m = 128, and X is a 128 × 128 matrix representing the image of a coin. The
blurring operation is represented as follows: Consider a 128 × 128 symmetric tridiagonal matrix B with
entries 1
B(i, i) = 1 − 2L, i = 1, 2, . . . n
B(i, i + 1) = L, i = 1, 2, . . . n − 1; B(i + 1, i) = L, i = 1, 2, . . . n − 1
where L = 0.45. Then the blurring operator is
A = B
10
The file ”hw2data.m” represents the 128 × 128 noisy data matrix Dˆ , obtained as
Dˆ = AX + ξ
where ξ is an unknown 128 × 128 noise matrix. To visualize Dˆ , in MATLAB you may execute:
load hw2data.m; D = hw2data; imagesc(D); colormap(gray)
Jour job: Implement the TSVD and the Tikhonov regularization procedures to reconstruct an approximation Xˆ of the unknown image X such that the year and the inscriptions on the coin can be read.
• (50 points) Implement a function Xˆ = regularize(A, D, method, p ˆ ) that takes as input the noisy
data matrix Dˆ, the blurring matrix A and the regularization parameter p, and returns Xˆ, the
regularization approximation to X. The input value method is used to distinguish between the
TSVD solution (5) in which case p is the truncation index and the Tikhonov solution (8-9) in which
case p is the regularization parameter λ. Alternatively, you may write a function for each method
tsvd and tikhonov, respectively.
• (20 points) Find values of k and λ such that from Xˆ the year and the inscriptions on the coin can
be read. Plot the l − curve.
Things to hand in: listing of the tsvd and tikhonov functions, value of index k (tsvd) and parameter
λ (Tikhonov), reconstructed image: year, inscriptions; a plot of the singular values σi on a log10 scale and
the filter factors fi
; a plot of the l − curve.
1
such matrix results from discretization of the 1-D heat equation ut − kuxx = 0 with L = k∆t/(∆x)
2
.
