# Robust approximate fitting

Section 6.4.2 Boyd & Vandenberghe "Convex Optimization" Original by Lieven Vandenberghe Adapted for Convex by Joelle Skaf - 10/03/05

Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14 Original CVX code and plots here: http://web.cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/fig6_15.html

Consider the least-squares problem: minimize $\|(A + tB)x - b\|_2$ where $t$ is an uncertain parameter in [-1,1] Three approximate solutions are found:

1. nominal optimal (that is, letting t=0)
2. stochastic robust approximation: minimize $\mathbb{E}\|(A+tB)x - b\|_2$ assuming $u$ is uniformly distributed on [-1,1]. (reduces to minimizing $\mathbb{E} \|(A+tB)x-b\|^2 = \|A*x-b\|^2 + x^TPx$ where $P = \mathbb{E}(t^2) B^TB = (1/3) B^TB$ )
3. worst-case robust approximation: minimize $\mathrm{sup}_{-1\leq u\leq 1} \|(A+tB)x - b\|_2$ (reduces to minimizing $\max\{\|(A-B)x - b\|_2, \|(A+B)x - b\|_2\}$ ).
using Convex, LinearAlgebra, SCS

Input Data

m = 20;
n = 10;
A = randn(m, n);
(U, S, V) = svd(A);
S = diagm(exp10.(range(-1, stop = 1, length = n)));
A = U[:, 1:n] * S * V';

B = randn(m, n);
B = B / norm(B);

b = randn(m, 1);
x = Variable(n)
Variable
size: (10, 1)
sign: real
vexity: affine
id: 731…529

## Case 1: nominal optimal solution

p = minimize(norm(A * x - b, 2))
solve!(p, SCS.Optimizer; silent = true)
Problem statistics
problem is DCP         : true
number of variables    : 1 (10 scalar elements)
number of constraints  : 0 (0 scalar elements)
number of coefficients : 220
number of atoms        : 3

Solution summary
termination status : OPTIMAL
primal status      : FEASIBLE_POINT
dual status        : FEASIBLE_POINT
objective value    : 2.7707

Expression graph
minimize
└─ norm2 (convex; positive)
└─ + (affine; real)
├─ * (affine; real)
│  ├─ …
│  └─ …
└─ 20×1 Matrix{Float64}

x_nom = evaluate(x)
10-element Vector{Float64}:
-6.147557441718144
-3.92517806828471
1.8241081203844198
0.1336847809167995
1.9598446027248864
3.17769814387929
0.1419527766662775
2.7995641000409432
3.5105221463085585
-3.629450477579609

## Case 2: stochastic robust approximation

P = 1 / 3 * B' * B;
p = minimize(square(pos(norm(A * x - b))) + quadform(x, Symmetric(P)))
solve!(p, SCS.Optimizer; silent = true)
Problem statistics
problem is DCP         : true
number of variables    : 1 (10 scalar elements)
number of constraints  : 0 (0 scalar elements)
number of coefficients : 324
number of atoms        : 10

Solution summary
termination status : OPTIMAL
primal status      : FEASIBLE_POINT
dual status        : FEASIBLE_POINT
objective value    : 8.9964

Expression graph
minimize
└─ + (convex; positive)
├─ qol_elem (convex; positive)
│  ├─ max (convex; positive)
│  │  ├─ …
│  │  └─ …
│  └─ [1.0;;]
└─ * (convex; positive)
├─ [1;;]
└─ qol_elem (convex; positive)
├─ …
└─ …

x_stoch = evaluate(x)
10-element Vector{Float64}:
-2.139035198954485
-1.5819018146437105
1.2606246446646314
1.4681316168515772
0.3594838777145165
1.04496829235293
0.13343142310339703
1.2455003891918222
0.9726954793491382
-0.946783630476293

## Case 3: worst-case robust approximation

p = minimize(max(norm((A - B) * x - b), norm((A + B) * x - b)))
solve!(p, SCS.Optimizer; silent = true)
Problem statistics
problem is DCP         : true
number of variables    : 1 (10 scalar elements)
number of constraints  : 0 (0 scalar elements)
number of coefficients : 440
number of atoms        : 7

Solution summary
termination status : OPTIMAL
primal status      : FEASIBLE_POINT
dual status        : FEASIBLE_POINT
objective value    : 3.1074

Expression graph
minimize
└─ max (convex; positive)
├─ norm2 (convex; positive)
│  └─ + (affine; real)
│     ├─ …
│     └─ …
└─ norm2 (convex; positive)
└─ + (affine; real)
├─ …
└─ …

x_wc = evaluate(x)
10-element Vector{Float64}:
-1.4077935065589873
-0.7708528545924042
0.434098606488757
0.7150650375788932
0.3549832867546875
0.5765731953566152
-0.14934048869524402
0.5790725012045513
0.3882743052999983
-0.30233119908541983

## Plots

Here we plot the residuals.

parvals = range(-2, stop = 2, length = 100);

errvals(x) = [norm((A + parvals[k] * B) * x - b) for k in eachindex(parvals)]
errvals_ls = errvals(x_nom)
errvals_stoch = errvals(x_stoch)
errvals_wc = errvals(x_wc)

using Plots
plot(parvals, errvals_ls, label = "Nominal problem")
plot!(parvals, errvals_stoch, label = "Stochastic Robust Approximation")
plot!(parvals, errvals_wc, label = "Worst-Case Robust Approximation")
plot!(
title = "Residual r(u) vs a parameter u for three approximate solutions",
xlabel = "u",
ylabel = "r(u) = ||A(u)x-b||_2",
)