Momentum-Net: Fast and convergent iterative neural network for inverse problems
Il Yong Chun Zhengyu Huang Hongki Lim Jeffrey A. Fessler
eess.IV cs.CV cs.LG math.OC









Iterative neural networks (INN) are rapidly gaining attention 
for solving inverse problems in imaging, image processing, and computer vision.
INNs combine regression NNs and an iterative model-based image reconstruction (MBIR) algorithm, 
often leading to both good generalization capability and 
outperforming reconstruction quality over existing MBIR optimization models.
This paper proposes the first fast and convergent INN architecture, Momentum-Net,
by generalizing a block-wise MBIR algorithm that uses momentum and majorizers with regression NNs.
For fast MBIR, Momentum-Net uses momentum terms in extrapolation modules,
and noniterative MBIR modules at each iteration by using majorizers,
where each iteration of Momentum-Net consists of three core modules:
image refining, extrapolation, and MBIR.
Momentum-Net guarantees convergence to a fixed-point
for general differentiable (non)convex MBIR functions (or data-fit terms) and convex feasible sets, 
under two asymptomatic conditions.
To consider data-fit variations across training and testing samples, 
we also propose a regularization parameter selection scheme 
based on the spectral spread of majorization matrices.
Numerical experiments for light-field photography using a focal stack and sparse-view computational tomography 
demonstrate that, given identical regression NN architectures,
Momentum-Net significantly improves MBIR speed and accuracy over several existing INNs;
it significantly improves reconstruction quality 
compared to a state-of-the-art MBIR method in each application.  




Iterative neural network, deep learning, 
model-based image reconstruction, inverse problems, 
block proximal extrapolated gradient method, block coordinate descent method, 
light-field photography, X-ray computational tomography.



Momentum-Net: Fast and convergent 
 
iterative neural network for inverse problems

    Il Yong Chun, Member, IEEE, Zhengyu Huang^∗, Hongki Lim^∗, Student Member, IEEE, 
 and Jeffrey A. Fessler, Fellow, IEEE





The authors indicated by asterisks (^∗) contributed equally to this work.
This work is supported in part by the Keck Foundation, NIH U01 EB018753, NIH R01 EB023618, NIH R01 EB022075, and NSF IIS 1838179.
Il Yong Chun was with the Department of Electrical Engineering and Computer Science, The University of Michigan, Ann Arbor, MI 48019 USA,
and is now with the Department of Electrical Engineering, the University of Hawai'i at Mānoa, Honolulu, HI 96822 USA (email: iychun@hawaii.edu).
Zhengyu Huang, Hongki Lim, and Jeffrey A. Fessler are with the Department of Electrical Engineering and Computer Science, The University of Michigan, Ann Arbor, MI 48019 USA (email: 
zyhuang@umich.edu; hongki@umich.edu; fessler@umich.edu).
This paper has appendices. The prefix “A” indicate the numbers in section, theorem, equation, figure, table, and footnote in the appendices.
This paper was presented in part at the Allerton Conference on Communication, Control, and Computing, Monticello, IL, USA, in Oct. 2018.


    
===============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



















§ INTRODUCTION



Deep regression neural network (NN) methods 
have been actively studied for solving diverse inverse problems,
due to their effectiveness at mapping noisy signals into clean signals. 
Examples include 
image denoising 
<cit.>,
image deconvolution 
<cit.>,
image super-resolution 
<cit.>,
magnetic resonance imaging (MRI) 
<cit.>,
X-ray computational tomography (CT) 
<cit.>,
and light-field (LF) photography <cit.>.
However, regression NNs with a greater mapping capability
have increased overfitting/hallucination risks <cit.>. 
An alternative approach to solving inverse problems
is an iterative NN (INN) that combines 
regression NNs – called refiners or denoisers – with an unrolled iterative model-based image reconstruction (MBIR) algorithm 
<cit.>.
This alternative approach can regulate overfitting of regression NNs,
by balancing physical data-fit of MBIR and prior information estimated by refining NNs <cit.>.
This soft-refiner approach has been successfully applied to several extreme imaging systems, 
e.g., highly undersampled MRI
<cit.>,
low-dose or sparse-view CT
<cit.>, 
and low-count emission tomography <cit.>.




 §.§ Notation
 


This section provides mathematical notations.
We use f(x;y) to denote a function f of x given y.
We use ·_p to denote the ℓ^p-norm and write ·· for the standard inner product on ^N.  
The weighted ℓ^2-norm with a Hermitian positive definite matrix A is denoted by ·_A =  A^1/2 (·) _2.
The Frobenius norm of a matrix is denoted by ·_F. 
( · )^T, ( · )^H, and ( · )^* indicate the transpose, complex conjugate transpose (Hermitian transpose), and complex conjugate, respectively. 
(·) denotes the conversion of a vector into a diagonal matrix or diagonal elements of a matrix into a vector.
For (self-adjoint) matrices A,B ∈^N × N, the notation B ≼ A denotes that A-B is a positive semi-definite matrix. 







 §.§ From block-wise optimization to INN
 

To recover signals x ∈^N from measurements y ∈^m, 
consider the following MBIR optimization problem:


_x ∈ F(x;y,z),   F(x;y,z) ≜f(x;y) + γ/2  x - z _2^2,
P0

where  is a set of feasible points,
f(x; y) is data-fit function,
γ is a regularization parameter, 
and z ∈^N is some high-quality approximation to the true unknown signal x.
The data-fit f(x; y) measures deviations of model-based predictions of x from data y, 
considering models of imaging physics (or image formation) and noise statistics in y.
In sys:mbir:z, the signal recovery accuracy increases as the quality of z improves <cit.>; 
however, it is difficult to obtain such z in practice. 
Alternatively, there has been a growing trend
in learning sparsifying regularizers (e.g., convolutional regularizers <cit.>) from training datasets
and applying the trained regularizers to the following block-wise MBIR problem:
_x ∈ f(x;y) + min_ζ r(x, ζ; ).
Here, a learned regularizer min_ζ r(x, ζ; )
quantifies consistency between x and refined sparse signal ζ via some learned operators .
Recently, we have constructed INNs
by generalizing the corresponding block-wise MBIR updates with regression NNs without convergence analysis
<cit.>.
In existing INNs, two major challenges exist: convergence and acceleration.








 §.§ Challenges in existing INNs: Convergence
 


Existing convergence analysis has some practical limitations.
The original form of plug-and-play (PnP <cit.>) is motivated by the alternating direction method of multipliers (ADMM <cit.>),
and its fixed-point convergence has been analyzed with consensus equilibrium perspectives <cit.>.
However, similar to ADMM, its practical convergence depends on how one selects ADMM penalty parameters.
For example, <cit.> reported unstable convergence behaviors
of PnP-ADMM with fixed ADMM parameters.
To moderate this problem,
<cit.> proposed a scheme that adaptively controls the ADMM parameters based on relative residuals.
Similar to the residual balancing technique <cit.>,
the scheme in <cit.> requires tuning initial parameters.
Regularization by Denoising (RED <cit.>)
is an alternative
that moderates some such limitations.
In particular, RED aims to make a clear connection between optimization and a denoiser , by defining its prior term by (scaled) x^T (x - (x)).
Nonetheless, <cit.> showed that many practical denoisers do not satisfy the Jacobian symmetry in <cit.>, and proposed a less restrictive method, score-matching by denoising.




The convergence analysis of the INN inspired by the relaxed projected gradient descent (RPGD) method in <cit.> has the least restrictive conditions on the regression NN among the existing INNs. 
This method replaces the projector of a projected gradient descent method with an image refining NN.
However, the RPGD-inspired INN directly applies an image refining NN to gradient descent updates of data-fit;
thus, this INN relies heavily on the mapping performance of a refining NN and can have overfitting risks, similar to non-MBIR regression NNs, e.g., FBPConvNet <cit.>.
In addition, it exploits the data-fit term only for the first few iterations <cit.>.
We refer the perspective used in RPGD-inspired INN and its related works <cit.> as hard-refiner: different from soft-refiners, these methods do not use a refining NN as a regularizer. 
More recently, 
<cit.> presented convergence analysis
for an INN inspired by a proximal gradient descent method.
However, their analysis is based on noiseless measurements,
which is typically impractical.




Broadly speaking, existing convergence analysis largely depends
on the (firmly) nonexpansive property
of image refining NNs <cit.>, <cit.>, <cit.>. 
However, except for a single-hidden layer convolutional NN (CNN),
it is yet unclear which analytical conditions guarantee the non-expansiveness of general refining NNs <cit.>.
To guarantee convergence of INNs even when using possibly expansive image refining NNs,
we proposed a method that normalizes the output signal of image refining NNs
by their Lipschitz constants
<cit.>.
However, if one uses expansive NNs that are identical across iterations,
it is difficult to obtain best image recovery
with that normalization scheme.
The spectral normalization based training <cit.> can ensure the non-expansiveness of refining NNs
by single-step power iteration. 
However, similar to the normalization method in <cit.>, 
refining NNs trained with the spectral normalization method <cit.> 
degraded the image reconstruction accuracy for an INN using iteration-wise refining NNs <cit.>.
In addition, there does not yet exist theoretical convergence results
when refining NNs change across iterations,
yet iteration-wise refining NNs are widely studied 
<cit.>. 
Finally, existing analysis considers only a narrow class of data-fit terms:
most analyses consider a quadratic function with a linear imaging model <cit.>
or more generally, a convex cost function <cit.> that can be minimized with a practical closed-form solution. 
No theoretical convergence results exist for general (non)convex data-fit terms, iteration-wise NN denoisers, and a general set of feasible points.








 §.§ Challenges in existing INNs: Acceleration



Compared to non-MBIR regression NNs that do not exploit the data-fit f(x;y) in sys:mbir:z,
INNs require more computation
because they consider the imaging physics.
Computation increases as the imaging system or image formation model becomes larger-scale,
e.g., LF photography from a focal stack, 3D CT, parallel MRI using many receive coils, and image super-resolution.
Thus, acceleration becomes crucial for INNs.

First,
consider the existing methods motivated by ADMM or block coordinate descent (BCD) method: 
examples include PnP-ADMM <cit.>, 
RED-ADMM <cit.>, 
MoDL <cit.>, 
BCD-Net <cit.>, etc.
These methods can require multiple inner iterations
to balance data-fit and prior information estimated by trained refining NNs,
increasing total MBIR time.
For example, in solving such problems, each outer iteration involves
x^(i+1)=_x F(x;y,z^(i+1)),
where F is given as in sys:mbir:z
and z^(i+1) is the output from the ith image refining NN.
For LF imaging system using a focal stack data <cit.>, 
solving the above problem requires multiple iterations,
and the total computational cost scale with the numbers of photosensors and sub-aperture images.
In addition,
nonconvexity of the data-fit term f(x;y) can break
convergence guarantees of these methods,
because in general, the proximal mapping _x f(x;y) + γ x - z^(i+1)_2^2 is no longer nonexpansive.

Second,
consider the existing works motivated by gradient descent methods
<cit.>.
These methods resolve the inner iteration issue; 
however, they lack a sophisticated step-size control or backtracking scheme that influences convergence guarantee and acceleration.
Accelerated proximal gradient (APG) methods using momentum terms can significantly accelerate convergence rates
for solving composite convex problems <cit.>,
so we expect that INN methods in the second class have yet to be maximally accelerated. 
The work in <cit.> applied PnP to the APG method <cit.>; 
<cit.> applied PnP to the primal-dual splitting (PDS) algorithm <cit.>.
However, similar to RPGD <cit.>, 
these are hard-refiner methods using some state-of-the-art denoisiers (e.g., BM3D <cit.>) but not trained NNs.
Those methods lack convergence analyses and guarantees may be limited to convex data-fit function.










 §.§ Contributions and organization of the paper
 

This paper proposes Momentum-Net, the first INN architecture that aims for fast and convergent MBIR.
The architecture of Momentum-Net is motivated by applying the Block Proximal Extrapolated Gradient method using a Majorizer (BPEG-M) <cit.> to MBIR using trainable convolutional autoencoders <cit.>.
Specifically, each iteration of Momentum-Net consists of three core modules: image refining, extrapolation, and MBIR. 
At each Momentum-Net iteration, an extrapolation module uses momentum from previous updates to amplify the changes in subsequent iterations and accelerate convergence,
and an MBIR module is noniterative.
In addition, Momentum-Net resolves the convergence issues mentioned in <ref>: for general differentiable (non)convex data-fit terms and convex feasible sets, it guarantees convergence to a point that satisfies fixed-point and critical point conditions, under some mild conditions and two asymptotic conditions, i.e., asymptotically nonexpansive paired refining NNs and asymptotically block-coordinate minimizer.



The remainder of this paper is organized as follows.
<ref> constructs the Momentum-Net architecture 
motivated by BPEG-M algorithm that solves MBIR problem using a learnable convolutional regularizer, 
describes its relation to existing works, 
analyzes its convergence, and summarizes the benefits of Momentum-Net over existing INNs.
<ref> provides details of training INNs, including image refining NN architectures, single-hidden layer or shallow CNN (sCNN)
and multi-hidden layer or deep CNN (dCNN), and training loss function,
and proposes a regularization parameter selection scheme to consider data-fit variations across training and testing samples.
<ref> considers two extreme imaging applications: sparse-view CT and LF photography using a focal stack. 
<ref> reports numerical experiments of applications where
the proposed Momentum-Net using extrapolation significantly improves MBIR speed and accuracy, 
over the existing INNs, 
BCD-Net <cit.>,
Momentum-Net using no extrapolation <cit.>,
ADMM-Net <cit.>,
and PnP-PDS <cit.> using refining NNs.
Furthermore, <ref> reports numerical experiments where 
Momentum-Net significantly improves reconstruction quality compared to a state-of-the-art MBIR method in each application.










§ MOMENTUM-NET: WHERE BPEG-M MEETS NNS FOR INVERSE PROBLEMS







 §.§ Motivation: BPEG-M algorithm for MBIR using learnable convolutional regularizer
 


This section motivates the proposed Momentum-Net architecture, based on our previous works <cit.>.
Consider the following approach for recovering signal x from measurements y (see the setup of block multi-(non)convex problems in <ref>):


    _x ∈   f(x; y) + γ( min_{ ζ_k } r(x, { ζ_k }; { h_k }) ),


    r(x, { ζ_k }; { h_k }) ≜∑_k=1^K 1/2  h_k x - ζ_k  _2^2 + β_k  ζ_k _1,

