add a user define option for threshold in bregman #3

ziyiyin97 · 2022-04-04T17:43:15Z

No description provided.

mloubout · 2022-04-04T17:52:57Z

src/bregman.jl

@@ -124,7 +126,7 @@ function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams,
        # Update z variable
        @. z = z + d
        # Get λ at first iteration
-        i == 1 && (λ = abs(T(quantile(abs.(z), options.quantile))))
+        i == 1 && (sol.λ = λ = isnothing(options.lambda) ? abs(T(quantile(abs.(z), options.quantile))) : abs(T(options.lambda)))


That's a lot of if in one line but I guess.

I'm still waiting on an answer as to why you would ever end up with a complex-valued x

So that will be easy if I want to solve AC'x=b where x is the complex curvelet coefficients. It is easier for me to take this form if I want to do weightings on x later

This is not the right way to use this as you are not solving the correct problem. You are always solving the L1-l2 elastic net on the curvelet coefficient so x is your image and is real. if you give z as the input and C' instead of C you are solving sparsity in the image domain.

If you want to do weighting you need to provide W*C instead of C

OK I see your point. However, I think they are different ways to form the optimization problem. Either way works because they should be mathematically equivalent to solve for either the image or its curvelet coefficients. So I still think a PR to fix the support of complex number is needed. Any thought?

Either way does not work. You are confusing mathematical equivalence and implementation equivalence you can't make a program magically understand a change of variable (which is what you are doing). The documentation fo this implementation is very explicit about what it solves.

I don't have an issue with complex number supports but you are not doing what you think you are doing and it has to be for the case where x is complex, the actual primal x.

The problem of provide W*C instead of C is that its adjoint is not inverse. Then during soft-thresholding we have to do a W inverse.

That's because you are trying to do a weighted l1 which is a different problem. You could add weighted l1 with the standard identity as default but again you need to be careful about what problem you are solving.

OK I agree what you are saying now but we don't have to add anything new: just solve AC'Wx = b and use the identity as the sparsity transform will do the work. The solution is C'WX in the end. That's much simpler than adding a weighted l1

That's not the same as you will computethe threshold on W*z not on z.

No we don't have to. The forward and adjoint of W can already bump the entries on support 2 times so these entries are much easier to pass the threshold

mloubout · 2022-04-07T17:10:08Z

src/bregman.jl

@@ -124,15 +127,29 @@ function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams,
        # Update z variable
        @. z = z + d
        # Get λ at first iteration
-        i == 1 && (λ = abs(T(quantile(abs.(z), options.quantile))))
+        if i == 1


This is very ugly, it takes more space than the actual algorithm, make a proper function with dispatch for this.

mloubout · 2022-04-07T17:10:42Z

src/bregman.jl

+            if length(λ) == 1
+                @printf("%10d %15.5e %15.5e %15.5e %15.5e \n",i, t, obj_fun, f, λ)
+            else
+                @printf("%10d %15.5e %15.5e %15.5e %5s \n",i, t, obj_fun, f, "vector")


that's not informative, choose a value like min/max/mean/ or tuple or something to give proper information

mloubout · 2022-04-17T00:14:18Z

examples/denoising.jl

@@ -20,11 +20,11 @@ imgn= img .+ .01f0*randn(Float32, size(img))
 b = A*vec(imgn)

 # setup bregamn
-opt = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true)
-opt2 = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true, spg=true)
+opt = bregman_options(maxIter=200, verbose=2, alpha=1, antichatter=true)


I really don't like this. Now there is options all over the place that completely defeat the purpose. And why is the TD moved

I remove TD so that it's optional (default is identity matrix)

What do you mean "Now there is options all over the place that completely defeat the purpose"?

like do you mean I should move these lambda, lambdafunc etc inside the bregmanparams?

Now there is the options and the additional kwargs. Like user have to set options then have to put new ones separately like lambda. I still think its a bit messy. Put back the old basic interface and use dispatch to properly set the other cases in a clean way not adding 50 kwargs

Can't dispatch over a quantile or a pre-set threshold. They are both a single number. Any thought?

99% of what you added can be easily done via simple dispatch rather than all these additional keyword arguments on top of the options, for exampel would have takien you one line to define

bregman(A, x, b) = bregman(A, LinearAlgebra.I, x::Array{T}, b)

And everything satays clean an easy to use. Try to make an effort making it user friendly not "works for me" friendly. All the rest is the same can be easily done via dispatch and/or options preprocessing and setup.

ok thanks for your suggestion on this LinearAlgebra.I thing. I agree this works in a more reasonable structure.

