Skip to content

Commit

Permalink
Tidy src/atoms/second_order_cone
Browse files Browse the repository at this point in the history
  • Loading branch information
odow committed Jan 10, 2024
1 parent 8dfc7f6 commit 42181be
Show file tree
Hide file tree
Showing 8 changed files with 74 additions and 129 deletions.
13 changes: 4 additions & 9 deletions src/atoms/second_order_cone/geomean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,28 +9,23 @@ mutable struct GeoMeanAtom <: AbstractExpr
error("geo mean must take arguments of the same size")
elseif any(x -> sign(x) == ComplexSign(), args)
error("The arguments should be real, not complex")
else
children = args
return new(children, sz)
end
return new(args, sz)
end
end

head(io::IO, ::GeoMeanAtom) = print(io, "geomean")

function Base.sign(q::GeoMeanAtom)
return Positive()
end
Base.sign(::GeoMeanAtom) = Positive()

function monotonicity(q::GeoMeanAtom)
return fill(Nondecreasing(), length(q.children))
end

function curvature(q::GeoMeanAtom)
return ConcaveVexity()
end
curvature(::GeoMeanAtom) = ConcaveVexity()

_geomean(scalar_args...) = prod(scalar_args)^(1 / length(scalar_args))

function evaluate(q::GeoMeanAtom)
n = length(q.children)
children = evaluate.(q.children)
Expand Down
23 changes: 7 additions & 16 deletions src/atoms/second_order_cone/huber.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,28 +5,21 @@ mutable struct HuberAtom <: AbstractExpr

function HuberAtom(x::AbstractExpr, M::Real)
if sign(x) == ComplexSign()
error("Arguemt must be real")
error("Argument must be real")

Check warning on line 8 in src/atoms/second_order_cone/huber.jl

View check run for this annotation

Codecov / codecov/patch

src/atoms/second_order_cone/huber.jl#L8

Added line #L8 was not covered by tests
elseif M <= 0
error("Huber parameter must by a positive scalar")
end
children = (x,)
return new(children, x.size, M)
return new((x,), x.size, M)
end
end

head(io::IO, ::HuberAtom) = print(io, "huber")

function Base.sign(x::HuberAtom)
return Positive()
end
Base.sign(::HuberAtom) = Positive()

function monotonicity(x::HuberAtom)
return (Nondecreasing() * sign(x.children[1]),)
end
monotonicity(x::HuberAtom) = (Nondecreasing() * sign(x.children[1]),)

function curvature(x::HuberAtom)
return ConvexVexity()
end
curvature(::HuberAtom) = ConvexVexity()

function evaluate(x::HuberAtom)
c = evaluate(x.children[1])
Expand All @@ -44,11 +37,9 @@ function new_conic_form!(context::Context, x::HuberAtom)
c = x.children[1]
s = Variable(c.size)
n = Variable(c.size)

# objective given by s.^2 + 2 * M * |n|
objective = conic_form!(context, square(s) + 2 * x.M * abs(n))
add_constraint!(context, c == s + n)
return objective
# objective given by s.^2 + 2 * M * |n|
return conic_form!(context, square(s) + 2 * x.M * abs(n))
end

huber(x::AbstractExpr, M::Real = 1.0) = HuberAtom(x, M)
8 changes: 5 additions & 3 deletions src/atoms/second_order_cone/norm.jl
Original file line number Diff line number Diff line change
@@ -1,20 +1,22 @@
# deprecate these soon
norm_inf(x::AbstractExpr) = maximum(abs(x))

norm_1(x::AbstractExpr) = sum(abs(x))

norm_fro(x::AbstractExpr) = LinearAlgebra.norm2(vec(x))