where  is a closed set, 
f(x; y) + γ r(x, {ζ_k }; { h_k }) is a (continuosly) differentiable (non)convex function in x, 
min_{ζ_k } r(x, {ζ_k } ; { h_k }) is a learnable convolutional regularizer <cit.>, 
{ζ_k : k = 1,…,K } is a set of sparse features that correspond to { h_k  x  }, 
{ h_k ∈^R : k = 1,…,K } is a set of trainable filters, 
and R and K denote the size and number of trained filters, respectively.








Problem sys:recov   caol can be viewed as a two-block optimization problem in terms of the image x and the features {ζ_k }.
We solve sys:recov   caol using the recent BPEG-M optimization framework <cit.>
that has attractive convergence guarantee and rapidly solved several block optimization problems <cit.>.
BPEG-M has the following key ideas for each block optimization problem
(see details in <ref>):




  * M_b-Lipschitz continuity for the gradient of the b block optimization problem, ∀ b:

[] 
A function g: ^n →^n is M-Lipschitz continuous on ^n if there exists a (symmetric) positive definite matrix M such that

g(u) - g(v)_M^-1 ≤u - v_M,   ∀u,v ∈^n.


Definition <ref> is a more general concept than the classical Lipschitz continuity.



  * A sharper majorization matrix M that gives a tighter bound in Definition <ref> 
leads to a tighter quadratic majorization bound in the following lemma: 

[Quadratic majorization via M-Lipschitz continuous gradients ] 
Let f(u) : ^n →. If ∇ f is M-Lipschitz continuous, then 





f(u) ≤f(v) + ∇_u f(v)u-v + 1/2 u - v_M^2,   ∀u,v ∈^n.



Having tighter majorization bounds, 
sharper majorization matrices tend to accelerate BPEG-M convergence.


  * The majorized block problems are proximable, i.e., proximal mapping of majorized function is easily computable depending on the properties of b block majorizer and regularizer, M_b and r_b, 
where the proximal mapping operator is defined by


Prox_r_b^M_b (z) ≜_u   1/2 u - z_M_b^2 + r_b(u),   ∀b.



  * Block-wise extrapolation and momentum terms to accelerate convergence.






Suppose that 1) gradient of f(x; y) + γ r(x, {ζ_k }; { h_k }) is M-Lipschitz continuous at an extrapolated point x́^(i+1), ∀ i; 2) filters in sys:recov   caol satisfy the tight-frame (TF) condition, 
∑_k=1^K  h_k  u _2^2 =  u _2^2, ∀ u,
for some boundary conditions <cit.>.
Applying the BPEG-M framework (see Algorithm <ref>) to solving sys:recov   caol leads to the following block updates:



9.5pt11.4pt

z^(i+1)    = ∑_k=1^K flip(h_k^*) _β_k ( h_k x^(i) ),



x́^(i+1)    = x^(i) + E^(i+1) ( x^(i) - x^(i-1) ),



x^(i+1)     = Prox__^M^(i+1)  
(  x́^(i+1) - ( M^(i+1) )^-1 ∇F (x́^(i+1); y, z^(i+1))   ),

where E^(i+1) is an extrapolation matrix that is given in up:Ex:cvx or up:Ex:ncvx below,
M^(i+1) is a (scaled) majorization matrix for ∇ F(x; y, z^(i+1)) that is given in up:MFtilde below, ∀ i,
the proximal operator Prox__^M^(i+1) (·) in eq:bpgm:x is given by eq:d:prox,
and _ (x) is the characteristic function of set  (i.e., _ equals to 0 if x ∈, and ∞ otherwise).


Proximal mapping update eq:bpgm:zk has a single-hidden layer convolutional autoencoder architecture 
that consists of encoding convolution, nonlinear thresholding, and decoding convolution,
where flip(·) flips a filter along each dimension,
and the soft-thresholding operator _α (u): ^N →^N is defined by 


( _α (u) )_n ≜{ [ u_n - α·(u_n),   | u_n | > α,;             0,              , ] .

for n = 1,…,N, in which (·) is the sign function.
See details of deriving BPEG-M updates eq:bpgm:zk–eq:bpgm:x in <ref>.
The BPEG-M updates in eq:bpgm:zk–eq:bpgm:x guarantee convergence to a critical point, 
when MBIR problem sys:recov   caol satisfies some mild conditions, e.g., lower-boundedness and existence of critical points; see
Assumption S.1 in <ref>. 

The following section generalizes the BPEG-M updates in eq:bpgm:zk–eq:bpgm:x and constructs the Momentum-Net architecture.












 §.§ Architecture
 


This section establishes the INN architecture of Momentum-Net by generalizing BPEG-M updates eq:bpgm:zk–eq:bpgm:x that solve sys:recov   caol.
Specifically, we replace the proximal mapping in eq:bpgm:zk with a general image refining NN _θ (·), 
where θ denotes the trainable parameters.
To effectively remove iteration-wise artifacts and give best signal estimates at each iteration,
we further generalize a refining NN _θ (·) to 
iteration-wise image refining NNs {_θ^(i+1) (·) : i = 0,…,N_iter-1 }, 
where θ^(i+1) denotes the parameters for the i iteration refining NN _θ^(i+1), 
and N_iter is the number of Momentum-Net iterations.
The iteration-wise NNs are particularly useful for reducing overfitting risks in regression,
because _θ^(i+1) is responsible for removing noise features only at the i iteration, 
and thus one does not need to greatly increase dimensions of its parameter θ^(i+1) <cit.>.
In low-dose CT reconstruction, for example, the refining NNs at the early and later iterations remove streak artifacts and Gaussian-like noise, respectively <cit.>.








Each iteration of Momentum-Net consists of 1) image refining, 2) extrapolation, and 3) MBIR modules, corresponding to the BPEG-M updates eq:bpgm:zk, eq:bpgm:xacute, and eq:bpgm:x, respectively. 
See the architecture of Momentum-Net in Fig. <ref>(a) and Algorithm <ref>.
At the i iteration, Momentum-Net performs the following three processes:





  * Refining: 
The i image refining module gives the refined image z^(i+1), by applying the i refining NN, _θ^(i+1), to an input image at the i iteration, x^(i) (i.e., image estimate from the (i-1) iteration).
Different from existing INNs, e.g., ADMM-Net <cit.>, PnP-ADMM <cit.>, RED <cit.>, 
MoDL <cit.>, BCD-Net <cit.> (see Fig. <ref>(b)), TNRD <cit.>, 
we apply ρ-relaxation with ρ∈ (0,1); see eq:momnet:map.
The parameter ρ controls the strength of inference from refining NNs, 
but does not affect the convergence guarantee of Momentum-Net.
Proper selection of ρ can improve MBIR accuracy (see <ref>).


  * Extrapolation: 
The i extrapolation module gives the  extrapolated point x́^(i+1), based on momentum terms x^(i) - x^(i-1); see eq:momnet:exp.
Intuitively speaking, momentum is information from previous updates to amplify the changes in subsequent iterations. 
Its effectiveness has been shown in diverse optimization literature, e.g., convex optimization <cit.> and block optimization <cit.>.


  * MBIR: 
Given a refined image z^(i+1) and a measurement vector y,
the i MBIR module eq:momnet:mbir applies the proximal operator Prox__^M^(i+1) (·)
to the extrapolated gradient update using a quadratic majorizer of F(x;y,z^(i+1)), 
where F is defined in sys:mbir:z.
Intuitively speaking, this step solves a majorized version of the following MBIR problem at the extrapolated point x́^(i+1):


min_x ∈ F(x;y,z^(i+1)),
P1

and gives a reconstructed image x^(i+1).
In Momentum-Net, we consider (non)convex differentiable MBIR cost functions F with M-Lipschitz continuous gradients,
and a convex and closed set .
For a wide range of large-scale inverse imaging problems, 
the majorized MBIR problem eq:momnet:mbir has a practical closed-form solution and thus, does not require an iterative solver,
depending on the properties of practically invertible majorization matrices M^(i+1) and constraints.
Examples of M^(i+1)- combinations that give a noniterative solution for eq:momnet:mbir include 
scaled identity and diagonal matrices with a box constraint and the non-negativity constraint,
and matrices decomposable by unitary transforms, e.g., a circulant matrix <cit.>, with = ^N.
The updated image x^(i+1) is the input to the next Momentum-Net iteration.











The followings are details of Momentum-Net in Algorithm <ref>.
A scaled majorization matrix is


M^(i+1) = λ·M^(i+1) ≻0,     λ≥1,

where M^(i+1)∈^N × N is 
a symmetric positive definite majorization matrix of ∇ F(x;y,z^(i+1))
in the sense of M-Lipschitz continuity (see Definition <ref>).
In up:MFtilde, λ = 1 and λ > 1 for convex and nonconvex F(x;y,z^(i+1)) (or convex and nonconvex f(x;y)), respectively.
We design the extrapolation matrices as follows:







    


    E^(i+1) =
δ^2 m^(i) ·( M^(i+1) )^-1/2 
( M^(i) )^1/2;
 


    


    E^(i+1) = 
δ^2 m^(i) ·λ-1/2 (λ+1) ·( M^(i+1) )^-1/2 
( M^(i) )^1/2, 


for some δ< 1 and { 0 ≤ m^(i)≤ 1 :∀ i }.
We update the momentum coefficients { m^(i+1):∀ i } by the following formula <cit.>:


m^(i+1) = θ^(i)  - 1/θ^(i+1),     θ^(i+1) = 1 + √(1 + 4 (θ^(i))^2)/2;

if F(x;y,z^(i+1)) has a sharp majorizer, i.e., ∇ F(x;y,z^(i+1)) has M^(i+1) 
such that the corresponding bound in Definition <ref> is tight, then we set m^(i+1) = 0, ∀ i.
<ref> lists parameters of Momentum-Net, 
and summarizes selection guidelines or gives default values.









 §.§ Relations to previous works
 

Several existing MBIR methods can be viewed as a special case of Momentum-Net:


(MBIR model sys:recov   caol using convolutional autoencoders that satisfy the TF condition <cit.>). 
The BPEG-M updates in eq:bpgm:zk–eq:bpgm:x are special cases of the modules in Momentum-Net (Algorithm <ref>), 
with {_θ^(i+1) ( · ) =∑_k=1^K flip(h_k^*) _β_k (h_k  (·) ) :∀ i } 
(i.e., single hidden-layer convolutional autoencoder <cit.>) and ρ≈ 1.
These give a clear mathematical connection between a denoiser eq:bpgm:zk and cost function sys:recov   caol.
One can find a similar relation between a multi-hidden layer CNN and a multi-layer convolutional regularizer <cit.>.





(INNs inspired by gradient descent method, e.g., TNRD <cit.>).
Removing extrapolation modules, i.e., setting { E^(i+1) = 0 : ∀ i } in eq:momnet:exp, and setting ρ≈ 1, Momentum-Net becomes the existing INN in <cit.>.



(BCD-Net for image denoising <cit.>).
To obtain a clean image x ∈^N from a noisy image y ∈^N corrupted by an additive white Gaussian noise (AWGN),
MBIR problem sys:mbir considers the data-fit f(x;y) = 1/2 y -  x _W^2 
with the inverse covariance matrix W = 1/σ^2 I, where σ^2 is a variance of AWGN, 
and the box constraint = [0, U]^N with an upper bound U > 0.
For this f(x;y), the MBIR module eq:momnet:mbir can use the exact majorizer 
{M^(i+1) = ( 1/σ^2 + γ) I  } and one does not need to use the extrapolation module eq:momnet:exp,
i.e., { E^(i+1) = 0  }.
Thus, Momentum-Net (with ρ≈ 1) becomes BCD-Net.




(BCD-Net for undersampled single-coil MRI <cit.>).
To obtain an object magnetization x ∈^N from a k-space data y ∈^m obtained by undersampling (e.g., compressed sensing <cit.>) MRI,
MBIR problem sys:mbir considers the data-fit f(x;y) = 1/2 y -  A x _W^2 
with an undersampling Fourier operator A
(disregarding relaxation effects and considering Cartesian k-space),
the inverse covariance matrix W = 1/σ^2 I, where σ^2 is a variance of complex AWGN <cit.>,
and = ^N.
For this f(x;y), the MBIR module eq:momnet:mbir can use the exact majorizer {M^(i+1) = F_^H ( 1/σ^2 P + γ I) F_} that is practically invertible,
where F_ is the discrete Fourier transform and 
P is a diagonal matrix with either 0 or 1 
(their positions correspond to sampling pattern in k-space), 
and the extrapolation module eq:momnet:exp uses the zero extrapolation matrices { E^(i+1) = 0 }.
Thus, Momentum-Net (with ρ≈ 1) becomes BCD-Net.



The following section analyzes the convergence of Momentum-Net.













 §.§ Convergence analysis
 


In practice, INNs, i.e., unrolled or PnP methods using refining NNs,
are trained and used with a specific number of iterations.
Nevertheless, similar to optimization algorithms, studying convergence properties of INNs with N_iter→∞ <cit.> is important; 
in particular, it is crucial to know if a given INN tends to converge as N_iter increases.
For INNs using iteration-wise refining NNs, e.g., BCD-Net <cit.> and proposed Momentum-Net, 
we expect that refiners converge,
i.e., their image refining capacity converges,
because information provided by data-fit function f(x;y) in MBIR (e.g., likelihood) reaches some bound after a certain number of iterations.
Fig. <ref> illustrates that dCNN parameters of Momentum-Net tend to converge for different applications.
(The similar behavior was reported for sCNN refiners in BCD-Net <cit.>.)
Although refiners do not completely converge, in practice,  
one could use a refining NN at a sufficiently large iteration number, e.g., N_iter= 100 in Momentum-Net,  
for the later iterations.


There are two key challenges in analyzing the convergence of Momentum-Net in Algorithm <ref>: 
both challenges relate to its image refining modules eq:momnet:map.
First, image refining NNs _θ^(i+1) change across iterations; 
even if they are identical across iterations, they are not necessarily nonexpansive operators <cit.> in practice.
Second, the iteration-wise refining NNs are not necessarily proximal mapping operators, 
i.e., they are not written explicitly in the form of eq:d:prox.
This section proposes two new asymptotic definitions to overcome these challenges, 
and then uses those conditions to analyze convergence properties of Momentum-Net in Algorithm <ref>.





  §.§.§ Preliminaries
 

To resolve the challenge of iteration-wise refining NNs and the practical difficulty in guaranteeing their non-expansiveness, 
we introduce the following generalized definition of the non-expansiveness <cit.>.


