$30
Exercises 129
The matrix A (with coefficients aij ) is available from the class webpage (see below), and
was constructed as follows. We take
aij = r−2
ij max{cos θij , 0},
where rij denotes the distance between lamp j and the midpoint of patch i, and θij denotes
the angle between the upward normal of patch i and the vector from the midpoint of patch
i to lamp j, as shown in the figure. This model takes into account “self-shading” (i.e., the
fact that a patch is illuminated only by lamps in the halfspace it faces) but not shading of
one patch caused by another. Of course we could use a more complex illumination model,
including shading and even reflections. This just changes the matrix relating the lamp
powers to the patch illumination levels.
The problem is to determine lamp powers that make the illumination levels Ii close to a
given desired level Ides. In other words, we want to choose the n-vector x such that
!n
j=1
aijxj ≈ Ides, i = 1, . . . , m,
but we also have to observe the power limits 0 ≤ xj ≤ 1. This is an example of a
constrained optimization problem. The objective is to achieve an illumination level that is
as uniform as possible; the constraint is that the components of x must satisfy 0 ≤ xj ≤ 1.
Finding the exact solution of this minimization problem requires specialized numerical
techniques for constrained optimization. However, we can solve it approximately using
least-squares.
In this problem we consider two approximate methods that are based on least-squares,
and compare them for the data generated using [A,Ides] = ch8ex8, with the MATLAB
ch8ex8.m. The elements of A are the coefficients aij . In this example we have m = 11,
n = 7 so A is 11 × 7, and Ides = 2.
(a) Saturate the least-squares solution. The first method is simply to ignore the bounds
on the lamp powers. We solve the least-squares problem
minimize !m
i=1
(
!n
j=1
aijxj − Ides)
2
ignoring the constraints 0 ≤ xj ≤ 1. If we are lucky, the solution will satisfy the
bounds 0 ≤ xj ≤ 1, for j = 1,...,n. If not, we replace xj with zero if xj < 0 and
with one if xj > 1.
Apply this method to the problem data generated by ch8ex8.m, and calculate the
resulting value of the cost function "m
i=1(Ii − Ides)
2.
(b) Weighted least-squares. The second method is to solve the problem
minimize !m
i=1
(
!n
j=1
aijxj − Ides)
2 + µ
!n
j=1
(xj − 0.5)2
,
where the constant µ ≥ 0 is used to attach a cost to the deviation of the powers from
the value 0.5, which lies in the middle of the power limits. For µ = 0, this is the
same least-squares problem as in part (a). If we take µ large enough, the solution
of this problem will satisfy 0 ≤ xj ≤ 1.
Formulate this problem as a least-squares problem in the variables x, and solve it
for µ = 1, µ = 2, µ = 3, etc., until you find a value of µ such that all components
of the solution
" x satisfy 0 ≤ xj ≤ 1. For that solution x, calculate the cost function m
i=1(Ii − Ides)
2 and compare with the value you obtained in part (a).
HW1Prob2.m
[A, Ides] = HW1Prob2
130 8 Linear least-squares
8.9 De-noising using least-squares. The figure shows a signal of length 1000, corrupted with
noise. We are asked to estimate the original signal. This is called signal reconstruction,
or de-noising, or smoothing. In this problem we apply a smoothing method based on
least-squares.
0 200 400 600 800 1000
−0.5
0
0.5
i
xcorr
We will represent the corrupted signal as a vector xcor of size 1000. (The values can
be obtained as xcor = ch8ex9 using the file ch8ex9.m.) The estimated signal (i.e., the
variable in the problem) will be represented as a vector ˆx of size 1000.
The idea of the method is as follows. We assume that the noise in the signal is the small
and rapidly varying component. To reconstruct the signal, we decompose xcor in two
parts
xcor = ˆx + v
where v is small and rapidly varying, and ˆx is close to xcor (ˆx ≈ xcor) and slowly varying
(ˆxi+1 ≈ xˆi). We can achieve such a decomposition by choosing ˆx as the solution of the
least-squares problem
minimize ∥x − xcor∥2 + µ
"999
i=1(xi+1 − xi)
2, (8.10)
where µ is a positive constant. The first term ∥x − xcor∥2 measures how much x deviates
from xcor. The second term, "999
i=1(xi+1 − xi)
2, penalizes rapid changes of the signal
between two samples. By minimizing a weighted sum of both terms, we obtain an estimate
xˆ that is close to xcor (i.e., has a small value of ∥xˆ − xcor∥2) and varies slowly (i.e., has a
small value of "999
i=1(ˆxi+1 − xˆi)
2). The parameter µ is used to adjust the relative weight
of both terms.
Problem (8.10) is a least-squares problem, because it can be expressed as
minimize ∥Ax − b∥2
where
A =
# I
õ D $
, b =
# xcor
0
$
,
and D is a 999 × 1000-matrix defined as
D =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
−1 1 00 ··· 0 0 00
0 −1 10 ··· 0 0 00
0 0 −1 1 ··· 0 0 00
.
.
. .
.
. .
.
. .
.
. .
.
. .
.
. .
.
. .
.
.
0 0 00 ··· −1 1 00
0 0 00 ··· 0 −1 10
0 0 00 ··· 0 0 −1 1
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
.
3)
xcor = HW1Prob3 HW1Prob3.
Exercises 131
The matrix A is quite large (1999 × 1000), but also very sparse, so we will solve the
least-squares problem using the Cholesky factorization method. You should verify that
the normal equations are given by
(I + µDT D)x = xcor. (8.11)
MATLAB and Octave provide special routines for solving sparse linear equations, and
they are used as follows. There are two types of matrices: full (or dense) and sparse. If
you define a matrix, it is considered full by default, unless you specify that it is sparse.
You can convert a full matrix to sparse format using the command A = sparse(A), and
a sparse matrix to full format using the command A = full(A).
When you type x = A\b where A is n×n, MATLAB chooses different algorithms depending on the type of A. If A is full it uses the standard method for general matrices (LU
or Cholesky factorization, depending on whether A is symmetric positive definite or not).
If A is sparse, it uses an LU or Cholesky factorization algorithm that takes advantage of
sparsity. In our application, the matrix I + µDT D is sparse (in fact tridiagonal), so if we
make sure to define it as a sparse matrix, the normal equations will be solved much more
quickly than if we ignore the sparsity.
The command to create a sparse zero matrix of dimension m×n is A = sparse(m,n). The
command A = speye(n) creates a sparse n × n-identity matrix. If you add or multiply
sparse matrices, the result is automatically considered sparse.
This means you can solve the normal equations (8.11) by the following MATLAB code
(assuming µ and xcor are defined):
D = sparse(999,1000);
D(:,1:999) = -speye(999);
D(:,2:1000) = D(:,2:1000) + speye(999);
xhat = (speye(1000) + mu*D’*D) \ xcor;
Solve the least-squares problem (8.10) with the vector xcor defined in ch8ex9.m, for three
values of µ: µ = 1, µ = 100, and µ = 10000. Plot the three reconstructed signals ˆx.
Discuss the effect of µ on the quality of the estimate ˆx.
8.10 Suppose you are given m (m ≥ 2) straight lines
Li = {pi + tiqi | ti ∈ R}, i = 1,...,m
in Rn. Each line is defined by two n-vectors pi, qi. The vector pi is a point on the line; the
vector qi specifies the direction. We assume that the vectors qi are normalized (∥qi∥ = 1)
and that at least two of them are linearly independent. (In other words, the vectors qi are
not all scalar multiples of the same vector, so the lines are not all parallel.) We denote by
di(y) = min ui∈Li
∥y − ui∥ = minti
∥y − pi − tiqi∥
the distance of a point y to the line Li.
Express the following problem as a linear least-squares problem. Find the point y ∈ Rn
that minimizes the sum of its squared distances to the m lines, i.e., find the solution of
the optimization problem
minimize "m
i=1
di(y)
2
with variable y. Express the least-squares problem in the standard form
minimize ∥Ax − b∥2
with a left-invertible (zero nullspace) matrix A.
(a) Clearly state what the variables x in the least-squares problem are and how A and
b are defined.
132 8 Linear least-squares
(b) Explain why A has a zero nullspace.
8.11 In this problem we use least-squares to fit a circle to given points (ui, vi) in a plane, as
shown in the figure.
We use (uc, vc) to denote the center of the circle and R for its radius. A point (u, v) is on
the circle if (u − uc)
2 + (v − vc)
2 = R2. We can therefore formulate the fitting problem as
minimize "m
i=1
+
(ui − uc)
2 + (vi − vc)
2 − R2,2
with variables uc, vc, R.
Show that this can be written as a linear least-squares problem if we make a change of
variables and use as variables uc, vc, and w = u2
c + v2
c − R2.
(a) Define A, b, and x in the equivalent linear least-squares formulation.
(b) Show that the optimal solution uc, vc, w of the least-squares problem satisfies u2
c +
v2
c − w ≥ 0. (This is necessary to compute R = √u2
c + v2
c − w from the result uc,
vc, w.)
Test your formulation on the problem data in the file ch8ex11.m on the course website.
The commands
[u,v] = ch8ex11
plot(u, v, ’o’);
axis square
will create a plot of the m = 50 points (ui, vi) in the figure. The following code plots the
50 points and the computed circle.
t = linspace(0, 2*pi, 1000);
plot(u, v, ’o’, R * cos(t) + uc, R * sin(t) + vc, ’-’);
axis square
(assuming your MATLAB variables are called uc, vc, and R).
The solution of a least-squares problem
8.12 Let A be a left-invertible m × n-matrix.
(a) Show that the (m + n) × (m + n) matrix
# I A
AT 0
$
is nonsingular.
4)
HW1Prob4.m
[u,v] = HW1Prob4