"""
norm(x::AbstractExpr, p::Real=2)
Computes the `p`-norm `‖x‖ₚ = (∑ᵢ |xᵢ|^p)^(1/p)` of a vector expression `x`. Matrices
are vectorized (i.e. `norm(x)` is the same as `norm(vec(x))`.)
Computes the `p`-norm `‖x‖ₚ = (∑ᵢ |xᵢ|^p)^(1/p)` of a vector expression `x`.
Matrices are vectorized (i.e., `norm(x)` is the same as `norm(vec(x))`.)
This function uses specialized methods for `p=1, 2, Inf`. For `p > 1` otherwise,
this function uses the procedure documented at
[`rational_to_socp.pdf`](https://github.com/jump-dev/Convex.jl/raw/master/docs/supplementary/rational_to_socp.pdf),
based on the paper "Second-order cone programming" by F. Alizadeh and D. Goldfarb,
Mathematical Programming, Series B, 95:3-51, 2001.
"""
function LinearAlgebra.norm(x::AbstractExpr, p::Real = 2)
if size(x, 2) > 1
Expand Down
49 changes: 11 additions & 38 deletions src/atoms/second_order_cone/norm2.jl
Original file line number Diff line number Diff line change
@@ -1,59 +1,32 @@
#############################################################################
# LinearAlgebra.norm2.jl
# Handles the euclidean norm (also called frobenius norm for matrices)
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

mutable struct EucNormAtom <: AbstractExpr
children::Tuple{AbstractExpr}
size::Tuple{Int,Int}

function EucNormAtom(x::AbstractExpr)
children = (x,)
return new(children, (1, 1))
end
EucNormAtom(x::AbstractExpr) = new((x,), (1, 1))
end

head(io::IO, ::EucNormAtom) = print(io, "norm2")

function Base.sign(x::EucNormAtom)
return Positive()
end
Base.sign(::EucNormAtom) = Positive()

function monotonicity(x::EucNormAtom)
return (sign(x.children[1]) * Nondecreasing(),)
end
monotonicity(x::EucNormAtom) = (sign(x.children[1]) * Nondecreasing(),)

function curvature(x::EucNormAtom)
return ConvexVexity()
end
curvature(::EucNormAtom) = ConvexVexity()

function evaluate(x::EucNormAtom)
return norm(vec(evaluate(x.children[1])))
end
evaluate(x::EucNormAtom) = norm(vec(evaluate(x.children[1])))

## Create a new variable euc_norm to represent the norm
## Additionally, create the second order conic constraint (euc_norm, x) in SOC
function new_conic_form!(context::Context{T}, A::EucNormAtom) where {T}
obj = conic_form!(context, only(AbstractTrees.children(A)))

x = only(AbstractTrees.children(A))
d = length(x)

t_obj = conic_form!(context, Variable())

f = operate(vcat, T, sign(A), t_obj, obj)
set = MOI.SecondOrderCone(d + 1)
MOI_add_constraint(context.model, f, set)

return t_obj
t = conic_form!(context, Variable())
obj = conic_form!(context, x)
f = operate(vcat, T, sign(A), t, obj)
MOI_add_constraint(context.model, f, MOI.SecondOrderCone(length(x) + 1))
return t
end

function LinearAlgebra.norm2(x::AbstractExpr)
if sign(x) == ComplexSign()
return EucNormAtom([real(x); imag(x)])
else
return EucNormAtom(x)
end
return EucNormAtom(x)
end
29 changes: 13 additions & 16 deletions src/atoms/second_order_cone/qol_elementwise.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,24 +8,19 @@ mutable struct QolElemAtom <: AbstractExpr
"elementwise quad over lin must take two arguments of the same size",
)
end
children = (x, y)
return new(children, x.size)
return new((x, y), x.size)
end
end

head(io::IO, ::QolElemAtom) = print(io, "qol_elem")

function Base.sign(q::QolElemAtom)
return Positive()
end
Base.sign(::QolElemAtom) = Positive()

function monotonicity(q::QolElemAtom)
return (sign(q.children[1]) * Nondecreasing(), Nonincreasing())
end

function curvature(q::QolElemAtom)
return ConvexVexity()
end
curvature(::QolElemAtom) = ConvexVexity()