[Asymptotically nonexpansive paired operators] 
A sequence of paired operators ( _θ^(i), _θ^(i+1) ) is asymptotically nonexpansive if there exist a summable nonnegative sequence {ϵ^(i+1)≥ 0 : ∑_i=0^∞ϵ^(i+1) < ∞} such that[
One could replace the bound in cond:eps with _θ^(i+1)(u) - _θ^(i) (v) _2^2 ≤ (1+ϵ^(i+1)) u - v_2^2 (and summable {ϵ^(i+1) : ∀ i }), and the proofs for our main arguments go through. 
]


 _θ^(i+1) (u) - _θ^(i) (v) _2^2 ≤u - v_2^2 + ϵ^(i+1),     ∀u,v, i.






When _θ^(i+1)=_θ and  ϵ^(i+1)= 0, ∀ i, 
Definition <ref> becomes the standard non-expansiveness of a mapping operator _θ.
If we replace the inequality (≤) with the strict inequality (<) in cond:eps, then we say that the sequence of paired operators ( _θ^(i+1), _θ^(i+1) ) is asymptotically contractive. (This stronger assumption is used to prove convergence of BCD-Net in Proposition <ref>.) 
Definition <ref> also implies that mapping operators _θ^(i+1) converge to some nonexpansive operator, 
if the corresponding parameters θ^(i+1) converge.




Definition <ref> incorporates a pairing property because Momentum-Net uses iteration-wise image refining NNs. 
Specifically, the pairing property helps prove convergence of Momentum-Net,  
by connecting image refining NNs at adjacent iterations.
Furthermore, the asymptotic property in Definition <ref> allows Momentum-Net to use expansive refining NNs (i.e., mapping operators having a Lipschitz constant larger than 1) for some iterations, while guaranteeing convergence; see Figs. <ref>(a3) and <ref>(b3).
Suppose that refining NNs are identical across iterations, i.e., _θ^(i+1)=_θ, ∀ i, 
similar to some existing INNs, e.g., PnP <cit.>, RED <cit.>, and other methods in <ref>.
In such cases, if _θ is expansive, Momentum-Net may diverge;
this property corresponds to the limitation of existing methods described in <ref>.
Momentum-Net moderates this issue by using iteration-wise refining NNs 
that satisfy the asymptotic paired non-expansiveness in Definition <ref>.


Because the sequence { z^(i+1) : ∀ i } in eq:momnet:map is not necessarily updated with a proximal mapping, 
we introduce a generalized definition of block-coordinate minimizers <cit.> for z^(i+1)-updates:


[Asymptotic block-coordinate minimizer]  
The update z^(i+1) is an asymptotic block-coordinate minimizer if there exists a summable nonnegative sequence {Δ^(i+1)≥ 0 : ∑_i=0^∞Δ^(i+1) < ∞} such that


 z^(i+1) - x^(i) _2^2 ≤ z^(i) - x^(i) _2^2 + Δ^(i+1),     ∀i.




Definition <ref> implies that as i →∞, the updates { z^(i+1) : i ≥ 0 } approach a block-coordinate minimizer trajectory that satisfies cond:delta with {Δ^(i+1)= 0 : i ≥ 0 }.
In particular, Δ^(i+1) quantifies how much the update z^(i+1) in eq:momnet:map perturbs a block-coordinate minimizer trajectory.
The bound z^(i+1)  - x^(i)_2^2 ≤ z^(i)  - x^(i)_2^2 always holds, ∀ i, when one uses the proximal mapping in eq:bpgm:zk within the BPEG-M framework.
Note that while applying trained Momentum-Net, cond:delta is easy to examine empirically, whereas cond:eps is harder to check.





 










  §.§.§ Assumptions
 

This section introduces and interprets the assumptions for convergence analysis of Momentum-Net in Algorithm <ref>:





  * Assumption 1) In MBIR problems sys:mbir,
(non)convex F(x;y,z^(i+1)) is (continuously) differentiable, proper, and lower-bounded in (F),[
F : ^n → (- ∞, + ∞] is proper if F ≠∅. 
F is lower bounded in (F) ≜{ u : F(u) < ∞} if inf_u ∈(F) F(u) > -∞.
]
∀ i, and  is convex and closed.
Algorithm <ref> has a fixed-point.



  * Assumption 2) ∇ F(x;y,z^(i+1)) is M^(i+1)-Lipschitz continuous with respect to x (see Definition <ref>),











where M^(i+1) is a iteration-wise majorization matrix that satisfies m_F,min I_N≼ M^(i+1)≼ m_F,max I_N with 0 < m_F,min≤ m_F,max < ∞, ∀ i.



  * Assumption 3) The extrapolation matrices E^(i+1)≽ 0 in up:Ex:cvx–up:Ex:ncvx satisfy the following conditions:




     


    ( E^(i+1) )^T M^(i+1) E^(i+1)
≼δ^2 ·M^(i),     δ< 1;



     


    ( E^(i+1) )^T M^(i+1) E^(i+1)
≼ δ^2 (λ-1)^2/4 (λ+1)^2 ·M^(i),     δ< 1.





  * Assumption 4)  
The sequence of paired operators ( _θ^(i+1), _θ^(i) ) is asymptotically nonexpansive with a summable sequence {ϵ^(i+i)≥ 0 };
the update z^(i+1) is an asymptotic block-coordinate minimizer with a summable sequence {Δ^(i+i)≥ 0 }.
The mapping functions {_θ^(i+1) : ∀ i } are continuous with respect to input points and the corresponding parameters {θ^(i+1) : ∀ i } are bounded.





Assumption 1 is a slight modification of Assumption S.1 of BPEG-M, and Assumptions 2–3 are identical to Assumptions S.2–S.3 of BPEG-M;
see Assumptions S.1–S.3 in <ref>.
The extrapolation matrix designs up:Ex:cvx and up:Ex:ncvx satisfy conditions cond:Ex:cvx and cond:Eu:ncvx in Assumption 3, respectively.


We provide empirical justifications for the first two conditions in Assumption 4.
First, Figs. <ref>(a2) and <ref>(a2) illustrate that 
paired refining NNs ( _θ^(i+1), _θ^(i) ) of Momentum-Net 
appear to be asymptotically nonexpansive 
in an application that has mild condition number variations across training data-fit majorization matrices.
Figs. <ref>(a3), <ref>(b3), <ref>(a3), and <ref>(b3)
illustrate for different applications that 
refining NNs {_θ^(i+1)} become nonexpansive:
their Lipschitz constants at the first several iterations are larger than 1, 
and their Lipschitz constants in later iterations become less than 1.
Alternatively, the asymptotic non-expansiveness of paired operators ( _θ^(i+1), _θ^(i) ) 
can be satisfied by a stronger assumption that the sequence {_θ^(i+1)} 
converges to some nonexpansive operator.
(Fig. <ref> illustrates that dCNN parameters of Momentum-Net appear to converge.)


Figs. <ref>(a3), <ref>(b3), <ref>(a3), and <ref>(b3)
illustrate for different applications that the z^(i+1)-updates are asymptotic block-coordinate minimizers.
Lemma <ref> and <ref> in the appendices provide 
a probabilistic justification for the asymptotic block-coordinate minimizer condition.













  §.§.§ Main convergence results
 


This section analyzes fixed-point and critical point convergence of Momentum-Net in Algorithm <ref>, under the assumptions in the previous section. 
We first show that differences between two consecutive iterates generated by Momentum-Net converge to zero: 



[Convergence properties]

Under Assumptions 1–4, let { x^(i+1), z^(i+1) : i ≥ 0 } be the sequence generated by Algorithm <ref>. 
Then, the sequence satisfies 


∑_i=0^∞  [ =1.5pt [ x^(i+1); z^(i+1) ] ] -  [ =1.5pt [ x^(i); z^(i) ] ] _2^2 < ∞,

and hence [ =1.5pt [ x^(i+1); z^(i+1) ]] -  [ =1.5pt [ x^(i); z^(i) ]] _2 → 0.


See <ref> in the appendices.





Using Proposition <ref>, our main theorem provides that any limit points of the sequence generated by Momentum-Net satisfy critical point and fixed-point conditions:


[A limit point satisfies both critical point and fixed-point conditions]

Under Assumptions 1–4 above, let { x^(i+1), z^(i+1) : i ≥ 0 } be the sequence generated by Algorithm <ref>.
Consider either a fixed majorization matrix with general structure, i.e., M^(i+1)= M for i ≥ 0, 
or a sequence of diagonal majorization matrices, i.e., { M^(i+1) : i ≥ 0 }.
Then, any limit point x̅ of { x^(i+1)} satisfies both the critical point condition:


∇F (x̅;y,z̅)  x - x̅  ≥0,     ∀x ∈,

where z̅ is a limit point of { z^(i+1)}, and the fixed-point condition:


[ [ x̅; x̅ ] ]  = __θ̅^M̅  ( [ [ x̅; x̅ ] ] ),

where [ [ x^(i+1);   x^(i) ]] =__θ^(i+1)^M^(i+1)( [ [   x^(i); x^(i-1) ]] ),
__θ^(i+1)^M^(i+1) (·) denotes performing the i updates in Algorithm <ref>, 
and θ̅ and M̅ is a limit point of {θ^(i+1)} and { M^(i+1)}, respectively.


See <ref> in the appendices.




Observe that, if = ^N or x̅ is an interior point of , 
t:momnet:convg:crit reduces to the first-order optimality condition 0 ∈∂ F(x̅;y,z̅),
where ∂ F(x) denotes the limiting subdifferential of F at x.
With additional isolation and boundedness assumptions for the points satisfying t:momnet:convg:crit and t:momnet:convg:fix,
we obtain whole sequence guarantees:





[Whole sequence convergence]

Consider the construction in Theorem <ref>.
Let  be the set of points satisfying the critical point condition in t:momnet:convg:crit and the fixed-point condition in t:momnet:convg:fix.
If { x^(i+1) : i ≥ 0 } is bounded, then dist (x^(i+1), ) → 0, where dist (u, ) ≜inf{ u - v  : v ∈} denotes the distance from u to , for any point u ∈^N and any subset ⊂^N.
If  contains uniformly isolated points, i.e., there exists η > 0 such that u - v ≥η for any distinct points u,v ∈, then { x^(i+1)} converges to a point in .


See <ref> in the appendices.






The boundedness assumption for { x^(i+1)} in Corollary <ref> is standard in block-wise optimization, e.g., <cit.>. The assumption can be satisfied if 
the set  is bounded (e.g., box constraints),
one chooses appropriate regularization parameters in Algorithm <ref> <cit.>, 
the function F(x;y,z) is coercive <cit.>, 
or the level set is bounded <cit.>.
However, for general F(x;y,z), it is hard to verify the isolation condition for the points in  in practice.
Instead, one may use Kurdyka-Łojasiewicz property <cit.> to analyze the whole sequence convergence with some appropriate modifications.



For simplicity, we focused our discussion to noniterative MBIR module eq:momnet:mbir.
However, Momentum-Net practically converges with any proximable MRIR function eq:momnet:mbir that may need an iterative solver,
if sufficient inner iterations are used. 
To maximize the computational benefit of Momentum-Net, 
one needs to make sure that majorized MBIR function eq:momnet:mbir is better proximable
over its original form sys:mbir.






 §.§ Benefits of Momentum-Net
 


Momentum-Net has several benefits over existing INNs:
[]


  * Benefits from refining module:
The image refining module eq:momnet:map can use iteration-wise image refining NNs {_θ^(i+1): i ≥ 0 }:
those are particularly useful to reduce overfitting risks by reducing dimensions of their parameters θ^(i+1) at each iteration <cit.>.
Iteration-wise refining NNs 
require less memory for training, compared to methods that use a single refining NN for all iterations, e.g., <cit.>.
Different from the existing methods mentioned in <ref>,
Momentum-Net does not require (firmly) nonexpansive mapping operators {_θ^(i+1)}
to guarantee convergence.
Instead, {_θ^(i+1)} in eq:momnet:map assumes a generalized notion of
the (firm) non-expansiveness condition assumed for convergence of the existing methods 
that use identical refining NNs across iterations, including 
PnP <cit.>, 
RED <cit.>, etc.
The generalized concept is the first practical condition to 
guarantee convergence of INNs using iteration-wise refining NNs; 
see Definition <ref>.



  * Benefits from extrapolation module:
The extrapolation module eq:momnet:exp uses the momentum terms x^(i) - x^(i-1)
that accelerate the convergence of Momentum-Net.
In particular, compared to the existing gradient-descent-inspired INNs,
e.g., TNRD <cit.>,
Momentum-Net converges faster.
(Note that the way the authors of <cit.> used momentum is less conventional.
The corresponding method, RED-APG <cit.>,
still can require multiple inner iterations to solve its quadratic MBIR problem, similar to BCD-Net-type methods.)



  * Benefits from MBIR module: 
The MBIR module eq:momnet:mbir does not require multiple inner iterations for a wide range of imaging problems and has both theoretical and practical benefits.
Note first that convergence analysis of INNs (including Momentum-Net) assumes that their MBIR operators are noniterative.
In other words, related convergence theory (e.g., Proposition <ref>) is inapplicable 
if iterative methods, particularly with insufficient number of iterations, are applied to MBIR modules.
Different from the existing BCD-Net-type methods <cit.> that can require iterative solvers for their MBIR modules,
MBIR module eq:momnet:mbir of Momentum-Net can have practical close-form solution (see examples in <ref>), 
and its corresponding convergence analysis (see <ref>) can hold stably for a wide range of imaging applications.
Second, combined with extrapolation module eq:momnet:exp, 
noniterative MBIR modules eq:momnet:mbir lead to faster MBIR, compared to the existing BCD-Net-type methods that can require multiple inner iterations for their MBIR modules for convergence.
Third, Momentum-Net guarantees convergence even for nonconvex MBIR cost function F(x;y,z) or nonconvex data-fit f(x;y) of which the gradient is M-Lipschitz continuous (see Definition <ref>), while existing INNs overlooked nonconvex F(x;y,z) or f(x;y).



Furthermore, <ref> analyzes the sequence convergence of BCD-Net <cit.>, 
and describes the convergence benefits of Momentum-Net over BCD-Net.








§ TRAINING INNS
 

This section describes training of all the INNs compared in this paper.





 §.§ Architecture of refining NNs and their training
 

For all INNs in this paper, we train the refining NN at each iteration
to remove artifacts from the input image x^(i) that is fed from the previous iteration. 
For the i iteration NN, we first consider the following sCNN architecture, residual single-hidden layer convolutional autoencoder:


_θ^(i+1) (u) = ∑_k=1^K d_k^(i+1) _exp(α_k^(i+1))  ( e_k^(i+1) u ) + u, 

