initial implementation for norms

examples/robust_uncertainty.jl

@@ -24,26 +24,16 @@ N = ceil((2+2log(2/𝛿))^2) + 1
 
 μhat = rand(d)
 M = rand(d,d)
-# Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')
-Σhat = (M-ones(d)*μhat')'*(M-ones(d)*μhat')
+Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')
 
 m = Model()
 
 @defVar(m, Σ[1:d,1:d], SDP)
 @defVar(m, u[1:d])
 @defVar(m, μ[1:d])
 
-@defVar(m, t1 >= 0)
-@defVar(m, L1[1:d])
-@addConstraint(m, L1 .== (μ-μhat))
-@addConstraint(m, sum{L1[i]^2, i=1:d} <= t1^2)
-@addConstraint(m, t1 <= Γ1(𝛿/2,N))
-
-@defVar(m, t2 >= 0)
-@defVar(m, L2[1:d,1:d])
-@addConstraint(m, L2 .== (Σ-Σhat))
-@addConstraint(m, sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2)
-@addConstraint(m, t2 <= Γ2(𝛿/2,N))
+@addConstraint(m, norm(μ-μhat) <= Γ1(𝛿/2,N))
+@addConstraint(m, vecnorm(Σ-Σhat) <= Γ2(𝛿/2,N))
 
 A = [(1-ɛ)/ɛ (u-μ)';
      (u-μ)     Σ   ]

src/JuMP.jl

@@ -20,8 +20,8 @@ using Compat
 
 export
 # Objects
-    Model, Variable, AffExpr, QuadExpr,
-    LinearConstraint, QuadConstraint, SDPConstraint,
+    Model, Variable, Norm, AffExpr, QuadExpr, NormExpr,
+    LinearConstraint, QuadConstraint, SDPConstraint, NormConstraint,
     ConstraintRef, LinConstrRef,
 # Functions
     # Model related
@@ -59,6 +59,7 @@ type Model
     linconstr#::Vector{LinearConstraint}
     quadconstr
     sosconstr
+    normconstr
     sdpconstr
 
     # Column data
@@ -117,7 +118,8 @@ function Model(;solver=UnsetSolver())
     if !isa(solver,MathProgBase.AbstractMathProgSolver)
         error("solver argument ($solver) must be an AbstractMathProgSolver")
     end
-    Model(zero(QuadExpr),:Min,LinearConstraint[],QuadConstraint[],SOSConstraint[],SDPConstraint[],
+    Model(zero(QuadExpr),:Min,LinearConstraint[],QuadConstraint[],
+          SOSConstraint[],NormConstraint[],SDPConstraint[],
           0,String[],String[],Float64[],Float64[],Symbol[],
           0,Float64[],Float64[],Float64[],nothing,solver,
           false,Any[],nothing,nothing,JuMPContainer[],
@@ -383,6 +385,7 @@ end
 typealias AffExpr GenericAffExpr{Float64,Variable}
 
 AffExpr() = zero(AffExpr)
+AffExpr(x::Union(Number,Variable)) = convert(AffExpr, x)
 
 Base.isempty(a::AffExpr) = (length(a.vars) == 0 && a.constant == 0.)
 Base.convert(::Type{AffExpr}, v::Variable) = AffExpr([v], [1.], 0.)
@@ -491,6 +494,64 @@ end
 
 getValue(arr::Array{QuadExpr}) = map(getValue, arr)
 
+##########################################################################
+# Norm
+# Container for √(∑ expr)
+type GenericNorm{Power,C,V}
+    terms::Vector{GenericAffExpr{C,V}}
+end
+
+typealias Norm{T} GenericNorm{T,Float64,Variable}
+
+Base.norm(x::Vector{Variable}) = vecnorm(x)
+Base.norm{C,V}(x::Array{GenericAffExpr{C,V}}) = vecnorm(x)
+Base.norm{T<:Union(Variable,GenericAffExpr)}(x::JuMPArray{T,1}) = vecnorm(x)
+function Base.norm(x::JuMPDict{Variable})
+    ndims(x) == 1 || error("Cannot use norm() on a multidimensional JuMPDict. Use vecnorm instead.")
+    arr = Array(Variable, length(x))
+    for (it,v) in enumerate(x)
+        arr[it] = v[end]
+    end
+    Norm{2}(arr)
+end
+
+Base.vecnorm(x::Array{Variable}) = Norm{2}(reshape(x,length(x)))
+Base.vecnorm{C,V}(x::Array{GenericAffExpr{C,V}}) = GenericNorm{2,C,V}(reshape(x,length(x)))
+Base.vecnorm{T<:Union(Variable,GenericAffExpr)}(x::JuMPArray{T}) =
+    Norm{2}(collect(x.innerArray))
+function Base.vecnorm(x::JuMPDict{Variable})
+    arr = Array(Variable, length(x))
+    for (it,v) in enumerate(x)
+        arr[it] = v[end]
+    end
+    Norm{2}(arr)
+end
+
+Base.copy{Power,C,V}(x::GenericNorm{Power,C,V}) = GenericNorm{Power,C,V}(copy(x.terms))
+
+Base.convert{Power,C,V}(::Type{GenericNorm{Power,C,V}}, x::Array) =
+    GenericNorm{Power,C,V}(convert(Vector{GenericAffExpr{C,V}}, vec(x)))
+
+##########################################################################
+# NormExpr
+# Container for expressions containing Norms and AffExprs
+type GenericNormExpr{Power,C,V}
+    norm::GenericNorm{Power,C,V}
+    coeff::C
+    aff::GenericAffExpr{C,V}
+end
+
+GenericNormExpr{Power,C,V}(norm::GenericNorm{Power,C,V}) =
+    GenericNormExpr{Power,C,V}(norm, one(C), zero(GenericAffExpr{C,V}))
+
+typealias NormExpr{Power} GenericNormExpr{Power,Float64,Variable}
+
+Base.copy{Power,C,V}(x::GenericNormExpr{Power,C,V}) =
+    GenericNormExpr{Power,C,V}(copy(x.norm), copy(x.coeff), copy(x.aff))
+
+Base.convert{Power,C,V}(::Type{GenericNormExpr{Power,C,V}}, x::GenericNorm{Power,C,V}) =
+    GenericNormExpr{Power,C,V}(x, one(C), zero(GenericAffExpr{C,V}))
+
 ##########################################################################
 # JuMPConstraint
 # abstract base for constraint types
@@ -703,6 +764,22 @@ function Base.copy(c::QuadConstraint, new_model::Model)
     return QuadConstraint(copy(c.terms, new_model), c.sense)
 end
 
+##########################################################################
+# SOCConstraint is a second-order cone constraint of the form
+# ||Ax-b||₂ + cᵀx + d ≤ 0
+
+type GenericNormConstraint{Power,T<:GenericNormExpr} <: JuMPConstraint
+    normexpr::T
+end
+
+typealias NormConstraint{Power} GenericNormConstraint{Power,NormExpr}
+
+function addConstraint{T<:GenericNormConstraint}(m::Model, c::T)
+    push!(m.normconstr,c)
+    m.internalModelLoaded = false
+    ConstraintRef{T}(m,length(m.normconstr))
+end
+
 ##########################################################################
 # ConstraintRef
 # Reference to a constraint for retrieving solution info
@@ -808,6 +885,8 @@ end
 
 ##########################################################################
 # Operator overloads
+typealias JuMPTypes Union(Variable,Norm,AffExpr,QuadExpr,NormExpr)
+
 include("operators.jl")
 if VERSION > v"0.4-"
     include(joinpath("v0.4","concatenation.jl"))
@@ -816,10 +895,10 @@ else
 end
 # Writers - we support MPS (MILP + QuadObj), LP (MILP)
 include("writers.jl")
-# Solvers
-include("solvers.jl")
 # Macros - @defVar, sum{}, etc.
 include("macros.jl")
+# Solvers
+include("solvers.jl")
 # Callbacks - lazy, cuts, ...
 include("callbacks.jl")
 # Pretty-printing of JuMP-defined types.

src/JuMPDict.jl

@@ -165,7 +165,6 @@ Base.size(x::JuMPArray)   = size(x.innerArray)
 Base.size(x::JuMPArray,k) = size(x.innerArray,k)
 Base.issym(x::JuMPArray) = issym(x.innerArray)
 Base.trace(x::OneIndexedArray) = trace(x.innerArray)
-Base.norm(x::OneIndexedArray) = norm(x.innerArray)
 Base.diag(x::OneIndexedArray) = diag(x.innerArray)
 Base.diagm{T}(x::JuMPArray{T,1,true}) = diagm(x.innerArray)
 

src/macros.jl

@@ -190,6 +190,19 @@ end
 
 _construct_constraint!(quad::QuadExpr, sense::Symbol) = QuadConstraint(quad, sense)
 
+function _construct_constraint!{Power}(normexpr::NormExpr{Power}, sense::Symbol)
+    # check that the constraint is SOC representable
+    if sense == :(<=)
+        normexpr.coeff >= 0 || error("Constraint $normexpr <= 0 is not representable. Check that it is convex.")
+        NormConstraint{Power}(normexpr)
+    elseif sense == :(>=)
+        normexpr.coeff <= 0 || error("Constraint $normexpr >= 0 is not representable. Check that it is convex.")
+        NormConstraint{Power}(-normexpr)
+    else
+        error("Invalid sense $sense is SOC constraint")
+    end
+end
+
 _construct_constraint!(x::Array, sense::Symbol) = map(c->_construct_constraint!(c,sense), x)
 
 _vectorize_like(x::Number, y::Array{AffExpr}) = fill(x, size(y))

src/operators.jl

@@ -22,6 +22,10 @@
 (+)(lhs::Number, rhs::Variable) = AffExpr([rhs],[+1.],convert(Float64,lhs))
 (-)(lhs::Number, rhs::Variable) = AffExpr([rhs],[-1.],convert(Float64,lhs))
 (*)(lhs::Number, rhs::Variable) = AffExpr([rhs],[convert(Float64,lhs)], 0.)
+# Number--Norm
+(+){Power,C,V}(lhs::Number, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs),  one(C), convert(C,lhs))
+(-){Power,C,V}(lhs::Number, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs), -one(C), convert(C,lhs))
+(*){Power,C,V}(lhs::Number, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs),     lhs, zero(GenericAffExpr{C,V}))
 # Number--GenericAffExpr
 (+)(lhs::Number, rhs::GenericAffExpr) = GenericAffExpr(copy(rhs.vars),copy(rhs.coeffs),lhs+rhs.constant)
 (-)(lhs::Number, rhs::GenericAffExpr) = GenericAffExpr(copy(rhs.vars),    -rhs.coeffs ,lhs-rhs.constant)
@@ -30,6 +34,10 @@
 (+)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2),copy(rhs.qcoeffs),lhs+rhs.aff)
 (-)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2),    -rhs.qcoeffs ,lhs-rhs.aff)
 (*)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2), lhs*rhs.qcoeffs ,lhs*rhs.aff)
