Closed AUK1939 closed 5 months ago
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
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
Thanks - I'm not sure how much performance gains are to be had by using fix! anyway
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...
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.
Related to https://github.com/jump-dev/Convex.jl/issues/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
The solution is
Then I create a new H matrix and fix! it and solve again
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?