where θ^(i+1)={ d_k^(i+1), α_k^(i+1), e_k^(i+1):∀ k  } is the parameter set of the i image refining NN, 
{ d_k^(i+1), e_k^(i+1)∈^R: k = 1,…,K } is a set of K decoding and encoding filters of size R, 
{exp( α_k^(i+1) ) : k = 1,…,K } is a set of K thresholding values, 
and _α (u) is the soft-thresholding operator with parameter α defined in eq:soft-threshold,
for i = 0,…,N_iter-1.
We use the exponential function exp(·) to prevent the thresholding parameters {α_k } from becoming negative during training.
We observed that the residual convolutional autoencoder in sys:auto:res gives better results compared to
the convolutional autoencoder, i.e., sys:auto:res without the second term <cit.>. 
This corresponds to the empirical result in <cit.> that 
having skip connections (e.g., the second term in sys:auto:res)
can improve generalization.
The sequence of paired sCNN refiners sys:auto:res can satisfy the asymptotic non-expansiveness, 
if its convergent refiner satisfies that 

σ_max (D^H D) ≤1/R,     σ_max (E^H E) ≤1/R,

where 
σ_max(·) is the largest eigenvalue of a matrix,
D≜ [d̅_1, …, d̅_K, δ_R], 
E≜ [e̅_1, …, e̅_K, δ_R], 
{d̅_k, e̅_k : ∀ k } are limit point filters, and δ_R is the Kronecker delta filter of size R.







For dCNN refiners, we use the following residual multi-hidden layer CNN architecture, 
a simplified DnCNN <cit.> using fewer layers, no pooling, and no batch normalization <cit.>
(we drop superscript indices (·)^(i) for simplicity):



_θ (u)    = u - ∑_k=1^K e_k^[L] u_k^[L-1],


u_k^[1]    = ReLU ( e_k^[1] u ),
  u_k^[l] = ReLU ( ∑_k'=1^K e_k,k'^[l] z_k'^[l-1] ), 

for k = 1,…, K and l = 2,…, L-1,
where θ={ e_k^[l], e_k,k'^[l]:∀ k,k',l } is the parameter set of each refining NN,
K is the number of feature maps,
L is the number of layers,
{ e_k^[l]∈^R : k = 1,…,K, l = 1,L } is a set of filters at the first and last dCNN layer,
{ e_k,k'^[l]∈^R : k,k' = 1,…,K, l = 2,…,L-1 } is a set of filters for remaining dCNN layer,

and the rectified linear unit activation function is defined by
ReLU (u) ≜max (0, u).






The training process of Momentum-Net requires S high-quality training images, { x_s: s = 1,…,S }, 
S training measurements simulated via imaging physics, { y_s: s = 1,…, S }, 
and S data-fits { f_s (x;y_s) : s = 1,…,S } 
and the corresponding majorization matrices { M_s^(i), M_s^(i): s = 1,…,S, i = 1,…,N_iter}.
Different from <cit.> that train convolutional autoencoders from the patch perspective, 
we train the image refining NNs in sys:auto:res–sys:dcnn from the convolution perspective (that does not store many overlapping patches, e.g., see <cit.>). 
From S training pairs ( x_s, x_s^(i) ), 
where { x_s^(i): s = 1,…,S } is a set of S reconstructed images at the (i-1) Momentum-Net iteration,
we train the i iteration image refining NN in sys:auto:res by solving the following optimization problem: 


θ^(i+1) 
= _ θ 1/2S ∑_s=1^S  x_s - _θ ( x_s^(i) ) _2^2,
P2

where θ^(i+1) is given as in sys:auto:res, for i = 0,…,N_iter- 1
(see some related properties in <ref>).
We solve the training optimization problems sys:auto:res:train by mini-batch stochastic optimization 
with the subdifferentials computed by the PyTorch  package.







 §.§ Regularization parameter selection based on spectral spread
 

When majorization matrices of training data-fits { f_s (x;y_s) : s = 1,…,S } have similar spectral properties, e.g., condition numbers,
the regularization parameter γ in sys:mbir is trainable by substituting eq:momnet:map
into eq:momnet:mbir and modifying the training cost sys:auto:res:train.
However, the condition numbers of data-fit majorizers can greatly differ due a variety of
imaging geometries or image formation systems, 
or noise levels in training measurements, etc.
See such examples in <ref>–<ref>.



To train Momentum-Net with diverse training data-fits,
we propose a parameter selection scheme based on the spectral spread of
their majorization matrices { M_f_s^(i)}.
For simplicity, consider majorization matrices of the form
M_s^(i)=M_s=λ ( M_f_s  + γ_s I ) ∀ i,
where the factor λ is selected by up:MFtilde and 
M_f_s is a symmetric positive semidefinite majorization matrix for f_s (x;y_s), ∀ s.
We select the regularization parameter γ_s for the s training sample as


γ_s = σ_spread (M_f_s)/χ,

where the spectral spread of a symmetric positive definite matrix is defined by σ_spread(·) ≜σ_max(·) - σ_min (·)
for σ_max (M_f_s) >σ_min (M_f_s) ≥ 0,
and σ_min (·) is the smallest eigenvalue of a matrix.
For the s training sample, 
a tunable factor χ controls γ_s in eq:reg:select
according to σ_spread (M_f_s), ∀ s.
The proposed parameter selection scheme also applies to testing Momentum-Net,
based on the tuned factor χ^⋆ in its training.
We observed that the proposed parameter selection scheme eq:reg:select
gives better MBIR accuracy than the condition number based selection scheme 
that is similarly used in selecting ADMM parameters <cit.> (for the two applications in <ref>).
One may further apply this scheme to iteration-wise majorization matrices M_s^(i)
and select iteration-wise regularization parameters γ_s^(i) accordingly.
For comparing different INNs, we apply eq:reg:select to all INNs.











§ EXPERIMENTAL RESULTS AND DISCUSSION
 

We investigated two extreme imaging applications: 
sparse-view CT and LF photography using a focal stack.
In particular, these two applications lack a practical closed-form solution 
for the MBIR modules of BCD-Net and ADMM-Net <cit.>, e.g., solving eq:bcdnet:recon.
For these applications, we compared the performances of the following five INNs: 
BCD-Net <cit.> 
(i.e., generalization of RED <cit.> and MoDL <cit.>), 
ADMM-Net <cit.>, i.e., PnP-ADMM <cit.> using iteration-wise refining NNs,
Momentum-Net without extrapolation 
(i.e., generalization of TNRD <cit.>), 
PDS-Net, i.e., PnP-PDS <cit.> using iteration-wise refining NNs,
and the proposed Momentum-Net using extrapolation.





 §.§ Experimental setup: Imaging
 



  §.§.§ Sparse-view CT
 


To reconstruct a linear attenuation coefficient image x ∈^N 
from post-log sinogram y ∈^m in sparse-view CT, 
the MBIR problem sys:mbir considers a data-fit f(x;y) = 1/2 y -  A x _W^2
and the non-negativity constraint = [0, ∞)^N,
where A ∈^m × N is an undersampled CT system matrix,
W ∈^m × m is a diagonal weighting matrix with elements 
{ W_m',m' = p_m'^2 / ( p_m' + σ^2 ) : ∀ m'  } based on 
a Poisson-Gaussian model <cit.> 
for the pre-log raw measurements p ∈^m 
with electronic readout noise variance σ^2.

We simulated 2D sparse-view sinograms of size m = 888 × 123 
– detectors or rays × regularly spaced projection views or angles, 
where 984 is the number of full views – with GE LightSpeed fan-beam geometry corresponding to a monoenergetic source with 
10^5 incident photons per ray and no background events, and electronic noise variance σ^2 = 5^2. 
We avoided an inverse crime in imaging simulation and 
reconstructed images of size N = 420 × 420 with a coarser grid Δ_x =Δ_y = 0.9766 mm; 
see details in <cit.>.





  §.§.§ LF photography using a focal stack



To reconstruct a LF x = [x_1^T, …, x_C' ]^T ∈^SN' 
that consists of C' sub-aperture images
from focal stack measurements y = [y_1^T, …, y_C^T]^T ∈^CN' 
that are collected by C photosensors, 
the MBIR problem sys:mbir considers a data-fit f(x;y) =1/2 y -  A x _2^2
and a box constraint = [0, U]^C'N' with U = 1 (or 255 without rescaling),
where A ∈^CN' × C'N' is a system matrix of LF imaging system using a focal stack
that is constructed blockwise with 
C · C' different convolution matrices {τ_c A'_c,c'∈^N' × N': c = 1,…,C, c' = 1,…,C' } <cit.>,
τ_c ∈ (0,1] is a transparency coefficient for the c detector,
and N' is the size of sub-aperture images, x_c', ∀ c'.[
Traditionally, one obtains focal stacks by physically moving imaging sensors and taking separate exposures across time.
Transparent photodetector arrays <cit.> allow one to collect focal stack data in a single exposure, 
making a practical LF camera using a focal stack.
If some photodetectors are not perfectly transparent, one can use τ_c < 1, for some c.
]
In general, a LF photography system using a focal stack is extremely under-determined, 
because C ≪ C'.



To avoid an inverse crime, our imaging simulation used higher-resolution synthetic LF dataset <cit.>
(we converted the original RGB sub-aperture images to grayscale ones 
by the  function in MATLAB,
for simplicity and smaller memory requirements in training).
We simulated C = 5 focal stack images of size N' = 255 × 255 
with 40 dB AWGN that models electronic noise at sensors, and setting transparency coefficients τ_c as 1, for c = 1,…,C.
The sensor positions were chosen such that five sensors focus at equally spaced depths;
specifically, the closest sensor to scenes and farthest sensor from scenes focus 
at two different depths that correspond to _min+ 0.2 and _max- 0.2, respectively,
where _max and _min are the approximate maximum and minimum disparity values specified in <cit.>.
We reconstructed 4D LFs that consist of S = 9 × 9 sub-aperture images of size N' = 255 × 255, 
with a coarser grid Δ_x =Δ_y = 0.13572 mm.













 §.§ Experimental setup: INNs
 




  §.§.§ Parameters of INNs
 


The parameters for the INNs compared in sparse-view CT experiments were defined as follows.
We considered two BCD-Nets (see Algorithm <ref>):
for one BCD-Net, we applied the APG method <cit.> with 10 inner iterations to eq:bcdnet:recon,
and set N_iter= 30; 
for the other BCD-Net, we applied the APG method with 3 inner iterations to eq:bcdnet:recon,
and set N_iter= 45.
For ADMM-Net,
we used the identical configurations as BCD-Net and set the ADMM penalty parameter to γ in sys:mbir, similar to <cit.>.
For Momentum-Net without extrapolation,
we chose N_iter= 100 and ρ= 1 - ε.
For the proposed Momentum-Net,
we chose N_iter= 100 and ρ= 0.5.
For PDS-Net, 
we set the first step size to γ_1 =γ^-1 and the second step size to γ_2 =γ_1^-1σ_max^-1 (M), 
per <cit.>.
For performance comparisons between different INNs,
all the INNs used sCNN refiners sys:auto:res with { R , K = 7^2 } to avoid the overfitting/hallucination risks.
For Momentum-Net using dCNN refiners, 
we chose L=4 layer dCNN sys:dcnn using R = 3^2 filters and K = 64 feature maps, following <cit.>.
(The chosen parameters gave lower RMSE values than { L= 6, R = 3^3 , K = 64 }, 
for identical regularization parameters.)
For comparing different MBIR methods,
Momentum used 
extrapolation, i.e., eq:momnet:exp with up:Ex:cvx and eq:mom_coeff, 
and { R = 7^2, K = 9^2 } for sys:auto:res.
We designed the majorization matrices as {M^(i+1)=( A^T W A 1 ) + γ I : i ≥ 0 },
using Lemma <ref> (A and W have nonnegative entries)
and setting λ= 1 by up:MFtilde.
We set an initial point of INNs, x^(0), to filtered-back projection (FBP) using a Hanning window.
The regularization parameters of all INNs
were selected by the scheme in <ref>
with χ^⋆= 167.64.
(This factor was estimated from the carefully chosen regularization parameter 
for sparse-view CT MBIR experiments using learned convolutional regularizers in <cit.>.)


The parameters for the INNs compared in experiments of LF photography using a focal stack were defined as follows.
We considered two BCD-Nets and two ADMM-Nets with the identical parameters listed above.
For Momentum-Net without extrapolation and the proposed Momentum-Net,
we set N_iter= 100 and ρ= 1 - ε.
For PDS-Net, we used the identical parameter setup described above.
For performance comparisons between different INNs,
all the INNs used sCNN refiners sys:auto:res with { R = 5^2, K = 32 }
(to avoid the overfitting risks)
in the epipolar domain.
For Momentum-Net using dCNN refiners, we chose L=6 layer dCNN sys:dcnn using R = 3^2 filters and K = 16 feature maps.
(The chosen parameters gave most accurate performances over the following setups, { L=4, R = 3^2, K = 16, 32, 64 }, 
given the identical regularization parameters.)
To generate _θ^(i+1) ( x^(i) ) in eq:momnet:map, 
we applied a sCNN sys:auto:res with { R = 5^2, K = 32 } or 
a dCNN sys:dcnn with { L=6, R = 3^2, K = 16 }
to two sets of horizontal and vertical epipolar plane images, 
and took the average of two LFs that were permuted back from 
refined horizontal and vertical epipolar plane image sets, ∀ i.[
Epipolar images are 2D slices of a 4D LF L_F ( c_x, c_y, c_u, c_v ),
where ( c_x, c_y ) and ( c_u, c_v ) are spatial and angular coordinates, respectively.
Specifically, each horizontal epipolar plane image are obtained by
fixing c_y and c_v,
and varying c_x and c_u;
and each vertical epipolar image are obtained by
fixing c_x and c_u,
and varying c_y and c_v.
]



We designed the majorization matrices as {M^(i+1)=( A^T A 1 ) + γ I : i ≥ 0 }, 
using Lemma <ref> and setting λ= 1 by up:MFtilde.
We set an initial point of INNs, x^(0), to A^T y rescaled in the interval [0, 1] (i.e., dividing by its max value).
The regularization parameters 
(i.e., γ in BCD-Net/Momentum-Net, the ADMM penalty parameter in ADMM-Net, and the first step size in PDS-Net)
were selected by the proposed scheme in <ref>
with χ^⋆= 1.5.
(We tuned the factor to achieve the best performances).

For different combinations of INNs and sCNN refiner sys:auto:res/dCNN refiner sys:dcnn, we use the following naming convention: the INN name-sCNN or dCNN.





  §.§.§ Training INNs
 

For sparse-view CT experiments,
we trained all the INNs from the chest CT dataset with 
{ x_s, y_s, f_s(x;y_s) =1/2 y_s - A x _W_s^2, M_s: s = 1,…,S, S = 142 }; we constructed the dataset by using XCAT phantom slices <cit.>. 
The CT experiment has mild data-fit variations across training samples: the standard deviation of the condition numbers (≜σ_max(·) / σ_min (·)) of { M_f_s=( A^T W_s A 1 ) :∀ s } is 1.1.
For experiments of LF photography using a focal stack,
we trained all the INNs from the LF photography dataset with
{ y_s, f_s(x;y_s) =1/2 y_s - A_s x _2^2, M_s: s = 1,…,S, S = 21 } and 
two sets of ground truth epipolar images, 
{ x_s,epi-h ,  x_s,epi-v: s = 1,…,S, S = 21 · (255 · 9) };
we constructed the dataset by excluding four unrealistic stratified scenes from the original LF dataset in <cit.> 
that consists of 28 LFs with diverse scene parameter and camera settings.
The LF experiment has large data-fit variations across training samples: the standard deviation of the condition numbers of { M_f_s=( A_s^T A_s 1 ) :∀ s } is 2245.5.


In training INNs for both the applications,
if not specified, we used identical training setups.
At each iteration of INNs, we solved sys:auto:res:train with the mini-batch version of Adam <cit.>
and trained iteration-wise sCNNs sys:auto:res or dCNNs sys:dcnn.
We selected the batch size and the number of epochs as follows:
for sparse-view CT reconstruction, we chose them as 20 & 300, and 20 & 200 for sCNN and dCNN refiners, respectively;
for LF photography using a focal stack, we chose them as 200 & 200, and 200 & 100, for sCNN and dCNN refiners, respectively.
We chose the learning rates for filters in sCNNs and dCNNs, 
and thresholding values {α_k^(i+1):∀ k, i } in sCNNs sys:auto:res,
as 10^-3 and 10^-1, respectively;
we reduced the learning rates by 10% every 10 epochs.
At the first iteration, we initialized filter coefficients with Kaiming uniform initialization <cit.>;
in the later iteration, i.e., at the i INN iteration, for i ≥ 2, we initialized filter coefficients 
from those learned from the previous iteration, i.e., (i-1) iteration 
(this also applies to initializing thresholding values).







  §.§.§ Testing trained INNs
 

In sparse-view CT reconstruction experiments, we tested trained INNs to two samples where 
ground truth images and the corresponding inverse covariance matrices (i.e., W in <ref>)
sufficiently differ from those in training samples
(i.e., they are a few cm away from training images).
We evaluated the reconstruction quality by the most conventional error metric in CT application,
RMSE (in HU), in a region of interest (ROI), 
where RMSE and HU stand for root-mean-square error and (modified) Hounsfield unit,
respectively,
and the ROI was a circular region that includes all the phantom tissues.
The  is defined by RMSE (x^⋆, x^true) ≜ ( ∑_j=1^N_ROI(x_j^⋆-x^true_j)^2/N_ROI )^1/2, 
where x^⋆ is a reconstructed image, 
x^true is a ground truth image, 
and N_ROI is the number of pixels in a ROI.
In addition, we compared the trained Momentum-Net (using extrapolation) 
to a standard MBIR method using a hand-crafted EP regularizer,
and an MBIR model using a learned convolutional regularizer <cit.> 
which is the state-of-the-art MBIR method within an unsupervised learning setup.
We finely tuned their regularization parameters to achieve the lowest RMSE.
See details of these two MBIR models in <ref>.



In experiments of LF photography using a focal stack, we tested trained INNs to three samples 
of which scene parameter and camera settings are different from those in training samples 
(all training and testing samples have different camera and scene parameters).
We evaluated the reconstruction quality by the most conventional error metric in LF photography application,
PSNR (in dB), where PSNR stands for peak signal-to-noise.
In addition, we compared the trained Momentum-Net (using extrapolation) 
to MBIR method using the state-of-the-art non-trained regularizer, 4D EP introduced in <cit.>.
(The low-rank plus sparse tensor decomposition model <cit.> failed 
when inverse crimes and measurement noise are considered.)
We finely tuned its regularization parameter to achieve the lowest RMSE values.
See details of this MBIR model in <ref>.
We further investigated impacts of the LF MBIR quality on a higher-level depth estimation application,
by applying the robust Spinning Parallelogram Operator (SPO) depth estimation method <cit.> to reconstructed LFs.


For comparing Momentum-Net with PDS-Net, we measured quality of refined images, z^(i+1) in eq:momnet:map,
because PDS-Net is a hard-refiner.


The imaging simulation and reconstruction experiments were based on 
the Michigan image reconstruction toolbox <cit.>,
and training INNs, i.e., solving sys:auto:res:train, was based on PyTorch 
(for sparse-view CT, we used ver. 1.2.0; for LF photography using a focal stack, we used ver. 0.3.1).
For sparse-view CT, 
single-precision MATLAB  and PyTorch implementations were tested 
on 2.6 GHz Intel Core i7 CPU with 16 GB RAM, 
and 1405 MHz Nvidia Titan Xp GPU with 12 GB RAM, respectively.
For LF photography using a focal stack, 
they were tested on 3.5 GHz AMD Threadripper 1920X CPU with 32 GB RAM,
and 1531 MHz Nvidia GTX 1080 Ti GPU with 11 GB RAM, respectively.














 §.§ Comparisons between different INNs
 


First, compare sCNN results in Figs. <ref>–<ref>
and Figs. <ref>–<ref>, for 
sparse-view CT and LF photography using a focal stack, respectively. For both applications, the
proposed Momentum-Net using extrapolation significantly improves MBIR speed and accuracy,
compared to the existing soft-refining INNs,
<cit.> 
that correspond to BCD-Net <cit.> or Momentum-Net using no extrapolation,
and ADMM-Net <cit.>,
and the existing hard-refining INN PDS-Net <cit.>.
(Note that 
BCD-Net and Momentum-Net require slightly less computational complexity per INN iteration,
compared to ADMM-Net and PDS-Net, respectively, due to having fewer modules.)
Fig. <ref> shows that to reach the 24 HU RMSE value in sparse-view CT reconstruction, 
the proposed Momentum-Net decreases MBIR time by 53.3% and 62.5%,
compared to Momentum-Net without extrapolation and BCD-Net using three inner iterations, respectively.
Fig. <ref> shows that to reach the 32 dB PSNR value in LF reconstruction from a focal stack, 
the proposed Momentum-Net decreases MBIR time by 36.5% and 61.5%,
compared to Momentum-Net without extrapolation and BCD-Net using three inner iterations, respectively.
In addition, Figs. <ref> and <ref> show that using extrapolation, 
i.e., eq:momnet:exp with up:Ex:cvx–eq:mom_coeff,
improves the performance of Momentum-Net versus iterations.



We conjecture that the larger performance gap between soft-refiner Momentum-Net and hard-refiner PDS-Net, 
in Fig. <ref>(a) compared to Fig. <ref>(a),
is because the LF problem needs stronger regularization, 
i.e., a smaller tuned factor χ^⋆ in eq:reg:select, than the CT problem.
Similarly, comparing Fig. <ref>(b) to Fig. <ref>(b) shows that 
the LF problem has small performance gaps between BCD-Net and ADMM-Net.

 









For both the applications,
using dCNN refiners sys:dcnn instead of  sCNN refiners sys:auto:res has a negligible effect on total run time of Momentum-Net,
because reconstruction time of MBIR modules eq:momnet:mbir (in CPUs) dominates inference time of 
image refining modules eq:momnet:map (in GPUs).
Compare results between Momentum-Net-sCNN and -dCNN in 
Figs. <ref> & <ref> and Tables <ref> & <ref>.




 




 






 







 §.§ Comparisons between different MBIR methods
 

In sparse-view CT using 12.5% of the full projection views,
Fig. <ref>(b)–(e) and Table <ref>(b)–(f) show that 
the proposed Momentum-Net 
achieves significantly better reconstruction quality compared to 
the conventional EP MBIR method 
and the state-of-the-art MBIR method within an unsupervised learning setup, 
MBIR model using a learned convolutional regularizer <cit.>.
In particular, Momentum-Net recovers 
both low- and high-contrast regions (e.g., soft tissues and bones, respectively)
more accurately than MBIR using a learned convolutional regularizer; see Fig. <ref>(c)–(e).
In addition, when their shallow convolutional autoencoders need identical computational complexities,
Momentum-Net achieves much faster MBIR compared to MBIR using a learned convolutional regularizer;
see Table <ref>(c)–(d). 

In LF photography using five focal sensors,
regardless of scene parameters and camera settings,
Momentum-Net consistently achieves significantly more accurate image recovery,
compared to MBIR model using the state-of-the-art non-trained regularizer, 4D EP <cit.>.
The effectiveness of Momentum-Net is more evident for a scene with less fine details.
See Fig. <ref>(b)–(d) and Table <ref>(b)–(d).
Regardless of the scene distances from LF imaging systems, the reconstructed LFs by Momentum-Net significantly improve the depth estimation accuracy over those reconstructed by the state-of-the-art non-trained regularizer, 4D EP <cit.>.
See Fig. <ref>(c)–(e) and Table <ref>(c)–(e).

In general,
Momentum-Net needs more computations per iteration than EP MBIR, 
because its refining NNs use more and larger filters than the small finite-difference filters in EP MBIR,
and EP MBIR algorithms can be often further accelerated by gradient approximations, 
e.g., ordered-subsets methods <cit.>.











§ CONCLUSIONS
 

Developing rapidly converging INNs is important, 
because 1) it leads to fast MBIR 
by reducing the computational complexity
in calculating data-fit gradients or 
applying refining NNs, 
and 2) training INNs with many iterations requires long training time
or it is challenging when refining NNs are fixed across INN iterations.
The proposed Momentum-Net framework is applicable for a wide range of inverse problems, while achieving fast and convergent MBIR.
To achieve fast MBIR, 
Momentum-Net uses momentum in extrapolation modules, 
and noniterative MBIR modules at each iteration via majorizers.
For sparse-view CT and LF photography using a focal stack,
Momentum-Net achieves faster and more accurate MBIR compared to 
the existing soft-refining INNs,
<cit.> 
that correspond to BCD-Net <cit.> or Momentum-Net using no extrapolation,
and ADMM-Net <cit.>,
and the existing hard-refining INN PDS-Net <cit.>.
When an application needs strong regularization strength, e.g., LF photography using limited detectors,
using dCNN refiners with moderate depth significantly improves the MBIR accuracy of Momentum-Net compared to sCNNs,
only marginally increasing total MBIR time.
In addition, Momentum-Net guarantees convergence to a fixed-point 
for general differentiable (non)convex MBIR functions (or data-fit terms) and convex feasible sets, 
under some mild conditions and two asymptotic conditions.
The proposed regularization parameter selection scheme uses the spectral spread of majorization matrices,
and is useful to consider data-fit variations across training/testing samples.