+# Number--NormExpr
+(+){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm),     rhs.coeff, lhs+rhs.aff)
+(-){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm),     rhs.coeff, lhs-rhs.aff)
+(*){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm), lhs*rhs.coeff, lhs*rhs.aff)
 
 # This is not well defined if variable types are different, but needed to avoid ambiguities
 (+)(lhs::GenericAffExpr, rhs::GenericAffExpr) = (+)(promote(lhs,rhs)...)
@@ -48,6 +56,9 @@
 (+)(lhs::Variable, rhs::Variable) = AffExpr([lhs,rhs], [1.,+1.], 0.)
 (-)(lhs::Variable, rhs::Variable) = AffExpr([lhs,rhs], [1.,-1.], 0.)
 (*)(lhs::Variable, rhs::Variable) = QuadExpr([lhs],[rhs],[1.],AffExpr(Variable[],Float64[],0.))
+# Variable--Norm
+(+){Power,C,V}(lhs::Variable, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs),  one(C), GenericAffExpr{C,V}(lhs))
+(-){Power,C,V}(lhs::Variable, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs), -one(C), GenericAffExpr{C,V}(lhs))
 # Variable--AffExpr
 (+){CoefType,VarType}(lhs::VarType, rhs::GenericAffExpr{CoefType,VarType}) =
     GenericAffExpr{CoefType,VarType}(vcat(rhs.vars,lhs),vcat( rhs.coeffs,one(CoefType)), rhs.constant)
