Sequentially solve a quadratic programming problem

[AUK1939]

Related to #398

I'm trying to solve a quadratic problem sequentially but only the first problem gets solved.

If I run the above code from 398

n=5
b = randn(n) 
x = Variable(n)
H=Semidefinite(n)
Hval=randn(n,n)
Hval.=Hval'*Hval+10*diagm(ones(n)) # symmetric positive definite
fix!(H,Hval)
problem = minimize(x'b +quadform(x,H), [x >= 0])
solve!(problem, SCS.Optimizer)

The solution is

 julia> x.value
 5×1 Matrix{Float64}:
 -1.1415298336189401e-5
 0.051775607935912626
 -2.338104123291265e-5
 -1.764813180548554e-5
 -1.8960736530572698e-5

Then I create a new H matrix and fix! it and solve again

Hval2=randn(n,n)
Hval2.=Hval2'*Hval2+10*diagm(ones(n)) # symmetric positive definite
fix!(H,Hval2)
solve!(problem, SCS.Optimizer)

But I get the exact same solution (below) as the first problem. Seems like the Hval2 value was ignored. I noticed this on versions 0.15.4 and the master branch (v0.16). Have I missed something or is this a bug?

 julia> x.value
5×1 Matrix{Float64}:
-1.1415298336189401e-5
0.051775607935912626
-2.338104123291265e-5
-1.764813180548554e-5
-1.8960736530572698e-5

[odow]

This looks like a bug:

julia> using Convex, LinearAlgebra, SCS

julia> function new_hval(n)
           Hval = randn(n, n)
           Hval .= Hval' * Hval + 10 * LinearAlgebra.I(n)
           return Hval
       end
new_hval (generic function with 1 method)

julia> function solve_hval(b, Hval)
           n = size(Hval, 1)
           x = Variable(n)
           H = Semidefinite(n)
           fix!(H, Hval)
           problem = minimize(x' * b + quadform(x, H), [x >= 0])
           solve!(problem, SCS.Optimizer; silent_solver = true)
           return x.value
       end
solve_hval (generic function with 2 methods)

julia> n = 2
2

julia> b = randn(n);

julia> Hval = new_hval(n)
2×2 Matrix{Float64}:
 10.4528      0.0800298
  0.0800298  10.0142

julia> Hval2 = new_hval(n)
2×2 Matrix{Float64}:
 10.6023    -0.095015
 -0.095015  10.0207

julia> solve_hval(b, Hval)
2×1 Matrix{Float64}:
 -2.224469395829347e-5
  0.05039542808465051

julia> solve_hval(b, Hval2)
2×1 Matrix{Float64}:
 -2.799539900451147e-5
  0.05036184395370259

julia> x = Variable(n);

julia> H = Semidefinite(n);

julia> fix!(H, Hval);

julia> problem = minimize(x' * b + quadform(x, H), [x >= 0]);

julia> solve!(problem, SCS.Optimizer; silent_solver = true)

julia> x.value
2×1 Matrix{Float64}:
 -2.224469395829347e-5
  0.05039542808465051

julia> fix!(H, Hval2);

julia> solve!(problem, SCS.Optimizer; silent_solver = true)

julia> x.value
2×1 Matrix{Float64}:
 -2.224469395829347e-5
  0.05039542808465051

[odow]

So the bad news is that this is never valid (specifically for quadform, because of the way it is implemented), and there's no simple fix #586 introduces an error message for this.

I'd do something like:

julia> using Convex

julia> import SCS

julia> n = 2;

julia> b = [-2.0, 1.0];

julia> A = [1.0 -0.5; -0.5 1.0];

julia> B = [2.0 -0.5; -0.5 2.0];

julia> x = Variable(n);

julia> constraints = [x >= 0];

julia> obj = x' * b;

julia> problem = minimize(obj + quadform(x, A), constraints);

julia> solve!(problem, SCS.Optimizer; silent_solver = true)

julia> x.value
2×1 Matrix{Float64}:
 1.0000045037700163
 8.936981015072835e-6

julia> problem = minimize(obj + quadform(x, B), constraints);

julia> solve!(problem, SCS.Optimizer; silent_solver = true)

julia> x.value
2×1 Matrix{Float64}:
  0.49999652105737485
 -1.0541790264797623e-5

[AUK1939]

Thanks - I'm not sure how much performance gains are to be had by using fix! anyway

[ericphanson]

I just want to apologize, this is a really bad bug! My "fix" (#399) made the error in #398 go away but did not actually make it work correctly, so we went from a error message to a silent correctness bug 😱. I hope no results were compromised by this...

[AUK1939]

Thanks, but honestly no need - there was no harm done. I was very much in the testing phase. Running the optimization (building the problem) in a loop was taking up a lot of memory so I thought I'd try the fix!(v, x) approach.

I was initially using JuMP.jl but preferred Convex.jl since my problem could be expressed much more succintly and warm starts were a breeze. Look forward to further developments. Keep up the good work.