18.06 Problem Set 4

Due Friday, March 4 at 11am.

Problem 1 (10+5)

Suppose, $x,y,a,b$ are $n$-component column vectors. Then we can compute the following scalar $\alpha$:

$$ \alpha = a^T (y x^T + x y^T) b $$

(a) Find a vector $c$ (made from $y,a,b$) such that we get the same $\alpha$ from: $$ \alpha = c^T x $$

(Fun fact: Any scalar that is a linear function of $x$ can be written as a dot product like this. The fancy name for this idea is the Riesz representation theorem.)

(b) Check your answer in Julia for some random vectors. Note that $x^T$ in Julia is written as x' for real column vectors and matrices.

a = [1,2,7,3]
b = [9,-2,4,5]
y = [0,2,3,1]
x = [-3,2,4,-5]

α = ???? # write aᵀ(yxᵀ+xyᵀ)b in Julia

c = ???? # some formula in terms of a, b, y

# this should output "true":
c' * x ≈ α

Problem 2 (10+10+5)

Let $f(x)$ be the function:

$$ f(x) = \Vert (B + xx^T)^{-1} b \Vert $$

where $B$ is an $n \times n$ matrix, and $b$ is an $n$-component column vector, and the input $x$ is an $n$-component column vector.

(a) Give a formula for $df = f(x+dx) - f(x)$, i.e. to first order in $dx$, as a linear operation on $dx$. (Your answer will have matrix inverses in it; that's okay.)

It might be helpful to break this into steps (a.k.a. a chain rule):

1. What is $d(\sqrt{\alpha})$ for a small $d\alpha$, i.e. to first order in the change in the scalar $\alpha$? (Just 18.01!)
2. In class, we saw $d(y^T y) = 2y^T dy$. So what is $d\Vert y\Vert = d\sqrt{y^T y}$ to first order in $dy$, the change for a small change $dy$ in $y$?
3. What is $d(xx^T)$ to first order in $dx$? (Product rule!)
4. In class, we saw that $d(A^{-1}) = -A^{-1} dA A^{-1}$. So what is $d\left[(B + xx^T)^{-1} b \right]$ to first order in $dx$?
5. Put it all together to write $df$ to first order in $dx$.

(b) Using your answer from problem 1, write $df = (\mbox{some vector})^T dx$, and hence give a formula for $\nabla f$.

(c) Check your $\nabla f$ formula by comparing it the numerically computed $f(x+dx) - f(x)$ for a small arbitrary $dx \in \mathbb{R}^4$.

You just need to fill in the "????" in the following code:

using LinearAlgebra

# some B matrix and b vector
B = [3   5   6  6
     2  -2  -5  8
     4   3   2  2
     5  -8  -1  3]
b = [1,-1,2,-2]

f(x) = norm((B + x*x') \ b) # = ‖(B+xxᵀ)⁻¹b‖

# random test inputs:
x = [1,-2,4,6]
dx = [-3,5,0,2] * 1e-8 # a small change in x

f(x + dx) - f(x)

∇f = ???? # your formula from (a)

# this should match f(x+dx)-f(x) above,
# at least to a few digits:
df = ∇f'*dx

Problem 3 (10+5+5+5)

Let $A_0= \begin{pmatrix} 3 & 1 & 2\\ 7 & 3 & -2\\ -2 & -5 & 4\\ \end{pmatrix}$, and consider the nonlinear function $g(x)$ mapping 3-component vector inputs $x$ to 3-component vector outputs: $$ g(x) = \underbrace{\begin{pmatrix} x_1 & & \\ & x_2 & \\ & & x_3 \end{pmatrix}}_{\mathrm{Diagonal}(x)\mbox{ in Julia}} A_0 x + x - \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} $$

Here, you will solve the nonlinear equations $g(x) = \vec{0}$ using Newton's method.

(a) Write down the $3\times 3$ Jacobian matrix $g'(x)$ as a function of $x$.

(b) Write a Julia function J(x) below that computes your Jacobian formula from (a).

(c) Check that g(x+dx) ≈ J(x)*dx for a small dx at some x.

(d) Plug your Jacobian, along with the provided g(x) function, into the Newton iteration below, and use the starting guess $x = (1,1,1)$. Give the resulting solution to at least 5 significant digits (in every component).

using LinearAlgebra

A₀ = [3 1 2; 7 3 -2; -2 -5 4]

g(x) = Diagonal(x) * (A₀*x) + x .- 1

(b) Jacobian:

J(x) = ???? # your answer from (b)

J([1,2,3]) # try it — this should give a 3x3 matrix

(c) random test:

x = randn(3)
dx = randn(3) * 1e-8

g(x+dx) - g(x)

J(x) * dx # this should match g(x+dx) - g(x) to a few digits

(d) Newton iteration

x = [1,1,1] # initial guess

# run 10 iterations of Newton's method
for i = 1:10
    # show some output:
    @show x
    @show g(x)
    println()

    x = x - J(x) \ g(x)
end

Problem 4 (10+5)

In problem 3 of pset 3, you considered $Ax = b$ for $$ A = \begin{pmatrix} 1 & 2 & 3 & 4 \\ 1 & 2 & 4 & 6 \\ 0 & 0 & 1 & 2 \end{pmatrix}, \qquad b = \begin{pmatrix} \alpha \\ 6 \\ 1 \end{pmatrix} $$ and found that it only has solutions for $\alpha = 5$.

(a) Find a basis for $N(A^T)$, the left nullspace.

(b) Because $N(A^T)^\perp = C(A)$, the equation $Ax = b$ only has solutions when $b$ is perpendicular to a basis for $N(A^T)$. Plug in your answer from (a) to check that it is orthogonal to $b$ only when $\alpha = 5$, consistent with your solution in pset 3.

Problem 5 (10)

Suppose $S$ is spanned by (1,2,3) and (4,5,6). Then $S^\perp$ is the nullspace of what matrix?