@@ -65,6 +76,32 @@ end
 # Variable--QuadExpr
 (+)(v::Variable, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),v+q.aff)
 (-)(v::Variable, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),    -q.qcoeffs ,v-q.aff)
+# Variable--NormExpr
+(+){T<:GenericNormExpr}(lhs::Variable, rhs::T) = T(copy(rhs.norm),  rhs.coeff, lhs+rhs.aff)
+(-){T<:GenericNormExpr}(lhs::Variable, rhs::T) = T(copy(rhs.norm), -rhs.coeff, lhs-rhs.aff)
+
+# Norm
+(+)(lhs::GenericNorm) = lhs
+(-){Power,C,V}(lhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(lhs), -one(C), zero(GenericAffExpr{C,V}))
+(*)(lhs::GenericNorm) = lhs
+# Norm--Number
+(+){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::Number) =
+    GenericNormExpr{Power,C,V}(copy(lhs), one(C), GenericAffExpr{C,V}( rhs))
+(-){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::Number) =
+    GenericNormExpr{Power,C,V}(copy(lhs), one(C), GenericAffExpr{C,V}(-rhs))
+(*){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::Number) =
+    GenericNormExpr{Power,C,V}(copy(lhs), copy(rhs), zero(GenericAffExpr{C,V}))
+(/){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::Number) =
+    GenericNormExpr{Power,C,V}(copy(lhs.terms), 1/rhs, zero(GenericAffExpr{C,V}))
+# Norm--Variable
+(+){Power}(lhs::Norm{Power},rhs::Variable) = NormExpr{Power}(copy(lhs), 1.0, AffExpr( rhs))
+(-){Power}(lhs::Norm{Power},rhs::Variable) = NormExpr{Power}(copy(lhs), 1.0, AffExpr(-rhs))
+# Norm--Norm
+# Norm--AffExpr
+(+){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::GenericAffExpr{C,V}) = GenericNormExpr{Power,C,V}(copy(lhs), one(C), copy(rhs))
+(-){Power,C,V}(lhs::GenericNorm{Power,C,V},rhs::GenericAffExpr{C,V}) = GenericNormExpr{Power,C,V}(copy(lhs), one(C),     -rhs)
+# Norm--QuadExpr
+# Norm--NormExpr
 
 # AffExpr (GenericAffExpr)
 (+)(lhs::GenericAffExpr) = lhs
