IndicatorConvexCone

[mfalt]

You seem to have defined IndicatorConvexCone <: IndicatorConvex but most cones are defined as immutable IndPSD <: IndicatorConvex end immutable IndFree <: IndicatorConvex end Is there a reason for this? I would be happy add a PR if not. It would also be nice for the package I'm writing (and for slight efficiency) if we had IndNonnegative (and IndNonpositive) as IndicatorConvexCone instead of IndNonnegative() = IndBox(0.0, +Inf). That way I would be able to restrict my algorithm to work with any IndicatorConvexCone in this package.

We might also consider adding a IndDualCone{T<:IndicatorConvexCone}. I already have an implementation for this in my code, and it should be easy to include here. It looks something like this:

type DualCone{T<:IndicatorConvex} <: IndicatorConvex
    C::T
end

dual{T<:IndicatorConvex}(C::T) = DualCone{T}(C)
dual(C::DualCone) = C.C

function prox!(cone::DualCone, x, y)
    prox!(cone.C, -x, y)
    for i = 1:length(x)
        y[i] = x[i] + y[i]
    end
end

and optionally we could add something like this

#Wrapper for dual prox to avoid double duals
function proxDual!(C::IndicatorConvex, x, y)
    prox!(C, -x, y)
    @simd for i = 1:length(x)
        y[i] = x[i] + y[i]
    end
end
proxDual!(C::DualCone, x, y) = prox!(C.C, x, y)

I'm happy to create a PR for this too if you are interested.

[lostella]

Regaring your first point: no, there's no reason for that. I guess that's just a mistake and IndPSD and IndFree should in fact be of type IndicatorConvexCone.

The same applies to IndNonpositive and IndNonnegative. I just made them as they are so as to have a shortcut for their implementation. It makes perfect sense (also for performance) to specialize their implementation. One could keep track of their correctness by comparing them with their IndBox-based implementation in tests.

As for your last point (dual cone), check out conjugate.jl, which already implements Moreau identity for computing the prox! of conjugate functions.

[lostella]

Also, keep in mind that we've gone for correctness first, leaving some optimizations aside. There's many @inbounds and @simd and so on that should be added here and there.

[lostella]

What could make sense would be to define

immutable Conjugate{T <: ProximableConvex} <: ProximableConvexConjugate
  f::T
end

(note that the type of T is different from what it currently is: I guess conjugacy mostly makes sense when the function is convex?) and then have your trick to avoid double conjugation (which yes, only makes sense in the case of closed, convex functions, in which case f** = f).

[mfalt]

Great, I will create a PR. I hadn't seen the conjugate file, nice work. I think we should think over how we handle the types though, ít's nice that we allow T<:Real but we should be careful about things like function prox!{R <: Real}(g::Conjugate, x::AbstractArray{R}, y::AbstractArray{R}, gamma::Real=1.0) for example

function g{T<:Real}(a::T, b::Real=1.0)
       return a*b
end

g(Float32(2))
#Gives 1.0 <: Float64

function g2{T<:Real}(a::T, b::Real=one(T))
       return a*b
end

g2(Float32(2))
#Gives 2.0f0 <: Float32

The same thing works for zero(T) but I'm not sure if there is a nice solution for Inf.

Did you see my comment on the issue you closed by the way?

[lostella]

Of course! I sent you an email about one hour ago but maybe I found the wrong address. Do you use skype? You can add me, ID: lorestella.

[mfalt]

I don't see why it should not be possible to have functions for both convex and non-convex functions? What about the following setup, where we don't have a specific Conjugate type (haven't thought it through completely)

immutable ConjugateFunction{T <: ProximableFunction} <: ProximableConvex
  f::T
end

Conjugate{T<:ProximableFunction}(x::T) = ConjugateFunction{T}(x)
Conjugate{T<:ProximableConvex}(x::ConjugateFunction{T}) = x.f

which would allow for retaining properties under conjugates (for example Cone), without having different names for the conjugate, for example

immutable PolarCone{T <: IndicatorConvexCone} <: IndicatorConvexCone
  f::T
end
Conjugate{T<:IndicatorConvexCone}(x::T) = PolarCone{T}(x)
Conjugate{T<:IndicatorConvexCone}(x::ConjugateCone{T}) = x.f

Which should work fine with something like

ConjugateType = Union{ConjugateFunction,PolarCone}

prox!{R <: Real}(g::ConjugateType, x::AbstractArray{R}, y::AbstractArray{R}) = ...

as well as the special case, including evaluating

prox!{R <: Real}(g::PolarCone, x::AbstractArray{R}, y::AbstractArray{R}) ) = ...