There are a number of avenues for future work.
First, we expect to further improve performances of Momentum-Net (e.g., MBIR time and accuracy) by using sharper majorizer designs.
Second, we expect to further reduce MBIR time of Momentum-Net with the stochastic gradient perspective (e.g., ordered subset <cit.>).
On the regularization parameter selection side, 
our future work is learning the factor χ in eq:reg:select from datasets while training refining NNs.





























IEEEtran




















 













23.423.4Momentum-Net: Fast and convergent 
 
iterative neural network for inverse problems 
 
(Appendices)







This supplementary material for <cit.> 
a) reviews the Block Proximal Extrapolated Gradient method using a Majorizer (BPEG-M) 
<cit.>, 
b) lists parameters of Momentum-Net, and summarizes selection guidelines or gives default values,
c) compares the convergence properties between Momentum-Net and BCD-Net, and
d) provides mathematical proofs or detailed descriptions to support several arguments in
the main manuscript. 
We use the prefix A for the numbers
in section, theorem, equation, figure, table, and footnote in the supplement.






§ BPEG-M: REVIEW



This section explains block multi-(non)convex optimization problems,
and summarizes the state-of-the-art method for block multi-(non)convex optimization, 
BPEG-M <cit.>, 
along with its convergence guarantees.




 §.§ Block multi-(non)convex optimization
 

In a block optimization problem, the variables of the underlying optimization problem are treated either as a single block or multiple disjoint blocks. 
In block multi-(non)convex optimization, we consider the following problem:


min_u   F(u_1,…,u_B) ≜f(u_1,…,u_B) + ∑_b=1^B r_b (u_b)


where variable u is decomposed into B blocks u_1,…,u_B ({ u_b ∈^n_b : b=1,…,B }), 
f is assumed to be (continuously) differentiable,
but functions { r_b : b = 1,…,B } are not necessarily differentiable. 
The function r_b can incorporate the constraint u_b ∈_b, by allowing r_b's to be extended-valued, e.g., r_b (u_b) = ∞ if u_b ∉_b, for b=1,…,B.
It is standard to assume that both f and { r_b } are proper and closed, and the sets {_b } are closed.
We consider either that sys:multiConvx has block-wise convexity (but sys:multiConvx is jointly nonconvex in general) <cit.>
or that f, { r_b }, or {_b } are not necessarily convex <cit.>.
Importantly, r_b can include (non)convex and nonsmooth ℓ^p (quasi-)norm, p ∈ [0, 1].
The next section introduces our optimization framework that solves sys:multiConvx.


The following sections review BPEG-M <cit.>, the state-of-the-art optimization framework for solving block multi-(non)convex problems, when used with sufficiently sharp majorizers.
BPEG-M uses block-wise extrapolation, majorization, and proximal mapping.
By using a more general Lipschitz continuity (see Definition <ref>) for block-wise gradients, BPEG-M is particularly useful for rapidly calculating majorizers involved with large-scale problems, and successfully applied to some large-scale machine learning and computational imaging problems; see <cit.> and references therein.





















 §.§ BPEG-M



This section summarizes the BPEG-M framework.
Using Definition <ref> and Lemma <ref>, the proposed method, BPEG-M, is given as follows.
To solve sys:multiConvx, we minimize majorizers of F cyclically over each block u_1,…,u_B, while fixing the remaining blocks at their previously updated variables. Let u_b^(i+1) be the value of u_b after its i update, and define


f_b^(i+1)(u_b) ≜f ( u_1^(i+1), …, u_b-1^(i+1), u_b, u_b+1^(i), …, u_B^(i) ),

for all b,i.
At the b block of the i iteration, we apply Lemma <ref> to functional f_b^(i+1)(u_b) with a M^(i+1)-Lipschitz continuous gradient at the extrapolated point ú_b^(i+1), and minimize a majorized function.
In other words, we consider the updates





u_b^(i+1) 
   = _ u_b     ∇f_b^(i+1) (ú_b^(i+1))  u_b - ú_b^(i+1)  
