I’m working on a set of optimization problems where I need to maximize something that looks a lot like a quadratic form but isn’t. The following example is one such problem:

Maximize f^TAf

where f = [1, x, x^2] and A is a 3x3 PSD matrix. I believe JuMP paired with a semidefinite solver like SCS will provide everything I need for the task, but I can’t seem to get Julia to accept my program. I feel like I should be able to do something like:

model = Model(SCS.Optimizer);
@variable(model, x)
f = [1, x, x^2]
@constraint(model, c1, x <= 1);
@constraint(model, c2, x >= -1);
@objective(model, Max, f'*A*f)
optimize!(model)

but an error is thrown when I define my objective. The particular error is:

ERROR: *(::GenericQuadExpr,::GenericQuadExpr) is not defined. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective. If you are using an `@NL` macro and you encountered this error message, it is because you are attempting to use another unsupported function which calls this method internally.

SCS supports only linear and quadratic functions. Your problem is not convex because it contains cubic and quartic terms.

In addition, even if you could write your f, if A is PSD, then you cannot maximize the objective because it is a convex function. Do you mean Min instead?

You could reformulate your problem to this:

using JuMP, SCS
A = rand(3, 3)
A += A'
model = Model(SCS.Optimizer)
@variable(model, x)
@variable(model, y)
@constraint(model, y >= x^2)
f = [1, x, y]
@objective(model, Min, f'*A*f)
optimize!(model)

but it might not be what you’re trying to achieve.

Apologies, I forgot some critical information (but thank you for your quick response!). x is bounded from above and below; -1<=x<=1. Still, the problem is nonconvex, but with a finite search space there should be a solver that can handle it right? Something like Gloptipoly in Matlab I’m thinking.

julia> using JuMP, Gurobi
julia> A = rand(3, 3)
3×3 Matrix{Float64}:
0.760386 0.407191 0.723833
0.604206 0.364107 0.623264
0.356807 0.60037 0.521208
julia> A += A'
3×3 Matrix{Float64}:
1.52077 1.0114 1.08064
1.0114 0.728214 1.22363
1.08064 1.22363 1.04242
julia> model = Model(Gurobi.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Gurobi
julia> set_attribute(model, "NonConvex", 2)
Set parameter NonConvex to value 2
julia> @variable(model, -1 <= x <= 1)
x
julia> @variable(model, y)
y
julia> @constraint(model, y == x^2)
-x² + y = 0
julia> f = [1, x, y]
3-element Vector{AffExpr}:
1
x
y
julia> @objective(model, Max, f' * A * f)
0.7282140530433079 x² + 2.447268930048988 x*y + 1.0424158567197273 y² + 2.0227956982742628 x + 2.16127995587574 y + 1.520771522270953
julia> optimize!(model)
Set parameter NonConvex to value 2
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])
CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0x6fc2d173
Model has 3 quadratic objective terms
Model has 1 quadratic constraint
Coefficient statistics:
Matrix range [0e+00, 0e+00]
QMatrix range [1e+00, 1e+00]
QLMatrix range [1e+00, 1e+00]
Objective range [2e+00, 2e+00]
QObjective range [1e+00, 5e+00]
Bounds range [1e+00, 1e+00]
RHS range [0e+00, 0e+00]
Continuous model is non-convex -- solving as a MIP
Presolve time: 0.00s
Presolved: 7 rows, 5 columns, 15 nonzeros
Presolved model has 3 bilinear constraint(s)
Variable types: 5 continuous, 0 integer (0 binary)
Root relaxation: objective 9.922746e+00, 2 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
* 0 0 0 9.9227460 9.92275 0.00% - 0s
Explored 1 nodes (2 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 1: 9.92275
Optimal solution found (tolerance 1.00e-04)
Best objective 9.922746016233e+00, best bound 9.922746016233e+00, gap 0.0000%
User-callback calls 100, time in user-callback 0.00 sec
julia> value(x)
1.0