Upper-Bounding the Regularization Constant for Convex Sparse Signal Reconstruction
Renliang Gu Aleksandar Dogandžić

Image Stitching by Line-guided Local Warping
with Global Similarity Constraint
    Tianzhu Xiang^1, Gui-Song Xia^1, Xiang Bai^2, Liangpei Zhang^1



^1State Key Lab. LIESMARS, Wuhan University, Wuhan, China.

^2Electronic Information School, Huazhong University of Science and Technology, China.

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






  Consider reconstructing a signal  by minimizing a weighted sum of a 
  convex differentiable NLL (data-fidelity) term and a convex 
  regularization term that imposes a convex-set constraint on  and 
  enforces its sparsity using ℓ_1-norm analysis regularization.  
  We compute upper bounds on the regularization tuning constant beyond 
  which the regularization term overwhelmingly dominates the NLL term 
  so that the set of minimum points of the objective function does not 
  change.  Necessary and sufficient conditions for irrelevance of sparse 
  signal regularization and a condition for the existence of finite upper 
  bounds are established.  We formulate an optimization problem for finding 
  these bounds when the regularization term can be globally minimized by a 
  feasible  and  also develop an ADMM type method for their 
  computation.  Simulation examples show that the derived and empirical 
  bounds match.


 
 
 

CPU



§ INTRODUCTION



Selection of the regularization tuning constant u > 0 in convex 
Tikhonov-type <cit.> penalized NLL minimization 
rCl
  
  f_u()=() +u r()

is a challenging problem critical for obtaining accurate estimates of the 
signal  <cit.>.  Too little regularization leads to 
unstable reconstructions with large noise and artifacts due to, for 
example, aliasing.  With too much regularization, the reconstructions  are 
too smooth and often degenerate to constant signals.  Finding bounds on the 
regularization constant u or finding conditions for the irrelevance of 
signal regularization has received little attention.  In this paper, we 
determine upper bounds on u beyond which the regularization term r() 
overwhelmingly dominates the NLL term () in 
(<ref>) so that the minima of the objective function 
f_u() do not change. For a linear measurement model with white 
Gaussian noise and ℓ_1-norm regularization, a closed-form expression 
for such a bound is determined in <cit.>; see also 
Example <ref>.    The obtained bounds can be used to design 
continuation procedures <cit.> that gradually 
decrease u from a large starting point down to the desired value, which 
improves the numerical stability and convergence speed of the resulting 
minimization algorithm by taking advantage of the fact that penalized 
NLL schemes converge faster for smoother problems with larger u 
<cit.>.  In some scenarios, users can monitor the 
reconstructions as u decreases and terminate when the result is 
satisfactory.

Consider a convex NLL () and a regularization term
c
  
  r() =  𝕀_C() +  Ψ^H _1   

that imposes a convex-set constraint on , ∈ C ⊆R^p, and sparsity of an appropriate linearly transformed 
, where Ψ∈C^p× p'  is a known 
sparsifying dictionary matrix.
Assume that the NLL () is differentiable and lower bounded
within the closed convex set C, and satisfies
c
  
  () ⊇C

which ensures that  () is computable for all ∈ C.
Define the convex sets of solutions to min_ f_u(),
min_ r(), and min_∈ Q():[The use of 
  “≤” in the definitions of Q and ^♢
 in (<ref>) and (<ref>) makes it easier to identify both 
 as convex sets.]

rCl
  
  _u
        
  
     | 
    f_u()=min_f_u()
  

  
  Q         |r() = min_r()  
  

     =   
  
    ∈C  |   Ψ^H _1 ≤min_∈C Ψ^H _1
  
  

  
  ^♢      
  
    ∈Q |() ≤min_∈Q ()
  
  ≠∅

where the existence of ^♢ is ensured by the assumption that 
() is lower bounded in C.


We review the notation: “^*”, “^T”, “^H”, “^+”, 
·_p, ·, ⊗, “≽”, “≼”, 
I_N, 1_N × 1, and 0_N×1 denote complex 
conjugation, transpose, Hermitian transpose, Moore-Penrose matrix inverse, 
ℓ_p-norm over the complex vector space C^N defined by 
_p^p = ∑_i=1^N |z_i|^p for = (z_i) ∈C^N, 
absolute value, Kronecker product, elementwise versions of “≥” and 
“≤”, the identity matrix of size N and the N × 1 vectors of 
ones and zeros, respectively (replaced by I,1, and 0 when the 
dimensions can be inferred).  𝕀_C() =  0,    ∈ C 

 +∞,    otherwise,
