Skip to content

Added Unit Simplex #150

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Added Unit Simplex #150

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 47 additions & 1 deletion src/functions/indSimplex.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# indicator of a simplex

export IndSimplex
export IndSimplex, IndUnitSimplex

"""
IndSimplex(a=1.0)
Expand Down Expand Up @@ -107,3 +107,49 @@ function prox_naive(f::IndSimplex, x, gamma)
end
return v, R(0)
end

"""
IndUnitSimplex(a=1.0)

Return the indicator of the unit simplex
```math
S = \\left\\{ x : x \\geq 0, \\sum_i x_i \\leq a \\right\\}.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm yeah, naming is tough here. Would this set be the convex hull of S u {0}, where S is that of IndSimplex? (i.e. with equality on the sum)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, the origin 0 \in R^k is the projection to R^k of e_{k+1} \in R^{k+1}. So I guess this set is the projection to R^k of the k-simplex (which lives in R^{k+1} originally). If so, this could suggest a better name maybe.

Am I making sense? Works in my mind for k=2 at least (hehe).

```

By default `a=1`, therefore ``S`` is the probability simplex of dimension n+1.
"""
struct IndUnitSimplex{R}
a::R
function IndUnitSimplex{R}(a::R) where R
if a <= 0
error("parameter a must be positive")
else
new(a)
end
end
end

is_convex(f::Type{<:IndUnitSimplex}) = true
is_set(f::Type{<:IndUnitSimplex}) = true

IndUnitSimplex(a::R=1) where R = IndUnitSimplex{R}(a)

function (f::IndUnitSimplex)(x)
R = eltype(x)
if all(x .>= 0) && sum(x) <= f.a + eps(f.a)
return R(0)
end
return R(Inf)
end

function prox!(y, f::IndUnitSimplex{R}, x, gamma) where {R}
fx = zero(R)
for i in eachindex(x)
y[i] = max(x[i], zero(R))
fx += y[i]
end
if fx > f.a
simplex_proj_condat!(y, f.a, x)
end
return eltype(x)(0)
end
2 changes: 2 additions & 0 deletions test/test_calls.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ test_cases_spec = [
"right" => [
( (IndSimplex(),), randn(Float32, 10) ),
( (IndSimplex(),), randn(Float64, 10) ),
( (IndUnitSimplex(),), randn(Float32, 10) ),
( (IndUnitSimplex(),), randn(Float64, 10) ),
( (IndNonnegative(), rand()), randn(Float64, 10) ),
( (IndZero(),), randn(Float64, 10) ),
( (IndBox(-1, 1),), randn(Float32, 10) ),
Expand Down
14 changes: 12 additions & 2 deletions test/test_equivalences.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,22 @@ for i = 1:N
r = 5*rand()
f = IndSimplex(r)
g = IndBallL1(r)
h = IndUnitSimplex(r)

y1, fy1 = prox(f, abs.(x))
y1 = sign.(x).*y1
y1_l1ball = sign.(x).*y1
y2, gy2 = prox(g, x)
y3, hy3 = prox(h, abs.(x))

@test y1 ≈ y2
@test y1_l1ball ≈ y2
@test y1 ≈ y3

x2 = abs.(x) * 0.5 * r
x2 ./= sum(x2)

y_probsimplex = prox(f, x2)
y_unit_simplex = prox(h, x2)
@test norm(y_probsimplex) >= norm(y_unit_simplex)
end

# projecting onto the simplex
Expand Down
Loading