Module src.fhs

Functions specific to the implementation of the Faceted HyperSARA algorithm [Thouvenin2021].

src.fhs.facetHyperSARA(y, epsilon, A, At, pU, G, W, param, Qx, Qy, K, wavelet, filter_length, nlevel, window_type, spectral_chunk, n_channels, overlap_size, M, N, oy, ox, warmstart_name, checkpoint_name, flag_dimensionality_reduction, Lambda, varargin)

Implementation of the Faceted HyperSARA imaging algorithm.

Implementation of the reweighting / PDFB algorithm underlying the Faceted HyperSARA imaging approach [Thouvenin2021]. At each iteration of the outer reweighting scheme, the PDFB solves a problem of the form

\[\min_{X \in \mathbb{R}_+^{N \times L}} \sum_{q=1}^Q \big( \overline{\mu}_q \| D_q \overline{S}_q X \|_{\overline{\omega}_q, *} + \mu \|\Psi_q^\dagger S_q X \|_{\omega_q, 2,1} \big) + \sum_{\ell, b} \iota_{ \{\| y_{\ell, b} - \Phi_{\ell, b}(\cdot) \|_2 \leq \varepsilon_{\ell, b}\} } (x_\ell).\]
Parameters

  • y (composite, cell) – Blocks of visibilities {L}{nblocks_l}.

  • epsilon (composite, cell) – Data fidelity constraints {L}{nblocks_l}.

  • A (function handle) – Measurement operator.

  • At (function handle) – Adjoint of the measurement operator.

  • pU (composite, cell) – Preconditioning matrices {L}{nblocks_l}.

  • G (composite, cell of sparse matrices) – Degridding matrices {L}{nblocks_l}.

  • W (composite, cell of int[:]) – Masks for selection of the blocks of visibilities.

  • param (struct) – Algorithm parameters.

  • Qx (int) – Number of spatial facets along spatial axis x.

  • Qy (int) – Number of spatial facets along spatial axis y.

  • K (int) – Number of data workers.

  • wavelet (cell of string) – List of wavelet basis (for the SARA prior). If the Dirac basis ('self') is considered, it needs to appear in last position.

  • filter_length (int[:]) – Length of each wavelet filter considered for the SARA dictionary (0 by convention for the Dirac basis).

  • nlevel (int) – Number of wavelet decomposition levels.

  • window_type (string ("triangular" by default, "hamming" or) – "pc" (piecewise-constant)) Type of apodization window considered for the faceted nuclear norm prior.

  • spectral_chunk (cell of int[:]) – List of indices of spectral channels handled by each of the K data workers

  • n_channels (int) – Total number of spectral channels.

  • overlap_size (int[2]) – Overlap size between consecutive facets along each axis (y and x) for the faceted low-rankness prior.

  • M (int) – Spatial image size along axis y.

  • N (int) – Spatial image size along axis x.

  • oy (int) – Fourier oversampling factor along axis y.

  • ox (int) – Fourier oversampling factor along axis x.

  • warmstart_name (string) – Name of a valid .mat file to initialize the solver (warm-start).

  • checkpoint_name (string) – Name defining the name for checkpoint files saved throughout the iterations.

  • flag_dimensionality_reduction (bool) – Flag to indicate data dimensiomality reduction is used.

  • Lambda (composite, cell of complex[:, :]) – Dimensionality reduction matrix. (?)

  • varargin{1} (double[:, :]) – Initial value for the primal variable [M*N, L].

  • varargin{2} (double[:, :]) – Ground truth wideband image (only when debugging synthetic data experiments) [M*N, L].

Returns

xsol – Reconstructed wideband image [M*N, L].

Return type

double[:, :]

Important

  • The param structure should contain the following fields.

param.verbose (string)

Print log or not

param.save_intermediate_results_mat (bool)

Flag to save all estimates in a .mat file after each re-weighted problem. The file can be used as a warmstart.

param.nu0 (double)

Norm of the splitting operator used to defined the faceted low-rankness dual variable (identity for rectangular window, thus fixed to 1).

param.nu1 (double)