+ 1/2  u_b - ú_b^(i+1) _M_b^(i+1)^2 
+ r_b (u_b)


   = Prox_r_b^M_b^(i+1)  (  ú_b^(i+1) - (  M_b^(i+1)  )^-1 ∇f_b^(i+1) (ú_b^(i+1)) _extrapolated gradient step using a majorizer of f_b^(i+1)  ), 


where 


ú_b^(i+1) = u_b^(i) + E_b^(i+1) ( u_b^(i) - u_b^(i-1) ),

the proximal operator is defined by eq:d:prox,
∇ f_b^(i+1) (ú_b^(i+1)) is the block-partial gradient of f at ú_b^(i+1), 
a scaled majorization matrix is given by


M_b^(i+1) = λ_b ·M_b^(i+1) ≻0,     λ_b ≥1,

and M_b^(i+1)∈^n_b × n_b is a symmetric positive definite majorization matrix of ∇ f_b^(i+1)(u_b).
In up:xacute, the ^n_b × n_b matrix E_b^(i+1)≽ 0 is an extrapolation matrix that accelerates convergence in solving block multi-convex problems <cit.>.
We design it to satisfy conditions cond:Eb:cvx or cond:Eb:ncvx below.
In up:Mtilde, {λ_b = 1 : ∀ b } and {λ_b > 1 : ∀ b }, for block multi-convex and block multi-nonconvex problems, respectively.



For some f_b^(i+1) having sharp majorizers, we expect that extrapolation up:xacute has no benefits in accelerating convergence, and use { E_b^(i+1) = 0 : ∀ i }.
Other than the blocks having sharp majorizers, one can apply some increasing momentum coefficient formula <cit.> to the corresponding extrapolation matrices.
The choice in <cit.> accelerated BPEG-M for some machine learning and data science applications.
Algorithm <ref> summarizes these updates.









 §.§ Convergence results
 

This section summarizes convergence results of Algorithm <ref> under the following assumptions:





  * Assumption S.1)
In sys:multiConvx, F is proper and lower bounded in (F) ≜{ u : F(u) < ∞}.
In addition,
[]


   for block multi-convex sys:multiConvx, f is differentiable and sys:multiConvx has a Nash point or block-coordinate minimizer[
Given a feasible set , a point u^⋆∈(F) ∪ is a critical point (or stationary point) of F if the directional derivative d^T ∇ F(u^⋆) ≥ 0 for any feasible direction d at u^⋆.
If u^⋆ is an interior point of , then the condition is equivalent to 0 ∈∂ F(u^⋆).
] (see its definition in <cit.>);


   for block multi-nonconvex sys:multiConvx, f is continuously differentiable, r_b is lower semicontinuous[
F is lower semicontinuous at point u_0 if lim inf_u → u_0 F(u) ≥ F(u_0).
], ∀ b, and sys:multiConvx has a critical point u^⋆ that satisfies 0 ∈∂ F(u^⋆).




  * Assumption S.2) ∇ f_b^(i+1) (u_b) is M-Lipschitz continuous with respect to u_b, i.e.,




 ∇f_b^(i+1) (u) - ∇f_b^(i+1) (v)  _(  M_b^(i+1)  )^-1 ≤u - v_M_b^(i+1),

for u,v ∈^n_b, where M_b^(i+1) is a bounded majorization matrix.



  * Assumption S.3) The extrapolation matrices E_b^(i+1)≽ 0 satisfy that




    
( E_b^(i+1) )^T M_b^(i+1) E_b^(i+1) ≼δ^2 ·M_b^(i);



    
( E_b^(i+1) )^T M_b^(i+1) E_b^(i+1) ≼δ^2 (λ_b - 1)^2/4 (λ_b + 1)^2 ·M_b^(i),

with δ < 1, ∀ b,i.




[Block multi-convex sys:multiConvx: A limit point is a Nash point ]

Under Assumptions S.1–S.3, let { u^(i+1) : i ≥ 0} be the sequence generated by Algorithm <ref>. Then any limit point of { u^(i+1) : i ≥ 0 } is a Nash point of sys:multiConvx.




[Block multi-nonconvex sys:multiConvx: A limit point is a critical point ]

Under Assumptions S.1–S.3, let { u^(i+1) : i ≥ 0} be the sequence generated by Algorithm <ref>.
Then any limit point of { u^(i+1) : i ≥ 0 } is a critical point of sys:multiConvx. 




Theorems <ref>–<ref> imply that, if there exists a critical point for sys:multiConvx, i.e., 0 ∈∂ F(u^⋆), then any limit point of { u^(i+1) : i ≥ 0 } is a critical point.
One can further show global convergence under some conditions: if { u^(i+1) : i ≥ 0 } is bounded and the critical points are isolated, then { u^(i+1) : i ≥ 0 } converges to a critical point <cit.>, <cit.>.









 §.§ Application of BPEG-M to solving block multi-(non)convex problem sys:recov   caol



For update eq:bpgm:zk, we do not use extrapolation, i.e., up:xacute, 
since the corresponding majorization matrices are sharp, so one obtains tight majorization bounds in Lemma <ref>. 
See, for example, <cit.>.
For updates eq:bpgm:zk and eq:bpgm:x, we rewrite ∑_k=1^K  h_k  x - ζ_k _2^2 as x - ∑_k=1^K flip(h_k^*)  ζ_k _2^2 
by using the TF condition in <ref> <cit.>, <cit.>.







§ EMPIRICAL MEASURES RELATED TO THE CONVERGENCE OF MOMENTUM-NET USING SCNN REFINERS
 

This section provides empirical measures related to Assumption 4 for Momentum-Net using single-hidden layer autoencoders sys:auto:res; 
see Fig. <ref> below. 
We estimated the sequence {ϵ^(i): i = 2,…,N_lyr} in Definition <ref>, 
 the sequence {Δ^(i): i = 2,…,N_lyr} in Definition <ref>, 
 and the Lipschitz constants {κ^(i): i = 1,…,N_lyr} of refining NNs {_θ^(i):∀ i }, 
 based on a hundred sets of randomly selected training samples related with the corresponding bounds of the measures, 
 e.g., u and v in cond:eps are training input to _θ^(i+1) and _θ^(i) in eq:momnet:map, respectively.

 







§ PROBABILISTIC JUSTIFICATION FOR THE ASYMPTOTIC BLOCK-COORDINATE MINIMIZER CONDITION IN ASSUMPTION 4
 



This section introduces a useful result for an asymptotic block-coordinate minimizer z^(i+1):
the following lemma provides a probabilistic bound for x^(i) - z^(i+1)_2^2 in cond:delta, 
given a subgaussian vector z^(i+1) - z^(i) with independent and zero-mean entries.




[Probabilistic bounds for x^(i) - z^(i+1)_2^2]

Assume that z^(i+1) - z^(i) is a zero-mean subgaussian vector of which entries are independent and zero-mean subgaussian variables.
Then, each bound in cond:delta holds with probability at least 

1 - exp (  - (  z^(i+1) - z^(i) _2^2 + Δ^(i+1) )^2 /8ρ·σ^(i+1) ·  _θ^(i+1) (x^(i)) - x^(i) _2^2 ),

where σ^(i+1) is a subgaussian parameter for z^(i+1) - z^(i), 
and a random variable is subgaussian with parameter σ if  { | · | ≥ t }≤ 2 exp( - t^2/2σ ) for t ≥ 0.


First, observe that

 x^(i) - z^(i+1) _2^2 
   =  x^(i) - z^(i) - ( z^(i+1) - z^(i) ) _2^2


   =   x^(i) - z^(i) _2^2 +  z^(i+1) - z^(i) _2^2 - 2 x^(i) - z^(i)z^(i+1) - z^(i)


   =   x^(i) - z^(i) _2^2 +  z^(i+1) - z^(i) _2^2 - 2  z^(i+1) - z^(i) + ρ(x^(i) - _θ^(i+1) (x^(i)) )  z^(i+1) - z^(i) 
 


   =   x^(i) - z^(i) _2^2 -  z^(i+1) - z^(i) _2^2 + 2 ρ _θ^(i+1) (x^(i)) - x^(i)  z^(i+1) - z^(i) 
 

where the inequality eq:l:delta:prob:1 holds by x^(i) = ρ x^(i) - ρ_θ^(i+1) + z^(i+1) via eq:momnet:map.
We now obtain a probablistic bound for the third quantity in eq:l:delta:prob:2 via a concentration inequality. 
The concentration inequality on the sum of independent zero-mean subgaussian variables (e.g., <cit.>) yields that for any t^(i+1)≥ 0

 
{   _θ^(i+1) (x^(i)) - x^(i) z^(i+1) - z^(i)  ≥t^(i+1) }
≤exp( -  ( t^(i+1) )^2 / 2 σ^(i+1)   _θ^(i+1) (x^(i)) - x^(i) _2^2  )

where σ^(i+1) is given as in Lemma <ref>.
Applying the result eq:l:delta:prob:3 with t^(i+1) = 1/2ρ (  z^(i+1) - z^(i)_2^2 + Δ^(i+1) ) to the bound eq:l:delta:prob:2 completes the proofs.



Lemma <ref> implies that, 
given sufficiently large Δ^(i+1), 
or sufficiently small σ^(i+1) (e.g., variance for a Gaussian random vector z^(i+1) - z^(i)) 
or _θ^(i+1) (x^(i)) - x^(i)_2^2, 
bound cond:delta is satisfied with high probability, for each i.
In particular, 
Δ^(i+1) can be large for the first several iterations;
if paired operators ( _θ^(i+1), _θ^(i) ) in eq:momnet:map map their input images to similar output images 
(e.g., the trained NNs _θ^(i+1) and _θ^(i) have good refining capabilities for x^(i) and x^(i-1)), 
then σ^(i+1) is small;
if the regularization parameter γ in eq:momnet:mbir is sufficiently large, then _θ^(i+1) (x^(i)) - x^(i)_2^2 is small.








§ PROOFS OF PROPOSITION <REF>



First, we show that ∑_i=0^∞ x^(i+1) - x^(i)_2^2 < ∞ for convex and nonconvex F(x;y,z^(i+1)) cases. 



  * Convex F(x;y,z^(i+1)) case:
Using Assumption 2 and {M^(i+1)= M^(i+1):∀ i } for the convex case via up:MFtilde, we obtain the following results for any :



     F ( x^(i); y,  z^(i) ) - F ( x^(i+1); y,  z^(i+1) ) + γΔ^(i+1)


   ≥F ( x^(i); y,  z^(i+1) ) - F ( x^(i+1); y,  z^(i+1) ) 



   ≥1/2  x^(i+1) - x́^(i+1)  _M^(i+1)^2 + ( x́^(i+1) - x^(i) )^T M^(i+1) ( x^(i+1) - x́^(i+1) )
 


    = 1/2  x^(i+1) -  x^(i) _M^(i+1)^2 - 1/2  E^(i+1) ( x^(i) -  x^(i-1) ) _M^(i+1)^2
 


   ≥1/2  x^(i+1) -  x^(i) _M^(i+1)^2 - δ^2/2  x^(i) -  x^(i-1) _M^(i)^2  



where the inequality eq:p:momnet-convg-cnvx:1 uses the condition cond:delta in Assumption 4,
the inequality eq:p:momnet-convg-cnvx:2 is obtained by using the results in <cit.>,
the equality eq:p:momnet-convg-cnvx:3 uses the extrapolation formula eq:momnet:exp and the symmetry of M^(i+1),
the inequality eq:p:momnet-convg-cnvx:4 holds by Assumption 3.


Summing the inequality of F(x^(i); y,z^(i)) - F(x^(i+1); y, z^(i+1)) + γΔ^(i+1) in eq:p:momnet-convg-cnvx:4 over i=0,…,N_lyr-1, we obtain




F ( x^(0); y, z^(0) )  - F ( x^(N_lyr); y, z^(N_lyr) ) 
+ γ∑_i=0^N_lyr-1 Δ^(i+1)
   ≥∑_i=0^N_lyr-1 1/2  x^(i+1) -  x^(i) _M^(i+1)^2 - δ^2/2  x^(i) -  x^(i-1) _M^(i)^2 


   ≥∑_i=0^N_lyr-1 1-δ^2/2   x^(i+1) -  x^(i) _M^(i+1)^2 


   ≥∑_i=0^N_lyr-1 m_F,min (1-δ^2)/2  x^(i+1) -  x^(i) _2^2 


where the inequality eq:p:momnet-convg-cnvx:7 holds by Assumption 2.
Due to the lower boundedness of F(x;y,z) in Assumption 1 and the summability of {Δ^(i+1)≥ 0 : ∀ i } in Assumption 4,
taking N_lyr→∞ gives


∑_i=0^∞  x^(i+1) - x^(i)_2^2 < ∞.





  * Nonconvex F(x;y,z^(i+1)) case:
Using Assumption 2, we obtain the following results without assuming that F(x;y,z^(i+1)) is convex:

      F ( x^(i); y,  z^(i) ) - F ( x^(i+1); y,  z^(i+1) ) + γΔ^(i+1)


   ≥F ( x^(i); y,  z^(i+1) ) - F ( x^(i+1); y,  z^(i+1) ) 
 


   ≥λ- 1/4  x^(i+1) - x^(i) _M^(i+1)^2 -  (λ+ 1)^2 /λ- 1  x^(i) - x́_b^(i+1) _M^(i+1)^2 
 


   = λ-1/4  x^(i+1) -  x^(i) _M^(i+1)^2 - (λ+1)^2/λ-1  E^(i+1) ( x^(i) -  x^(i-1) ) _M^(i+1)^2
 


   ≥λ-1/4 (  x^(i+1) -  x^(i) _M^(i+1)^2 - δ^2  x^(i) -  x^(i-1) _M^(i)^2 )
 

where the inequality eq:p:momnet-convg-ncnvx:0 uses the condition cond:delta in Assumption 4,
the inequality eq:p:momnet-convg-ncnvx:1 use the results in <cit.>,
the equality eq:p:momnet-convg-ncnvx:2 holds by eq:momnet:exp,
the inequality eq:p:momnet-convg-ncnvx:3 is obtained by Assumption 3.


Summing the inequality of F(x^(i); y,z^(i)) - F(x^(i+1); y, z^(i+1)) + γΔ^(i+1) in eq:p:momnet-convg-ncnvx:3 over i=0,…,N_lyr-1, we obtain




