Multidimensional arrays

[1oly]

Hello,

Thank you for developing and maintaining such a useful package :)

I'm trying to figure out the best way to deal with multidimensional arrays. I've touched this subject before in https://github.com/kul-forbes/AbstractOperators.jl/pull/13, which helped but my method is still sub-optimal. I've tried to make a MWE to illustrate my workflow:

using ImageFiltering, FFTW, AbstractOperators, ProximalAlgorithms, ProximalOperators

n = 25
p1 = Kernel.gaussian(3)
p2 = Kernel.gaussian(6)

xtrue = zeros(n,n,2)
xtrue[13,13,:] .= 1
b1 = imfilter(xtrue[:,:,1], p1)
b2 = imfilter(xtrue[:,:,2], p2)

b = zeros(n,n,2)
b[:,:,1] = b1
b[:,:,2] = b2

sp = round.(Int,(size(p2) .- size(p1) )./2)
p1_pad = parent(padarray(p1, Fill(0,sp,sp)))

kern = zeros(n,n,2)
kern[:,:,1] = p1_pad
kern[:,:,2] = p2

center = round.(Int,size(p2)./2).+1
Fps = fft(circshift(kern,center),(1,2))
ifft_F = IDFT(Complex{Float64}, (n,n))
fft_F = DFT((n,n))

solver = ProximalAlgorithms.ForwardBackward(tol = 1e-12;maxit=20000,fast=true,adaptive=true)
g = IndNonnegative()
bc = complex(b)
xp1 = similar(kern)
for i = 1:size(kern,3)
    f = Translate(SqrNormL2(), -bc[:,:,i])
    psf = DiagOp(Fps[:,:,i])
    A = Compose(ifft_F,Compose(psf,fft_F))
    xp1[:,:,i],_ = solver(zeros(n,n), f=f, A=A, g=g)
end

all(isapprox.(xp1[13,13,:],1.0;rtol=1e-4)) # true

##### 3D CASE: #######

ifft_F = IDFT(Complex{Float64}, (n,n,2),(1,2))
fft_F = DFT((n,n,2),(1,2))

f = Translate(SqrNormL2(), -bc)
psf = DiagOp(Fps)
A = Compose(ifft_F,Compose(psf,fft_F))
xp2 = similar(kern)
xp2, _= solver(zeros(n,n,2), f=f, A=A, g=g)

all(isapprox.(xp2[13,13,:],1.0;rtol=1e-4)) # false

Basically, I'm doing image deblurring in acoustic imaging, where I have many 2D 'images' for different frequency bins. So in my setup I'm doing a lot of processing on 3D arrays. In the above example I therefore need to loop over the 3rd dimension of kern/psf and b and construct f and A in each loop.

I see two optimization strategies here:

1. Implement broadcasting for operators such as Translate, DiagOp and Compose to allow in-place operations .=

2. Modify solvers to accept 3D arrays (it actually already does, see MWE) and calculate gamma from a slice along a certain dimension (this is currently the issue I'm facing).

I might have misunderstood the solver interface, so forgive me if I'm mistaken. I would appreciate any comments on how to improve the performance.

Thanks! /Oliver

[nantonel]

It looks like you're calling the solvers correctly.

Regarding point 1. these are functions that are not exported by ProximalAlgorithms to this package but rather to AbstractOperators and ProximalOperators. So this issue is not related to this package.

Are you saying that having 1. and 2. will get you to all(isapprox.(xp2[13,13,:],1.0;rtol=1e-4)) = true?

[1oly]

Hi, and sorry if I wasn't clear. I'm trying to reduce allocations arising from redefining f and A at every iteration. This might be possible by simply using a different approach, i.e., without changing anything in ProximalAlgorithms or AbstractOperators. Or possibly it requires either 1. or 2. Broadcasting would help because f and A could be overwritten in-place.

The issue in the 3D case is that the value xp2[13,13,:] is underestimated possibly due to the computation of gamma with a 3D operator A (I haven't investigated how gamma is computed exactly).

[nantonel]

I'm not sure that having gamma set to the optimal one will get your condition satisfied. Having optimal gamma might only get you to the optimal solution with less iterations. My impression is that your condition may simply be unfeasible. However, I'm not familiar with your problem.

Concerning the allocations of f and A, this is related to ProximalOperators and AbstractOperator respectively. Just open an issue in the right package and we can discuss about it!

[1oly]

Thanks! Will look into the details and open a new issue if the issue persist.