function evaluate(q::QolElemAtom)
return (evaluate(q.children[1]) .^ 2) ./ evaluate(q.children[2])
Expand All @@ -36,7 +31,6 @@ function new_conic_form!(context::Context{T}, q::QolElemAtom) where {T}
t = Variable(sz[1], sz[2])
t_obj = conic_form!(context, t)
x, y = q.children

for i in 1:length(q.children[1])
add_constraint!(
context,
Expand All @@ -50,11 +44,12 @@ end
qol_elementwise(x::AbstractExpr, y::AbstractExpr) = QolElemAtom(x, y)

function Base.Broadcast.broadcasted(::typeof(^), x::AbstractExpr, k::Int)
return k == 2 ? QolElemAtom(x, constant(ones(x.size[1], x.size[2]))) :
error("raising variables to powers other than 2 is not implemented")
if k != 2
error("raising variables to powers other than 2 is not implemented")

Check warning on line 48 in src/atoms/second_order_cone/qol_elementwise.jl

View check run for this annotation

Codecov / codecov/patch

src/atoms/second_order_cone/qol_elementwise.jl#L48

Added line #L48 was not covered by tests
end
return QolElemAtom(x, constant(ones(x.size[1], x.size[2])))
end

# handle literal case
function Base.Broadcast.broadcasted(
::typeof(Base.literal_pow),
::typeof(^),
Expand All @@ -65,13 +60,16 @@ function Base.Broadcast.broadcasted(
end

invpos(x::AbstractExpr) = QolElemAtom(constant(ones(x.size[1], x.size[2])), x)

function Base.Broadcast.broadcasted(::typeof(/), x::Value, y::AbstractExpr)
return dotmultiply(constant(x), invpos(y))
end

function Base.:/(x::Value, y::AbstractExpr)
return size(y) == (1, 1) ? MultiplyAtom(constant(x), invpos(y)) :
error("cannot divide by a variable of size $(size(y))")
if size(y) != (1, 1)
error("cannot divide by a variable of size $(size(y))")

Check warning on line 70 in src/atoms/second_order_cone/qol_elementwise.jl

View check run for this annotation

Codecov / codecov/patch

src/atoms/second_order_cone/qol_elementwise.jl#L70

Added line #L70 was not covered by tests
end
return MultiplyAtom(constant(x), invpos(y))
end

sumsquares(x::AbstractExpr) = square(norm2(x))
Expand All @@ -81,7 +79,6 @@ function square(x::AbstractExpr)
error(
"Square of complex number is not DCP. Did you mean square_modulus?",
)
else
QolElemAtom(x, constant(ones(x.size[1], x.size[2])))
end
return QolElemAtom(x, constant(ones(x.size[1], x.size[2])))
end
43 changes: 22 additions & 21 deletions src/atoms/second_order_cone/quadform.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@ end
"""
is_psd(A; tol)
Check whether `A` is positive semi-definite by computing a LDLᵀ factorization of `A + tol*I`
Check whether `A` is positive semi-definite by computing a LDLᵀ factorization of
`A + tol*I`
"""
function is_psd(A; tol = sqrt(eps(float(real(eltype(A))))))
T = eltype(A)
Expand All @@ -14,16 +15,17 @@ function is_psd(A; tol = sqrt(eps(float(real(eltype(A))))))
# * dense fallack otherwise
if T <: AbstractFloat || T <: LinearAlgebra.BlasFloat
return is_psd(SparseArrays.sparse(A); tol = tol)
else
return is_psd(Matrix(A); tol = tol)
end
return is_psd(Matrix(A); tol = tol)
end

function is_psd(
A::SparseArrays.SparseMatrixCSC{Complex{T}};
tol::T = sqrt(eps(T)),
) where {T<:LinearAlgebra.BlasReal}
return LinearAlgebra.isposdef(A + tol * LinearAlgebra.I)
end

function is_psd(
A::SparseArrays.SparseMatrixCSC{T};
tol::T = sqrt(eps(T)),
Expand All @@ -39,33 +41,30 @@ function is_psd(
return minimum(d) >= 0
catch err
# If the matrix could not be factorized, then it is not PSD
isa(err, LDLFactorizations.SQDException) && return false
rethrow() # Something else happened
if err isa LDLFactorizations.SQDException
return false

Check warning on line 45 in src/atoms/second_order_cone/quadform.jl

View check run for this annotation

Codecov / codecov/patch

src/atoms/second_order_cone/quadform.jl#L44-L45

Added lines #L44 - L45 were not covered by tests
end
rethrow(err)

Check warning on line 47 in src/atoms/second_order_cone/quadform.jl

View check run for this annotation

Codecov / codecov/patch

src/atoms/second_order_cone/quadform.jl#L47

Added line #L47 was not covered by tests
end
end

function is_psd(A::Matrix; tol = sqrt(eps(float(real(eltype(A))))))
return LinearAlgebra.isposdef(A + tol * LinearAlgebra.I)
end

function quadform(x::AbstractExpr, A::Value; assume_psd = false)
if length(size(A)) != 2 || size(A, 1) != size(A, 2)
error("Quadratic form only takes square matrices")
end
if !LinearAlgebra.ishermitian(A)
elseif !LinearAlgebra.ishermitian(A)
error("Quadratic form only defined for Hermitian matrices")
end
if assume_psd
factor = 1
factor = if assume_psd || is_psd(A)
1
elseif is_psd(-A)
-1
else
if is_psd(A)
factor = 1
elseif is_psd(-A)
factor = -1
else
error("Quadratic forms supported only for semidefinite matrices")
end
error("Quadratic forms supported only for semidefinite matrices")
end

P = sqrt(LinearAlgebra.Hermitian(factor * A))
return factor * square(norm2(P * x))
end
Expand All @@ -75,10 +74,12 @@ end
Represents `x' * A * x` where either:
* `x` is a vector-valued variable and `A` is a positive semidefinite or negative semidefinite matrix (and in particular Hermitian or real symmetric). If `assume_psd=true`, then
`A` will be assumed to be positive semidefinite. Otherwise, `Convex.is_psd` will be used to check if `A` is positive semidefinite or negative semidefinite.
* or `A` is a matrix-valued variable and `x` is a vector.
* `x` is a vector-valued variable and `A` is a positive semidefinite or
negative semidefinite matrix (and in particular Hermitian or real symmetric).
If `assume_psd=true`, then `A` will be assumed to be positive semidefinite.
Otherwise, `Convex.is_psd` will be used to check if `A` is positive
semidefinite or negative semidefinite.
* or `A` is a matrix-valued variable and `x` is a vector.
"""
function quadform(x::AbstractExpr, A::AbstractExpr; assume_psd = false)
if vexity(x) == ConstVexity()
Expand Down
11 changes: 3 additions & 8 deletions src/atoms/second_order_cone/quadoverlin.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,24 +6,19 @@ mutable struct QuadOverLinAtom <: AbstractExpr
if x.size[2] != 1 && y.size != (1, 1)
error("quad over lin arguments must be a vector and a scalar")
end
children = (x, y)
return new(children, (1, 1))
return new((x, y), (1, 1))
end
end

head(io::IO, ::QuadOverLinAtom) = print(io, "qol")

function Base.sign(q::QuadOverLinAtom)
return Positive()
end
Base.sign(::QuadOverLinAtom) = Positive()

function monotonicity(q::QuadOverLinAtom)
return (sign(q.children[1]) * Nondecreasing(), Nonincreasing())
end

function curvature(q::QuadOverLinAtom)
return ConvexVexity()
end
curvature(::QuadOverLinAtom) = ConvexVexity()

function evaluate(q::QuadOverLinAtom)
x = evaluate(q.children[1])
Expand Down
Loading

0 comments on commit 42181be

Please sign in to comment.