F ( x^(0); y, z^(0) )  - F ( x^(N_lyr); y, z^(N_lyr) ) 
+ γ·∑_i=0^N_lyr-1 Δ^(i+1) 
   ≥∑_i=0^N_lyr-1 λ-1/4 (  x^(i+1) -  x^(i) _M^(i+1)^2 - δ^2  x^(i) -  x^(i-1) _M^(i)^2 )


   ≥∑_i=0^N_lyr-1 (λ-1)(1-δ^2)/2  x^(i+1) -  x^(i) _M^(i+1)^2


   ≥∑_i=0^N_lyr-1 m_F,min (λ-1) (1-δ^2)/2  x^(i+1) - x^(i) _2^2,

where we follow the arguments in obtaining eq:p:momnet-convg-cnvx:7 above.
Again, using the lower boundedness of F(x;y,z) and the summability of {Δ^(i+1)≥ 0 : ∀ i },
taking N_lyr→∞ gives the result eq:p:momnet:convg:sqSum:u for nonconvex F(x;y,z^(i+1)). 



Second, we show that ∑_i=0^∞ z^(i+1) - z^(i)_2^2 < ∞. Observe

z^(i+1) - z^(i)_2^2
   =  (1-ρ) ( x^(i) - x^(i-1) ) + ρ( _θ^(i+1) (x^(i)) - _θ^(i) (x^(i-1)) ) _2^2


   ≤(1-ρ)  x^(i) - x^(i-1) _2^2 + ρ _θ^(i+1) (x^(i)) - _θ^(i) (x^(i-1)) _2^2


   ≤ x^(i) - x^(i-1) _2^2 + ρϵ^(i+1)


where the first equality uses the image mapping formula in eq:momnet:map,
the first inequality holds by applying Jensen's inequality to the (convex) squared ℓ^2-norm,
the second inequality is obtained by using the asymptotically non-expansiveness of the paired operators ( _θ^(i+1),  _θ^(i) ) in Assumption 4. 
Summing the inequality of z^(i+1) - z^(i)_2^2 in eq:p:momnet-convg:mapdiff over i=0,…,N_lyr-1, we obtain


∑_i=0^N_lyr-1  z^(i+1) - z^(i)_2^2 
   ≤∑_i=0^N_lyr-2  x^(i+1) - x^(i) _2^2 + ρ∑_i=0^N_lyr-1 ϵ^(i+1),

where we used x^(0) = x^(-1) as given in Algorithm <ref>.
By taking N_lyr→∞ in eq:p:momnet:sqSum:u0, using result eq:p:momnet:convg:sqSum:u, and the summability of the sequence {ϵ^(i+1) : i ≥ 0 }, we obtain


∑_i=0^∞  z^(i+1) - z^(i)_2^2 < ∞.

Combining the results in eq:p:momnet:convg:sqSum:u and eq:p:momnet:convg:sqSum:u0 completes the proofs.








§ PROOFS OF THEOREM <REF>



Let x̅ be a limit point of { x^(i) : i ≥ 0 } and { x^(i_j)} be the subsequence converging to x̅.
Let z̅ be a limit point of { z^(i) : i ≥ 0 } and { z^(i_j)} be the subsequence converging to z̅.
The closedness of  implies that x̅∈.
Using the results in Proposition <ref>, 
{ x^(i_j+1)} and { z^(i_j+1)} also converge to x̅ and z̅, respectively.
Taking another subsequence if necessary, the subsequence { M^(i_j+1)} converges to some M̅, since M^(i+1) is bounded by Assumption 2.
The subsequences {θ^(i_j+1)} converge to some θ̅, since x^(i_j+1)→x̅, z^(i_j+1)→z̅, and {θ^(i+1)} is bounded via Assumption 4. 




Next, we show that the convex proximal minimization eq:t:momnet:convg:cvx:1 below is continuous in the sense that the output point x^(i_j+1) continuously depends on the input points x́^(i_j+1) and z^(i_j+1), and majorization matrix M^(i_j+1):

x^(i_j+1)    = _ x ∈   ∇F(x́^(i_j+1); y, z^(i_j+1))x - x́^(i_j+1) + 1/2 x - x́^(i_j+1)_M^(i_j+1)^2 



   = Prox__^M^(i_j+1)  ( x́^(i_j+1) - ( M^(i_j+1) )^-1 ∇F (x́^(i_j+1); y, z^(i_j+1)) ).


where the proximal mapping operator Prox__^M^(i_j+1) (·) is given as in eq:d:prox.
We consider the two cases of majorization matrices { M^(i+1)} given in Theorem <ref>:


  * For a sequence of diagonal majorization matrices, i.e., { M^(i+1) : i ≥ 0 },
one can obtain the continuity of the convex proximal minimization eq:t:momnet:convg:cvx:1 with respect to x́^(i_j+1), z^(i_j+1), and M^(i_j+1), 
by extending the existing results in <cit.>, <cit.> with the separability of eq:t:momnet:convg:cvx:1 to element-wise optimization problems.



  * For a fixed general majorization matrix, i.e., M = M^(i+1), ∀ i, we obtain that the convex proximal minimization eq:t:momnet:convg:cvx:M:1 below is continuous with respect to the input points x́^(i_j+1) and z^(i_j+1):

x^(i_j+1) 
   = _ x ∈   ∇F(x́^(i_j+1); y, z^(i_j+1))x - x́^(i_j+1) + 1/2 x - x́^(i_j+1)_M^2
 

 
   = Prox__^M  ( x́^(i_j+1) - M^-1 ∇F (x́^(i_j+1); y, z^(i_j+1)) )


   = ( Id + M^-1 ∂̂ _ )^-1 ( x́^(i_j+1) - M^-1 ∇F (x́^(i_j+1); y, z^(i_j+1)) )


where ∂̂ f(x) is the subdifferential of f at x and Id denotes the identity operator, 
and the proximal mapping of _ relative to ·_M is uniquely determined by 
the resolvent of the operator M^-1∂̂_ in eq:t:momnet:convg:cvx:M:2.


First, we obtain that the operator M^-1∂̂_ is monotone.
For a convex extended-valued function f_e: ^N →∪{∞}, observe that M^-1∂̂ f_e is a monotone operator:

 M^-1 ∂̂ f_e (u) - M^-1 ∂̂ f_e (v)  u - v 
=   M^-1 M _= I ∂̂ f_e (M ũ) -  M^-1 M _=I ∂̂ f_e (M ṽ)  M ũ - M ṽ  
≥0,     ∀u, v,


where the equality uses the variable change { u = Mũ, v = Mṽ},  
a chain rule of the subdifferential of a composition of a convex extended-valued function 
and an affine mapping <cit.>, and the symmetry of M,
and the inequality holds because the subdifferential of convex extended-valued function is a monotone operator <cit.>.
Because characteristic function of a convex set is extended-valued function, 
the result in eq:t:momnet:convg:monotone implies that the operator M^-1∂̂_ is monotone.
Second, note that the resolvent of a monotone operator M^-1∂̂_ (with a parameter 1), i.e., 
( Id + M^-1∂̂_ )^-1 in eq:t:momnet:convg:cvx:M:2, 
is nonexpansive <cit.> and thus continuous.
We now obtain that the convex proximal minimization eq:t:momnet:convg:cvx:M:1 is continuous with respect to the input points x́^(i_j+1) and z^(i_j+1), 
because the proximal mapping operator ( Id + M^-1∂̂_ )^-1 in eq:t:momnet:convg:cvx:M:2, the affine mapping M^-1, and ∇ F(x;y,z) are continuous with respect to their input points.








For the two cases above, 
using the fact that x^(i_j+1)→x̅, x́^(i_j+1)→x̅, z^(i_j+1)→z̅, and M^(i_j+1)→M̅ (or M̅ = M for the { M^(i+1) = M } case) as j →∞, eq:t:momnet:convg:cvx:1 becomes


x̅ = _ x ∈   ∇F(x̅; y,z̅)x - x̅ + 1/2 x - x̅_M̅^2 .

Thus, x̅ satisfies the first-order optimality condition of min_x ∈ F(x;y,z̅):

∇F (x̅;y, z̅)  x - x̅  ≥0,     ,

and this completes the proof of the first result. 


Next, note that the result in Proposition <ref> imply


 __θ^(i+1)^M^(i+1)  ( [ =1.5pt [   x^(i); x^(i-1) ] ] ) - [ =1.5pt [   x^(i); x^(i-1) ] ] _2 →0.

Additionally, note that a function __θ^(i+1)^M^(i+1) - I is continuous.
To see this, observe that the convex proximal mapping in eq:momnet:mbir is continuous (see the obtained results above), and _θ^(i+1) is continuous (see Assumption 4).
Combining t:momnet:convg:cvx:convg_to0, 
the convergence of { M^(i_j+1), _θ^(i_j+1)}, 
and the continuity of __θ^(i+1)^M^(i+1) - I, 
we obtain [x̅^T , x̅^T]^T =__θ̅^M̅ ( [x̅^T , x̅^T]^T ), and this completes the proofs of the second result.








§ PROOFS OF COROLLARY <REF>




To prove the first result, we use proof by contradiction. 
Suppose that dist (x^(i), ) ↛ 0.
Then there exists ϵ > 0 and a subsequence { x^(i_j)} such that dist (x^(i_j), ) ≥ϵ, ∀ j.
However, the boundedness assumption of { x^(i_j)} in Corollary <ref> implies that there must exist a limit point x̅∈ via Theorem <ref>.
This is a contradiction, and gives the first result (via the result in Proposition <ref>).
Under the isolation point assumption in Corollary <ref>, using the obtained results, x^(i+1) - x^(i)_2 → 0 (via Proposition <ref>) and dist (x^(i+1), ) → 0, and the following the proofs in <cit.>, we obtain the second result.







§ MOMENTUM-NET VS. BCD-NET



This section compares the convergence properties of Momentum-Net (Algorithm <ref>)  
and BCD-Net (Algorithm <ref>).
We first show that for convex f(x;y) and , 
the sequence of reconstructed images generated by BCD-Net converges:

[Sequence convergence]

In Algorithm <ref>, let f(x;y) be convex and subdifferentiable, and  be convex.
Assume that the paired operators ( _θ^(i+1),  _θ^(i) ) are asymptotically contractive, i.e.,

   _θ^(i+1) (u) -  _θ^(i)  (v) _2 <  u - v_2 + ϵ^(i+1),

with ∑_i=0^∞ϵ^(i+1) < ∞ and {ϵ^(i+i)∈ [0, ∞) : ∀ i }, ∀ u,v, i.
Then, the sequence { x^(i+1) : i ≥ 0 } generated by Algorithm <ref> is convergent.


We rewrite the updates in Algorithm <ref> as follows:

x^(i+1) 
    = _x ∈   f(x;y) + γ/2  x -  _θ^(i+1) ( x^(i) ) _2^2 
= Prox_f + _^γI  (  _θ^(i+1) ( x^(i) ) )


   =  ( Id + γ^-1 ∂̂ ( f(x;y) + _ ) )^-1 (  _θ^(i+1) ( x^(i) ) )


    =: ^(i+1) (x^(i)).

We first show that the paired operators {^(i+1), ^(i)} is asymptotically contractive:

     ^(i+1) (u) - ^(i) (v) _2


   =  ( Id + γ^-1 ∂̂ ( f(x;y) + _ ) )^-1 (  _θ^(i+1) ( u ) ) - ( Id + γ^-1 ∂̂ ( f(x;y) + _ ) )^-1 (  _θ^(i) ( v ) ) _2 


   ≤  _θ^(i+1) ( u ) -  _θ^(i) ( v ) _2



   ≤L'  u - v _2 + ϵ^(i+1)  u - v _2,