function (f::PolarCone){T <: Union{Real, Complex}}(x::AbstractArray{T})
  #We can evaluate this
end

[mfalt]

I think we could even allow ?

immutable PolarCone{T <: IndicatorConvex} <: IndicatorConvexCone
  f::T
end
Conjugate{T<:IndicatorConvex}(x::T) = PolarCone{T}(x)
Conjugate{T<:IndicatorConvexCone}(x::ConjugateCone{T}) = x.f

[lostella]

Well, it makes perfect sense to consider the conjugate of nonconvex functions. What doesn't work in general, though, is Moreau's decomposition for computing prox!: that works everywhere only for convex functions, see for example Bauschke, Combettes, Thm 14.3. I'm not aware of extensions of that identity for nonconvex functions (maybe not being valid everywhere), but they might be out there somewhere.

[mfalt]

Yeah, you are right, I wasn't thinking properly. However, I still think there is some merit to allowing the dual of a cone remain a cone. But we might consider restricting duals to be of a convex function.

I guess some conjugates can be created immediately too, like the IndFree -> IndZero, IndNonnegative -> IndNonpositive, and NormL1 -> IndBallLinf without wrapping them. These are obviously minor improvements that can be added later.

[lostella]

Anyway, yes, I like that structure.

As for your last example, though, the conjugate of a general indicator function of a convex set wouldn't be a cone necessarily.

[mfalt]

Are you sure? I'll look up the necessary conditions, or reference tomorrow.

[lostella]

Yes, I'm sure about that. In general, the conjugate of ind_A, the indicator of set A, is sup_A, the support function of A: this is not even an indicator function in general. For example, the conjugate of the indicator of an infinity-norm ball, is a (weighted) 1-norm.

[mfalt]

You are of course right again, I even wrote NormL1 -> IndBallLinf in the previous comment, I was confusing it with the polar/dual cone of a convex set.

[mfalt]

I think we might want to rethink some parts here when it comes to dual cones. For example IndExpDual() = Precompose(Conjugate(IndExpPrimal()), -1.0) I see no clear way of making this return something of type IndicatorConvexCone. The only way I see would be to have both Precompose and Conjugate be type aliases as in ConjugateType = Union{ConjugateFunction,PolarCone} but this seems like it could result in an explosion in the number of types if we want to retain properties under these operations. And when it comes to other types like precompose and postcompose, I'm not sure that properties like Cone and Convex is always retained, independent of the additional arguments. Maybe we should consider a way of verifying these preperties instead of coding them into the types? Like iscone, isaffine and so on? I imagine that this might result in some overhead though. This is more general than just cones, for example Precompose{ProximableConvex} is not ProximableConvex, Precompose{IndicatorConvex} is not IndicatorConvex and so on. The specific dual cones could however still be solved by

type DualCone{T<:IndicatorConvexCone} <: IndicatorConvexCone
    C::T
end

[lostella]

Yes, I perfectly understand your point. I'm open to suggestions on this, I had already realized that relying too much on the type system for distinguishing between "types" of functions was a bit restrictive and didn't mix well, for example, with calculus rules as you said.

[lostella]

@mfalt: the thing you said about Real types, I think we should wait until 0.6 where we'll be able to do nice things like

MaybeComplex{R} = Union{R, Complex{R}} # this is 0.6 syntax

function prox!{R <: Real, C <: MaybeComplex{R}}(y::AbstractArray{C}, f::SomeProximableType, x::AbstractArray{C}, gamma::R)
  ...
end

I mean, I think it makes sense to rethink the types of arguments in this way: fix a Real type, have all arguments based on that Real type. Currently I believe we can't do something like the snippet above, so maybe we should wait until we decide to drop 0.5 (it doesn't really make sense to support all these pre-1.0 versions).

[lostella]

Ok, I'm closing this and opening a new one with the issues that came out of the discussion!