Cleaning variable coefficients #324
Description
It has been requested that I post the following hack for a clean routine I wrote for polynomials with tunable coefficients in a SOS SDP program (minimal example at the end of this post). In particular, take a polynomial with tunable coefficients, which you would like to optimize the problem over. If the polynomial has coefficients which are monomial in the tunable quantities, the cleanup routine is trivial. However, if the polynomials has coefficients which are polynomial in the tunable quantities, the cleanup routine is not as straightforward due to the structure of the @variable macro. I missed this the first time I wrote the clean-up routine, which was quite important; my problem was able to be solved much more efficiently (I say this loosely) once I discovered this extra clean-up caveat.
This is likely not the most efficient routine to use in this setting, but it would be particularly helpful to either have a note in the documentation that one should be careful how they do this or a reference to a new command for cleanup that does this under the hood. @blegat attempted to do the latter back in March but it didn't get pushed, see here for details.
Thanks, Elizabeth
(The basis is not taken from any real problem but rather is randomly chosen to make terms overlap in some way, which can occur when manipulating multiple functions with tunable coefficients, with constraints on their sum.)
using JuMP
using SumOfSquares
using DynamicPolynomials
using Mosek
using MosekTools
@polyvar x y
basis = vcat(x^2+y^2, x*y, x^3+y^3, x^4, y^4, x+y, x-y, x^2-y^2);
model = SOSModel(Mosek.Optimizer);
MTB_coeffs = @variable(model, [1:length(basis)])
C = vcat(1.2e-9, 1.2, 2.3e-6, 2.2e-10, 3.e-11, 2.1, 4.3, 1.1)
MTB = sum(MTB_coeffs.*C.*basis);
tol = 1e-8;
for i in 1:length(MTB.x)
if length(MTB.a[i].terms.vals) == 1
MTB.a[i].terms.vals = ifelse.((MTB.a[i].terms.vals .< -tol ) .| (MTB.a[i].terms.vals .> tol), MTB.a[i].terms.vals, 0.0)
elseif length(MTB.a[i].terms.vals) >= 2
for ii in 1:length(MTB.a[i].terms.vals)
MTB.a[i].terms.vals[ii] = ifelse.((MTB.a[i].terms.vals[ii] .< -tol ) .| (MTB.a[i].terms.vals[ii] .> tol ), MTB.a[i].terms.vals[ii], 0.0)
end
else
MTB.a[i].constant = ifelse.((MTB.a[i].constant .< -tol ) .| (MTB.a[i].constant .> tol ), MTB.a[i].constant, 0.0);
end
end
MTB = sum(MTB.a.*MTB.x);