[1]:

#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
#
#  Authored by:    Tommi Heikkilä (LUT)
#  Reviewed by:    Margaret Duff (STFC-UKRI), Jakob Jørgensen (DTU)

[2]:

import cil
print(f"Using CIL version {cil.__version__}")

Using CIL version 24.2.0

Controlled Wavelet Sparsity using callbacks (3D data)#

This notebook showcases an automated method for tuning the regularization parameter towards a good value during the optimization process. It is organized as follows.

Set up the data and the appropriate geometry. The projection images need to be converted to sinogram.
Brief introduction to a representation system called wavelets and explaining what wavelet regularization should do.
- Traditional run with fixed regularization parameter using FISTA
Brief explanation of the controlled wavelet sparsity method including the original source and how the CIL implementation is done.
Some comments on the limitations of the method

3D laboratory micro-CT, fan-beam data of a seashell#

This example uses the seashell dataset 20211124_seashell.zip from https://doi.org/10.5281/zenodo.6983008 :

https://zenodo.org/records/6983008/files/20211124_seashell.zip?download=1

If running locally please download the data, unzip the folder and update the ‘path’ variable below.

[3]:

path = './data'

[4]:

# remove some annoying warnings
import logging
logger = logging.getLogger('dxchange')
logger.setLevel(logging.ERROR)

First import all of the modules we will need:

[5]:

import os
import numpy as np

# cil imports
from cil.framework import ImageData, ImageGeometry
from cil.framework import AcquisitionGeometry, AcquisitionData

from cil.optimisation.functions import LeastSquares, L1Sparsity, ScaledFunction
from cil.optimisation.operators import WaveletOperator
from cil.optimisation.algorithms import FISTA

from cil.plugins.astra import ProjectionOperator

from cil.utilities.display import show2D, show_geometry
from cil.utilities.jupyter import islicer

from cil.io import TIFFStackReader

import scipy.io

Load the seashell projections images (TIFFs). There are 721 projections in total, with 0.5 degree interval. We do not need all of them so let’s create data with sparser set of projections.

[6]:

