Midterm

You may use whatever reading materials you want – the course notes should be sufficient, but you may also look online or at books. Make sure you cite any sources beyond the course notes. For this exam, you may ask questions of the course staff, but you should not work with others, whether inside or outside the class.

md"""

# Midterm

You may use whatever reading materials you want -- the course notes should be

sufficient, but you may also look online or at books. Make sure you cite any

sources beyond the course notes. For this exam, you may ask questions of the

course staff, but you should not work with others, whether inside or outside the

class.

"""

212 μs

using LinearAlgebra

273 μs

1: Efficient operations (5 points)

Rewrite each expression to have the indicated complexity.

md"""

#### 1: Efficient operations (5 points)

Rewrite each expression to have the indicated complexity.

"""

163 μs

p1d (generic function with 1 method)

begin

# 1 point: Rewrite to be O(1) time

function p1a(x, k)

ek = zeros(length(x))

ek[k] = 1.0

ek'*x

end

# 1 point: Rewrite to be O(n) time

function p1b(d, u, v, x)

(Diagonal(d)+u*v')*x

end

# 1 point: Rewrite to be O(n^2) time

function p1c(A, u, v, x)

A*u*v'*A*x

end

# 2 points: Rewrite this to take O(mn^2) time rather than O((m+n)^3)

function p1d(A, b)

m, n = size(A)

M = [I A; A' zeros(n,n)]

rx = M\[b; zeros(n)]

end

1.2 ms

2: Floating point fiddling (5 points)

Rewrite each expression or better floating point error.

md"""

#### 2: Floating point fiddling (5 points)

Rewrite each expression or better floating point error.

"""

168 μs

p2c (generic function with 1 method)

begin

# 1 point: Rewrite for good accuracy when z is near zero

function p2a(z)

1.0-cos(z)

end

# 1 point: Rewrite for good accuracy when z is large

function p2b(z)

1.0/(1.0+z)+1.0/(1.0-z)

end

# 3 points: Rewrite to avoid overflow and to maintain good accuracy for small U;

# you may assume U is upper triangular and all non-negative. One point here

# is for complexity -- you should be able to do this in O(n) time.

# Hint: you may want to look up log1p (the inverse of expm1 from HW)

function p2c(U)

log(det(I+U))

end

684 μs

3: Normwise nonsense (5 points)

(1 point) Argue that $∥ A x ∥_{\infty} \leq (max_{i j} | a_{i j} |) ∥ x ∥_{1}$ .
(2 points) If $P A = L U$ is computed using the usual partial pivoting strategy, argue that $∥ A ∥_{1} \leq n ∥ U ∥_{1}$ .
(2 points) The matrix exponential can be defined formally as $\exp (A) = \sum_{k = 0}^{\infty} \frac{1}{k!} A^{k}$ . Argue that for any operator norm, $∥ \exp (A) ∥ \leq \exp (∥ A ∥)$ .

md"""

#### 3: Normwise nonsense (5 points)

1. (1 point) Argue that $\|Ax\|_\infty \leq \left( \max_{ij} |a_{ij}| \right) \|x\|_1$.

2. (2 points) If $PA = LU$ is computed using the usual partial pivoting strategy, argue that $\|A\|_1 \leq n \|U\|_1$.

3. (2 points) The matrix exponential can be defined formally as $\exp(A) = \sum_{k=0}^\infty \frac{1}{k!} A^k$. Argue that for any operator norm, $\|\exp(A)\| \leq \exp(\|A\|)$.

"""

1.5 ms

4: Delightful derivatives (4 points)

Let $A \in R^{m \times n}$ be well conditioned with $m > n$ , and consider $x (s) = (A + s E)^{†} b$ for some given $E \in R^{m \times n}$ . Complete the following code for computing $x (0)$ and $d x (0) / d s$ (you may use the normal equations approach).

md"""

#### 4: Delightful derivatives (4 points)

Let $A \in \mathbb{R}^{m \times n}$ be well conditioned with $m > n$, and consider

$x(s) = (A+sE)^\dagger b$ for some given $E \in \mathbb{R}^{m \times n}$.

Complete the following code for computing $x(0)$ and $dx(0)/ds$ (you may use the normal equations approach).

"""

196 μs

deriv_ls (generic function with 1 method)

function deriv_ls(A, E, b)

# Setup: your code should include these two lines

m, n = size(A)

F = cholesky(A'*A)

# TO DO: Replace the following two lines with your code

x = A\b # Rewrite for O(mn) time using F (1 point)

dx = zeros(n) # Rewrite to be correct and take O(mn) time (3 points)

x, dx

end

710 μs

check_deriv_ls (generic function with 1 method)

function check_deriv_ls()

# Finite difference check on derivative

h = 1e-6

A = rand(10,3)

E = rand(10,3)

b = rand(10)

xp = (A+h*E)\b

xm = (A-h*E)\b

dx_fd = (xp-xm)/(2h)

x, dx = deriv_ls(A, E, b)

norm(dx-dx_fd)/norm(dx)

end

711 μs

Inf

check_deriv_ls()

102 μs

5: Extending Cholesky (3 points)

Consider the block matrix

$M = [\begin{matrix} A & b \\ b^{T} & d \end{matrix}]$

Complete the following code to extend a Cholesky factorization $A = R^{T} R$ to a Cholesky factor of $M$ . Your code should take $O (n^{2})$ time where $A \in R^{n \times n}$

md"""

#### 5: Extending Cholesky (3 points)

Consider the block matrix

$$M = \begin{bmatrix} A & b \\ b^T & d \end{bmatrix}$$

Complete the following code to extend a Cholesky factorization $A = R^T R$ to a Cholesky factor of $M$. Your code should take $O(n^2)$ time where $A \in \mathbb{R}^{n \times n}$

"""

200 μs

extend_cholesky (generic function with 1 method)

function extend_cholesky(R, b, d)

# Placeholder: does the job, but costs O(n^3)

# Replace with your own code (3 points)

M = [R'*R b; b' d]

F = cholesky(M)

end

384 μs

You may use the following code to sanity check your solution.

md"""

You may use the following code to sanity check your solution.

"""

128 μs

test_extend_cholesky (generic function with 1 method)

function test_extend_cholesky()

W = rand(10,6)

M = W'*W

A = M[1:5,1:5]

b = M[1:5,end]

d = M[end,end]

FA = cholesky(A)

R = extend_cholesky(FA.U, b, d)

norm(R'*R-M)/norm(M)

end

686 μs

5.581511538241275e-17

test_extend_cholesky()

147 ms

6: Least squares limbo (5 points)

Consider the least squares problem of minimizing

$ϕ (x, y) = ∥ A x + B y - d ∥^{2}$

and suppose you are given a code solveAls such that solveAls(d) computes $A^{†} d$ . In this problem, we will solve the least squares problem for $x$ and $y$ using this solver as a building block.

(1 point) Write the normal equations for the least squares problem as a block 2-by-2 system.
(1 point) Do block Gaussian elimination to write a Schur complement system that can be solved for $y$
(1 point) Rewrite the Schur complement system so that $A$ is not used directly – only use the combinations $A^{T} B$ and $A^{†}$ appear (and anything involving $B$ on its own is fine).
(2 points) Complete the following Julia code to solve the least squares problem.

md"""

#### 6: Least squares limbo (5 points)

Consider the least squares problem of minimizing

$$\phi(x, y) = \|Ax + By - d\|^2$$

and suppose you are given a code `solveAls` such that `solveAls(d)` computes $A^\dagger d$. In this problem, we will solve the least squares problem for $x$ and $y$ using this solver as a building block.

1. (1 point) Write the normal equations for the least squares problem as a block 2-by-2 system.

2. (1 point) Do block Gaussian elimination to write a Schur complement system that can be solved for $y$

3. (1 point) Rewrite the Schur complement system so that $A$ is not used directly -- only use the combinations $A^T B$ and $A^\dagger$ appear (and anything involving $B$ on its own is fine).

4. (2 points) Complete the following Julia code to solve the least squares problem.

"""

441 μs

solveABls (generic function with 1 method)

function solveABls(solveAls, ATB, B, d)

# Placeholder: Replace these lines!

x = solveAls(d)

y = B\d

# Return the solution.

x, y

end

330 μs

You may use the following test to sanity check your solution.

md"""

You may use the following test to sanity check your solution.

"""

119 μs

test_solveABls (generic function with 1 method)

function test_solveABls()

A = rand(10,4)

B = rand(10,3)

d = rand(10)

x, y = solveABls((d) -> A\d, A'*B, B, d)

xyref = [A B]\d

norm(xyref - [x; y])/norm(xyref)

end

900 μs

1.0249950910118462

test_solveABls()

51.2 μs

7: Your turn

(1 point) Share one thing in the class you think is working well.
(1 point) Share one thing you think could work better (concrete recommendations are great!).
(1 point) What do you consider the most difficult concept from the first part of the course?

md"""

#### 7: Your turn

1. (1 point) Share one thing in the class you think is working well.

2. (1 point) Share one thing you think could work better (concrete recommendations are great!).

3. (1 point) What do you consider the most difficult concept from the first part of the course?

"""

274 μs