∀ u,v, 
where the inequality eq:p:bcdnet:nonexp holds because 
the subdifferential of the convex extended-valued function f(x;y) + _
(the characteristic function of a convex set , _, is convex, and the sum of the two convex functions, f(x;y) + _, is convex) 
is a monotone operator <cit.>, 
and the resolvent of a monotone relation with a positive parameter, 
i.e., ( Id + γ^-1∂̂ ( f(x;y) + _ ) )^-1 with γ^-1 > 0,
is nonexpansive <cit.>,
and the inequality eq:p:bcdnet:lipschitz holds by L' < 1 via the contractiveness of the paired operators ( _θ^(i+1),  _θ^(i) ), ∀ i.
Note that the inequality eq:p:bcdnet:nonexp does not hold for nonconvex f(x;y) and/or .
Considering that L' < 1, we show that the sequence { x^(i+1) : i ≥ 0 } is Cauchy sequence:

 x^(i+l) - x^(i) _2 
   =  ( x^(i+l) - x^(i+l-1)  ) + …+ ( x^(i+1) - x^(i) ) _2


   ≤ x^(i+l) - x^(i+l-1)   + …+  x^(i+1) - x^(i) _2


   ≤( L'^l-1 + …+ 1 )  x^(i+1) - x^(i)  _2 + ( ϵ^(i+l) + …+ ϵ^(i+1) )


   ≤1/1 - L'  x^(i+1) - x^(i)  _2 + ∑_i'=1^l ϵ^(i+i')

where the second inequality uses the result in eq:p:bcdnet:lipschitz. 
Since the sequence { x^(i+1) : i ≥ 0 } is Cauchy sequence, { x^(i+1) : i ≥ 0 } is convergent, and this completes the proofs.



In terms of guaranteeing convergence, 
BCD-Net has three theoretical or practical limitations compared to Momentum-Net:






  * Different from Momentum-Net, BCD-Net assumes the asymptotic contractive condition for the paired operators {_θ^(i+1),  _θ^(i)}. 
When image mapping operators in eq:bcdnet:map are identical across iterations, i.e., {_θ = _θ^(i+1) : i ≥ 0 }, then _θ is assumed to be contractive.
On the other hand, a mapping operator (identical across iterations) of Momentum-Net only needs to be nonexpansive.
Note, however, that when f(x;y) = 1/2 y - Ax _W^2 with A^H W A ≻ 0 (e.g., Example <ref>), 
BCD-Net can guarantee the sequence convergence with the asymptotically nonexpansive paired operators ( _θ^(i+1),  _θ^(i) ) 
(see Definition <ref>) <cit.>.


  * When one applies an iterative solver to eq:bcdnet:recon, 
there always exist some numerical errors and these obstruct the sequence convergence guarantee in Proposition <ref>. To guarantee sufficiently small numerical errors from iterative methods solving eq:bcdnet:recon (so that one can find a critical point solution for the MBIR problem eq:bcdnet:recon), one needs to use sufficiently many inner iterations that can substantially slow down entire MBIR.


  * BCD-Net does not guarantee the sequence convergence for nonconvex data-fit f(x;y), whereas Momentum-Net guarantees convergence to a fixed-point for both convex f(x;y) and nonconvex f(x;y).








§ FOR THE SCNN ARCHITECTURE SYS:AUTO:RES, CONNECTION BETWEEN CONVOLUTIONAL TRAINING LOSS SYS:AUTO:RES:TRAIN AND ITS PATCH-BASED TRAINING LOSS
 

This section shows that given the sCNN architecture sys:auto:res,
the convolutional training loss in sys:auto:res:train has three advantages over 
the patch-based training loss in <cit.>
that may use all the extracted overlapping patches of size R:


  * The corresponding patch-based loss does not model the patch aggregation process that is inherently modeled in sys:auto:res.

  * It is an upper bound of the convolutional loss sys:auto:res:train.

  * It requires about R times more memory than sys:auto:res:train.

We prove the benefits of sys:auto:res:train using the following lemma.




The loss function sys:auto:res:train for training the residual convolutional autoencoder in sys:auto:res is bounded by the patch-based loss function:


1/2L ∑_s=1^S  x̂_s^(i) -  1/R ∑_k=1^K d_k ⊛_α_k ( e_k ⊛x_s^(i) ) _2^2 
≤1/2LR ∑_s=1^S  X_s^(i)  - D _α̃ ( E X_s^(i) ) _F^2,

where the residual is defined by x̂_s^(i)≜ x_s - x_s^(i), 
{ x_s, x_s^(i)} are given as in sys:auto:res:train, 
X_s∈^R × V_s and X_s∈^R × V_s are 
the l training data matrices whose columns are V_s vectorized patches extracted from 
the images x̂_s and x_s (with the circulant boundary condition and the stride parameter 1), respectively,
D ≜ [d_1, …, d_K] ∈^R × K is a decoding filter matrix,
and E ≜ [e_1^*, …, e_K^*]^H ∈^K × R is an encoding filter matrix.
Here, the definition of soft-thresholding operator in eq:soft-threshold is generalized by


( _α̃ (u) )_k ≜{ [ u_k - α_k ·(u_k),    | u_k | > α_k,;                0,                 , ] .

for K = 1,…,K, where α̃ = [α_1,…,α_K]^T.
See other related notations in sys:auto:res.



First, we have the following reformulation <cit.>:


[ [ e_1 u;     ⋮; e_K u ] ] =  P' [ [ E P_1;     ⋮; E P_N ] ] _ u,   ∀u,

where P' ∈^KN × KN is a permutation matrix, E is defined in Lemma l:train, and P_n ∈^R × N is the n patch extraction operator for n = 1,…,N.
Considering that 1/R∑_k=1^K 𝖿𝗅𝗂𝗉(e_k^*) ⊛ ( e_k ⊛ u ) = E^H E u via the definition of E in eq:l:train:reform (see also the reformulation technique in <cit.>), 
we obtain the following reformulation result:


1/R ∑_k=1^K 𝖿𝗅𝗂𝗉(e_k^*) ⊛_α_k ( e_k ⊛x_s^(i) ) =  1/R ∑_n=1^N P_n^H E^H _α ( E P_n x_s^(i) ) 

where the soft-thresholding operators {_α_k(·) : ∀ k } and _α̃(·) are defined in eq:soft-threshold:gen and
we use the permutation invariance of the thresholding operator _α (·), i.e., _α ( P (·) ) =  P ·_α ( · ) for any α.
Finally, we obtain the result in eq:l:train as follows:


1/2L ∑_s=1^S  x̂_s^(i) -  1/R ∑_k=1^K d_k ⊛_α_k ( e_k ⊛x_s^(i) ) _2^2 
   = 1/2L  ∑_s=1^S  x̂_s^(i) -  1/R ∑_n=1^N P_n^H D _α̃ ( E P_n x_s^(i) )  _2^2


 
    = 1/2L  ∑_s=1^S  1/R ∑_n=1^N P_n^H P_n x̂_s^(i) -  1/R ∑_n=1^N P_n^H D _α̃ ( E P_n x_s^(i) )  _2^2


    =  1/2LR^2 ∑_s=1^S  ∑_n=1^N P_n^H ( x̂_l,n^(i) - D _α̃ ( E x_l,n^(i) ) ) _2^2


 
    ≤1/2LR ∑_s=1^S ∑_n=1^N  x̂_l,n^(i) - D _α̃ ( E x_l,n^(i) ) _2^2  


   =  1/2LR ∑_s=1^S  X_s^(i) - D _α̃ ( E X_s^(i) )  _F^2,

where D is defined in Lemma <ref>,
{x̂_l,n^(i) = P_n x̂_s^(i)∈^R, x_l,n^(i) = P_n x_s^(i)∈^R: n = 1,…,N } is a set of extracted patches, 
the training matrices {X_s^(i), X_s^(i)} are defined by X_s^(i)≜ [ x̂_l,n^(i), …, x̂_l,N^(i) ] 
and X_s^(i)≜ [ x_l,1^(i), …, x_l,N^(i) ]. 
Here, the equality eq:reform:dist:1 uses the result in eq:reform:idt, 
the equality eq:reform:dist:2 holds by ∑_n=1^N P_n^H P_n = R · I (for the circulant boundary condition in Lemma <ref>), 
and the inequality eq:reform:dist:3 holds by PP^H ≼ R · I with P≜ [P_1^H ⋯ P_N^H]^H.



Lemma <ref> reveals that
when the patch-based training approach extract all the R-size overlapping patches,
1) the corresponding patch-based loss is an upper bound of the convolutional loss sys:auto:res:train;
2) it requires about R-times larger memory than sys:auto:res:train 
because V_s≈ R N_s for x ∈^N_s and the boundary condition described in Lemma <ref>, ∀ l;
and 3) it misses modeling the patch aggregation process that is inherently modeled in sys:auto:res
–  see that the patch aggregation operator ∑_n=1^N P_n^H (·)_n is removed in the inequality eq:reform:dist:3 in the proof of Lemma <ref>.
In addition, different from the patch-based training approach <cit.>, 
i.e., training with the function on the right-hand side in eq:l:train, 
one can use different sizes of filters { e_k, d_k : ∀ k } 
in the convolutional training loss,
i.e., the function on the left-hand side in eq:l:train.








§ DETAILS OF EXPERIMENTAL SETUP




 §.§ Majorization matrix designs for quadratic data-fit


For (real-valued) quadratic data-fit f(x;y) in the form of 1/2 y -  A x _W^2, 
if a majorization matrix M exists such that A^H W A ≼ M, 
it is straightforward to verify that the gradient of quadratic data-fit f(x;y) satisfies the M-Lipchitz continuity in Definition <ref>, i.e.,

 ∇f (u;y) - ∇f (v;y) _M^-1 
=   A^H W A u - A^H W A v _M^-1
≤ u - v _M^2,   ∀u,v ∈^N.

because the assumption A^T W A ≼ M ⇔ M^-1/2 A^T W A M^-1/2≼ I implies that 
the eigenspectrum of M^-1/2 A^T W A M^-1/2 lies in the interval [0, 1], 
and gives the following result:

( M^-1/2 A^T W A M^-1/2 )^2 ≼I 
⇔( A^T W A  ) M^-1 ( A^T W A ) ≼M.


Next, we review a useful lemma in designing majorization matrices for a wide class of quadratic data-fit f(x;y):

[<cit.>]

For a (possibly complex-valued) matrix A and a diagonal matrix W with non-negative entries, A^H W A  ≼( | A^H | W | A | 1 ), where | A | denotes the matrix consisting of the absolute values of the elements of A.







 §.§ Parameters for MBIR optimization models: Sparse-view CT reconstruction
 


For MBIR model using EP regularization, 
we combined a EP regularizer ∑_n =1^N∑_n' ∈_nι_nι_n'φ(x_n - x_n')
and the data-fit f(x;y) in <ref>, 
where _n is the set of indices of the neighborhood, 
ι_n and ι_n' are parameters that encourage uniform noise <cit.>,
and φ(·) is the Lange penalty function, i.e., φ(t) = δ^2 (  | t/δ | - log(1+| t/δ |) ), with δ= 10 in HU.
We chose the regularization parameter (e.g., γ in sys:mbir:z) as 2^15.5.
We ran the relaxed linearized augmented Lagrangian method <cit.>
with 100 iterations and 12 ordered-subsets,
and initialized the EP MBIR algorithms with a conventional FBP method using a Hanning window.



For MBIR model using a learned convolutional regularizer <cit.>, 
we trained convolutional regularizer with filters of { h_k ∈^R : R = K = 7^2 } 
via CAOL <cit.> in an unsupervised training manner;
see training details in <cit.>.
The regularization parameters (e.g., γ in sys:recov   caol) were selected
by applying the spectral spread based selection scheme in <ref>
with the tuned factor χ^⋆= 167.64.
We selected the spatial-strength-controlling hard-thresholding parameter 
(i.e., α' in <cit.>) as follows:
for Test samples #1–2, we chose it is as 10^-10 and 6^-11, respectively.
We initialized the MBIR model using a learned regularizer with the EP MBIR results obtained above.
We terminated the iterations if the relative error stopping criterion (e.g., <cit.>) is met before reaching the maximum number of iterations.
We set the tolerance value as 10^-13 and the maximum number of iterations to 4 × 10^3.







 §.§ Parameters for MBIR optimization models: LF photography using a focal stack
 

For MBIR model using 4D EP regularization <cit.>, 
we combined a 4D EP regularizer ∑_n =1^N∑_n' ∈_nφ(x_n - x_n')
and the data-fit f(x;y) in <ref>, 
where _n is the set of indices of the 4D neighborhood,
and φ(·) is the hyperbola penalty function, i.e., φ(t) = δ^2 ( √(1+ | t/δ |^2) - 1 ).
We selected the hyperbola function parameter δ and 
regularization parameter (e.g., γ in sys:mbir:z) as follows: 
for Test samples #1–3, we chose them as 
{δ= 10^-4, γ= 10^3 }, 
{δ= 10^-1, γ= 10^7 }, 
and {δ= 10^-1, γ= 5 × 10^3 }, respectively.
We ran the conjugate gradient method with 100 iterations,
and initialized the 4D EP MBIR algorithms with A^T y rescaled in the interval [0, 1].






 §.§ Reconstruction accuracy and depth estimation accuracy of different MBIR methods
 

Tables <ref>–<ref> below provide reconstruction accuracy numerics of different MBIR methods in sparse-view CT reconstruction 
and LF photography using a focal stack,
and reports the SPO depth estimation <cit.> accuracy numerics on reconstructed LFs from different MBIR methods.




















 §.§ Reconstructed images and estimated depths with noniterative analytical methods
 


This section provides reconstructed images by an analytical back-projection method in sparse-view CT reconstruction 
and LF photography using a focal stack (see the first two columns in Fig. <ref>),
and estimated depths from reconstructed LFs via the SPO depth estimation method <cit.> (see the third column in Fig. <ref>(c)).
Results in Fig. <ref> below are supplementary to Fig. <ref>, Fig. <ref>, and Fig. <ref>,
and the first two columns visualize initial input images to INN methods.


 







§ HOW TO CHOOSE PARAMETERS OF IMAGE REFINING MODULES IN SOFT-REFINING INNS?




In soft-refining INNs using iterative-wise refining NNs, 
one does not need to greatly increase parameter dimensions of refining NNs <cit.>.
The natural question then arises, How one can choose between sCNN sys:auto:res and dCNN sys:dcnn refiners, 
and select their parameters (R, K, and L)?
The first answer to this question depends on some understanding of data-fit f(x;y) in MBIR problem sys:mbir,
e.g., the regularization strength γ and the condition number variations across training data-fit majorizers.
(An additional criteria could be general understandings between sample size/diversity and parameter dimension of NNs.)


For example, the sparse-view CT system in <ref> needs moderate regularization strength (χ^⋆= 167.64) 
and the majorization matrices of its training data-fits have mild condition number variations (the standard deviation is 1.1).
training data-fits have mild parameter variations across samples.
Comparing results between Momentum-Net-sCNN and -dCNN in Fig. <ref> and Table <ref>
demonstrates that sCNN sys:auto:res seems suffice.
Table <ref>(d)–(e) shows that one can further improve the refining accuracy of sCNN sys:auto:res by increasing its width, i.e., K.
The LF photography system using a limited focal stack in <ref> needs a large γ value (χ^⋆= 1.5),
and the majorization matrices of its training data-fits have large condition number variations (the standard deviation is 2245.5).
Comparing results between Momentum-Net-sCNN and -dCNN in Fig. <ref> and Table <ref>
demonstrates that dCNN sys:dcnn yields higher PSNR than sCNN sys:auto:res.
For dCNN sys:dcnn, we observed increasing its depth, i.e., L, up to a certain number 
is more effective than increasing its width, i.e., K, as briefly discussed in <ref>.

For choosing the relaxation parameter ρ in eq:momnet:map, 
we also suggest considering the regularization strength in eq:momnet:mbir.
For an application that needs moderate regularization strength,
e.g., sparse-view CT in <ref>,
we suggest setting ρ to 0.5 so as to mix information between input and output of refining NNs,
rather than 1 - ε that does not mix input and output.
For an application that needs strong regularization,
e.g., LF photography using a limited focal stack in <ref>,
we suggest using ρ= 1 - ε than ρ= 0.5.
Results in the next section validate this suggestion.




 §.§ Performance of Momentum-Net with different relaxation parameters ρ in eq:momnet:map
 


Fig. <ref> below compares the performances of Momentum-Net-sCNN with different ρ values.
The results in Fig. <ref> support the ρ selection guideline in <ref>.
One can maximize the MBIR accuracy of Momentum-Net by properly selecting ρ.


Note that ρ∈ (0,1) controls strength of inference from refining NNs in eq:momnet:map, but does not affect the convergence guarantee of Momentum-Net. 
Fig. <ref> illustrates that Momentum-Net appears to converge regardless of ρ values.










§ PARAMETERS OF MOMENTUM-NET
 


Table <ref> below lists parameters of Momentum-Net, and summarizes selection guidelines or default values.
Similar to BCD-Net/ADMM-Net, 
the main tuning jobs to maximize the performance of Momentum-Net
include selecting architectures of refining NNs {_θ^(i):∀ i } in eq:momnet:map,
and choosing a regularization parameter γ in eq:momnet:mbir by tuning χ in <ref>.
One can simplify the tuning process by using the selection guidelines in <ref> for selecting architectures of {_θ^(i):∀ i }, and
training χ in <ref>.
Note that one designs majorization matrices { M^(i):∀ i } rather than tuning them:
majorization matrices can be analytically designed, e.g., Lemma <ref> as used in <ref>;
one can algorithmically design them <cit.>.
Tighter majorization matrices are expected to further accelerate the convergence of Momentum-Net <cit.>.


















IEEEtran
referencesSupp_Bobby






