##
using Plots, LinearAlgebra, Random, StatsBase, SparseArrays
## Least-squares fitting with Convex.jl
# Convex.jl is a convex programming language that allows you
# to specify an objective function easily and solve the
# resulting optimization problem with a variety of different solvers.
# To install Convex.jl, use these lines
# using Pkg
# Pkg.add("Convex")
# Pkg.add("SCS") # one of the solvers that works with Convex.jl
##
using Convex, SCS
## Step 1: Create a set of data with a linear model
m=40
n=2
A = randn(m,n)
xex = [5;1] # "b" = 5 in ax+b and a=1
pts = -10.0 .+ 20*rand(m,1)
A = [ones(m,1) pts]
b = A*xex + .5*randn(m,1)
40×1 Matrix{Float64}:
7.112218015637595
5.008976064454254
10.349890144049471
2.035473353820581
8.02977700056039
4.983462204016269
1.64056258539329
11.94147733642886
9.179368412202574
-3.7486628508050055
12.66784944128789
0.5872959164766413
5.734967709857549
⋮
13.642056131866449
-1.913458342661348
0.45649134439013667
9.39043879443364
6.514474624992099
-1.315791038770222
11.205161377719158
10.401999295562684
13.325416283262845
6.635651950909954
14.569780393084162
6.6127055785568665
## Show the linear model
scatter(pts,b, dpi=300)
xlabel!("x")
ylabel!("y")
##
## The convex.jl problem
x = Variable(n)
problem = minimize(sumsquares(b - A*x))
solve!(problem, SCS.Optimizer)
xls = x.value
@show x
x = Variable
size: (2, 1)
sign: real
vexity: affine
id: 263…089
value: [4.742037141032825, 0.9962500528062521]
------------------------------------------------------------------
SCS v3.2.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 5, constraints m: 46
cones: z: primal zero / dual free vars: 1
l: linear vars: 1
q: soc vars: 44, qsize: 2
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 86, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.71e+01 1.00e+00 1.62e+01 -8.03e+00 1.00e-01 1.35e-04
225| 4.67e-08 6.75e-09 3.52e-07 6.85e+00 3.39e-01 1.30e-03
------------------------------------------------------------------
status: solved
timings: total: 1.30e-03s = setup: 1.04e-04s + solve: 1.19e-03s
lin-sys: 1.23e-04s, cones: 3.80e-05s, accel: 9.42e-04s
------------------------------------------------------------------
objective = 6.853359
------------------------------------------------------------------
Variable size: (2, 1) sign: real vexity: affine id: 263…089 value: [4.742037141032825, 0.9962500528062521]
## Show the least squares fit
scatter(pts,b;label="data")
plot!([-11; 11], [1 -11; 1 11]*xls;label="fit")
#xaxis!([-11 11])
title!("Least-square fit")
xlabel!("x")
ylabel!("y")
## Now we add outliers
outliers = [-9.5; 9]
outvals = [20; -15]
A = [A; ones(length(outliers),1) outliers]
b = [b; outvals]
m = size(A,1)
pts = [pts;outliers]
42×1 Matrix{Float64}:
2.595368847148814
-0.19560042566858904
5.773594728758546
-2.85234985277879
2.7540327514596648
0.44745280288393374
-3.68306773175898
7.393619809972648
3.6693954526433803
-7.879210489562574
8.834756140905004
-4.908972989599915
0.2897349839968921
⋮
-4.480217401691943
4.6875999501597505
1.130579839115608
-6.458210362243415
6.473198498664164
5.876328798249141
8.924650006723724
2.4593716746040943
9.833485820768349
2.172775504499773
-9.5
9.0
## Show the new data
scatter(pts,b)
xlabel!("x")
ylabel!("y")
## Look at the LS fit
x = Variable(n)
problem = minimize(sumsquares(b - A*x))
solve!(problem, SCS.Optimizer)
xls = x.value
plot!([-11; 11], [1 -11; 1 11]*xls;label="fit")
------------------------------------------------------------------
SCS v3.2.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 5, constraints m: 48
cones: z: primal zero / dual free vars: 1
l: linear vars: 1
q: soc vars: 46, qsize: 2
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 90, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 3.57e+01 1.00e+00 2.75e+01 -2.36e+00 1.00e-01 5.66e-05
250| 7.85e+00 2.32e-02 7.14e-01 9.50e+02 2.28e-02 3.51e-04
500| 6.95e-04 1.02e-05 1.60e-03 1.27e+03 2.28e-02 6.08e-04
------------------------------------------------------------------
status: solved
timings: total: 6.09e-04s = setup: 4.66e-05s + solve: 5.62e-04s
lin-sys: 2.37e-04s, cones: 6.44e-05s, accel: 1.12e-04s
------------------------------------------------------------------
objective = 1271.046141
------------------------------------------------------------------
## Solve the Huber problem and look at the fit
x = Variable(n)
problem = minimize(sum(huber(b - A*x)))
solve!(problem, SCS.Optimizer)
xr = x.value
plot!([-11; 11], [1 -11; 1 11]*xr;label="fit_huber")
------------------------------------------------------------------
SCS v3.2.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 171, constraints m: 295
cones: z: primal zero / dual free vars: 43
l: linear vars: 126
q: soc vars: 126, qsize: 42
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 547, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 5.25e+01 1.48e+00 3.12e+03 -1.34e+03 1.00e-01 2.34e-04
125| 4.05e-04 5.46e-05 4.33e-04 1.11e+02 1.00e-01 1.42e-03
------------------------------------------------------------------
status: solved
timings: total: 1.42e-03s = setup: 2.07e-04s + solve: 1.21e-03s
lin-sys: 6.06e-04s, cones: 1.51e-04s, accel: 2.07e-04s
------------------------------------------------------------------
objective = 111.435693
------------------------------------------------------------------