Upper bound on the norm of the SARA operator \(\Psi^\dagger\).

param.nu2 (double)

Upper bound on the norm of the measurement operator \(\Phi\)

param.gamma (double)

Regularization parameter (sparsity prior).

param.gamma0 (double[:)

Regularization parameter (faceted low-rankness prior).

param.pdfb_min_iter (int)

Minimum number of iterations (PDFB).

param.pdfb_max_iter (int)

Maximum number of iterations (PDFB).

param.pdfb_rel_var (double)

Relative variation tolerance (PDFB).

param.pdfb_fidelity_tolerance (double)

Tolerance to check data constraints are satisfied (PDFB).

param.reweighting_max_iter (int)

Maximum number of reweighting steps.

param.reweighting_min_iter (int)

Minimum number of reweighting steps (to reach “noise” level).

param.reweighting_rel_var (double)

Tolerance relative variation (reweighting).

param.reweighting_alpha (double)

Starting reweighting parameter.

param.reweighting_sig (double)

Noise level (in wavelet space)

param.reweighting_sig_bar (double)

Noise level (singular value space)

param.elipse_proj_max_iter (int)

Maximum number of iterations for the FB algo that implements the projection onto the ellipsoid.

param.elipse_proj_min_iter (int)

Minimum number of iterations for the FB algo that implements the projection onto the ellipsoid.

parma.elipse_proj_eps (double)

Stopping criterion for the projection.

  • The following fileds are added to param in the course of the algorithm to be able to restart it from a previous state

param.reweighting_alpha (double)

Current state of the reweighting parameter.

param.init_reweight_step_count (int)

Iteration index of the current reweighting step when the checkpoint file is written.

param.init_reweight_last_iter_step (int)

Global iteration index (unrolled over all pdfb iterations performed in the reweighting scheme) from which to restart the algorithm.

param.init_t_start (int)

Global iteration index (unrolled over all pdfb iterations performed in the reweighting scheme) from which to restart the algorithm (param.init_t_start = param.init_reweight_last_iter_step + 1).

  • The following fileds are involved in the adaptive estimation of the noise level

param.adjust_flag_noise (bool)

Flag to activate the adaptive procedure to estimate the noise level.

param.adjust_noise_min_iter (int)

Minimum number of iterations to enable the adjustement of the estimate of the noise level.

param.adjust_noise_rel_var (double)

Tolerance relative variation (reweighting) to to enable the adjustement of the estimate of the noise level.

param.adjust_noise_start_iter (int)

Number of iterations to force triggering noise adjustement.

  • The following fileds are added to param in the course of the algorithm to be able to restart it from a previous state

param.reweighting_alpha (double)

Current state of the reweighting parameter.

param.init_reweight_step_count (int)

Iteration index of the current reweighting step when the checkpoint file is written.

param.init_reweight_last_iter_step (int)

Global iteration index (unrolled over all pdfb iterations performed in the reweighting scheme) from which to restart the algorithm.

param.init_t_start (int)

Global iteration index (unrolled over all pdfb iterations performed in the reweighting scheme) from which to restart the algorithm (param.init_t_start = param.init_reweight_last_iter_step + 1).

Note

The checkpoint file saved throughout the iterations is composed of the following variables (to be verified).

param (struct)

Current state of the parameter structure (updated with auxiliary parameters describing the state of the solver when the checkpoint file has been written).

xsol (double[:, :])

Current state of the image estimate [M*N, L].

res (double[:, :, :])

Current state of the wideband residual image.

g (double[:, :])

Auxiliary variable involved in the update of the primal variable.

epsilon (cell of cell of double)

(Updated) Value of the l2-ball radii.

t_block (cell of cell of int)

Index of the last iteration where the weigths have been updated.

proj (cell of cell of complex[:])

Auxiliary variable involved in the update of the data fidelity dual variables.

norm_res (cell of cell of double)

Norm of the residual image (per channel per block).

v0 (cell of double[:])

Dual variables associated with the low-rankness prior {Q}.

v1 (cell of double[:, :])

Dual variables associated with the sparsity prior {Q}.

v2 (cell of cell complex[:])

Dual variables associated with the data fidelity terms (each cell corresponding to a data block) {K}.

weights0 (cell of double[:])

Weights associated with the faceted low-rankness prior {Q}.

weights1 (cell of double[:])

Weights associated with the sparsity prior {Q}.

end_iter (int)

Last iteration (unrolled over all pdfb iterations).

t_facet (double)

Time to update all the variables assoacited with a spatial facet..

t_data (double)

Time to update the data fidelity dual variables.

rel_val (double)

Relative variation of the solution across the iterations.

src.fhs.fhs_compute_facet_prior(x_facet, Iq, offset, status_q, nlevel, wavelet, Ncoefs_q, dims_overlap_ref_q, offsetLq, offsetRq, crop_sparsity, crop_low_rank, spatial_weights, size_v1)

Compute the value of the faceted prior (\(\ell_{2,1}\) + nuclear norm).

Compute the value of the faceted prior (\(\ell_{2,1}\) and nuclear norm) for a single facet. This version includes a spatial correction of the faceted nuclear norm (tapering window).

Parameters

  • x_facet (double[:, :, :]) – Overlapping image facet [M, N, L].

  • Iq (int[:]) – Starting index of the non-overlapping facet [1, 2].

  • offset (int[:]) – offset to be used from one dictionary to another (different overlap needed for each dictionary -> cropping) {nDictionaries}.

  • status_q (int[:]) – Status of the current facet (last or first facet along vert. / hrz. direction).

  • nlevel (int) – Depth of the wavelet decompositions.

  • wavelet (cell (strings)) – Name of the wavelet dictionaries.

  • Ncoefs_q (int[:, :]) – Size of the wavelet decompositions at each scale.

  • dims_overlap_ref_q (int[:]) – Dimensions of the facet [1, 2].

  • offsetLq (int[:]) – Amount of zero-padding from the “left” [1, 2].

  • offsetRq (int[:]) – Amount of zero-padding from the “right” [1, 2].

  • crop_sparsity (int[:]) – Relative cropping necessary for the facet \(\ell_{2,1}\) norm [1, 2].

  • crop_low_rank (int[:]) – Relative cropping necessary for the facet nuclear norm [1, 2].

  • spatial_weights (double[:, :]) – Spatial weights (apodization window) applied to the facet nuclear norm.

  • size_v1 (int[:]) – Size of the dual variable associated with the facet \(\ell_{2,1}\) norm [1, 2].

Returns

  • l21_norm (double) – Facet \(\ell_{2,1}\) norm.

  • nuclear_norm (double) – Facet nuclear norm.

  • ..note:: – The l21 and nuclear norms do not act on the same facet, hence the cropping step.

src.fhs.fhs_initialize_data_worker(y)

Initialize the variables stored on the data workers.

Parameters

y (cell) – Input data, stored as 2-layer cell, the first for channels, the second for data blocks within each channel.

Returns

  • v2 (cell) – Data fidelity dual variable. Same structure as y.

  • norm_res (cell) – Norm of the residual. Same structure as y.

  • t_block (cell) – Last iteration at which each block has been updated {L}{nblocks}. Same structure as y.

  • proj (cell) – Result of the projection step {L}{nblocks}[M, 1]. Same structure as y.

src.fhs.fhs_initialize_dual_and_weights(x_facet, I, offset, status, nlevel, wavelet, Ncoefs, dims_o, n_channels, dims_overlap_ref, offsetL, offsetR, reweight_alpha, crop_sparsity, crop_low_rank, apodization_window, sig, sig_bar)

Initialize dual variables and weights associated with the faceted low-rankness and average joint-sparsity prior.

Parameters

  • x_facet (double[:, :, :]) – Overlapping image facet [M, N, L].

  • I (int[:]) – Starting index of the non-overlapping facet [1, 2].

  • offset (int[:]) – Offset to be used from one dictionary to another (different overlap needed for each dictionary -> cropping) {nDictionaries}.

  • status (int[:]) – Status of the current facet (last or first facet along vert. / hrz. direction) [1, 2].

  • nlevel (int) – Depth of the wavelet decompositions.

  • wavelet (cell, string) – Name of the wavelet transforms involved in the faceted average joint-sparsity prior.

  • Ncoefs (int[:]) – Size of the wavelet decompositions at each scale.

  • dims_o (int[:]) – Dimension of a facet (with overlap) [1, 2].

  • n_channels (int) – Number of spectral channels.

  • dims_overlap_ref (int[:]) – Dimension of the facet [1, 2].

  • offsetL (int[:]) – Amount of zero-pading from the “left” [1, 2].

  • offsetR (int[:]) – Amount of zero-padding from the “right” [1, 2].

  • reweight_alpha (double) – Reweighting parameter.

  • crop_sparsity (int[:]) – Amount of pixels to be cropped from the facet along each dimension to retrieve the pixels over which the current facet’s sparstiy prior is acting [1, 2].

  • crop_low_rank (int[:]) – Amount of pixels to be cropped from the facet along each dimension to retrieve the pixels over which the current facet’s sparstiy prior is acting [1, 2].

  • apodization_window (double[:, :]) – Apodization window used in the faceted low-rankness prior.

  • sig (double) – Noise level for the weights (joint-sparsity prior).

  • sig_bar (double) – Noise level for the weights (low-rankness prior).

Returns

  • v0 (double[:, :]) – Dual variable associated with the low-rankness prior.

  • v1 (double[:, :]) – Dual variable associated with the average joint-sparsity prior.

  • weights0 (double[:]) – Weigths associated with the low-rankness prior.

  • weights1 (double[:]) – Weigths associated with the average joint-sparsity prior.

src.fhs.fhs_initialize_facet_dual(Ncoefs, dims_o, n_channels, nlevel)

Initialize dual variables (constant overlap).

Initialize all the dual variables for a given facet (nuclear and \(\ell_{2, 1}\) norms).

Parameters

  • Ncoefs (int[:, :]) – Number of wavelet coefficients at each scale.

  • dims_o (int[:]) – Dimension of a facet (with overlap) [1, 2].

  • n_channels (int) – Number of spectral channels.

  • nlevel (int) – Depth of the wavelet decompositions.

Returns

  • v0 (double[:, :]) – Dual variable associated with the nuclear norm.

  • v1 (double[:, :]) – Dual variable associated with the \(\ell_{2,1}\) norm.

  • weights0 (double[:]) – Weigths associated with the nuclear norm.

  • weights1 (double[:]) – Weigths ssociated with the \(\ell_{2,1}\) norm.

src.fhs.fhs_setup_priors(Qx, Qy, I, dims, dims_o, dims_overlap_ref, I_overlap, dims_overlap, status, offsetL, offsetR, Ncoefs, temLIdxs, temRIdxs, window_type, overlap_size)

Setup all the composite arrays for the faceted low-rankness and average joint sparsity priors involed in Faceted HyperSARA.

Parameters

  • Qx (int) – Number of facets (axis x).

  • Qy (int) – Number of facets (axis y).

  • I (int[:]) – Starting index of the non-overlapping tile [1, 2].

  • dims (int[:]) – Size of the non-overlapping tile [1, 2].

  • dims_o (int[:]) – Dimension of a facet (with overlap) [1, 2].

  • dims_overlap_ref (int[:]) – Dimension of the facet [1, 2].

  • I_overlap (int[:]) – Starting index of the overlapping facets [Q, 1].

  • dims_overlap (int[:, :]) – Size of the overlapping facets [Q, 2].

  • status (int[:]) – Status of the current facet (last or first facet along vert. / hrz. direction) [1, 2].

  • offsetL (int[:, :]) – Amount of zero-pading from the “left” [Q, 2].

  • offsetR (int[:, :]) – Amount of zero-padding from the “right” [Q, 2].

  • Ncoefs (cell of int[:, :]) – Size of the wavelet decomposition across each scale.

  • temLIdxs (int[:, :]) – Amount of cropping from the “left” [Q, 2].

  • temRIdxs (int[:, :]) – Amount of cropping from the “right” [Q, 2].

  • window_type (string) – Type of apodization window considered for the faceted low-rankness prior.

  • overlap_size (int[:]) – Number of overlapping pixels in the current facet along each direction [1, 2].

Returns

  • Iq (int[:]) – Starting index of the non-overlapping base facet [1, 2].

  • dims_q (int[:]) – Dimensions of the non-overlapping base facet [1, 2].

  • dims_oq (int[:]) – Diemnsions of the overlapping facet [1, 2].

  • dims_overlap_ref_q (int[:]) – Dimension of the facet [1, 2].

  • I_overlap_q (int[:]) – Starting index of the facet [1, 2].

  • dims_overlap_q (int[:]) – Dimensions of the facet [1, 2].

  • status_q (int[:]) – Status of the current facet (last or first facet along vert. or hrz. direction) [ndict, 2].

  • offsetLq (int[:]) – Amount of zero-pading from the “left” [1, 2].

  • offsetRq (int[:]) – Amount of zero-pading from the “right” [1, 2].

  • Ncoefs_q (int[:, :]) – Size of the wavelet decompositions at each scale.

  • temLIdxs_q (int[:]) – Amount of cropping from the “left” [1, 2].

  • temRIdxs_q (int[:]) – Amount of cropping from the “right” [1, 2].

  • overlap_g_south (int[:]) – Size of the overlap for the south neighbour [2,1].

  • overlap_g_east (int[:]) – Size of the overlap for the east neighbour [2,1].

  • overlap_g_south_east (int[:]) – Size of the overlap for the south-east neighbour [2,1].

  • overlap (int[:]) – Number of overlapping pixels along each direction.

  • w (double[:, :]) – Apodization window considered for the faceted low-rankness prior.

  • crop_low_rank (int[:]) – [Relative cropping necessary for the faceted low-rankness prior [1, 2].

  • crop_sparsity (int[:]) – Relative cropping necessary for the faceted joint-sparsity prior [1, 2].

src.fhs.fhs_update_dual_lowrankness(v0, xhat, apodization_window, weights0, beta0)

Update the dual variable related to the weighted nuclear norm prior.

Update step for the facet dual variable related to the nuclear norm prior over a given facet. This version includes a spatial correction for the faceted nuclear norm (tapering window).

Parameters

  • v0 (double[:]) – Dual variable associated with the nuclear norm prior [min(M*N, L), 1].

  • xhat (double[:, :, :]) – Auxiliary variable related to the wideband image [M, N, L].

  • apodization_window (double[:, :]) – Spatial weights applied to the facet nuclear norm [M, N, L].

  • weights0 (double[:]) – Weights for the reweighting [min(M*N, L), 1].

  • beta0 (double) – Thresholding parameter (gamma0/sigma0).

Returns

  • v0 (double[:, :]) – Dual variable associated with the nuclear norm prior [min(M*N, L), 1].

  • g0 (double[:, :, :]) – Auxiliary variable for the update of the primal variable [M, N, L].

src.fhs.fhs_update_dual_sparsity(v1, x_facet, apodization_window, beta1, Iq, dims_q, I_overlap_q, dims_overlap_q, offset, status_q, nlevel, wavelet, Ncoefs_q, temLIdxs_q, temRIdxs_q, offsetLq, offsetRq, dims_overlap_ref_q)

Update the dual variable associated with the facet l21-norm prior.

Update a facet dual variable associated with the l21-norm prior.

Parameters

  • v1 (double[:, :]) – Dual variable associated with the \(\ell_{2,1}\) prior [s, L].

  • x_facet (double[:, :, :]) – Overlapping image facet [M, N, L].

  • apodization_window (double[:, :]) – Weights to mitigate tessellation aretefacts due to the overlap between facets [M, N].

  • beta1 (double) – Ratio between regularization and convergence parameter (gamma1 / sigma1).

  • Iq (int[:]) – Starting index of the non-overlapping base facet [1, 2].

  • dims_q (int[:]) – Dimensions of the non-overlapping base facet [1, 2].

  • I_overlap_q (int[:]) – Starting index of the facet [1, 2].

  • dims_overlap_q (int[:]) – Dimensions of the facet [1, 2].

  • offset (cell of int[:]) – Offset to be used from one dictionary to another (different overlap needed for each dictionary -> cropping) {nDictionaries}.

  • status_q (int[:]) – Status of the current facet (last or first facet along vert. or hrz. direction) [ndict, 2].

  • nlevel (int) – Depth of the wavelet decompositions.

  • wavelet (cell (string)) – Name of the wavelet dictionaries {ndict}.

  • Ncoefs_q (int[:]) – Size of the wavelet decompositions at each scale.

  • temLIdxs_q (int[:]) – Amount of cropping from the “left” [1, 2].

  • temRIdxs_q (int[:]) – Amount of cropping from the “right” [1, 2].

  • offsetLq (int[:]) – Amount of zero-pading from the “left” [1, 2].

  • offsetRq (int[:]) – Amount of zero-pading from the “right” [1, 2].

  • dims_overlap_ref_q (int[:]) – Dimension of the facet [1, 2].

Returns

  • v1 (double[:, :]) – Dual variable associated with the \(\ell_{2,1}\) prior [s, L].

  • g (double[:, :, :]) – Auxiliary variable for the update of the primal variable [M, N, L].

src.fhs.fhs_update_primal(xsol, g)

Update the wideband facet (primal variable).

Update step for the primal variable in the preconditioned primal-dual algorithm.

Parameters

  • xsol (double[:, :, :]) – Wideband facet [M, N, L].

  • g (double[:, :, :]) – Auxiliary variable related to the wideband facet [M, N, L].

Returns

  • xsol (double[:, :, :]) – Udated wideband facet [M, N, L].

  • xhat (double[:, :, :]) – Auxiliary variable related to the wideband facet [M, N, L].

  • rel_x (double) – Squared relative variation between two consecutive iterates.

  • norm_x (double) – Squared Euclidean norm of xsol.

src.fhs.fhs_update_weights(x_facet, size_v1, I, offset, status, nlevel, wavelet, Ncoefs, dims_overlap_ref, offsetL, offsetR, reweight_alpha, crop_sparsity, crop_low_rank, apodization_window, sig, sig_bar)

Update the weights of the per facet priors.

Update the weights involved in the reweighting procedure applied to the faceted l21 and nuclear norms. This version includes a spatial weighting correction for the nuclear norm (tapering window).

Parameters

  • x_facet (double[:, :, :]) – Overlapping image facet [M, N, L].

  • size_v1 (int[:]) – Size of the dual variable associated with the facet \(\ell_{2,1}\) norm [1, 2].

  • I (int[:]) – Starting index of the non-overlapping facet [1, 2].

  • offset (cell of int[:]) – Offset to be used from one dictionary to another (different overlap needed for each dictionary -> cropping) {nDictionaries}.

  • status (int[:]) – Status of the current facet (last or first facet along vert. / hrz. direction) [1, 2].

  • nlevel (int) – Depth of the wavelet decompositions.

  • wavelet (cell of string) – Name of the wavelet dictionaries.

  • Ncoefs (int[:]) – Size of the wavelet decompositions at each scale.

  • dims_overlap_ref (int[:]) – Dimension of the facet [1, 2].

  • offsetL (int[:]) – Amount of zero-pading from the “left” [1, 2].

  • offsetR (int[:]) – Amount of zero-padding from the “right” [1, 2].

  • reweight_alpha (double) – Reweighting parameter.

  • crop_sparsity (int[:]) – Relative cropping necessary for the facet l21-norm [1, 2].

  • crop_low_rank (int[:]) – Relative cropping necessary for the facet nuclear norm [1, 2].

  • apodization_window (double[:, :]) – Spatial weights applied to the facet nuclear norm (same size as x_facet after cropping by crop_low_rank).

  • sig (double) – Noise level (wavelet space).

  • sig_bar (double) – Noise level (singular value space).

Returns

  • weights1 (double[:, :]) – Weights for the reweighting of the facet \(\ell_{2,1}\)-norm [size_v1(1), size_v1(2)].

  • weights0 (double[:, :]) – Weights for the reweighting of the nuclear norm [min(M*N, L), 1].