FISTA for tomographic reconstruction

[ndjurabe]

Hello, I am trying to perform FISTA reconstruction but the f_prox = f.proximal(gamma) line in proximal_gradient_solvers.py throws me an error because (I guess) my regularizer functional f falls into the "functional(operator)" class (in functional.py) which is not implemented. Am I not defining something correctly?

My data discrepancy functional g (||A(x) - sino||^2_2) is defined as follows (based on some other odl examples):

l2_norm_sq = 0.5 * odl.solvers.L2NormSquared(forward_operator.range).translated(sinogram)
data_discr= l2_norm_sq*forward_operator

and the regularizer f (TV(x)) as follows:

grad = odl.Gradient(reco_space)
reg_param = 0.05
regularizer = reg_param * odl.solvers.L2NormSquared(grad.range)*grad 

(I also tried:

l1_norm = odl.solvers.GroupL1Norm(grad.range)
regularizer = odl.solvers.MoreauEnvelope(l1_norm, sigma=reg_param)*grad

but that gave the same result)

and then calling the solver with:

rec = reco_space.zero()
odl.solvers.accelerated_proximal_gradient(rec, f=regularizer, g=data_discr, niter=50, gamma=step_size, callback=callback)

[JevgenijaAksjonova]

It might be a bit confusing, but one must define f to be a regularizer and g to be data discrepancy, see
https://odlgroup.github.io/odl/generated/odl.solvers.nonsmooth.proximal_gradient_solvers.accelerated_proximal_gradient.html

and also this example

https://github.com/odlgroup/odl/blob/master/examples/solvers/proximal_gradient_wavelet_tomography.py

[ndjurabe]

It might be a bit confusing, but one must define f to be a regularizer and g to be data discrepancy, see
https://odlgroup.github.io/odl/generated/odl.solvers.nonsmooth.proximal_gradient_solvers.accelerated_proximal_gradient.html

and also this example

https://github.com/odlgroup/odl/blob/master/examples/solvers/proximal_gradient_wavelet_tomography.py

Thanks for getting back to me! Soon after posting my question I realized this and made the edit accordingly (not sure if you can see the edit). I have looked at the proximal gradient wavelet tomography example and ran it but I don't seem to be able to get essentially the same thing to work without the wavelet operator (so just || Ax-f ||_2^2 + Tv(x), where A is the ray_trafo operator defined using Astra). Could you let me know if I am defining my data discrepancy/ TV operators wrong somehow?

[JevgenijaAksjonova]

Why not using

regularizer = reg_param * odl.solvers.L1Norm(grad.range)*grad

?

[ndjurabe]

Why not using

regularizer = reg_param * odl.solvers.L1Norm(grad.range)*grad

?

I tried that but it still gives me the same "not implemented" error:

Traceback (most recent call last):
File "/usr/local/anaconda3/envs/tomosipo38/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3441, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "", line 1, in
regularizer.proximal(gamma)
File "/usr/local/anaconda3/envs/tomosipo38/lib/python3.8/site-packages/odl/solvers/functional/functional.py", line 157, in proximal
raise NotImplementedError(
NotImplementedError: no proximal operator implemented for functional FunctionalComp(FunctionalLeftScalarMult(L1Norm(ProductSpace(uniform_discr([-5., -5.], [ 5., 5.], (480, 480), dtype='float32'), 2)), 0.05), Gradient(
uniform_discr([-5., -5.], [ 5., 5.], (480, 480), dtype='float32')
))

[JevgenijaAksjonova]

Right, the proximal is not implemented for the composition of l1 and gradient. This is because there is no closed form formula for that. If you want to do TV regularization, you have to apply another optimization algorithm, for instance, ADMM:

https://github.com/odlgroup/odl/blob/master/examples/solvers/admm_tomography.py

[ndjurabe]

ok, so then there's no proximal for TV I guess - thank you for the help! :)