ProximalOperators.jl/src/functions/elasticNet.jl at 4e0f0054a56b9f48379e538687777aecba678fd3 · hakkelt/ProximalOperators.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
# elastic-net regularization

export ElasticNet

"""
    ElasticNet(μ=1, λ=1)

Return the function
```math
f(x) = μ\\|x\\|_1 + (λ/2)\\|x\\|^2,
```
for nonnegative parameters `μ` and `λ`.
"""
struct ElasticNet{R, S}
    mu::R
    lambda::S
    function ElasticNet{R, S}(mu::R, lambda::S) where {R, S}
        if lambda < 0 || mu < 0
            error("parameters `μ` and `λ` must be nonnegative")
        else
            new(mu, lambda)
        end
    end
end

is_separable(f::Type{<:ElasticNet}) = true
is_proximable(f::Type{<:ElasticNet}) = true
is_convex(f::Type{<:ElasticNet}) = true

ElasticNet(mu::R=1, lambda::S=1) where {R, S} = ElasticNet{R, S}(mu, lambda)

function (f::ElasticNet)(x)
    R = real(eltype(x))
    return f.mu * norm(x, 1) + f.lambda / R(2) * norm(x, 2)^2
end

function prox!(y, f::ElasticNet, x, gamma)
    R = real(eltype(x))
    sqnorm2x = R(0)
    norm1x = R(0)
    gm = gamma * f.mu
    gl = gamma * f.lambda
    for i in eachindex(x)
        y[i] = (x[i] + (x[i] <= -gm ? gm : (x[i] >= gm ? -gm : -x[i])))/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y, f::ElasticNet, x, gamma::AbstractArray)
    R = real(eltype(x))
    sqnorm2x = R(0)
    norm1x = R(0)
    for i in eachindex(x)
        gm = gamma[i] * f.mu
        gl = gamma[i] * f.lambda
        y[i] = (x[i] + (x[i] <= -gm ? gm : (x[i] >= gm ? -gm : -x[i])))/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y, f::ElasticNet, x::AbstractArray{<:Complex}, gamma)
    R = real(eltype(x))
    sqnorm2x = R(0)
    norm1x = R(0)
    gm = gamma * f.mu
    gl = gamma * f.lambda
    for i in eachindex(x)
        y[i] = sign(x[i]) * max(0, abs(x[i]) - gm)/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y, f::ElasticNet, x::AbstractArray{<:Complex}, gamma::AbstractArray)
    R = real(eltype(x))
    sqnorm2x = R(0)
    norm1x = R(0)
    for i in eachindex(x)
        gm = gamma[i] * f.mu
        gl = gamma[i] * f.lambda
        y[i] = sign(x[i]) * max(0, abs(x[i]) - gm)/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function gradient!(y, f::ElasticNet, x)
    R = real(eltype(x))
    # Gradient of 1 norm
    y .= f.mu .* sign.(x)
    # Gradient of 2 norm
    y .+= f.lambda .* x
    return f.mu * norm(x, 1) + f.lambda / R(2) * norm(x, 2)^2
end

function prox_naive(f::ElasticNet, x, gamma)
    R = real(eltype(x))
    uz = max.(0, abs.(x) .- gamma .* f.mu)./(1 .+ f.lambda .* gamma)
    return sign.(x) .* uz, f.mu * norm(uz, 1) + f.lambda / R(2) * norm(uz)^2
end