MIRAG.optim package

MIRAG.optim.sparse_coding

Problem solving module for GPR using the ADMM, for solving the following problem:

\(\left\{\begin{array}{ll}\mathbf{minimize}\quad &\frac{1}{2}||\sum_{k}\mathbf{H}_k\star \mathbf{C}_k - \mathbf{Y} ||_2^2 +\epsilon \sum_k{||\mathbf{S}_k||_1} \\ \text{t.q.} \quad&\mathbf{C}_k =\mathbf{S}_k\end{array}\right.\)

Variant with optimization by Huber norm instead of norm 2

References

4(1,2)

‘Convolutional dictionary learning: A comparative review and new algorithms’,Garcia-Cardona, C., & Wohlberg, B. (2018). Avalaible at: https://arxiv.org/pdf/1709.02893.pdf

Warning

To do

• Complete class documentation (for all attributes)

• Add logging for tracking

• Check dictionary strucure such as: \(D^H= \mathrm{conjugate}(D)\)

• Add over-relaxation parameter (see source_separation and 3)

MIRAG.optim.sparse_coding.addm_iteration_norm2(lambdaS, c, DHY, DF, DF_H, U, rho, S, penalty='l1', m=- 1)

function for computing an iteration of ADMM with the 2 norm Based on 4

Parameters

• lambdaS (float) – parsimony parameter

• c (float) – pre-computation of the term for Sherman_MorrisonF

• DHY (complex) – fft from dictionary * original image

• DF (complex) – fft of the dictionary

• DF_H (complex) – conjugate of the dictionary fft

• U (float) – dual variable

• rho (float) – penalty parameter

• S (float) – auxiliary variable \(\mathbf{S}_k\)

• penalty (str{"l1", "l0", "FirmThresholding", "l*"}, optional) – data attachment penalty, basic \(\sum_k{||\mathbf{S}_k||_1\)

• m (int{-1}, optional) – Number of workers (cores) used for the fft

Returns

• primal_new (float) – variable primal \(\mathbf{C}_k^i\) (Nx,Ny,K)

• auxiliary_new (float) – dual variable \(\mathbf{S}_k^i\) (Nx,Ny,K)

• Dal (float) – convolution product dictionary + map coeff \(\sum_k{\mathbf{C}_k^i\star\mathbf{H}_k^i}\) (Nx,Ny,1)

MIRAG.optim.sparse_coding.addm_iteration_normH(lambdaS, DF, DFH, U, rho, S, primal_tilde, YF, grad_iter_max=50, beta=0.001, thresh=250, penalty='l1', m=- 1)

function for computing an iteration of ADMM with the HUber norm

Parameters

• lambdaS (float) – parsimony parameter

• DF (complex) – fft of the dictionary (Nx,Ny,K)

• DF_H (complex) – conjugate of the dictionary fft

• U (float) – dual variable (Nx,Ny,K)

• rho (float) – penalty parameter

• S (float) – auxiliary variable \(\mathbf{S}_k\) (Nx,Ny,K)

• primal_tilde (float) – fft of the primal \(\mathbf{C}_k\) (Nx,Ny,K)

• YF (float) – fft of the original image (Nx,Ny,1)

• grad_iter_max{50} (int, optional) – maximum iteration of the gradient descent

• beta (float{0.001}, optional) – learning rate parameter of the gradient

• thresh (int, optional) – [description], by default 250

• penalty (str{"l1", "l0", "FirmThresholding", "l*"}, optional) – penalty of attachment to data, basic \(\sum_k{||\mathbf{S}_k||_1\)

• m (int{-1}, optional) – Number of workers (cores) used for the fft

Returns

• primal_new (float) – variable primal \(\mathbf{C}_k^i\) (Nx,Ny,K)

• primal_tilde (float) – fft of the primal variable \(\mathbf{C}_k^i\) (Nx,Ny,K)

• auxiliary_new (float) – dual variable \(\mathbf{S}_k^i\) (Nx,Ny,K)

• Dal (float) – convolution product dictionary + map coeff \(\sum_k{\mathbf{C}_k^i\star\mathbf{H}_k^i}\) (Nx,Ny,1)

MIRAG.optim.sparse_coding.cost_function_addm1(Y, rho, old_auxiliary_, old_dal, dal, primal, auxiliary)

ADMM cost and error calculation function

Parameters

• Y (float) – original image (Nx,Ny,1)

• rho (float) – penalty parameter

• old_auxiliary (float) – auxiliary variable of the previous iteration \(\mathbf{S}_k^{i-1}\)

• old_dal (float) – reconstruction of the previous iteration \(\sum_k{\mathbf{C}_k^{i-1}\star\mathbf{H}_k}^{i-1}\)

• dal (float) – convolution product dictionary + coeff map \(\sum_k{\mathbf{C}_k^i\star\mathbf{H}_k^i}\)

• primal (float) – variable primal \(\mathbf{C}_k^i\)

• auxiliary (float) – auxiliary variable \(\mathbf{S}_k^{i-1}\)

Returns

out – variation of primal, reconstruction error, error of primal and error of dual

Return type

list of float

MIRAG.optim.source_separation

Problem solving module for GPR using the ADMM, for solving the following problem:

\(\left\{\begin{array}{ll}\mathbf{minimize} \quad &||\mathbf{L}||_* + \epsilon \sum_k{||\mathbf{S}_k||_1}\\ \mathbf{t.q.}\quad& \mathbf{Y} = \mathbf{L} +\sum_{k}\mathbf{H}_k\star \mathbf{C}_k \\ &\quad \mathbf{S}_k=\mathbf{C}_k \end{array}\right.\)

References

5

‘Sparse Decomposition of the GPR Useful Signal from Hyperbola Dictionary’, Guillaume Terasse, Jean-Marie Nicolas, Emmanuel Trouvé and Émeline Drouet Avalaible at: https://hal.archives-ouvertes.fr/hal-01351242

Warning

To do

• Add logging for tracking

• Add stop condition on the variation of parameters in addition to the number of iterations

MIRAG.optim.source_separation.addm_iteration_norm2(lambdaS, c, Y, Dal, DF, DF_H, Us, Uy, rhoS, rhoL, S, dim1, over_relax, penalty='l1', m=- 1)

function for computing an iteration of ADMM with norm 2 and the hollow matrix Based on 4 and 5

Parameters

• lambdaS (float) – parsimony parameter

• c (float) – pre-calculation of the term for Sherman_MorrisonF

• Y (float) – original image

• Dal (float) – convolution product dictionary + coeff map (Nx,Ny,1)

• DF (complex) – fft of the dictionary (Nx,Ny,K)

• DF_H (complex) – conjugate of the dictionary fft (Nx,Ny,K)

• Us (float) – dual variable of S

• Uy (float) – dual variable of L

• rhoS (float) – penalty parameter on S

• rhoL (float) – penalty parameter on L

• S (float) – auxiliary variable \(\mathbf{S}_k\)

• dim1 (int) – dimensions of the image (ex: [256,256,1])

• over_relax (float) – over-relaxation parameter (improves convergence for \(\alpha\sim 1.6\))

• penalty (str{"l1", "l0", "FirmThresholding"}, optional) – data attachment penalty, basic \(\sum_k{||\mathbf{S}_k||_1\)

• m (int{-1}, optional) – Number of workers (cores) used for the fft

Returns

• primal_new (float) – variable primal \(\mathbf{C}_k^i\) (Nx,Ny,K)

• auxiliary_new (float) – dual variable \(\mathbf{S}_k^i\) (Nx,Ny,K)

• L (float) – hollow matrix (Nx,Ny,1)

• Dal (float) – convolution product dictionary + map coeff ?:math:sum_k{mathbf{C}_k^istarmathbf{H}_k^i}? (Nx,Ny,1)

MIRAG.optim.source_separation.cost_function_addm2(Y, rhoS, rhoL, old_auxiliary, old_dal, old_L, L, dal, primal, auxiliary)

ADMM cost and error calculation function with 2 constraints.

Parameters

• Y (float) – original image (Nx,Ny,1)

• rhoS (float) – penalty parameter on variable S

• rhoL (float) – penalty parameter on variable L

• old_auxiliary (float) – auxiliary variable of the previous iteration \(\mathbf{S}_k^{i-1}\)

• old_dal (float) – reconstruction of the previous iteration \(\sum_k{\mathbf{C}_k^{i-1}\star\mathbf{H}_k}^{i-1}\)

• old_L (float) – variable of the hollow matrix of the previous iteration \(\mathbf{L}^{i-1}\)

• L (float) – variable of the hollow matrix \(\mathbf{L}\)

• dal (float) – reconstruction

• primal (float) – primal variable \(\mathbf{C}_k^i\)

• auxiliary (float) – auxiliary variable \(\mathbf{S}_k^{i-1}\)

Returns

• var_prim_ (float) – variation of the primal

• error_rec_ (float) – reconstruction error

• error_primal_ (float) – primal error

• error_dual_S (float) – error of dual of S \(\mathbf{U_S}\)

• error_dual_L (float) – error from dual of L \(\mathbf{U_L}\)

MIRAG.optim.huber_source_separation

MIRAG.optim.huber_source_separation.ADMMSourceSepHUB

alias of <MagicMock spec=’str’ id=’140306494904816’>

MIRAG.optim.admm_func

Complementary module for the ADMM convolutional

References

1(1,2,3)

‘Multichannel sparse recovery of complex-valued signals using Huber’s criterion’,Esa Ollila Avalaible at: https://arxiv.org/pdf/1504.04184.pdf

2

‘Robust Principal Component Analysis?’, Candes & al. Avalaible at: http://www.columbia.edu/~jw2966/papers/CLMW11-JACM.pdf

3(1,2)

‘Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers p.23’, Stephen Boyd Avalaible at: https://web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf

MIRAG.optim.admm_func.SVD_gpr(ref, rank)

Perform a SVD on the reference image and dump the first n rank singular values

Parameters

• ref (ndarray) – reference image

• rank (int) – rank of the SVD

Returns

A_remake – Reconstructed reference image without the rank-n singular values

Return type

ndarray

MIRAG.optim.admm_func.Sherman_MorrisonF(DF_H, b, c, rho)

Solve a diagonal block linear system with a scaled identity term using the Sherman-Morrison equation

The solution is obtained by independently solving a set of linear systems of the form (see wohlberg-2015-efficient)

\((a\cdot a^H +\rho I)x = b\)

In this equation inner products and matrix products are taken along the 3rd dimension of the corresponding multi-dimensional arrays; the solutions are independent over the 1st and 2nd (and 4th, if non-singleton) dimensions.

Parameters

• DF_H (complex) – Conjugate of Multi-dimensional array containing \(a^H\)

• b (complex) – Multi-dimensional array containing b

• c (complex) – Multi-dimensional array containing pre-computed quantities \(a^H/(a^H\cdot a +\\rho)\)

• rho (float) – Scalar rho

Returns

x – Multi-dimensional array containing linear system solution

Return type

complex

Notes

Adapted from matlab code : Brendt Wohlberg <brendt@lanl.gov> Modified: 2014-12-18

MIRAG.optim.admm_func.diag_thresh(m, n, s)

Return a diagonal matrix of dimension mxn from a vector s Use after np.linalg.svd to get a matrix from S

Parameters

• m (int) – vertical dimension of the desired matrix

• n (int) – horizontal dimension of the desired matrix

• s (float) – vector to diagonalize

Returns

out – Diagonal matrix (from S)

Return type

float

MIRAG.optim.admm_func.gradient_huber(S_tilde, X_m_tilde, D_m_tilde, D_m_H_tilde, Z_m_tilde, rho, delta)

Function to compute the gradient of the Huber function of the following minimization problem:

\(\underset{\mathbf{\hat{c}}}{\operatorname{argmin}} \sum_N\mathcal{H}_\delta(\mathbf{\hat{D}} \cdot \mathbf{\hat{c}}-\mathbf{\hat{y}})+ \frac{\rho}{2}\Big|\Big| {\mathbf{\hat{c}}-\mathbf{\hat{z}}}\Big|\Big|_2^2\)

The gradient is computed as follows:

\(\nabla_\mathbf{\hat{c}}f=\{\mathbf{\tilde{DH}_m}\}_k\odot\Psi_\delta \left(\sum_{k}\mathbf{\tilde{D}_m}_k\odot\mathbf{\tilde{X}_m}_k -\mathbf{\tilde{S}}\right)+\rho({\mathbf{\tilde{X}_m}-\mathbf{\tilde{Z}_m}})\)

Parameters

• S_tilde (complex) – fft of the original signal (Nx,Ny)

• X_m_tilde (complex) – fft of the coefficients maps \(\mathbf{\hat{C}_k}\) (Nx,Ny,K)

• D_m_tilde (complex) – fft of the dictionary \(\mathbf{\hat{D}_k}\) (Nx,Ny,K)

• D_m_H_tilde (complex) – fft hermitian of the dictionary \(\mathbf{\hat{D}_k}^H\) (Nx,Ny,K)

• Z_m_tilde (complex) – fft (auxiliary variable- dual variable) (Nx,Ny,K)

• rho (float) – regularization parameter

• delta (float) – threshold parameter

Returns

out – gradient of the Huber function (NX,NY,K)

Return type

complex

MIRAG.optim.admm_func.huber_complex(x, delta)

Huber function applied to a complex array as proposed in 1

\(\rho_{H, \delta}(x)=\left\{\begin{array}{ll}|x|^{2}, & \text { for }|x| \leq \delta \\2 \delta|x|-\delta^{2}, & \text { for }|x|>\delta\end{array}\right.\)

Parameters

• x (complex) – array to apply the Huber function

• delta (float) – threshold parameter

Returns

out – array with the Huber function applied

Return type

complex

MIRAG.optim.admm_func.loss_derivative_function_complex(x, delta)

Derivative of the Huber function, based on 1

\(\psi_{H, \delta}(x)=\left\{\begin{array}{ll} x, & \text { for }|x| \leq \delta \\ \delta \operatorname{sign}(x), & \text { for }|x|>\delta \end{array}\right.\)

Parameters

• x (complex) – array to apply the Huber function

• delta (float) – threshold parameter

Returns

out – array with the derivative of the Huber function applied

Return type

complex

MIRAG.optim.admm_func.proxH(x, seuil, rho)

proximal operator of the Huber function

Parameters

• x (ndarray) – input signal

• seuil (float) – threshold

• rho (float) – regularization parameter

Returns

result of the proximal operator

Return type

ndarray

MIRAG.optim.admm_func.roll_fft(alpha, t, x)

Correction of the position of the coefficients of the maps according to the of the central position of the hyperbolas + Summation of the C_k

Parameters

• alpha (float) – Complex tensor (M x N x K) of the coefficients maps

• t (int) – central position (pixel) of the used hyperbolas (ordinate).

• x (int) – central position (pixel) of the used hyperbolas (abscissa).

Returns

out – Corrected and reduced matrix (M x N)

Return type

float

MIRAG.optim.admm_func.sign_complex_array(x)

Determine the sign of a complex array, based on 1

\(\operatorname{sign}(e)=\left\{\begin{array}{ll} e /|e|, & \text { for } e \neq 0 \\ 0, & \text { for } e=0 \end{array}\right.\)

Parameters

x (complex) – array to determine the sign

Returns

out – array with the sign of the input array

Return type

int

MIRAG.optim.admm_func.thresholding_nuclear(M, i, lambdaS, rho)

Singular value thresholding function It is the resultant of the proximal operator associated to the nuclear norm.

\(\underset{\mathbf{{L}}}{\operatorname{argmin}} \lambda||\mathbf{L}||_* + \frac{\rho_L}{2}\Big|\Big|\mathbf{X}-\mathbf{L}\Big|\Big|_2^2 = \mathrm{prox}_{||.||_*,\lambda/\rho_L}(\mathbf{X}) \\ \ \mathrm{avec}\ \ \mathrm{prox}_{||.||_*,\lambda/\rho_L}(x)= \mathcal{T}_{\lambda/\rho_L}\left(x\right)\)

Based on 2

Parameters

• M (float) – tensor of dimension Nx,Ny,K

• i (int) – i th layer to threshold (0<i<K)

• lambdaS (float) – Parsimony parameter (if existing >0)

• rho (float) – Penalty parameter >0

Returns

L – Result of the minimization for a layer i

Return type

float

MIRAG.optim.admm_func.update_rhoLS_adp(er_prim, er_dual, rhoLS)

Adaptation function of the penalty parameters in the the ADMM from the errors made on the primal and dual. Based on 3

Parameters

• er_prim (float) – error of the primal

• er_dual (float) – error of the dual

• rhoLS (float) – value of the parameter to be updated

Returns

• rhoLS (float) – updated penalty parameter

• k (float) – expansion/contraction factor used