Description
The goal of this assignment is to illustrate issues related with inversion of a diffusion process in the
presence of noisy data and the need for regularization. You are encouraged to use any software available in
your favorite programming language. In particular, Matlab provides a friendly syntax for implementation.
Mathematical Background:
Consider a linear system of equations
Ax = d (1)
where x ∈ R
n
is an unknown vector (transmitted data), d ∈ R
n
is a given vector (observed/received 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,
xˆ = A−1dˆ = A−1
(d + ξ) = x + A−1
ξ
| {z }
amplif ied noise
(3)
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 + ξ (4)
The direct inversion approach provides
Xb = A−1Db (5)
The data file clown1.m contains a 200 × 320 matrix X (”the truth”) that represents the picture
clown1.jpg. Alternatively, in Matlab use ”load clown” since ”clown.mat” should be available in the Matlab
built in data files. The image is loaded by default in a matrix X. Then use ”A=X; imagesc(A)”, ”colormap(gray)”.
Blurring (data transmission process) is represented by a diffusion process as follows: Consider a 200×200
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.1. Then the blurring operator is B to the power k
A = B
k
such that higher values of k increase the amount of diffusion. We will use a range of values k = 1, 2, . . . 20.
The noise ξ in received data is represented as a matrix with entries random numbers taken from a
normal distribution with mean 0 and standard deviation 0.01,
ξ = 0.01 ∗ randn(200, 320)
• Task 1 (10 points) Provide the plot of the condition number of the matrix A = Bk
for k = 1 : 20.
• Task 2 (20 points) Generate the noisy receive data image Db as in (4) and the reconstructed image
Xb given by the direct inversion approach (5). Provide the plots of the absolute error and the relative
error in the reconstructed image Xˆ for k = 1 : 20,
kX − Xˆ k,
kX − Xˆ k
kXk
using the Frobenius (Euclidean) matrix norm.
• Task 3 (10 points) Provide the reconstructed image for k = 1, k = 5 and k = 20
• Task 4 (10 points) Repeat tasks 2 and 3 for a noise matrix with standard deviation 0.1,
ξ = 0.1 ∗ randn(200, 320)
