Multiple regularizers

[hakkelt]

Hi, I wanted to try multiple regularizers, but I couldn't manage to get it to work. I checked the source code, and it seemed to support this feature. I tried with the following (among many other variations, but this seemed to most logical way based on the code):

nx, ny = 10, 10
traj = MRIReco.CartesianTrajectory(Float32, ny, nx)
acq = MRIReco.AcquisitionData(traj, reshape([ones(ComplexF32, nx*ny, 1)], 1, 1, 1))
params = Dict{Symbol,Any}()
params[:reco] = "standard"
params[:reconSize] = (nx, ny)
params[:reg] = [MRIReco.L1Regularization(0.05), MRIReco.TVRegularization(0.2, shape = (nx, ny))]
params[:sparseTrafo] = ["Wavelet", nothing]
reconstruction(acq, params);

I always get this error: AssertionError: reg and regTrafo must have the same length
Does MRIReco support this (assuming ADMM as the algorithm)?

[nHackel]

@hakkelt sorry for the very late reply, this issue must have slipped past me.

Yes we support multiple regularizers. The solver that support multiple proximal term is often applied applied separately from the transformation in which the proximal term is applied. Here you can see the documentation for ADMM.

In your case, you'd need to supply an additional regTrafo command of equal length to your reg vector:

# for no transformation you can supply an opEye with the appropriate eltype and size(AHA, 1)
params[:regTrafo] = [opEye(...), opEye(...))]

# in your case it would probably be:
params[:regTrafo] = [WaveletOp(ComplexF32, shape = (nx, ny)), opEye(...)]

[nHackel]

A better user-experience would be if we transform the sparseTrafo to regTrafo is applicable and a user didn't supply any

[hakkelt]

Thanks, @nHackel, I managed to make my code work based on your suggestion. For the record, the fixed code snippet is the following:

nx, ny = 10, 10
traj = MRIReco.CartesianTrajectory(Float32, ny, nx)
acq = MRIReco.AcquisitionData(traj, reshape([ones(ComplexF32, nx*ny, 1)], 1, 1, 1))
params = Dict{Symbol,Any}()
params[:reco] = "standard"
params[:reconSize] = (nx, ny)
params[:reg] = [MRIReco.L1Regularization(0.05), MRIReco.TVRegularization(0.2, shape = (nx, ny))]
params[:regTrafo] = [WaveletOp(ComplexF32, shape = (nx, ny)), opEye(ComplexF32, nx*ny, nx*ny)]
reconstruction(acq, params);

Btw, I find sparseTrafo and regTrafo quite confusing. What is the difference between the two? I could not find any mentions of regTrafo in the documentation.

[nHackel]

They are essentially an implementation detail of the underlying solvers of RegularizedLeastSquares.jl. SparseTrafo result in a regularization term of form:

sparseOp = WaveletOp(...)
core = L1Regularization(...)
reg = TransformedRegularization(core, sparseOp) # i.e. apply core term in sparseOp domain

while regTrafo is a feature of certain solver of RegularizedLeastSquares.jl such as ADMM and SplitBregman. Those solver have a built in transformation to another domain where regularization terms are applied.

MRIReco is only aware of sparseTrafo at the moment and the rest is just passed on to the solver. The better way would be for MRIReco to always accept sparseTrafo and either built a TransformedRegularization term or to translate to regTrafo depending on the solver used