HW 6 for CS 4220

You may (and probably should) talk about problems with the each other, with the TAs, and with me, providing attribution for any good ideas you might get. Your final write-up should be your own.

md"""

# HW 6 for CS 4220

You may (and probably should) talk about problems with the each other, with the TAs, and with me, providing attribution for any good ideas you might get. Your final write-up should be your own.

"""

2.4 ms

using LinearAlgebra

318 μs

using Plots

3.0 s

1. Non-negative least squares

Consider the non-negative least squares problem

$$\mbox{minimize } \frac{1}{2} \|Ax-b\|^2 \mbox{ s.t. } \forall i, x_i \geq 0$$

Write a projected gradient iteration to solve this problem for the instance below (I used 5000 iterations with a fixed step size of $10^{-4}$).
Based on your approximate solution, do a "polishing" step where you assume the constraint is active for components that are sufficiently small. This should allow you to compute the solution to machine tolerance.
Write a function to compute a residual for the KKT conditions. You should check all parts of the KKT conditions: stationarity, primal feasibility, dual feasibility, and complementary slackness. What are the residual norms before and after polishing?

md"""

## 1. Non-negative least squares

Consider the non-negative least squares problem

$$\mbox{minimize } \frac{1}{2} \|Ax-b\|^2 \mbox{ s.t. } \forall i, x_i \geq 0$$

- Write a projected gradient iteration to solve this problem for the instance below (I used 5000 iterations with a fixed step size of $10^{-4}$).

- Based on your approximate solution, do a "polishing" step where you assume the constraint is active for components that are sufficiently small. This should allow you to compute the solution to machine tolerance.

- Write a function to compute a residual for the KKT conditions. You should check all parts of the KKT conditions: stationarity, primal feasibility, dual feasibility, and complementary slackness. What are the residual norms before and after polishing?

"""

368 μs

let

A = [9.0 4.0 7.0 5.0;

1.0 5.0 1.0 4.0;

9.0 6.0 9.0 0.0;

8.0 6.0 7.0 7.0;

3.0 1.0 1.0 6.0]

b = [29.2; 41.0; 22.4; 32.3; 20.9]

function kkt_resid(A, x, b)

# TODO: Return KKT residual norms (don't worry about making them relative)

end

# Projected gradient iteration

x = max.(A\b, 0.0) # Initial guess

for k = 1:5000

# TODO: Projected gradient steps

end

# TODO: Polish the solution

# TODO: Print the solutions and KKT residual norms before/after polishing

end

3.8 ms

2. Continuation

Write a continuation routine to trace out the solution to

$$x^2-y^2 + (x^7+y^7)/8 + (x^4+y^4)/4 = 0.5$$

from about $x = -3$ to about $x = 3$. Use pseudoarclength continuation, and space the points out by about $h = 0.01$.

md"""

## 2. Continuation

Write a continuation routine to trace out the solution to

$$x^2-y^2 + (x^7+y^7)/8 + (x^4+y^4)/4 = 0.5$$

from about $x = -3$ to about $x = 3$. Use pseudoarclength continuation, and space the points out by about $h = 0.01$.

"""

18.0 ms

let

f(x, y) = x^2 - y^2 + (x^7 + y^7)/8 + (x^4 + y^4)/4 - 0.5

df(x, y) = [ 2*x + 7*x^6/8 + x^3; -2*y + 7*y^6/8 + y^3 ]

# Set up storage for the (x,y) pairs along the curve

xs = []

ys = []

# Find initial point at x = -3.0

# TODO

# Pseudo-arclength continuation steps up to x > 3.0

# TODO

# Plot the computed curve together with the zero set from the contour plotter

plot(xs, ys, linewidth=2, linecolor=:black, legend=false)

contour!(range(-3.0,3.0, length=100), range(-3.0,3.0, length=100), f,

levels=[0.0], linewidth=2, linecolor=:red, linestyle=:dash)

end

340 ms