skip = 6
roi = {'axis_0': (0, -1, skip), 'axis_1': -1, 'axis_2': -1}
reader = TIFFStackReader(file_name=path + '/20211124_seashell', transpose=False, roi=roi)
angles = np.linspace(0., 360., 720//skip, endpoint=False)

Create geometry#

We have to manually set the values for the correct geometry. These can be found in the documentation of the data: 20211124_seashell_.txt.

[7]:

# Acquisition geometry parameters
SOD = 210.66 # Source-Origin-distance
print(f'SOD: {SOD:.3f} mm')
SDD = 553.74 # Source-Detector-distance
print(f'SDD: {SDD:.3f} mm')
ODD = SDD - SOD # Origin-detector-distance
print(f'ODD: {ODD:.3f} mm')
pixelSize = 0.05
print(f'Detector pixel size: {pixelSize:.3f} mm')

num_pixels = (2240, 2368) # Projection image width and height

SOD: 210.660 mm
SDD: 553.740 mm
ODD: 343.080 mm
Detector pixel size: 0.050 mm

[8]:

ag_big = AcquisitionGeometry.create_Cone3D([0,-SOD, 0], [0,ODD, 0], units='mm', detector_direction_x=[1, 0, 0], rotation_axis_direction=[0, 0, -1])
ag_big.set_panel(num_pixels, pixel_size=(pixelSize, pixelSize), origin='bottom-left')
ag_big.set_angles(angles)
ag_big.set_labels(('angle', 'vertical', 'horizontal'))
print(ag_big)

3D Cone-beam tomography
System configuration:
        Source position: [   0.  , -210.66,    0.  ]
        Rotation axis position: [0., 0., 0.]
        Rotation axis direction: [ 0.,  0., -1.]
        Detector position: [  0.  , 343.08,   0.  ]
        Detector direction x: [1., 0., 0.]
        Detector direction y: [0., 0., 1.]
Panel configuration:
        Number of pixels: [2240 2368]
        Pixel size: [0.05 0.05]
        Pixel origin: bottom-left
Channel configuration:
        Number of channels: 1
Acquisition description:
        Number of positions: 120
        Angles 0-9 in degrees: [ 0.,  3.,  6.,  9., 12., 15., 18., 21., 24., 27.]
        Angles 110-119 in degrees: [330., 333., 336., 339., 342., 345., 348., 351., 354., 357.]
        Full angular array can be accessed with acquisition_data.geometry.angles
Distances in units: mm

[9]:

show_geometry(ag_big);

../../_images/demos_controlledWaveletSparsity_3d_13_0.png

Read data and add geometry#

[10]:

big_raw_data = reader.read_as_AcquisitionData(ag_big)

[11]:

print(big_raw_data)
show2D(big_raw_data);

Number of dimensions: 3
Shape: (120, 2368, 2240)
Axis labels: (<AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.VERTICAL: 'vertical'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)

../../_images/demos_controlledWaveletSparsity_3d_16_1.png

Process the data#

We want to bin the data to smaller resolution and convert from transmission to absorption images

[12]:

from cil.processors import Binner, TransmissionAbsorptionConverter


binning = 10
binner_processor = Binner(roi={'horizontal': (220, 2020, binning), 'vertical': (84, 2284, binning)})
binner_processor.set_input(big_raw_data)
raw_data = binner_processor.get_output()

# Get the raw intensity from the mean of the k largest values
k = 100
max_k = np.partition(raw_data.as_array(), -k, axis=None)[-k:].mean()
transmission_processor = TransmissionAbsorptionConverter(white_level=max_k)
transmission_processor.set_input(raw_data)
data = transmission_processor.get_output()

show2D(data)
print(data)

../../_images/demos_controlledWaveletSparsity_3d_18_0.png

Number of dimensions: 3
Shape: (120, 220, 180)
Axis labels: (<AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.VERTICAL: 'vertical'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)

Image geometry and data#

[13]:

# Setup image geometry
(_, nz, nx) = data.shape
ig = data.geometry.get_ImageGeometry()
ig.voxel_num_x = nx
ig.voxel_num_y = nx
ig.voxel_num_z = nz
print(ig)

Number of channels: 1
channel_spacing: 1.0
voxel_num : x180,y180,z220
voxel_size : x0.1902156246613934,y0.1902156246613934,z0.1902156246613934
center : x0,y0,z0

Reconstructions#

[14]:

data.reorder('astra')
ag = data.geometry
A = ProjectionOperator(ig, ag, 'gpu')

x_bp = A.adjoint(data) # Simple backprojection
show2D(x_bp, size=(10,10), title="Naive backprojection");

../../_images/demos_controlledWaveletSparsity_3d_22_0.png

[15]:

islicer(x_bp)

[15]:

Wavelet based regularization#

Classical approach in imaging problems (such as CT, denoising, etc.) is to try represent the key information about the solution using a wavelet basis. By construction a relatively sparse solution (i.e. only few nonzero components) should suffice, as opposed to pixel-based representation.

In short, we have the least squares data mismatch term of form

\[\frac{1}{2} \| A x - m_\delta \|_2^2,\]

where \(A\) is the forward operator, \(m_\delta\) is the noisy data and \(x\) is the unknown solution. Then we add a regularization term of the form

\[\alpha \| W(x) \|_1 = \alpha \sum_{j,k} | c(x)_{j,k} |,\]

where \(\alpha > 0\) is the regularization term (balacing our trust in data and need in regularization) and \(W\) is the Discrete Wavelet Transform operator which decomposes any input \(x\) (signal, image or volume) into a sequence of coefficients \(c(x)_{j,k}\) of varying scale (\(j \in \mathbb{N}\)) and location (\(k \in \mathbb{Z}^d\)), such that

\[x = \sum_{j,k} c(x)_{j,k} \psi_{j,k}.\]

Wavelet transform#

Classical choice is one of the orthogonal wavelets of Ingrid Daubechies, usually distinguished by their smoothness properties. For example the Daubechies-2 wavelets always give continuous approximations.

[16]:

wname = 'db2'
level = 4

W = WaveletOperator(ig, level=level, wname=wname)

Fast Iterative Soft Thresholding Algorithm (FISTA)#

Equally classic choice of minimization algrorithm in the (Fast) Iterative Soft Thresholding Algorithm (ISTA), originally developed for wavelet-based denoising in
> Daubechies, I, Defrise, M, & De Mol, C..
An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.
Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 57.11 (2004): 1413-1457.

Although we will run the accelerated version.

[17]:

f = LeastSquares(A, data)
g = L1Sparsity(W)

tau = 1e-5  # Regularization parameter

fista = FISTA(initial=ig.allocate(0), f=f, g=tau*g)

Let’s run arbitrary number of iterations (200)

[18]:

fista.run(200)

[19]:

x_fista = fista.get_output()
show2D(x_fista, size=(10,10), fix_range=(0,0.4), title=f"FISTA with {tau = }", num_cols=3, slice_list=[(0, 90), (1, 90), (2, 90)]);

../../_images/demos_controlledWaveletSparsity_3d_31_0.png

[20]:

islicer(x_fista, minmax=(0,0.4), title=f"FISTA with {tau = },")

[20]:

Looks okay but how do we know the choice of regularization parameter was good? Could we get a better result with another value? Do we have the patience to try out dozens or even hundreds of values?

Automated Controlled Sparsity#

We follow the automated tuning strategy as introduced in
> Purisha, Z., Rimpeläinen, J., Bubba, T., & Siltanen, S. (2017).
Controlled wavelet domain sparsity for x-ray tomography.
Measurement Science and Technology, 29(1), 014002.
DOI:10.1088/1361-6501/aa9260

We’ll try to implement it as a StepSizeRule and callback. However we are NOT actually changing the step_size (which is based on the Lipschitz constant of the data mistmatch term) but the regularization parameter tau. However, we can update it every iterate by using StepSizeRule (in fact we could pass it an actual StepSizeRule to perform but let’s just stick with constant value for now).

[21]:

from SparsityControl import ControlledSparsity, DesiredSparsity

[22]:

reg_weight = 1e-2
# Automated tuning of regularization parameter, smaller step_weight makes the tuning slightly slower -> not as severe overshoot in reg. param. value
sparsity_ctrl = ControlledSparsity(desired_sparsity=0.3, step_weight=0.1)
fista_controlled = FISTA(initial=ig.allocate(0), f=f, g=reg_weight*g, step_size=sparsity_ctrl, max_iteration=500)

[23]:

# Sparsity based stopping rule
sparsity_stop = DesiredSparsity(print_stuff=2, rel_change_tol=5e-4)

[24]:

fista_controlled.run(500, callbacks=[sparsity_stop])

Iteration : 5 | Current sparsity: 0.751 | Regularization param.: 1.02e-02 | Desired sparsity: 0.300 | Difference: 0.451 | Beta: 2.55e-04 | Relative change: 1.27e-01
Iteration : 10 | Current sparsity: 0.747 | Regularization param.: 1.07e-02 | Desired sparsity: 0.300 | Difference: 0.447 | Beta: 2.55e-04 | Relative change: 6.41e-02
Iteration : 15 | Current sparsity: 0.739 | Regularization param.: 1.13e-02 | Desired sparsity: 0.300 | Difference: 0.439 | Beta: 2.55e-04 | Relative change: 3.75e-02
Iteration : 20 | Current sparsity: 0.731 | Regularization param.: 1.18e-02 | Desired sparsity: 0.300 | Difference: 0.431 | Beta: 2.55e-04 | Relative change: 2.38e-02
Iteration : 25 | Current sparsity: 0.722 | Regularization param.: 1.24e-02 | Desired sparsity: 0.300 | Difference: 0.422 | Beta: 2.55e-04 | Relative change: 1.58e-02
Iteration : 30 | Current sparsity: 0.713 | Regularization param.: 1.29e-02 | Desired sparsity: 0.300 | Difference: 0.413 | Beta: 2.55e-04 | Relative change: 1.09e-02
Iteration : 35 | Current sparsity: 0.703 | Regularization param.: 1.34e-02 | Desired sparsity: 0.300 | Difference: 0.403 | Beta: 2.55e-04 | Relative change: 7.90e-03
Iteration : 40 | Current sparsity: 0.692 | Regularization param.: 1.39e-02 | Desired sparsity: 0.300 | Difference: 0.392 | Beta: 2.55e-04 | Relative change: 6.20e-03
Iteration : 45 | Current sparsity: 0.682 | Regularization param.: 1.44e-02 | Desired sparsity: 0.300 | Difference: 0.382 | Beta: 2.55e-04 | Relative change: 5.25e-03
Iteration : 50 | Current sparsity: 0.671 | Regularization param.: 1.49e-02 | Desired sparsity: 0.300 | Difference: 0.371 | Beta: 2.55e-04 | Relative change: 4.65e-03
Iteration : 55 | Current sparsity: 0.660 | Regularization param.: 1.54e-02 | Desired sparsity: 0.300 | Difference: 0.360 | Beta: 2.55e-04 | Relative change: 4.17e-03
Iteration : 60 | Current sparsity: 0.649 | Regularization param.: 1.58e-02 | Desired sparsity: 0.300 | Difference: 0.349 | Beta: 2.55e-04 | Relative change: 3.74e-03
Iteration : 65 | Current sparsity: 0.637 | Regularization param.: 1.63e-02 | Desired sparsity: 0.300 | Difference: 0.337 | Beta: 2.55e-04 | Relative change: 3.37e-03
Iteration : 70 | Current sparsity: 0.626 | Regularization param.: 1.67e-02 | Desired sparsity: 0.300 | Difference: 0.326 | Beta: 2.55e-04 | Relative change: 3.07e-03
Iteration : 75 | Current sparsity: 0.614 | Regularization param.: 1.71e-02 | Desired sparsity: 0.300 | Difference: 0.314 | Beta: 2.55e-04 | Relative change: 2.84e-03
Iteration : 80 | Current sparsity: 0.602 | Regularization param.: 1.75e-02 | Desired sparsity: 0.300 | Difference: 0.302 | Beta: 2.55e-04 | Relative change: 2.66e-03
Iteration : 85 | Current sparsity: 0.590 | Regularization param.: 1.79e-02 | Desired sparsity: 0.300 | Difference: 0.290 | Beta: 2.55e-04 | Relative change: 2.51e-03
Iteration : 90 | Current sparsity: 0.578 | Regularization param.: 1.82e-02 | Desired sparsity: 0.300 | Difference: 0.278 | Beta: 2.55e-04 | Relative change: 2.38e-03
Iteration : 95 | Current sparsity: 0.566 | Regularization param.: 1.86e-02 | Desired sparsity: 0.300 | Difference: 0.266 | Beta: 2.55e-04 | Relative change: 2.26e-03
Iteration : 100 | Current sparsity: 0.554 | Regularization param.: 1.89e-02 | Desired sparsity: 0.300 | Difference: 0.254 | Beta: 2.55e-04 | Relative change: 2.15e-03
Iteration : 105 | Current sparsity: 0.542 | Regularization param.: 1.92e-02 | Desired sparsity: 0.300 | Difference: 0.242 | Beta: 2.55e-04 | Relative change: 2.06e-03
Iteration : 110 | Current sparsity: 0.530 | Regularization param.: 1.95e-02 | Desired sparsity: 0.300 | Difference: 0.230 | Beta: 2.55e-04 | Relative change: 1.97e-03
Iteration : 115 | Current sparsity: 0.518 | Regularization param.: 1.98e-02 | Desired sparsity: 0.300 | Difference: 0.218 | Beta: 2.55e-04 | Relative change: 1.89e-03
Iteration : 120 | Current sparsity: 0.506 | Regularization param.: 2.01e-02 | Desired sparsity: 0.300 | Difference: 0.206 | Beta: 2.55e-04 | Relative change: 1.82e-03
Iteration : 125 | Current sparsity: 0.495 | Regularization param.: 2.03e-02 | Desired sparsity: 0.300 | Difference: 0.195 | Beta: 2.55e-04 | Relative change: 1.75e-03
Iteration : 130 | Current sparsity: 0.484 | Regularization param.: 2.05e-02 | Desired sparsity: 0.300 | Difference: 0.184 | Beta: 2.55e-04 | Relative change: 1.68e-03
Iteration : 135 | Current sparsity: 0.473 | Regularization param.: 2.08e-02 | Desired sparsity: 0.300 | Difference: 0.173 | Beta: 2.55e-04 | Relative change: 1.62e-03
Iteration : 140 | Current sparsity: 0.463 | Regularization param.: 2.10e-02 | Desired sparsity: 0.300 | Difference: 0.163 | Beta: 2.55e-04 | Relative change: 1.56e-03
Iteration : 145 | Current sparsity: 0.453 | Regularization param.: 2.12e-02 | Desired sparsity: 0.300 | Difference: 0.153 | Beta: 2.55e-04 | Relative change: 1.50e-03
Iteration : 150 | Current sparsity: 0.444 | Regularization param.: 2.14e-02 | Desired sparsity: 0.300 | Difference: 0.144 | Beta: 2.55e-04 | Relative change: 1.45e-03
Iteration : 155 | Current sparsity: 0.435 | Regularization param.: 2.15e-02 | Desired sparsity: 0.300 | Difference: 0.135 | Beta: 2.55e-04 | Relative change: 1.40e-03
Iteration : 160 | Current sparsity: 0.426 | Regularization param.: 2.17e-02 | Desired sparsity: 0.300 | Difference: 0.126 | Beta: 2.55e-04 | Relative change: 1.35e-03
Iteration : 165 | Current sparsity: 0.418 | Regularization param.: 2.19e-02 | Desired sparsity: 0.300 | Difference: 0.118 | Beta: 2.55e-04 | Relative change: 1.31e-03
Iteration : 170 | Current sparsity: 0.410 | Regularization param.: 2.20e-02 | Desired sparsity: 0.300 | Difference: 0.110 | Beta: 2.55e-04 | Relative change: 1.27e-03
Iteration : 175 | Current sparsity: 0.403 | Regularization param.: 2.21e-02 | Desired sparsity: 0.300 | Difference: 0.103 | Beta: 2.55e-04 | Relative change: 1.23e-03
Iteration : 180 | Current sparsity: 0.396 | Regularization param.: 2.23e-02 | Desired sparsity: 0.300 | Difference: 0.096 | Beta: 2.55e-04 | Relative change: 1.19e-03
Iteration : 185 | Current sparsity: 0.390 | Regularization param.: 2.24e-02 | Desired sparsity: 0.300 | Difference: 0.090 | Beta: 2.55e-04 | Relative change: 1.15e-03
Iteration : 190 | Current sparsity: 0.384 | Regularization param.: 2.25e-02 | Desired sparsity: 0.300 | Difference: 0.084 | Beta: 2.55e-04 | Relative change: 1.12e-03
Iteration : 195 | Current sparsity: 0.378 | Regularization param.: 2.26e-02 | Desired sparsity: 0.300 | Difference: 0.078 | Beta: 2.55e-04 | Relative change: 1.09e-03
Iteration : 200 | Current sparsity: 0.373 | Regularization param.: 2.27e-02 | Desired sparsity: 0.300 | Difference: 0.073 | Beta: 2.55e-04 | Relative change: 1.05e-03
Iteration : 205 | Current sparsity: 0.368 | Regularization param.: 2.28e-02 | Desired sparsity: 0.300 | Difference: 0.068 | Beta: 2.55e-04 | Relative change: 1.02e-03
Iteration : 210 | Current sparsity: 0.364 | Regularization param.: 2.29e-02 | Desired sparsity: 0.300 | Difference: 0.064 | Beta: 2.55e-04 | Relative change: 9.95e-04
Iteration : 215 | Current sparsity: 0.359 | Regularization param.: 2.29e-02 | Desired sparsity: 0.300 | Difference: 0.059 | Beta: 2.55e-04 | Relative change: 9.67e-04
Iteration : 220 | Current sparsity: 0.356 | Regularization param.: 2.30e-02 | Desired sparsity: 0.300 | Difference: 0.056 | Beta: 2.55e-04 | Relative change: 9.40e-04
Iteration : 225 | Current sparsity: 0.352 | Regularization param.: 2.31e-02 | Desired sparsity: 0.300 | Difference: 0.052 | Beta: 2.55e-04 | Relative change: 9.13e-04
Iteration : 230 | Current sparsity: 0.348 | Regularization param.: 2.32e-02 | Desired sparsity: 0.300 | Difference: 0.048 | Beta: 2.55e-04 | Relative change: 8.88e-04
Iteration : 235 | Current sparsity: 0.345 | Regularization param.: 2.32e-02 | Desired sparsity: 0.300 | Difference: 0.045 | Beta: 2.55e-04 | Relative change: 8.64e-04
Iteration : 240 | Current sparsity: 0.342 | Regularization param.: 2.33e-02 | Desired sparsity: 0.300 | Difference: 0.042 | Beta: 2.55e-04 | Relative change: 8.40e-04
Iteration : 245 | Current sparsity: 0.340 | Regularization param.: 2.33e-02 | Desired sparsity: 0.300 | Difference: 0.040 | Beta: 2.55e-04 | Relative change: 8.18e-04
Iteration : 250 | Current sparsity: 0.337 | Regularization param.: 2.34e-02 | Desired sparsity: 0.300 | Difference: 0.037 | Beta: 2.55e-04 | Relative change: 7.96e-04
Iteration : 255 | Current sparsity: 0.335 | Regularization param.: 2.34e-02 | Desired sparsity: 0.300 | Difference: 0.035 | Beta: 2.55e-04 | Relative change: 7.74e-04
Iteration : 260 | Current sparsity: 0.332 | Regularization param.: 2.35e-02 | Desired sparsity: 0.300 | Difference: 0.032 | Beta: 2.55e-04 | Relative change: 7.54e-04
Iteration : 265 | Current sparsity: 0.330 | Regularization param.: 2.35e-02 | Desired sparsity: 0.300 | Difference: 0.030 | Beta: 2.55e-04 | Relative change: 7.34e-04
Iteration : 270 | Current sparsity: 0.328 | Regularization param.: 2.35e-02 | Desired sparsity: 0.300 | Difference: 0.028 | Beta: 2.55e-04 | Relative change: 7.15e-04
Iteration : 275 | Current sparsity: 0.326 | Regularization param.: 2.36e-02 | Desired sparsity: 0.300 | Difference: 0.026 | Beta: 2.55e-04 | Relative change: 6.96e-04
Iteration : 280 | Current sparsity: 0.325 | Regularization param.: 2.36e-02 | Desired sparsity: 0.300 | Difference: 0.025 | Beta: 2.55e-04 | Relative change: 6.78e-04
Iteration : 285 | Current sparsity: 0.323 | Regularization param.: 2.36e-02 | Desired sparsity: 0.300 | Difference: 0.023 | Beta: 2.55e-04 | Relative change: 6.61e-04
Iteration : 290 | Current sparsity: 0.322 | Regularization param.: 2.37e-02 | Desired sparsity: 0.300 | Difference: 0.022 | Beta: 2.55e-04 | Relative change: 6.44e-04
Iteration : 295 | Current sparsity: 0.320 | Regularization param.: 2.37e-02 | Desired sparsity: 0.300 | Difference: 0.020 | Beta: 2.55e-04 | Relative change: 6.28e-04
Iteration : 300 | Current sparsity: 0.319 | Regularization param.: 2.37e-02 | Desired sparsity: 0.300 | Difference: 0.019 | Beta: 2.55e-04 | Relative change: 6.12e-04
Iteration : 305 | Current sparsity: 0.317 | Regularization param.: 2.37e-02 | Desired sparsity: 0.300 | Difference: 0.017 | Beta: 2.55e-04 | Relative change: 5.97e-04
Iteration : 310 | Current sparsity: 0.316 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.016 | Beta: 2.55e-04 | Relative change: 5.82e-04
Iteration : 315 | Current sparsity: 0.315 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.015 | Beta: 2.55e-04 | Relative change: 5.68e-04
Iteration : 320 | Current sparsity: 0.314 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.014 | Beta: 2.55e-04 | Relative change: 5.54e-04
Iteration : 325 | Current sparsity: 0.313 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.013 | Beta: 2.55e-04 | Relative change: 5.41e-04
Iteration : 330 | Current sparsity: 0.312 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.012 | Beta: 2.55e-04 | Relative change: 5.28e-04
Iteration : 335 | Current sparsity: 0.311 | Regularization param.: 2.38e-02 | Desired sparsity: 0.300 | Difference: 0.011 | Beta: 2.55e-04 | Relative change: 5.15e-04
Iteration : 340 | Current sparsity: 0.311 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.011 | Beta: 2.55e-04 | Relative change: 5.03e-04
Iteration : 345 | Current sparsity: 0.310 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.010 | Beta: 2.55e-04 | Relative change: 4.91e-04
Iteration : 350 | Current sparsity: 0.309 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.009 | Beta: 2.55e-04 | Relative change: 4.80e-04
Iteration : 355 | Current sparsity: 0.308 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.008 | Beta: 2.55e-04 | Relative change: 4.69e-04
Iteration : 360 | Current sparsity: 0.308 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.008 | Beta: 2.55e-04 | Relative change: 4.58e-04
Iteration : 365 | Current sparsity: 0.307 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.007 | Beta: 2.55e-04 | Relative change: 4.48e-04
Iteration : 370 | Current sparsity: 0.306 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.006 | Beta: 2.55e-04 | Relative change: 4.38e-04
Iteration : 375 | Current sparsity: 0.306 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.006 | Beta: 2.55e-04 | Relative change: 4.28e-04
Iteration : 380 | Current sparsity: 0.305 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.005 | Beta: 2.55e-04 | Relative change: 4.19e-04
Iteration : 385 | Current sparsity: 0.305 | Regularization param.: 2.39e-02 | Desired sparsity: 0.300 | Difference: 0.005 | Beta: 2.55e-04 | Relative change: 4.10e-04
Desired sparsity reached after 385 iterations!

[25]:

x_controlled = fista_controlled.solution
show2D(x_fista, title="Fixed reg. param.", fix_range=(0,0.4), size=(10,20), num_cols=3, slice_list=[(0, 90), (1, 90), (2, 90)]);
show2D(x_controlled, title="Automated reg. param.", fix_range=(0,0.4), size=(10,20), num_cols=3, slice_list=[(0, 90), (1, 90), (2, 90)]);

../../_images/demos_controlledWaveletSparsity_3d_39_0.png

../../_images/demos_controlledWaveletSparsity_3d_39_1.png

We can use islicer to view the automated controlled solution:

[26]:

islicer(x_controlled, minmax=(0,0.4), title=f"Automated parameter choice")

[26]:

[27]:

import matplotlib.pyplot as plt

N = fista_controlled.iterations[-1]
fig, ax = plt.subplots(2)
ax[0].plot(sparsity_ctrl.sparsity_vals, label='Current sparsity')
ax[0].hlines(sparsity_ctrl.desired_sparsity, xmin=0, xmax=N, colors='r', label='Desired sparsity')
ax[0].set_ylabel("Sparsity level")
ax[0].legend()

ax[1].semilogy(sparsity_ctrl.reg_param_vals)
ax[1].set_ylabel("Regularization parameter")
ax[1].set_xlabel("Iteration")
for a in ax.flat:
    a.label_outer()

../../_images/demos_controlledWaveletSparsity_3d_42_0.png

Here we can see how increasing the regularization parameter (lower plot) drives the sparsity level down towards the desired level (top plot). The iteration can be stopped once the suitable level has been reached (and the relative change is small enough).

Notice also how the tuning is adaptive, in the sense that as we get closer to the desired sparsity level, the incremental changes in the regularization parameter get finer.

Limitations of the method#

The method does not fully automate the choice of regularization parameter as the desired sparsity level and some of the tuning parameters still need to be chosen. However, chosing a good level of sparsity is much easier, the values are clearly bounded between 0 and 1 (corresponding to 0% and 100% of coefficients being “meaningful”) and the tuning parameters mostly affect the speed on the method. Moreover, with bit of experience, it is relatively easy to tell early whether the choice will be good and restart the algorithm if needed.

Table of Contents