What's your take on this threshold thing? Following your suggestion, I should add it in the bregmanparams, right? Since a pre-set threshold and a quantile are both single value. Can't dispatch.

Options are all kwargs so you can filter then and do so pre-setup based on input

ziyiyin97 · 2022-04-20T02:06:06Z

Now previous functionality still works

mloubout · 2022-04-20T02:31:41Z

examples/denoising.jl

@@ -20,11 +20,11 @@ imgn= img .+ .01f0*randn(Float32, size(img))
 b = A*vec(imgn)

 # setup bregamn
-opt = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true)
-opt2 = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true, spg=true)
+opt = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true, TD=W)


at that point let's just put A and x and b in options too........

I think the bregman function itself should show A and b (or the function to calculate function value and gradient) since we are solving a linear system. In terms of the sparsifying transform ... I slightly prefer putting it into the options (since we have different "options" for sparsifying transform)

mloubout

I don't mind moving TD to the options but this still needs cleanup

mloubout · 2022-04-20T02:51:41Z

src/bregman.jl

-bregman_options(;verbose=1, progTol=1e-8, maxIter=20, store_trace=false, antichatter=true, quantile=.95, alpha=.5, spg=false) =
-                BregmanParams(verbose, progTol, maxIter, store_trace, antichatter, quantile, alpha, spg)
+bregman_options(;verbose=1, progTol=1e-8, maxIter=20, store_trace=false, antichatter=true, alpha=.5, spg=false, TD=LinearAlgebra.I, quantile=.95, λ=nothing, λfunc=nothing) =
+                BregmanParams(verbose, progTol, maxIter, store_trace, antichatter, alpha, spg, TD, quantile, λ, λfunc)


This is still messy

Some of these options are not compatible

You wouldn't need to check all these if options.x below if you just properly filter the input options to setup the problem

mloubout · 2022-04-20T02:52:30Z

src/bregman.jl

+            λfunc = z->quantile(abs.(z), options.quantile)
+        end
+    end
+    return bregman(funobj, x, options, λfunc)


now you have to call the algo with λfunc both in options and as an input which you wuldn't with properly setup options

mloubout · 2022-04-20T02:53:49Z

src/bregman.jl

+end
+
+function bregman(funobj::Function, TD, x::AbstractArray{T}, options=bregman_options()) where {T}
+    options.TD = TD


Mutating structure and objects is in general very bad practice coding-wise, so this needs a proper deprecation warning.

mloubout · 2022-04-20T02:54:13Z

src/bregman.jl

@@ -124,16 +160,14 @@ function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams,
        # Update z variable
        @. z = z + d
        # Get λ at first iteration
-        i == 1 && (λ = abs(T(quantile(abs.(z), options.quantile))))
+        (i == 1) && (sol.λ = λ = abs.(T.(λfunc(z))))


or you can just use sol.λ everywhere

mloubout · 2022-04-20T02:55:01Z

src/bregman.jl

 """
-function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams, TD=nothing) where {T}
+


No blank line between docstring and function

mloubout · 2022-04-20T02:55:32Z

src/bregman.jl

 """
-function bregman(A, TD, x::Array{T}, b, options) where {T}
+


No blank line between docstring and code

mloubout · 2022-04-20T02:56:09Z

src/bregman.jl


 """
-    bregman(A, TD, x, b, options)
+    bregman(A, x, b; options)


TD isn't defined anymore in problem def line 47

mloubout · 2022-04-20T02:57:04Z

src/bregman.jl

 end

 """
    bregman_options(;verbose=1, optTol=1e-6, progTol=1e-8, maxIter=20
-                    store_trace=false, linesearch=false, alpha=.25, spg=false)
+                    store_trace=false, λ=.2, alpha=.25, spg=false)


where does .2 comes from

just to inform user that you can do something like this

mloubout · 2022-04-20T02:58:21Z

test/test_bregman.jl

@@ -5,9 +5,9 @@ using LinearAlgebra

 N1 = 100
 N2 = div(N1, 2) + 5
-A = randn(N1, N2)
+A = randn(ComplexF32, N1, N2)


If you want to add a case add an extra one do not modify an existing one especially if you are disabling options

mloubout · 2022-04-20T02:59:27Z

test/test_bregman.jl

@@ -20,8 +20,8 @@ function obj(x)
    return fun, grad
 end

-opt = bregman_options(maxIter=200, progTol=0, verbose=2)
-sol = bregman(obj, 1 .+ randn(N2), opt)
+opt = bregman_options(maxIter=200, progTol=0, verbose=2, antichatter=false) # anti chatter now only works with real number


This needs to raise a proper error if attempted. Why isn't it supported?