@@ -85,6 +122,9 @@ end
 (-){CoefType,VarType}(lhs::GenericAffExpr{CoefType,VarType}, rhs::VarType) = GenericAffExpr{CoefType,VarType}(vcat(lhs.vars,rhs),vcat(lhs.coeffs,-one(CoefType)),lhs.constant)
 (*)(lhs::AffExpr, rhs::Variable) = (*)(rhs,lhs)
 (/)(lhs::AffExpr, rhs::Variable) = error("Cannot divide affine expression by a variable")
+# AffExpr--Norm
+(+){Power,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs),  one(C), copy(lhs))
+(-){Power,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNorm{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs), -one(C), copy(lhs))
 # AffExpr--AffExpr
 (+){T<:GenericAffExpr}(lhs::T, rhs::T) = GenericAffExpr(vcat(lhs.vars,rhs.vars),vcat(lhs.coeffs, rhs.coeffs),lhs.constant+rhs.constant)
 (-){T<:GenericAffExpr}(lhs::T, rhs::T) = GenericAffExpr(vcat(lhs.vars,rhs.vars),vcat(lhs.coeffs,-rhs.coeffs),lhs.constant-rhs.constant)
@@ -145,7 +185,9 @@ end
 # AffExpr--QuadExpr
 (+)(a::AffExpr, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),a+q.aff)
 (-)(a::AffExpr, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),    -q.qcoeffs ,a-q.aff)
