max constraints in JuMP

Similar to this question, how can I formulate a constraint involving max of the form \sum_{i=1}^m [z_i + t]_+ \leq c where [x]_+ = \max(x,0), both t,c \in \mathbb{R}, and z = My \in \mathbb{R}^m for variable y \in \mathbb{R}^n with data matrix M \in \mathbb{R}^{m \times n}. My attempt is

m, n = size(M)
@variable(model, y[1:n])
@variable(model, t)
zz = [M[i,:]'*y + t for i=1:m]
z = [max(zz[i], 0.0) for i=1:m]
@constraint(model, sum(z) <= c)

but I got the following error:

MethodError: no method matching isless(::Float64, ::JuMP.GenericAffExpr{Float64,JuMP.Variable})

I tried adding a non-negative variable s and doing the following

m, n = size(M)
@variable(model, y[1:n])
@variable(model, t)
@variable(model, s >= 0.0)
z = [M[i,:]'*y + t for i=1:m]
@constraint(model, z >= s)
@constraint(model, sum(z) <= c)

but wasnâ€™t sure if this is the standard way of adding this type of constraint. For example, is this what Convex.jl does, or where could I check?

This makes no sense to me. yy is a bunch of expressions, and the â€śconstraintâ€ť is not really constraining anything. Also data[i] should probably be data[:,i]. In the question you referred to, the OP was trying to minimize the sum of terms each of which is of the form max(expr, 0), so the standard way is to define z and constrain z >= expr and z >= 0, then minimize z instead, which due to the minimization behavior of the optimization algorithm, z will eventually land on either expr or 0, and in the process expr will be pushed as low as possible until it gets to or below 0, then z will not push it down anymore.

1 Like

Thanks, my original question was offâ€¦ Edited to reflect my issue better

@variable(model, s[1:m] >= 0.0)

@constraint(model, [i=1:m], s[i] >= z[i])

@constraint(model, sum(s[i] for i=1:m) <= c)

Note that under this transformation, there is no guarantee that s[i] will be exactly equal to max(z[i], 0) (may or may not be the case depending on your objective). s[i] can take any value greater than or equal to max(z[i], 0) as long as sum(s) is less than or equal to c. But satisfying the transformed constraints clearly guarantees satisfying your original constraints. That is for every feasible solution (y_f,t_f,s_f) to the transformed constraints, the solution (y_f, t_f) is guaranteed to satisfy your original constraints. And conversely, for every feasible solution (y_f,t_f) to the original constraints, there exists at least 1 feasible solution (y_f,t_f,s_f) to the transformed constraints by setting s_f[i] = max(z[i], 0) for example, so you are not â€śmissingâ€ť on any feasible solution to the original constraints by replacing them with the transformed set of constraints.

2 Likes

Thanks â€“ the feasibility argument makes sense, but why did you constrain the sum of s[i] and not z?. Is there a way to doe the max(z[i] + t, 0) constraint precisely? In Convex.jl, I was able to solve the toy problem:

using Convex
using Mosek
m = 10
n = 3
t = Variable(1)
y = Variable(n)
M = randn(m,n)
c = 100
r = randn(n)
problem = minimize(r'*y)
problem.constraints += [t >= 0]
z = max.([M[i,:]'*y + t for i=1:m], 0.0)
problem.constraints += [sum(z) <= c]
solve!(problem, MosekSolver(),verbose=true)
y.value

but donâ€™t entirely know how Convex is treating that max constraint.

The above transformation will do precisely what you want, just ignore s in the final solution, it doesnâ€™t add any semantics to your model. As for Convex.jl, I havenâ€™t read the src but one possibility is that Convex.jl is doing the same transformation under the hood and is hiding away the s details from you, so you only see what you want to see. Or another possibility is that the max expression is processed as a non-smooth function using sub-gradients, I have no idea which one it is using.

1 Like