# Primal-Dual Hybrid Gradient for L1 regularization

**URL:** https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201
**Category:** Optimization (Mathematical)
**Tags:** optimization
**Created:** [July 18, 2024, 9:05pm UTC](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201 "2024-07-18T21:05:35Z")
**Posts on this page:** 4
**Page:** 2

<div class="post-metadata">

### Author: ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)
#### Post date: [July 21, 2024, 1:57pm UTC](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201/21 "2024-07-21T13:57:08Z")

</div>

Replace

```julia
G1(f) = (Kf - b)'*(Kf-b) + alpha * sum(abs.(f))

```

with

```julia
G2(f) = (Kf - b)'*(Kf-b) + alpha * sL1(f,1.e-7)

```

and send G2 to your favorite smooth optimization package. Most of them can handle the nonnegativity constraint with ease.

I have a recent paper that uses this idea for a very similar problem, but it may be more mathematical than you need.

> **[Min-Max Optimization for Robust Nonlinear Least Squares Problems](https://arxiv.org/abs/2402.12679)**
>
> This paper considers min-max optimization for a class of robust nonlinear least squares problems. We show via an example that solving the first order optimality conditions defined by gradients of the objective function can lead to incorrect solutions...

---

<div class="post-metadata">

### Author: ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)
#### Post date: [July 21, 2024, 4:47pm UTC](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201/22 "2024-07-21T16:47:50Z")

</div>

Understood, here is an implementation:

> **Smooth L1 norm using OptimizationOptimJL.BFGS**
>
> ```julia
> using Optimization, OptimizationOptimJL, LinearAlgebra
> 
> K = rand(10, 10)
> g = rand(10)
> α = 1
> 
> function sL1(x::AbstractVector{T}, mu=1.e-7) where {T}
> snrm = sum(Cabs.(x, mu))
> return snrm
> end
> 
> function Cabs(x, mu)
> p = 4.0 * mu * mu
> Cabs = sqrt(x * x + p)
> return Cabs
> end
> 
> function G(f, p)
> 
> # u = f, p = [K,s,α]
> k = p[1]
> s = p[2]
> α = p[3]
> 
> r = K * f - s
> return r' * r + α * sL1(f, 1.e-7)
> 
> end
> 
> optf = Optimization.OptimizationFunction(G, Optimization.AutoForwardDiff())
> prob = Optimization.OptimizationProblem(optf, ones(size(K, 2)),
> (K, g, α),
> lb=zeros(size(K, 2)), ub=[Inf for i in 1:size(K, 2)])
> f = OptimizationOptimJL.solve(prob, OptimizationOptimJL.BFGS())
> 
> ```

It does converge to something, but it’s pretty slow compared to the other methods. Perhaps the way I’ve written it is far from optimal.  
I assume we can remove the r allocations from the cost function, but I’m not sure what would be the best way to do it.

---

<div class="post-metadata">

### Author: ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)
#### Post date: [July 21, 2024, 5:22pm UTC](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201/23 "2024-07-21T17:22:21Z")

</div>

It seems you have precomputed K and b, so the allocations are not coming from that.

If something else is doing what you want and is fast enough, you should probably use that. You can invest a lot of time in the Optimization.jl docs before you figure this problem out.

---

<div class="post-metadata">

### Author: ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)
#### Post date: [July 21, 2024, 5:34pm UTC](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201/24 "2024-07-21T17:34:33Z")

</div>

> [@aris](#):
>
> > [@kellertuer](#):
> >
> > You can of course also just be lucky that f \geq 0f≥0f \geq 0 is in your cases automatically fulfilled.
> 
> No luck unfortunately! I’ll try to dig into the documentation of ProximalAlgorithms a bit more to see if there’s an obvious way to do it.

I am not aware of algorithms that do constrained PDHG, so you would have to check with a cost that is pus infinity for negative input and check what it’s pros (and grad) are in the PDHG setting are; one rephrasing could be to check out to work on the set of positive numbers (as a manifold), but I am super biased there, since the Chambolle-Pock I linked can (among others) do exactly that (and I implemented that). You would still need special probes and gradients then.

[Previous page](https://discourse.julialang.org/t/primal-dual-hybrid-gradient-for-l1-regularization/117201.md?page=1)