-
+# AffExpr--NormExpr
+(+){Power,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNormExpr{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs.norm),  one(C), lhs+rhs.aff)
+(-){Power,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNormExpr{Power,C,V}) = GenericNormExpr{Power,C,V}(copy(rhs.norm), -one(C), lhs-rhs.aff)
 
 # QuadExpr
 (+)(lhs::QuadExpr) = lhs
@@ -172,6 +214,27 @@ end
 (-)(q1::QuadExpr, q2::QuadExpr) = QuadExpr( vcat(q1.qvars1, q2.qvars1),     vcat(q1.qvars2, q2.qvars2),
                                             vcat(q1.qcoeffs, -q2.qcoeffs),  q1.aff - q2.aff)
 
+# NormExpr
+(+)(lhs::GenericNormExpr) = lhs
+(-){T<:GenericNormExpr}(lhs::T) = T(copy(lhs.norm), -lhs.coeff, -lhs.aff)
+(*)(lhs::GenericNormExpr) = lhs
+# NormExpr--Number
+(+){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
+(-){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
+(*){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff*rhs, lhs.aff*rhs)
+(/){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff/rhs, lhs.aff/rhs)
+# NormExpr--Variable
+(+)(lhs::NormExpr,rhs::Variable) = NormExpr(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
+(-)(lhs::NormExpr,rhs::Variable) = NormExpr(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
+# NormExpr--Norm
+# NormExpr--AffExpr
+(+){Power,C,V}(lhs::GenericNormExpr{Power,C,V},rhs::GenericAffExpr{C,V}) =
+    GenericNormExpr{Power,C,V}(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
+(-){Power,C,V}(lhs::GenericNormExpr{Power,C,V},rhs::GenericAffExpr{C,V}) =
+    GenericNormExpr{Power,C,V}(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
+# NormExpr--QuadExpr
+# NormExpr--NormExpr
+
 # (==)(lhs::AffExpr,rhs::AffExpr) = (lhs.vars == rhs.vars) && (lhs.coeffs == rhs.coeffs) && (lhs.constant == rhs.constant)
 # (==)(lhs::QuadExpr,rhs::QuadExpr) = (lhs.qvars1 == rhs.qvars1) && (lhs.qvars2 == rhs.qvars2) && (lhs.qcoeffs == rhs.qcoeffs) && (lhs.aff == rhs.aff)
 
@@ -231,8 +294,6 @@ end
 #  - dot
 #############################################################################
 
-typealias JuMPTypes Union(Variable,AffExpr,QuadExpr)
-
 Base.sum(j::JuMPArray) = sum(j.innerArray)
 Base.sum(j::JuMPDict)  = sum(values(j.tupledict))
 Base.sum(j::JuMPArray{Variable}) = AffExpr(vec(j.innerArray), ones(length(j.innerArray)), 0.0)
@@ -286,8 +347,6 @@ end
 # A bunch of operator junk to make matrix multiplication and friends act
 # reasonably sane with JuMP types
 
-typealias JuMPTypes Union(Variable,AffExpr,QuadExpr)
-
 Base.promote_rule{R<:Real}(::Type{Variable},::Type{R}       ) = AffExpr
 Base.promote_rule(         ::Type{Variable},::Type{AffExpr} ) = AffExpr
 Base.promote_rule(         ::Type{Variable},::Type{QuadExpr}) = QuadExpr