All of the examples can be found in Jupyter notebook form here.

# 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 (i.e. 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
if VERSION < v"1.2.0-DEV.0"
LinearAlgebra.diagm(v::AbstractVector) = diagm(0 => v)
end

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: 451…878

Case 1: Nominal optimal solution

p = minimize(norm(A * x - b, 2))
solve!(p, () -> SCS.Optimizer(verbose=0))
x_nom = evaluate(x)
10-element Array{Float64,1}:
0.14011971595651504
-0.6267576851011559
0.38677870155089944
-2.1196481215776486
2.608178491683313
2.7863929513651997
-1.4229869874573207
3.4250281029277514
1.1881563248418205
0.6209660589666043 

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(verbose=0))
x_stoch = evaluate(x)
10-element Array{Float64,1}:
0.3021775477239739
-0.41551897138372085
0.08738385390693054
-1.4965058790960735
0.6935662856586672
1.562174850658455
-0.9742998235164324
2.5436069628043234
0.509474152335715
0.16229007524299716

Case 3: Worst-case robust approximation

p = minimize(max(norm((A - B) * x - b), norm((A + B) * x - b)))
solve!(p, () -> SCS.Optimizer(verbose=0))
x_wc = evaluate(x)
10-element Array{Float64,1}:
0.4503856705879205
0.12295021498367938
-0.3543861169777943
-1.048771359486388
-0.40928225562595466
0.7864578523377699
-0.5947431035116082
1.4542911733298822
0.5357692245624266
-0.24795958641619295

Plot residuals:

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

errvals(x) = [ norm((A + parvals[k] * B) * x - b) for k = 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")