Description
In this assignment, you will implement a Matlab function to decompose a
matrix A into the product of two matrices A = QR, where Q has orthonormal
columns and R is upper triangular. You will then use your decomposition to
solve a shape fitting problem.
You will need the files back_sub.m for Part 4 and generate_data.m and
visualize.m for Part 5. These can be found in the zip file provided on the
Moodle assignment page.
1 Submission Guidelines
You will submit a zip file that contains the following .m files on Moodle:
• ortho_decomp.m
• my_qr.m
• least_squares.m
• my_pack.m
• my_unpack.m
• design_matrix.m
• affine_fit.m
as well as any helper functions that are needed by your implementation.
As with previous assignments, please note the following rules:
1. Each file must be named as specified above.
2. Each function must return the specified value(s).
3. You may use Matlab’s array manipulation syntax, e.g. size(A), A(i,:),
and [A,B], and basic operations like addition, multiplication, and transpose
of matrices and vectors. High-level linear algebra functions such as inv,
qr, and A\b are not allowed except where specified. Please contact the
instructor with further questions.
1
4. This is an individual assignment, and no collaboration is allowed. Any
in-person or online discussion should stop before you start discussing or
designing a solution.
2 Orthogonal decomposition
Suppose you have an orthonormal basis {u1, . . . , un} for a subspace H ⊆ R
m.
Any vector v ∈ R
m can be decomposed into the sum of two orthogonal vectors,
vˆ ∈ H and vˆ
⊥ ∈ H⊥. Furthermore, since vˆ ∈ H, it can be expressed as
vˆ = c1u1 + · · · + cnun for some weights c1, . . . , cn. Implement a function to
perform this decomposition.
Specification:
function [c, v_perp] = ortho_decomp(U, v)
Input: an m × n matrix U =
u1 · · · un
with orthonormal columns, and a
vector v.
Output:
c: a vector c =
c1
.
.
.
cn
such that v = c1u1 + · · · + cnun + vˆ
⊥
v_perp: the residual vector v
⊥ such that ui
· vˆ
⊥ = 0 for all i.
Implementation notes:
Here vˆ is simply the orthogonal projection of v onto the subspace H. Since we
have an orthogonal basis of H, this should be easy to compute. In fact, since
the basis is orthonormal, it is possible to compute the projection in just one line
of code without any for loops. Hint: What are the entries of UT v?
Test cases:
1. U =
1 0
0 0
0 1
0 0
, v =
2
3
4
5
→ c =
2
4
, vˆ
⊥ =
0
3
0
5
2. U =
0.6
0.8
0
, v =
1
1
1
→ c =
1.4
, vˆ
⊥ =
0.16
−0.12
1
2
3. Generate a random orthonormal set of 5 vectors in R
10 using the command [U,~] = qr(randn(10,5),0), and generate another random integer
vector v = randi([-2,2], 10,1). Compute your orthogonal decomposition, [c, v_perp] = ortho_decomp(U, v). Verify that U*c + v_perp is
nearly the same as v, and U’*v_perp is nearly zero (both to within 10−12).
3 QR decomposition
Suppose you have an m × n matrix A which is “tall and thin”, i.e. with m > n,
and the columns of A are linearly independent. The Gram-Schmidt process
corresponds to a factorization A = QR, where Q is an m × n matrix with
orthonormal columns, and R is an n × n upper triangular matrix.1
Specification:
function [Q, R] = my_qr(A)
Input: an m × n matrix A with m > n.
Output: an m × n matrix Q with orthonormal columns, and an n × n upper
triangular matrix R, such that A = QR.
Implementation notes:
I recommend starting your implementation by first having your function compute
Q correctly. After that, you can modify your code to also compute R.
The columns of Q are essentially obtained by performing the Gram-Schmidt
process with normalization on the columns of A:
q1 = a1/ka1k,
aˆ
⊥
2 = a2 − projH1
a2 where H1 = span{q1},
q2 = aˆ
⊥
2 /kaˆ
⊥
2 k,
aˆ
⊥
3 = a3 − projH2
a3 where H2 = span{q1, q2},
q3 = aˆ
⊥
3 /kaˆ
⊥
3 k,
.
.
.
You can carry this out using your orthogonal decomposition function. For each
column i = 1, . . . , n, call ortho_decomp with an appropriate set of inputs, and
1Technically, what we are computing here is known as the “thin” or “reduced” QR decomposition. In the full QR decomposition, Q is square and orthogonal, and R is an m × n upper
triangular matrix whose lower (m − n) rows are all zero. In practice, it is also not computed
with the Gram-Schmidt process, but with other methods that are more numerically stable.
3
use the resulting vˆ
⊥ to fill in the ith column of Q. In the end, the columns of
Q should be an orthonormal set of vectors {q1, . . . , qn} which span Col A.
Once you can compute Q correctly, modify your function to compute R as well.
In principle, since QT Q = I, you could obtain R simply via QT A = QT QR = R.
However, this is more work than necessary, because you can actually fill in the
entries of R as you compute each column of Q. When you compute the ith
column of Q, the ortho_decomp call gives you both c and vˆ
⊥; can you use these
to fill in the i nonzero entries in the ith column of R? Hint: Consider the case
when A is already an upper triangular matrix, say A =
1 2 3
0 4 5
0 0 6
, and work
out by hand what c and vˆ
⊥ will be for each column.
Test cases:
1. A =
0 3 0
0 0 −4
−2 0 0
→ Q =
0 1 0
0 0 −1
−1 0 0
, R =
2 0 0
0 3 0
0 0 4
2. A =
1 1 0
1 0 0
1 1 1
1 0 1
→ Q =
1/2 1/2 −1/2
1/2 −1/2 −1/2
1/2 1/2 1/2
1/2 −1/2 1/2
, R =
2 1 1
0 1 0
0 0 1
3. Generate a random 10 × 5 integer matrix A = randi([-2,2], 10,5)
and compute your QR decomposition, [Q, R] = my_qr(A). Verify that
istriu(R) is true, Q’*Q is nearly the identity matrix, and Q*R is nearly
the same as A (both to within 10−12).
4
4 Least-squares problems
Use the QR decomposition to solve the least-squares problem Ax ≈ b.
Specification:
function x = least_squares(A, b)
Input: an m × n matrix A and a vector b ∈ R
m.
Output: a vector x ∈ R
n which is the least-squares solution of the overdetermined system Ax ≈ b, i.e. such that Ax − b ∈ (Col A)
⊥.
Implementation notes:
Do not solve the normal equations AT Ax = AT b directly, as this is a very
numerically unstable approach. Instead, compute the QR factorization of A and
use it to solve the least-squares problem with only a matrix-vector multiplication
and a back-substitution, as described in the textbook. Use the back_sub function
provided with this assignment to perform the back-substitution.
If you cannot get your my_qr function to work, you may use Matlab’s built-in qr
function here so you can still attempt Part 5. Call qr(A,0) to get a rectangular
Q and square R as desired.
Test cases:
1. A =
1
1
0
, b =
2
3
4
→ x =
2.5
2. A =
1 0 0
1 1 1
1 2 4
1 3 9
1 4 16
, b =
0
3
4
3
0
→ x =
0
4
−1
3. Generate a random 10 × 5 integer matrix A = randi([-2,2], 10,5)
and a random integer vector b = randi([-2,2], 10,1). Compute your
least-squares solution, x = least_squares(A, b), and obtain the residual
r = b – A*x. Verify that A’*r is nearly zero (to within 10−12).
5
5 Best-fitting transformations
An affine transformation is a combination of a linear transformation and a
translation (i.e. displacement) while retaining parallelism, i.e., parallel lines
remain parallel after the transformation. For example, a point
x
y
on A (orange
square) in Figure 1 can be transformed to
x
y
on B (red parallelogram) via a
linear transform M =
m11 m12
m21 m22
, i.e.,
x
y
=
m11 m12
m21 m22 x
y
.
Note that the parallel lines stay parallel. This transform can be combined with
a translation, moving the red parallelogram to blue parallelogram (C). This
composite transformation can be written as:
xe
ye
=
x
y
+
tx
ty
=
m11 m12
m21 m22 x
y
+
tx
ty
= M
x
y
+ t, (1)
where t =
tx
ty
is the translation vector.
11 12
2 21 22 1
2 m m 1
m
x
y m y
x
=
11 12
21 22
1
1
x
y
m m t
m y t
x
y m
x
= +
A
B
C
Figure 1: Affine transform.
In sum, Equation (1) can be written as
xe
ye
=
m11x + m12y + tx
m22x + m22y + ty
. (2)
6
5.1 Linear Equation from a Single Correspondence
Your task is given many correspondences2
(xi
, yi) ↔ (xei
, yei) where i is the
index for the correspondence, to compute the best affine transform parameters,
m11, m12, m21, m22, tx, ty (unknowns): We can rewrite Equation (2) by arranging
it with respect to the unknowns for, say, the first correspondence:
xe1
ye1
| {z }
pe1
=
x1 y1 0 0 1 0
0 0 x1 y1 0 1
| {z }
A1
m11
m12
m21
m22
tx
ty
| {z }
β
. (3)
or simply,
pe1 = A1β. (4)
Note that A1 depends on the original position of the point, p1 =
x1
y1
. The
correspondence produces 2 linear equations while the number of unknowns is
6, so at least 3 correspondences are needed to uniquely determine the affine
transform parameters.
Implement a function beta = my_pack(M, t) to obtain β from M =
m11 m12
m21 m22
and t =
tx
ty
, and its inverse function [M, t] = my_unpack(beta).
function beta = my_pack(M, t)
Input: a 2 × 2 matrix M and a vector t.
Output: a vector β ∈ R
6
containing the entries of M and t.
function [M, t] = my_unpack(beta)
Input: a vector β ∈ R
6
.
Output: a 2 × 2 matrix M and a vector t ∈ R
2
such that my_pack(M, t) = β.
Finally, construct the matrix Ai given pi
.
function Ai = design_matrix(pi)
Input: a vector pi =
xi
yi
∈ R
2
.
Output: the matrix Ai given by Equation (3).
2A point (x, y) in A is transformed to a point (x, e ye) in C. This pair of points forms a
correspondence.
7
5.2 Solving Linear System from n Correspondences
Equation (3) can be extended to include n correspondences as follow:
xe1
ye1
.
.
.
xen
yen
=
x1 y1 0 0 1 0
0 0 x1 y1 0 1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
xn yn 0 0 1 0
0 0 xn yn 0 1
m11
m12
m21
m22
tx
ty
, (5)
or again simply,
pe1
.
.
.
pen
=
A1
.
.
.
An
β, (6)
If all correspondences are noise-free, Equation (6) will be always satisfied. Due
to correspondence noise in practice, it cannot be satisfied, and therefore,
pe1
.
.
.
pen
≈
A1
.
.
.
An
β, (7)
Now we will compute the best affine transform parameters using Equation (7).
function [M, t] = affine_fit(P, P_tilde)
Input: 2×k correspondence matrices P =
p1 · · · pk
and Pe =
pe1 · · · pek
containing the reference points and transformed points as columns.
Output: a 2 × 2 matrix M and a vector t ∈ R
2
such that the affine transformation Mx + t best matches the given data. To compute this, you will have
to set up and solve the least-squares problem in Equation (7) to obtain β, then
unpack it to get M and t out.
Test cases:
For the first two tests, pick an arbitrary M and t, say M =
1 2
3 4
and t =
5
6
.
1. Check that [M_new, t_new] = my_unpack(my_pack(M, t)) recovers the
original values of M and t.
2. Verify that when x =
0
0
, design_matrix(x)*my_pack(M,t) gives back
t. Similarly, x =
1
0
and x =
0
1
should give you m1 + t and m2 + t
respectively, where m1 and m2 are the columns of M.
8
3. Applying affine_fit to P =
0 1 1 0
0 0 1 1
, Pe =
0.4 1 1.8 1.2
0.8 0 0.6 1.4
should give M =
0.6 0.8
−0.8 0.6
and t =
0.4
0.8
: a rotation and a translation.
p1
pe1
p2 = pe2
p3
pe3
p4
pe4
4. We have provided a function [P, P_tilde, M, t] = generate_data()
that produces some random test data. It does so by filling random values
in M and t, choosing random points pi
, then setting each pei to Mpi + t
plus a small amount of random noise. You can visualize the data by calling
the provided function visualize(P, P_tilde).
Call affine_fit(P, P_tilde) to obtain your own estimates of M and t.
The result should be close to the “ground truth” transformation returned
by generate_data, although due to the added noise it will not be exactly
identical. Call visualize(P, P_tilde, M, t) to see what your fit looks
like.
Note: If you want to test on a smaller data set, just use the first few columns
of P and Pe, for example P = P(:, 1:10) and P_tilde = P_tilde(:, 1:10).
9