Anti-chatter needs to record number of -1/1 when the gradient fluctuates. If things are complex then there is no -1/1

Yes there is no sign in complex which is why there is a different thresholding function for this case so the anti-chatter should follow it and use the angle instead of sign but I think it should follow

hmm in anti-chatter there is a sum of -1/1. How do we do it with angle? A sum of unit vectors with different angles? And then do inner product to calculate the scaling of gradient?

Probably something like that. I don't mind it leaving it as a todo and see if can make something like that work, would assume it would by summing the angles. But would adda warning or error if someone tries it so that it's readable not a lengthy error trace

gotcha make sense

mloubout · 2022-05-03T22:19:00Z

src/bregman.jl

 """
-function bregman(A, TD, x::Array{T}, b, options) where {T}
+function bregman(A, x::AbstractArray{T}, b::AbstractArray{T}, options::BregmanParams) where {T}


No! There is defaults for a reason

mloubout

Mostly there

mloubout · 2022-05-04T13:13:41Z

src/bregman.jl

+        if ~isnothing(λ) 
+            λfunc = z->λ
+        else
+            λfunc = z->SlimOptim.quantile(abs.(z), quantile)


Statistics.

mloubout · 2022-05-04T13:14:43Z

src/bregman.jl

 """
-function bregman(A, TD, x::Array{T}, b, options) where {T}
+function bregman(A, x::AbstractArray{T}, b::AbstractArray{T}, options::BregmanParams=bregman_options()) where {T}


AbstractVector no? May wanna add it to obj below as well

mloubout · 2022-05-04T13:16:09Z

src/bregman.jl

@@ -97,7 +116,7 @@ function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams,
    # Result structure
    sol = breglog(x, z)
    # Initialize λ
-    λ = abs(T(0))
+    sol.λ = abs(T(0))


Can't this be done in options init now?

mloubout · 2022-05-04T13:17:23Z

test/test_bregman.jl

@@ -36,3 +36,35 @@ part_nz = i -> norm(sol.x[i], 1)/N2
 @test part_nz(inds) < 1f-1
 @test part_n(ninds) < 1f-1
 @test sol.residual/sol.r_trace[1] < 1f-1
+
+# test complex
+A = randn(ComplexF32, N1, N2)


For T in [float, complex]...

mloubout · 2022-05-04T13:17:41Z

src/bregman.jl

        norm(x - sol.x) < options.progTol && (@printf("Step size below progTol\n"); break;)
        update!(sol; iter=i, ϕ=obj_fun, residual=f, x=x, z=z, g=g, store_trace=options.store_trace)
    end
    return sol
 end

+function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams, TD) where {T}
+    @warn "deprecation warning: TD should be put in BregmanParams when version >= 0.1.8; now overwritting TD in BregmanParams"


Add new syntax explicitly to message

add a user define option for threshold in bregman

0514449

mloubout reviewed Apr 4, 2022

View reviewed changes

have a custom thresholding function for first iter

32ff641

mloubout reviewed Apr 7, 2022

View reviewed changes

ziyiyin97 added 2 commits April 16, 2022 19:59

update options for lambda function

aa9dcc6

sorry about slightly changed api

24f4ec8

mloubout reviewed Apr 17, 2022

View reviewed changes

ziyiyin97 added 2 commits April 19, 2022 21:53

dispatch, make kwargs into options

b9dc130

doesnt change API

5604792

mloubout reviewed Apr 20, 2022

View reviewed changes

ziyiyin97 added 2 commits April 30, 2022 01:08

pre-process at bregmanparams

6c0f7d7

clean up tests and documentation

899325e

ziyiyin97 requested a review from mloubout May 2, 2022 15:32

TD is at the end for funobj bregman

649bebc

mloubout requested changes May 3, 2022

View reviewed changes

do defaults

0fad8aa

mloubout reviewed May 4, 2022

View reviewed changes

ziyiyin97 added 2 commits May 4, 2022 21:05

fixed all

52bdfec

don't need to be in same type

25f8817

ziyiyin97 requested a review from mloubout May 6, 2022 20:17

mloubout merged commit c43c82b into master May 6, 2022

mloubout deleted the bregman branch May 6, 2022 22:41

		"""
		function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams, TD=nothing) where {T}

		"""
		function bregman(A, TD, x::Array{T}, b, options) where {T}

add a user define option for threshold in bregman #3

add a user define option for threshold in bregman #3

Conversation

ziyiyin97 commented Apr 4, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout Apr 4, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ziyiyin97 commented Apr 20, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout Apr 4, 2022 •

edited

Loading