C = min_∈ C-_2^2,
and exp_∘
denote the indicator function, projection onto C, and the elementwise 
exponential function: exp_∘_i= exp a_i.

  Denote by 𝒩(A) and ℛ(A) the null space and range 
  (column space)  of a matrix A. These vector spaces are real or complex 
  depending on whether A is a real- or complex-valued matrix.
For a set S of complex vectors of size p, define S 
∈R^p
  |+ ∈ S  for some ∈R^p and S∩R^p∈R^p |+0∈ S, where = 
√(-1).
For A ∈C^M × N, 
c'c
  
  𝒩(A^H) ∩R^M =
  𝒩 A^T   ,    
   ℛ(A) =  ℛ( A )
 
are the real null space and range of A^T and 
A, respectively, where
c
  
  A A    A  ∈R^M 
  ×2 N.

If A in (<ref>)
 has full row rank, we can define
c
  A^A^H (A A^H)^-1
  

which reduces to A^+ for real-valued A.
The following are equivalent:
(ℛΨ)=R^p,
𝒩(Ψ^H) ∩R^p={0},
and
d=p, where
 c
    
    d ℛ(Ψ) ≤min(p,2p').
  
We can decompose Ψ as
  c
    Ψ= F Z
    
  
where F ∈R^p × d and  Z ∈C^d 
× p' with F = d and Z = d; 
Z =  Z   Z∈R^d 
× 2 p', consistent with the notation in (<ref>).  
Here, ℛF denotes the real range of the real-valued 
matrix F.
Clearly, d≥ 1 is of interest; otherwise Ψ=0.
Observe that  (see (<ref>))

 rCl
  
  (ΨZ^)   =   F

  
  ℛF   =   (ℛΨ).







The subdifferential of the indicator function 
N_C()=∂𝕀_C() is the normal cone to C at 



<cit.>
and, by the definition of a cone, satisfies
c"s
  
  N_C()=aN_C(),     for any a>0.

Define
c
  
  G(s)
    s/|s|,     s≠0

    w∈C ||w|≤1,     s=0
  

and its elementwise extension G() for vector arguments , which 
can be interpreted as twice the Wirtinger subdifferential of _1 
with respect to  <cit.>.  Note that 
^HG()=_1, and, when  is a real vector, 
G() is the subdifferential of _1 with 
respect to  <cit.>.





























  
  For Ψ∈C^p× p' and ∈R^p, the 
  subdifferential of Ψ^H_1 with respect to  is  
  c
    
    ∂_Ψ^H _1 = ΨGΨ^H.
  


  (<ref>) follows from
  rCl
    
    ∂_|_j^H|
    =
    _j G(_j^H)
  
  where _j is the jth column of Ψ.  We obtain 
  (<ref>) by replacing the linear transform matrix in 
  <cit.> with _j  _j^T.



We now use Lemma <ref>
to formulate the necessary and sufficient conditions for 
∈_u:

rCl
  
  0    ∈   
  u ΨGΨ^H 
  + ∇( ) + N_C( )

and ∈ Q:
rCl
  
  0   ∈    ΨGΨ^H
  +N_C
       
    
    
    
    
    
 


respectively.








When the signal vector =  X corresponds to an image X ∈R^J× K, its isotropic and anisotropic TV 
regularizations correspond to
<cit.>

  
  rCl"s"
    
    Ψ   =   Ψ_v+Ψ_h∈C^JK×JK     
    (isotropic)
    

    Ψ   =   [Ψ_v  Ψ_h]∈R^JK×2JK
        (anisotropic)
    
  

respectively, where
  Ψ_v=I_K⊗ D^T(J) and
  Ψ_h=D^T(K)⊗ I_J
are the vertical and horizontal difference matrices (similar to those in  
<cit.>),  and
c
  
  D(L)[r]
    1      -1                       

            1       -1              

                ⋱   ⋱     

                             1    -1

    0       0   ⋯    0     0
    ∈R^L×L

obtained by appending an all-zero row from below to
the (L-1) × L upper-trapezoidal matrix with first row 
1,-1,0,…,0;  note that D(1)=0.
Here, d = JK-1 and
c
  
  𝒩(Ψ^H)=ℛ(1)

for both the isotropic and anisotropic TV regularizations. 



The scenario where
rCl
  
  𝒩(Ψ^H) ∩C ≠∅
holds is of practical interest: then Q = 𝒩(Ψ^H) ∩ C and 
^♢∈^♢ globally minimize the regularization term: 
r^♢=0.  If (<ref>) holds and ^♢∈^♢, then G(Ψ^H^♢)=H, where
c
  
  H ∈C^p' ×1 |_∞≤1.

If, in addition to (<ref>),

  
  * d=p,
    then  ^♢=Q={0};
  
  * 𝒩(Ψ^H) ∩R^p=ℛ(1), 
    then   Q= ℛ(1) ∩ C and
    ^♢∈^♢ are constant signals of the 
    form ^♢ = 1 x_0^♢,   x_0^♢∈R.
   
   
   




In Section <ref>, we define and explain an upper bound U on 
useful regularization constants u and establish conditions under which signal 
sparsity regularization is irrelevant and finite U does not 
exist.
We then present an optimization problem for finding U when 
(<ref>) holds (Section <ref>), develop a general 
numerical method for computing bounds U (Section <ref>), present 
numerical examples (Section <ref>), and make concluding remarks 
(Section <ref>).





§ UPPER BOUND DEFINITION AND PROPERTIES



Define
c
  
  U infu ≥0  _u ∩Q ≠∅ .

  If _u ∩ Q = ∅ for all u, then finite U does not 
  exist, which we denote by U=+∞.


We now show that, if u ≥ U, then the
the set of minimum points _u of the objective function does not 
change.


  

  
  * For any u, _u∩ Q=^♢ if and only if _u∩ 
    Q≠∅. 
  
  * 
    Assuming _U∩ Q ≠∅ for some U ≥ 0, _u = 
    ^♢ for u > U.





  We first prove <ref>.
  Necessity follows by the existence of ^♢; see  
  (<ref>).
  We argue sufficiency by contradiction.  Consider any 
  _u∈_u∩ Q; i.e., _u minimizes both f_u() and 
  r().  If _u ∉^♢, there exists a ∈^♢ with ()<(_u) that, by the definition of 
  ^♢, also minimizes r().  Therefore, 
  f_u()=()+ur()<f_u(_u), which contradicts the assumption 
  _u∈_u.  Therefore, _u∩ Q⊆^♢.  If 
  there exists a ∈^♢⊆ Q such that 
  ∉_u, then f_u()>f_u(_u) which, since both  and 
  _u are in Q, implies that ()>(_u) and contradicts 
  the definition of ^♢.  Therefore, 
  ^♢⊆_u.

  We now prove <ref>. By <ref>, _U∩ 
  Q=^♢, which confirms <ref> for u=U.
  Consider
  now u>U, a ∈_U∩ Q = ^♢, and any ∈_u. Then,
  
    
    rCl
      ()+Ur()   ≥   ()+Ur()

      ()+ur()   ≥   ()+ur().
    
  
  By summing the two inequalities in (<ref>) and 
  rearranging, we obtain r() ≥ r().  Since ∈ Q,  is 
  also in Q; i.e., _u⊆ Q,
  which implies  _u=^♢ by <ref>.






As u increases, _u moves gradually towards Q and, according to the 
definition (<ref>), _u and Q do not intersect when u<U.  
Once u=U, the intersection of the two sets is ^♢, and, by
Remark <ref><ref>,
_u = ^♢ for all u > U.





































































 §.§ Irrelevant Signal Sparsity Regularization




  
  The following claims are equivalent:
  
    
  * 
      ^♢∩_0≠∅; i.e., there exists an 
      ^♢∈^♢ such that
      rCl
        
        0∈∇^♢+N_C^♢;
      
    
  * 
      
      ^♢⊆_0;
      and
    
  * 
      U=0; i.e., _0∩ Q≠∅.
  





  <ref> follows from <ref> because 
  ^♢⊆ Q.  <ref> follows from <ref> 
  by applying Remark <ref><ref> to obtain _0∩ 
  Q=^♢, which implies <ref>.  
  Finally, <ref> implies <ref>.









Having ∇(^♢)=0 for at least one 
^♢∈^♢ implies (<ref>) and is therefore 
a stronger condition than (<ref>).


  
  Consider ()=_2^2 and C=∈R^2 |-1_2× 1_2≤1.
  (Here, () could correspond to the Gaussian measurement model with 
  measurements equal to zero.)
  Since C is a circle within R_+^2, the objective functions 
  for the identity (Ψ=I_2) and 1D TV 
  
  
  sparsifying transforms are
  
    rCl's
      f_u()   =   x_1^2+x_2^2 + u(x_1+x_2)+𝕀_C(),     (identity)  
      

      f_u()   =   x_1^2+x_2^2 + u|x_1-x_2|+𝕀_C(),     (1D TV)
  
  
  respectively, where 𝒳_u = 𝒳^♢=Q={^♢}
  and
  ^♢ = 1-√(2)/21.
  Here, 
  ∇^♢=2-√(2)1_2× 1 
  and N_C^♢=a1| a≤0, which 
  confirms that (<ref>) holds.







 §.§ Condition for Infinite U and Guarantees for Finite U




  
  If there exists ^♢∈^♢ such that
  rCl
    
    ∇^♢+N_C^♢
    ∩(ℛ Ψ)=∅.
  
  then U=+∞.  When (<ref>)
  holds, the reverse is also true with a stronger claim: U=+∞ 
  implies (<ref>) for all ^♢∈^♢.


  First, we prove sufficiency by contradiction.  If a finite U exists, 
  then ^♢⊆_u for all u≥ U.  Therefore, 
  (<ref>) holds with  being any 
  ^♢∈^♢, which contradicts (<ref>).

  In the case where (<ref>) holds, we prove the necessity by 
  contradiction.
  If (<ref>) does not hold for all 
  ^♢∈^♢, there exist ∈ N_C(^♢) and 
  ∈C^p' such that
  c
    
    0 = ∇(^♢)+Ψ+.
  
  Since (<ref>) holds, Ψ^H ^♢ = 0 and 
  G(Ψ^H^♢)=H; see (<ref>).  When 
  u≥_∞, ∈ uH and Ψ∈
  uΨ GΨ^H^♢.   Therefore, 
  (<ref>) holds at =^♢ for all 
  u≥_∞, which contradicts U=+∞.





















  
  Consider ()=x_1+𝕀_R_+(x_1), Ψ=I_2, 
  and
  C=∈R^2 |-1_2× 1_2≤1.
  (Here, () could correspond to the
  (x_1) measurement model with measurement equal to zero.)
  Since C is a circle within R_+^2, the objective function 
  is
  rCl
    f_u()=(1+u)x_1+ux_2+𝕀_C()
  
  with 𝒳_u = {_u },  𝒳^♢ = Q = {^♢},
and

  rCl
    
    _u    =    1_2 ×1-1/√(2+2/u+1/u^2)
    [ 1+1/u;     1 ]

    
    ^♢   =    1-√(2)/21_2×1
  

which implies U=+∞, consistent with
the observation that _u ∩ Q = ∅.  Here, (<ref>) is 
not satisfied:  (<ref>) is only a sufficient condition for 
U=+∞ and does not hold in this example.   


 





  
  Consider ()=^2_2, 1D TV sparsifying transform 
  with Ψ=D^T(2), and
  C=∈R^2 |-2,   0^T_2^2≤2.
  Since C is a circle with x_1-x_2≥0, the objective function is
  
    rCl
      f_u()   =   _2^2+u|x_1-x_2|+𝕀_C()
      

         =   
      -12[u  -u]^T_2^2
      -u^2/2+𝕀_C()
    
  
  with 𝒳_u = 
      2-(1+4/u)/q(u),  
      1/q(u)
  ^T, q(u) √(1+4/u+8/u^2),
  and 𝒳^♢ = Q = {1_2× 1}, which implies 
  U=+∞.  Since (<ref>) holds in this example, 
  (<ref>) is necessary and sufficient for U=+∞.  Since  
  -1^T∇^♢=-4 and 
  N_C^♢=(-a,a)^T | a≥0,  
  (<ref>) holds.















































































































  §.§.§ Two cases of finite U



If  d=p and (<ref>) holds, then U must be finite:  in this 
case,  condition (<ref>) in Remark <ref> cannot 
hold, which is easy to confirm by substituting 
(ℛΨ)=R^p into 
(<ref>).






























U must also be finite if
c
    
  
  ^♢∩C≠∅.

Indeed, (<ref>)
implies (<ref>) and that for ^♢∈^♢∩ C,

  rCl
    
    N_C^♢    =    {0}

     
  ∇(^♢)    ∈    ℛ(Ψ) 


and hence (<ref>) cannot hold upon substituting 
(<ref>) and (<ref>).  Here, 
(<ref>) follows from 
0∈∇(^♢)+N_Q(^♢), the condition for 
optimality of the optimization problem min_∈ Q() that 
defines ^♢, by using the fact that N_Q(^♢) = 
ℛ(Ψ) when ^♢∈^♢∩ C.


























If (<ref>) holds then,  by Remark <ref>, U=0 if 
and only if ∇^♢=0. 




§ BOUNDS WHEN (<REF>) HOLDS




















































































































































































We now present an optimization problem for finding U when (<ref>) 
holds.



  
  Assume that (<ref>) holds and that the convex NLL 
  () is differentiable within ^♢.  Consider the 
  following optimization problem:
  
    
    s'c+l
      
      (P_0):     U_0(^♢) =     -15pt
      min_∈R^p,   ∈C^p'
      (^♢,,)_∞

             
      subject to 
      10pt  ∈N_C^♢

      
             
      48pt ∇^♢+∈ℛF
      
  
  with
  c
    
    (,,)
    + 
    Z^F^+∇+-(Z).
  
  Then,
  U_0(^♢)=U for all ^♢∈^♢ and U  in 
  (<ref>).

Here, 




U=+∞ if and only if the constraints in  (<ref>) and 
(<ref>) cannot be satisfied for any , which is 
equivalent to ^♢∈^♢ satisfying 
(<ref>) in Remark <ref>.



  
  Observe that G(Ψ^H^♢)=H for all  
  ^♢∈^♢ and
  c
    
    Ψ(,,)
    =∇+.
  
  due to (<ref>) and (<ref>), respectively.

  We first prove that ^♢⊆_u if u≥ 
  U_0(^♢).






  Consider any ^♢∈^♢ and denote by 
  , a pair , that 
  solves the minimization problem (P_0).  Since u≥ 
  U_0(^♢), there exists an ∈ H such that 
  ^♢,,+u=0.  
  Using (<ref>), we obtain
  c
    0=
      Ψ^♢,,
        +u
      
     
    =uΨ+∇(^♢)+
    

  
  which implies ^♢∈_u according to 
  (<ref>).

  Second, we prove that if u< U_0(^♢) for any 
  ^♢∈^♢, then ^♢∩_u=∅. We 
  employ proof by contradiction.  Suppose 
  ^♢∩_u≠∅; then, there exists an 
  ^♢∈^♢∩_u. According to (<ref>), 
  there exist an ∈ H and an ∈ 
  N_C(^♢) such that 
  0=uΨ+∇(^♢)+.  
  Using (<ref>), we have
  rCl
    
    0
       =   
    
      Ψu+^♢,,-u
    .
  
  Note that
  rCl
    u+^♢,,-u
       =   
    Z^F^+∇^♢++u(Z).  
    

    
  
  Inserting (<ref>) into (<ref>) and using 
  (<ref>) and the fact that F has full column rank leads to 
  0=F^+∇^♢++u(Z); 
  thus
  rCl
    0
       =   
    u+^♢,,-u.
  
  Now, 
  
  
  
  
  rearrange and use the fact that _∞≤1  (see 
  (<ref>)) to obtain
  c
    
    ^♢,,-u_∞= u-_∞≤u
    < U_0(^♢)
  
  which contradicts (<ref>), where U_0(^♢) is the 
  minimum.

  Finally, we prove by contradiction that U_0(^♢) is invariant 
  within ^♢ if ^♢ has more than one element.  
  Assume that there exist ^♢_1, ^♢_2 ∈^♢ 
  and u such that U_0(^♢_1)≤ u<U_0(^♢_2).  We 
  obtain contradictory results: ^♢_1∈_u and 
  ^♢∩_u≠∅ because u≥ U_0(^♢_1) 
  and u<U_0(^♢_2), respectively.  Therefore, 
  U=U_0(^♢) is invarant to ^♢∈^♢.
  

 The constraints on  in (<ref>) and 
 (<ref>) are equivalent to stating that 
 (<ref>) does not hold for any ^♢∈𝒳^♢; see also (<ref>).  If an  does not 
 exist that satisfies these constraints, (<ref>) holds and 
 U=+∞ according to Remark <ref>.


We make a few observations: (P_0) is a linear programming problem with 
linear constraints and can be solved using CVX <cit.> and Matlab's 
optimization toolbox upon identifying N_C(^♢) and 
ℛ(F) in (<ref>) and (<ref>), 
respectively.  Theorem <ref> requires differentiability of the 
NLL only at =^♢∈^♢.  If Ψ is real, 
then Z is real as well, the optimal  in (P_0) has zero imaginary 
component and the corresponding simplified version of Theorem <ref> 
follows and requires optimization in (P_0) with respect to real-valued 
∈R^p'.





If Ψ is real and d = p', then we can select Z=I, which leads to 
Z^=I and cancellation of the variable  in 
(<ref>) and simplification of (P_0).

We now specialize Theorem <ref> to two cases with finite U.


  
  If d=p and if (<ref>) holds, then U in (<ref>) can be 
  computed as
  c
    
    U = min_∈N_C0,   ∈C^p'
    
      +
      Ψ^∇0+-(Ψ)  
    _∞.  


  Theorem <ref> applies, 𝒳^♢={0}, and U 
  must be finite.  Setting F=I in (<ref>) leads to 
  (<ref>).

If C=R_+^p, then N_C(0)=R_-^p and the 
condition ∈  N_C0 reduces to ≼0.























  
  If (<ref>) holds, then U in (<ref>) can be computed 
  as
  rCl
    
    U   =   min_∈C^d
     + Z^F^+ ∇^♢ - (Z )   
    _∞

with any ^♢∈^♢∩ C.

 

  Thanks to (<ref>), (<ref>) and 
  (<ref>)–(<ref>) are satisfied, 
  Theorem <ref> applies,
  U must be finite, and =0 (by (<ref>)). By using 
  these facts,  we simplify (<ref>) to obtain (<ref>).


















If d=p and 0∈ C, then both Corollaries <ref>
and <ref> apply and the upper bound U can be obtained by 
setting =0  and N_C(0)={0} in
(<ref>) or
by setting  ^♢ = 0 and  F=I in  (<ref>).


  
  Consider a real invertible Ψ∈R^p× p.
  
    
  * If C=R_+^p,
      Corollary <ref> applies and (<ref>) becomes
  
  c
    U=
    min_≼0
    Ψ^-1 ∇0 + _∞.
    
  
  In this case, U=0 and signal sparsity regularization is irrelevant if 
  ∇(0) ≽0, which follows by inspection from  
  (<ref>), as well as from (<ref>) in 
  Remark <ref>. If Ψ=I, (<ref>)
  further reduces to U=-min0,min_i[∇(0)]_i.
  
  * 
  If 0∈ C, Corollaries <ref>
  and <ref> apply and the bound U simplifies to
  c
    
    U = Ψ^-1 ∇0  _∞.
  














    For Ψ=I and a linear measurement model with white Gaussian noise, 
    (<ref>) reduces to the expressions
    in <cit.> and <cit.>, 
    used in <cit.> to design its continuation scheme;  
    <cit.> and <cit.> also assume 
    C=R^p.
 




























[One-dimensional TV regularization]
  
  Consider 1D TV regularization with Ψ=D^T(p) ∈R^p × p obtained by setting K=1,J=p in 
  (<ref>); note that d=p-1.   Consider a constant signal 
  ^♢ = 1 x_0^♢∈^♢. Then 
  Theorem <ref>
  applies and yields
  
  rCl
    
    U=  min_∈N_C( 1 x_0^♢)  max_1 ≤j < p
    ∑_i=1^j    ∇1 x_0^♢+_i
  
  where we have used the factorization (<ref>)
  with F obtained by the block partitioning Ψ = F  0_p×1, Z=I_p-1 0_(p-1)×1, and 
  the fact that F^+ is equal to the (p-1)× p lower-triangular 
  matrix of ones.  When  (<ref>) holds, 1x_0^♢∈^♢∩ C,
  Corollary <ref> applies,
  = 0 (see (<ref>)), and (<ref>) reduces to:
  rCl
    
    U=max_1 ≤j < p
    ∑_i=1^j ∇ 1 x_0^♢_i.
  



  
  
  





  

































































































The bounds obtained by solving (P_0) are often simple but restricted 
to the scenario where (<ref>) holds.  In the following section, we 
remove assumption  (<ref>) and develop a general numerical method 
for finding U in (<ref>).



§ ADMM ALGORITHM FOR COMPUTING U



We focus on the nontrivial scenario where (<ref>)
does not hold and assume u>0.
We also assume that an ^♢∈𝒳^♢ is 
available, which will be sufficient to obtain the U in (<ref>).
We use the duality of norms <cit.>:
c
  
  Ψ^H_1 = max__∞≤1 
  ^HΨ^H

to rewrite the minimization of (<ref>) as the following 
min-max problem (see also (<ref>)):
rCl
  
  min_ max_ ()+u ^HΨ^H
  +𝕀_C()-𝕀_H().
  
Since the objective function in (<ref>) is convex with respect to 
 and concave with respect to , the optimal (,)= 
_u,_u is at the saddle point of
(<ref>) and satisfies

  
  rCl
    
    0   ∈   ∇_u
    +uΨ_u
    +N_C_u
    

    
    
    _u   ∈   G(Ψ^H_u).
  




















Now, select U as the smallest u for which 
(<ref>)–(<ref>) hold with _u=^♢:
c
  
  U=  1/v^♢
 ∇^♢_2
 
where v^♢, ^♢, ^♢ is the solution to 
the following constrained linear programming problem:
s'l+l+x*
  
  
   (P_1):     v,,minimize
     
  -v
  +𝕀_G(Ψ^H^♢)()
  +𝕀_N_C^♢() 

      subject to
      20pt
  
  v+Ψ+ =0
 
 
 
 
 
 

obtained from (<ref>)–(<ref>)
with _u and _u replaced by ^♢ and .  Here, 
c
  ∇^♢ / 
  ∇^♢_2
  

is the normalized gradient (for numerical stability) of the NLL at 
^♢;
∇^♢≠0 because (<ref>) 
does not hold.  Due to (<ref>), v=0 is a feasible point 
that satisfies the constraints (<ref>), which implies that 
v^♢≥0.  When (<ref>) holds, v has to be zero, 
implying U=+∞.


To solve (P_1) and find v^♢, we apply an iterative algorithm 
based on ADMM <cit.>


rCl
  
  ^(i+1)    =   -10pt min_
    ∈GΨ^H^♢
    
    v^(i)+
    Ψ +^(i) +^(i)
  _2^2
  

  v^(i+1)   =   ρ-^TΨ^(i+1)+^(i)+^(i)
  

  
  ^(i+1)   =   P_N_C^♢
    -v^(i+1)-Ψ^(i+1)-^(i)
  
  

  ^(i+1)   =   ^(i)+
    Ψ^(i+1)
    +v^(i+1)+^(i+1)
  


where ρ>0 is a tuning parameter for the ADMM iteration and we 
solve (<ref>) using the
 Broyden-Fletcher-Goldfarb-Shanno
optimization algorithm  with box constraints
<cit.> and PNPG algorithm  <cit.> for 
real and complex Ψ, respectively.  We initialize the iteration 
(<ref>) with v^(0)=1, ^(0)=0,  ^(0)=0, 
and ρ=1, where ρ is adaptively adjusted
thereafter using the scheme in <cit.>.









In special cases, (<ref>) simplifies.   If (<ref>) holds, 
then Ψ^H^♢ = 0 and the constraint in 
(<ref>) simplifies to _∞≤1; see 
(<ref>).  If (ΨΨ^H)=cI,   c>0, and 
Ψ∈R^p× p or Ψ∈C^p× p/2,























(<ref>) has the following analytical solution:
c
  
  ^(i+1) =GΨ^H^♢
  -1/cΨ^Hv^(i)+^(i)+^(i).

When (<ref>) holds, (<ref>) reduces to ^(i)=0 
for all i, thanks to (<ref>).

When Ψ is real, the constraints imposed by 
𝕀_G(Ψ^H^♢)() become linear and (P_1)  
becomes a linear programming problem with linear constraints.




retain-zero-exponent=true




§ NUMERICAL EXAMPLES



Matlab implementations of the presented examples are available at 
<https://github.com/isucsp/imgRecSrc/uBoundEx>.


In all numerical examples,
the empirical upper bounds U were obtained by a grid search over u with 
𝒳_u ={_u} obtained using the PNPG method 
<cit.>.




 §.§ Signal reconstruction for Gaussian linear model





We adopt the linear measurement model with white Gaussian noise and
scaled NLL  () = 0.5 -Φ_2^2,
where the elements of the sensing matrix Φ∈R^N × 
p are iid and drawn from the uniform distribution on a unit sphere.  
We reconstruct the nonnegative “skyline” signal _true∈R^1024 × 1  in 
<cit.> from noisy linear 
measurements  using the DWT and
1D TV regularizations, where
the DWT matrix Ψ is orthogonal  (ΨΨ^T = Ψ^T Ψ = 
I),
constructed using the Daubechies-4 wavelet with three decomposition levels.
Define the SNR as

    () = 
      10log_10Φ_true_2^2/Nσ^2

where σ^2 is the variance of the Gaussian noise added to Φ_true to create the noisy measurement vector .


For C=R_+^p and C=R^p with DWT 
regularization, 𝒳^♢={0} and 
Example <ref> applies and yields the upper bounds 
(<ref>) 





and (<ref>), respectively.

For TV regularization, we apply the result in Example <ref>.  
For C=R^p and C=R_+^p, we have   
𝒳^♢={1 x_0} and 𝒳^♢={1max(x_0,0) }, respectively, where
rCl
  
  x_0 min_x ∈R  (1 x)  = 
  1^TΦ^T/ Φ1^2_2.

If 1 x_0 ∈ C, which holds when C=R^p or 
when C=R_+^p and x_0 > 0, then the bound U is given by  
(<ref>).
 For  C=R_+^p and if x_0 ≤ 0, then
𝒳^♢={0} and (<ref>) applies.   In this 
case, U=0 if [∇(0)]_i ≥ 0 for i=1,…,p-1, which 
occurs only when [∇(0)]_i=0 for all i.









Table <ref> shows the theoretical and empirical bounds for 
DWT and TV regularizations and C=R_+^p and
C=R^p; we decrease the SNR from 
30-30 with independent noise realizations for 
different SNR.  The theoretical bounds in 
Sections <ref>
and <ref> coincide.








For DWT regularization, ^♢ is the same for both convex 
sets C and thus the upper bound U for C=R_+^p is always 
smaller than its counterpart for C=R^p, thanks to being 
optimized over variable  in (<ref>).
For TV regularization, when x_0>0, the upper bounds U coincide 
for both C because, in this case,
^♢ is the same for both C and ^♢∈ C.  
In the last row of Table <ref> we show the case where 
x_0≤0; then, ^♢ differs for the two convex sets C,
and the upper bound U for C=R_+^p is smaller than its 
counterpart for C=R^p, thanks to
being optimized over variable
  in (<ref>): compare (<ref>) with (<ref>).







































 §.§ PET image reconstruction from Poisson measurements




Consider PET reconstruction of the 128×128 concentration map 
_true in <cit.>, 
which represents simulated radiotracer activity in a human chest, from 
independent noisy Poisson-distributed measurements =(y_n) with means 
[Φ_true+]_n.  The choices of parameters in the PET 
system setup and concentration map _true have been taken from 
the IRT <cit.>.  Here,

c
  ()=1^T Φ+-
  + ∑_n,y_n≠0y_nlny_n/Φ+_n
  
  
and
c
  Φ= w (exp_∘(-S+))S
  ∈R_+^N ×p
  


is the known sensing matrix;  is the density map
needed to model the attenuation of the gamma rays 
<cit.>; =(b_i) is the known intercept term 
accounting for background radiation, scattering effect, and accidental 
coincidence;[The elements of the intercept term have been set to a 
  constant equal to 10 of the sample mean of Φ_true: =[1^TΦ_true/(10N)]1.]
 is a known vector that models the detector efficiency variation; and 
w>0 is a known scaling constant, which we use to control the expected 
total number of detected  photons due to electron-positron annihilation, 
1^T(-) = 1^TΦ_true, an SNR 
measure.  We collect the photons from 90 equally spaced directions 
over  180, with 128 radial samples at each direction.  
Here,   we adopt the parallel strip-integral matrix S 
<cit.> and use its implementation in the 
IRT <cit.>.


We now consider the nonnegative convex set C=R_+^p,
which ensures that (<ref>) holds, and 2D isotropic and anisotropic 
TV and DWT regularizations, where the 2D DWT matrix 
Ψ is constructed using the Daubechies-6 wavelet with six decomposition 
levels.


For TV regularizations, 𝒳^♢={1max(0,x_0) 
}, where x_0 = min_x ∈R(1 x), computed 
using the bisection method that finds the zero of
∂(1x)/∂ x,
which is an increasing function of x∈R_+.  Here, no search 
for x_0 is needed when [0]∂(1x) / ∂ x_x=0  
> 0, because in this case x_0<0. 
















We computed the theoretical bounds using the ADMM-type algorithm in 
Section <ref>.

Table <ref> shows the theoretical and empirical bounds for 
DWT and TV regularizations and the SNR 
1^TΦ_true
varying from e1 to e9,  with independent measurement 
realizations for different SNR.


















Denote the isotropic and anisotropic 2D TV bounds by U_iso 
and U_ani, respectively. Then, it is easy to show that when 
(<ref>) holds,
U_ani≤ U_iso≤√(2)U_ani, which follows by 
using the inequalities 
√(2)√(a^2+b^2)≥a+b≥√(a^2+b^2) and is 
confirmed in Table <ref>.


































§ CONCLUDING REMARKS



Future work will include obtaining simple expressions for upper bounds U 
for isotropic 2D TV regularization, based on Theorem <ref>.



