We address the problem of optimal scale-dependent parameter learning in total variation image denoising. Such problems are formulated as bilevel optimization instances with total variation denoising problems as lower-level constraints. For the bilevel problem, we are able to derive M-stationarity conditions, after characterizing the corresponding Mordukhovich generalized normal cone and verifying suitable constraint qualification conditions. We also derive B-stationarity conditions, after investigating the Lipschitz continuity and directional differentiability of the lower-level solution operator. A characterization of the Bouligand subdifferential of the solution mapping, by means of a properly defined linear system, is provided as well. Based on this characterization, we propose a two-phase non-smooth trust-region algorithm for the numerical solution of the bilevel problem and test it computationally for two particular experimental settings.

Keywords: Bilevel optimization, machine learning, variational models, total variation.

AMS: 49K99, 90C33, 68U10, 68T99, 65K10

§ INTRODUCTION

In this work we address the bilevel optimal selection of parameters in total variation image denoising models, from a variational and nonsmooth analysis perspectives. Bilevel techniques were proposed for optimal parameter selection of variational models in the seminal work <cit.>, provided a training dataset to learn from. The variational models considered there were based on the total variation (TV) seminorm and different noise statistics models were taken into account. Thereafter, apart of variational denoising problems, the bilevel learning framework has been considered for other imaging applications such as image restoration <cit.>, blind image deconvolution <cit.>, image segmentation <cit.>, nonlocal models <cit.> and kernel parameter estimation for support vector machines <cit.>.
The bilevel methodology was also extended by Cao et al. <cit.>, Hintermüller et al. <cit.> and Strong et al. <cit.> for learning different optimal scale-dependent total variation (TV) and total generalized variation (TGV) denoising parameters. Either using a training set or noise statistics of the image, scale-dependent problems face the risk of leading to overfitted optimal parameters (see, e.g., <cit.>). To gain generalization, alternative parameter functions may be considered that lie between a scalar and a fully scale-dependent weight. Examples of those intermediate functions are patch-dependent parameters, dictionary weights or basis functions coefficients.

One of the most challenging aspects in bilevel imaging learning, either with scalar or scale-dependent parameters, is the derivation of necessary optimality conditions that characterize the optimal solutions sharply and that may be used for numerical purposes. The approach pursued in <cit.> was based on a local regularization of the nonsmooth terms and an asymptotic analysis thereafter, yielding a C-stationarity system that does not fully characterize optimal parameters. A related approach, based on a dual reformulation of the lower level instance, was studied in <cit.>. Alternatively, a direct nonsmooth approach was considered in <cit.> and <cit.> to learn point spread functions in blind deconvolution models and the weight in front of the fidelity term in denoising models, respectively. In both such cases, the parameter affects the fidelity term and does not alter the structure of the primal-dual reformulation of the lower-level problem. Based on variational analysis tools, an M-stationarity system was then obtained.

Our aim in this paper consists in deriving sharp optimality conditions of Mordukhovich (M-) and Bouligand (B-) type for total variation bilevel learning problems, when the scale-dependent parameter appears within the regularizer. In such cases, differently from <cit.>, the nonnegativity constraint of the parameter alters the structure of the primal-dual characterization of the variational denoising model, leading to the appearance of not only a biactive, but also a triactive set. This involved nonsmoothness significantly complicates the derivation of the corresponding Fréchet and Mordukhovich normal cones, as well as the nonsmooth analysis of the lower-level solution mapping.

Firstly, after introducing the bilevel problem, we are able to reformulate it as a generalized mathematical program with equilibrium constraints (GMPEC). We rigorously characterize the corresponding Mordukhovich generalized normal cone (<Ref>) and verify a constraint qualification condition. This leads to the derivation of an M-stationarity system for the characterization of local minima (<Ref>).
Thereafter, we focus on B-stationarity conditions, which require the directional differentiability of the solution mapping. In that respect, we are able to prove the Lipschitz continuity and directional differentiability of the solution operator, and characterize the directional derivative by means of a variational inequality (<Ref>). In case of strict complementarity, the directional derivative becomes the Fréchet one and is characterized by a linear system. Also the Bouligand subdifferential of the solution operator is characterized by means of a general linear problem, based on a splitting of the biactive set and a properly defined subspace (<Ref>).

A further challenge is related with the numerical solution of the bilevel learning problems. Based on the properties of a local regularization technique, a numerical optimization algorithm of second-order type was considered in <cit.> and further explored in <cit.>, where dynamic sampling techniques were taken into account to accelerate the algorithm in the case of large datasets. A new approach handling the bilevel problem using a non-smooth setting can be found in the work of Ochs et. al. <cit.>, where the authors make use of algorithmic differentiation to find optimal parameters in variational models. An additional contribution of this paper is related to the solution of the bilevel instances by means of a nonsmooth trust-region algorithm. Specifically, based on the results in <cit.> and the Bouligand subdifferential characterization, we devise a two-phase trust-region algorithm for the solution of the problems (<Ref>). In the first phase, the Bouligand subdifferential is used in the quadratic trust-region subproblem, while in the second one, when the trust-region radius becomes small enough, a regularized gradient is utilized. We verify the performance of the proposed algorithm by means of scalar and scale-dependent experiments.

The paper is organized as follows. We will describe the bilevel learning problem in <Ref> along with the standing assumptions. In <Ref> we will characterize the Mordhukovich generalized normal cone and verify an appropriate constraint qualification condition, that will enable us to derive an M-stationary optimality system for the bilevel learning problem. In <Ref> we study the Lipschitz continuity and directional differentiability of the solution operator, with a proper characterization of the directional derivative, to obtain B-stationarity conditions for the bilevel problem. Thereafter, in <Ref> we characterize the Bouligand subdifferential of the solution mapping by means of a generalized linear system. This characterization of the Bouligand subdifferential will be of use in <Ref>, where we propose a non-smooth trust-region solution algorithm.  Finally, detailed numerical experiments will be provided in <Ref>.

§ PROBLEM FORMULATION

Variational image denoising is a known ill-posed inverse problem. Indeed, by considering an image as a $m_1\times m_2$ pixel matrix and mapping this matrix into a vector of length $m=m_1m_2$, a variational denoising model can be formulated as follows:

\begin{equation}
\min_{u\in\R^m} \;\Ecal(u) = \phi(u,f) + \mathcal{R}(u,\alpha),
\end{equation}

where $f$ is the noise contaminated image, $\alpha$ is a parameter that balances a data fidelity term $\phi$ and a regularization term $\mathcal{R}$.

This paper will focus on the case where the regularization term corresponds to the Isotropic Total Variation Semi-norm (TV). To characterize the TV regularizer, let us introduce the discrete gradient operator $\Kbb : \R^m\to\R^{n\times 2}$

\begin{equation*}
\Kbb u = (\Kbb_x u, \Kbb_y u)\in \R^n\times\R^n.
\end{equation*}

The matrix $\Kbb_x$ computes the difference in all neighboring pixel intensities in the $x$-direction and $\Kbb_y$ in the $y$-direction. Taking $u\in\R^m$ as a given image, the isotropic discrete total variation is then defined by $\|\Kbb u\|_{2,1} := \sum_{j=1}^n \|(\Kbb u)_j\|,$
where $\|\cdot\|$ is the euclidean norm and $(\Kbb u)_j$ corresponds to the $j$-th row of $\Kbb u$.
This regularizer induces sparsity on the gradients of the image, which favors piecewise constant images with sparse edges <cit.>.

Now, regarding the balance parameter $\alpha$, it can be considered as a scalar value in the case we assume a uniform noise distribution across the entire image. However, such assumption is known to be unrealistic for most applications. Therefore, we will explore different types of balance parameters, starting with a scale-dependent parameter $\alpha\in\R^n_+$ and the following form of the regularization term:

\begin{equation}\label{eq:data_fidelity_term}
\mathcal{R}(u,\alpha) = \sum_{j=1}^n\alpha_j\|(\Kbb u)_j\|.
\end{equation}

With these considerations in mind, and assuming a $C^2$ and strongly convex data fidelity term $\phi:\R^m\to\R$, our variational denoising model reads as follows:

\begin{equation}\label{eq:sd_denoising_model}
\min_{u\in\R^n} \;\mathcal{E}(u) := \phi(u,f) + \sum_{j=1}^n\alpha_j\|(\Kbb u)_j\|.
\end{equation}

In certain circumstances, such scale-dependent parameter may lead to an overfitted solution, which is why an intermediate approach is often desirable. One of such approaches considers a parameter $\alpha\in\R^p$ with $p<<n$ as piecewise constant in different patches of the image. Along the paper, we will derive the theoretical results only for the scale-dependent model, since they easily extend to the  patch-based and scalar models.

We will focus on bilevel strategies for finding optimal regularization parameters for the image denoising model <ref>. That is, given a training set containing $N$ pairs of damaged and “noise-free” pairs of images, we will find an optimal regularization parameter $\alpha\in\R^n_+$ that minimizes a continuously differentiable, proper and strongly convex loss function $J:\R^m\to\R$, that measures the quality of the reconstruction with respect to the ground truth images in the training set.

Using training set pairs $(\bar{u}_i,f_i)$, $i=1,\dots,N$, with noisy images $f_i$ and ground truth images $\bar{u}_i$, our bilevel parameter learning problem reads as follows:
\beginmini!
α∈^n, α≥0∑_i=1^N J(u_i,u̅_i)

u_i ∈_u∈^m ϕ(u,f_i) + ∑_j=1^nα_j(u)_j ,    i=1,…,N
\endmini!

Regarding the spatially dependent TV denoising problem <ref>, we can guarantee the existence of a unique solution as well as characterize it through a necessary and sufficient optimality condition.

Problem <ref> has a unique solution $u^*\in\R^m$. Moreover, a necessary and sufficient condition is given by the following variational inequality of the second kind

\begin{equation}\label{eq:nec_suff_condition_vi}
\scalar{\phi'(u^*)}{v-u^*} + \sum_{j=1}^n\alpha_j\|(\Kbb v)_j\| - \sum_{j=1}^n\alpha_j\|(\Kbb u^*)_j\|\ge 0,\;\forall v \in\R^m.
\end{equation}

The strong convexity of $\phi$, along with the convexity of the total variation semi-norm, yields the strong convexity of the lower level optimization problem. Consequently this problem has a unique optimizer $u^*$. The variational inequality then follows by standard arguments (see, e.g., <cit.>)

Using duality techniques <cit.>, the variational inequality of the second kind (<ref>) can be equivalently written in primal-dual form, yielding the following reformulation of the lower-level problem:

\begin{align}
\phi'(u) + \Kbb^\top q &= 0,\\
\scalar{q_j}{(\Kbb u)_j} - \alpha_j\|(\Kbb u)_j\|&=0,\;&\forall j=1,\dots,n,\\
\|q_j\|-\alpha_j&\le 0,\;&\forall j=1,\dots,n,
\end{align}

Consequently, the bilevel parameter learning problem, for a single training pair, can be written as:
\beginmini!
ϕ'(u) + ^⊤q = 0
q_j(u)_j - α_j(u)_j= 0,    ∀j=1,…,n
q_j-α_j≤0,    ∀j=1,…,n
α_j≥0,    ∀j=1,…,n.
\endmini!
For the sake of clarity in the exposition, we restrict hereafter the analysis to the case of a single training pair. The results are, however, easily extendable to larger training sets.

§ MORDUKHOVICH STATIONARITY

In this section we address stationarity conditions for the bilevel problem <ref> using variational analysis tools. Specifically, we reformulate the bilevel problem as a generalized mathematical program with equilibrium constraints (GMPEC) and verify a constraint qualification condition based on the variational geometry of the solution set.

A generalized mathematical program with equilibrium constraints may be formulated in the following general form:

0∈F_1(x,y) + Q(F_2(x,y))

where $f$ is locally Lipschitz on $\R^n\times\R^m$, $F_1: \R^n\times\R^m \to \R^m$ and $F_2: \R^n\times\R^m \to \R^l$ are continuously differentiable, $\omega\subset \R^n\times\R^m$ is non-empty and closed and $Q:\R^l \rightrightarrows \R^m$ is a multifunction with closed graph. A constraint qualification condition for GMPEC that guarantee the existence of KKT multipliers and its corresponding stationarity system is privede in the next result.

Let $(x^*,y^*)$ be a local solution of <ref> and the following constraint qualification

\begin{multline}\label{eq:outrata_cq_conditions}
\left.
\begin{aligned}
\begin{bmatrix}
(\nabla_x F_2(x^*,y^*))^\top & -(\nabla_x F_1(x^*,y^*))^\top\\
(\nabla_y F_2(x^*,y^*))^\top & -(\nabla_y F_1(x^*,y^*))^\top
\end{bmatrix}
\begin{bmatrix}
\end{bmatrix}\in -N^M_\omega(x^*,y^*),\\
(w,z)\in N^M_{\gph Q}(F_2(x^*,y^*),-F_1(x^*.y^*))
\end{aligned}\right\}\text{implies}\begin{cases}
\end{cases}
\end{multline}

holds true. Then there exist a pair $(\xi,\eta)\in\partial f(x^*,y^*)$, a pair $(\gamma,\delta)\in N^M_\omega(x^*,y^*)$, and a KKT pair $(w^*,z^*)\in N^M_{\gph Q}(F_2(x^*,y^*),-F_1(x^*,y^*))$ such that

\begin{align*}
0 = \xi + (\nabla_x F_2(x^*,y^*))^\top w^* - (\nabla_x F_1(x^*,y^*))^\top z^* + \gamma,\\
0 = \eta + (\nabla_y F_2(x^*,y^*))^\top w^* - (\nabla_y F_1(x^*,y^*))^\top z^* + \delta.
\end{align*}

where $N^M_{\gph Q}$ stands for Mordukhovich generalized normal cone to the graph of $Q$.

In our case, using the primal-dual formulation <ref>, the constraints may be written as

\begin{equation}\label{eq:generalized_equation}
0\in \phi'(u) + Q(\alpha,u),
\end{equation}

where $Q:\R^n_+\times\R^m \rightrightarrows \R^m$ is the set-valued operator associated to the subdifferential of the Euclidean norm, i.e.,

\begin{equation}\label{eq:set_valued_Q}
Q(\alpha,u) := \left\{\Kbb^\top q:q\in \R^{n\times 2} , \begin{cases}
q_j = \alpha_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|},\;&\text{ if } \|(\Kbb u)_j\| \neq 0,\alpha_j \ge 0,\\
% q_j = 0,\;&\text{ if }\alpha_j=0,\\
\|q_j\|\le \alpha_j,\;&\text{ if }\|(\Kbb u)_j\|=0,\alpha_j \ge 0.
\end{cases}\right\}
\end{equation}

Equivalently, by making use of the definition of the graph of the multifunction $Q$ we can rewrite <ref> as follows

\begin{align}
\phi'(u) + \Kbb^\top q &= 0,\\
(\alpha,u,\Kbb^\top q)&\in \gph\;Q,\\
\end{align}

where $\gph Q:=\{(\alpha,u,\Kbb^\top q)\in\R^n_+\times\R^m\times\R^m: \Kbb^\top q \in Q(\alpha,u)\}$.

The constraint qualification in <ref> guarantees the existence of multipliers that allow the derivation of a stationarity system. In <cit.> such constraint qualification condition was circumvented by using Robinson strong regularity of the set valued equation. In contrast, the structure of the multifunction <ref>, depending both on the regularization parameter $\alpha$ and the image $u$, prevent us from using the same shortcut.

Using the structure of the set valued operator $Q$ presented in <ref>, let us introduce the following notation for the inactive, strongly active, biactive, zero-inactive and triactive sets, respectively:

\begin{align*}
\Ical(\alpha,u) &:=\{j\in\{1,\dots,n\}:(\Kbb u)_j\neq 0,\; \alpha_j>0\},\\
\Acal_s(\alpha,u) &:=\{j\in\{1,\dots,n\}:\|q_j\|<\alpha_j\},\\
\Bcal(\alpha,u) &:= \{j\in\{1,\dots,n\}:\|q_j\|=\alpha_j,\;(\Kbb u)_j = 0,\;\alpha_j > 0\},\\
\Ical_0(\alpha,u) &:=\{j\in\{1,\dots,n\}:(\Kbb u)_j\neq 0,\; \alpha_j=0\},\\
\Tcal(\alpha,u) &:= \{j\in\{1,\dots,n\}:\|q_j\|=\alpha_j,\;(\Kbb u)_j = 0,\;\alpha_j = 0\}.
\end{align*}

We will omit the arguments in the set notation whenever they can be inferred from the context. Let us note that the condition in the strongly active set $\Acal_s$ implies a strict positive parameter $\alpha_j>0$ for this index set.

The constraint qualification condition <ref> makes use of key fundamental definitions from Mordukhovich's generalized calculus. In particular, the Mordukhovich normal cone to the graph of the mutifunction $Q$. For the reader's convenience we will lay some basic concepts in variational geometry that will be used through this section. For a more rigorous review refer to, e.g., <cit.>

Let $C\subset\R^n$ and $\bar{x}\in cl C$. The Bouligand Tangent Cone to $C$ at $x$ is defined as

\begin{equation}
T_C(\bar{x}):=\{d\in\R^n:\exists t_k \to 0,\exists x_k\in C:\frac{x_k-\bar{x}}{t_k}\to d\}.
\end{equation}

Let $C\subset\R^n$ and $\bar{x}\in cl C$. The Fréchet Normal Cone is defined as the polar to $T_C(\bar{x})$, i.e.,

\begin{equation}
N^F_C(\bar{x})=\left\{d\in\R^n : \limsup_{x\overset{C}{\to}\bar{x}}\frac{\scalar{d}{x-\bar{x}}}{\|x-\bar{x}\|}\le 0\right\}=[T_C(\bar{x})]^\circ.
\end{equation}

Let $C\subset\R^n$ and $\bar{x}\in cl C$. The Mordukhovich Normal Cone is defined as

\begin{equation}
N^M_C(\bar{x}):=\{d\in\R^n:\exists x_k\to\bar{x},d_k\in N^F_C(x_k)\text{ s.t. }d_k\to d\}.
\end{equation}

If the set $C$ in <ref> is convex, $N^M_C(\bar{x})$ amounts to the standard normal cone to $C$ at $\bar{x}$. Otherwise, this cone is in general non-convex.

In the following lemmata we obtain the Bouligand tangent cone, the Fréchet normal cone and the Mordukhovich normal cone to the graph of the multifunction $Q$.

The Bouligand tangent cone to the graph of $Q$, described in <ref>, is given by

\begin{multline*}\label{eq:contingent_cone}
T_{\gph Q}(\alpha,u,\Kbb^\top q) =\\ \left\{ (\delta_\alpha,\delta_u,\Kbb^\top\delta_q):
\left \{
\begin{aligned}
& (\delta_q)_j - (\delta_\alpha)_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|} - \alpha_jT_j(\Kbb \delta_u)_j=0,&&\text{if }j\in\Ical,\\
& (\Kbb\delta_u)_j = 0, &&\text{if }j\in\Acal_s,\\
& \left.
\begin{aligned}
& (\Kbb\delta_u)_j = 0,\scalar{(\delta_q)_j}{q_j}\le \alpha_j(\delta_\alpha)_j\;\vee \\
& (\Kbb\delta_u)_j = \tilde{c}q_j (\tilde{c}\ge 0),\scalar{(\delta_q)_j}{q_j}=\alpha_j(\delta_\alpha)_j
\end{aligned} \right\}  &&\text{if }j\in\Bcal,\\
& (\delta_q)_j - (\delta_\alpha)_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|}=0,\;(\delta_\alpha)_j\ge 0&&\text{if }j\in\Ical_0,\\
& \left.
\begin{aligned}
& (\delta_\alpha)_j \ge 0,(\Kbb\delta_u)_j\in\R^2\backslash\{0\},(\delta_q)_j - (\delta_\alpha)_j\frac{(\Kbb \delta_u)_j}{\|(\Kbb \delta_u)_j\|}=0\;\vee \\
& (\delta_\alpha)_j \ge 0, (\Kbb \delta_u)_j = 0, \|(\delta_q)_j\| - (\delta_\alpha)_j \le 0\\
\end{aligned} \right\}  &&\text{if }j\in\Tcal,
\end{aligned} \right.
\right\}
\end{multline*}

\begin{equation*}
T_j(\Kbb v)_j = \frac{(\Kbb v)_j}{\|(\Kbb u)_j\|} - \frac{(\Kbb u)_j(\Kbb u)_j^\top(\Kbb v)_j}{\|(\Kbb u)_j\|^3}, \text{ for } v\in\R^n.
\end{equation*}

The tangent cone to the graph of the multifunction $Q$ is defined as

\begin{multline*}
T_{\gph Q}(\alpha,u,\Kbb^\top q) = \{(\delta_\alpha,\delta_u,\Kbb^\top \delta_q)\in\R^n\times\R^m\times\R^m: \exists t_k\to 0, (\alpha_k,u_k,\Kbb^\top q_k)\in \gph Q :\\\frac{1}{t_k}((\alpha_k,u_k,\Kbb^\top q_k)-(\alpha,u,\Kbb^\top q))\to (\delta_\alpha,\delta_u,\Kbb^\top \delta_q)\}.
\end{multline*}

In order to calculate this cone, we split the analysis into different cases, according to the definition of the multifunction $Q$.

Case 1: $\bm{j\in\Ical.}$
In this index set the dual variable can be uniquely characterized. According to <ref>, the following equation is fulfilled:

\begin{equation}\label{eq:lyusternik_manifold}
h_j(\alpha,u,\Kbb^\top q) = q_j-\alpha_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|}=0.
\end{equation}

Using Lyusternik's theorem <cit.>, the $j-$th component of the tangent direction then satisfies

\begin{equation*}
(\delta_q)_j-(\delta_\alpha)_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|} - \alpha_jT_j(\Kbb\delta_u)_j=0.
\end{equation*}

Case 2: $\bm{j\in\Acal_s.}$
In this index set we know $\|q_j\|<\alpha_j$. Therefore, this point can only be approximated by taking sequences in the strongly active set. For $n$ sufficiently large we then take sequences such that $(\Kbb u_k)_j=0$, $\|(q_k)_j\|<(\alpha_k)_j$. Taking the limit in $(\Kbb u_k)_j=0$, as $k \to \infty$, yields $(\Kbb \delta_u)_j = 0$. For the dual variable, let us take the sequence $(q_k)_j = q_j + t_k d$ with arbitrary $d\in\R^2$. It then follows that

\begin{equation*}
(\delta_q)_j = \lim_{t_k\to 0} \frac{(q_k)_j-q_j}{t_k} = d.
\end{equation*}

Since we took an arbitrary direction $d$ it yields $(\delta_q)_j\in\R^2$. Similarly, we get that $(\delta_\alpha)_j\in\R$.

Case 3: $\bm{j\in\Bcal.}$
There are three possible approximations to a point in this index set, via inactive, strongly active or biactive sequences. Since $\alpha_j>0$ in all these three cases, we can approximate it by sequences $(\alpha_k)_j<\alpha_j$ or $(\alpha_k)_j>\alpha_j$. Consequently, we get $(\delta_\alpha)_j\in\R$.

Now, if we approach a biactive point using a sequence in the inactive set, this sequence satisfies $(\Kbb u_k)_j\neq 0$. Taking in particular $(\Kbb u_k)_j = (\Kbb u)_j + t_k (\Kbb\delta_{u_k})_j = t_k (\Kbb\delta_{u_k})_j$, then the sequence of dual variables has the following form

\begin{equation}\label{eq:tan_cone_b_3}
(q_k)_j = (\alpha_k)_j\frac{(\Kbb u_k)_j}{\|(\Kbb u_k)_j\|} = (\alpha_k)_j\frac{(\Kbb\delta_{u_k})_j}{\|(\Kbb\delta_{u_k})_j\|}.
\end{equation}

Furthermore, considering the following product

\begin{equation*}
\scalar{(q_k)_j}{(\Kbb\delta_{u_k})_j} =\bigscalar{(\alpha_k)_j\frac{(\Kbb\delta_{u_k})_j}{\|(\Kbb\delta_{u_k})_j\|}}{(\Kbb\delta_{u_k})_j} = (\alpha_k)_j \|(\Kbb\delta_{u_k})_j\|,
\end{equation*}

and taking the limit as $k \to \infty$ we get $\scalar{q_j}{(\Kbb\delta_u)_j} = \alpha\|(\Kbb\delta_u)_j\|$. Recalling that in this index set $\|q_j\|=\alpha_j > 0$, we know both vectors are collinear, i.e.,

\begin{equation}\label{eq:tan_cone_b_0}
(\Kbb\delta_u)_j=\tilde{c} q_j,\;\text{ for some }\tilde{c}\ge 0.
\end{equation}

Using that $\|q_j\|=\alpha_j$, the following product holds

\begin{align*}
\bigscalar{\frac{(q_k)_j-q_j}{t_k}}{q_j} &= \frac{1}{t_k}(\scalar{(q_k)_j}{q_j}-\scalar{q_j}{q_j})= \frac{1}{t_k}(\scalar{(q_k)_j}{q_j} - \scalar{(q_k)_j}{(q_k)_j} + \scalar{(q_k)_j}{(q_k)_j} - \alpha_j^2),\\
&= \bigscalar{(q_k)_j}{\frac{q_j-(q_k)_j}{t_k}} + \frac{(\alpha_k)_j^2-\alpha_j^2}{t_k}.
\end{align*}

The last equality holds since from <ref>, we get the product $\scalar{(q_k)_j}{(q_k)_j} = (\alpha_k)_j^2$. Taking the limit as $t_k\to 0$ we get

\begin{equation}\label{eq:tan_cone_b_1}
\scalar{(\delta_q)_j}{q_j} = \alpha_j(\delta_\alpha)_j.
\end{equation}

Now, if the approximation is done through a sequence of strongly active points, we know the sequence satisfies $(\Kbb u_k)_j=0$ and $\|q_k\|<(\alpha_k)_j$. In this case we know $(\Kbb \delta_u)_j = 0$ and, using the Cauchy-Schwarz inequality, we get

\begin{equation*}
\bigscalar{\frac{(q_k)_j-q_j}{t_k}}{q_j} \le \frac{\alpha_j}{t_k}(\|(q_k)_j\|-\|q_j\|) < \alpha_j\left(\frac{(\alpha_k)_j-\alpha_j}{t_k}\right),
\end{equation*}

which implies that $\scalar{(\delta_q)_j}{q_j} \le \alpha_j(\delta_\alpha)_j$.

Finally, approximating through a biactive set sequence, we know again $(\Kbb\delta_u)_j=0$ and, estimating the product

\begin{equation*}
\bigscalar{\frac{(q_k)_j-q_j}{t_k}}{q_j} \le \frac{\alpha_j}{t_k}(\|(q_k)_j\|-\|q_j\|) = \alpha_j\left(\frac{(\alpha_k)_j-\alpha_j}{t_k}\right),
\end{equation*}

we get

\begin{equation}\label{eq:tan_cone_b_2}
\scalar{(\delta_q)_j}{q_j} \le \alpha_j(\delta_\alpha)_j.
\end{equation}

Case 4: $\bm{j\in\Ical_0.}$
We can approximate a point in the zero-inactive set by sequences in the inactive set. Therefore, considering a sequence $(\Kbb u_k)_j = (\Kbb u)_j + t_k v$ with $v\in\R^2$ arbitrary we have

\begin{equation}\label{eq:kdu_zero_inactive}
(\Kbb\delta_u)_j = \lim_{t_k\to 0}\frac{(\Kbb u_k)_j-(\Kbb u)_j}{t_k} = v.
\end{equation}

Since we took $v\in\R^2$ arbitrary, it follows that $(\Kbb\delta_u)\in\R^2$. Furthermore, since $q_j = \alpha_j((\Kbb u)_j/\|(\Kbb u)_j\|)$ and $\alpha_j=0$, then $q_j=0$ in this index set. Considering the sequence $(\alpha_k)_j= \alpha_j + t_k h$ we obtain

\begin{equation}\label{eq:tan_cone_I0_2}
(\delta_q)_j = \lim_{t_k\to 0}\frac{(q_k)_j}{t_k} = \lim_{t_k\to 0}h\frac{(\Kbb u_k)_j}{\|(\Kbb u_k)_j\|} = (\delta_\alpha)_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|}.
\end{equation}

Since $\alpha_j=0$, the only valid approximations are the ones coming from positive elements, thus $(\delta_\alpha)_j\ge 0$. Another possible approximation can be done through zero-inactive points, meaning $(q_k)_j=0$ and $(\alpha_k)_j=0$. This case implies $(\delta_q)_j=0$, $(\delta_\alpha)_j=0$ and $(\Kbb\delta_u)_j\in\R^2$, which is a particular case of <ref>.

Case 5: $\bm{j\in\Tcal.}$
There are four ways to approach a triactive point. As in the zero-inactive case, all valid approximations come from $(\alpha_k)_j\ge 0$, which again implies $(\delta_\alpha)_j\ge 0$. Similarly to the zero-inactive case <ref>, we also get $(\Kbb\delta_u)_j\in\R^2$.

Approximating through a sequence in the inactive set, we obtain

\begin{equation}\label{eq:tan_cone_I0_1}
(\delta_q)_j = \lim_{t_k\to 0}\frac{1}{t_k}\left((\alpha_k)_j\frac{(\Kbb u_k)_j}{\|(\Kbb u_k)_j\|}\right) = (\delta_\alpha)_j\frac{(\Kbb\delta_u)_j}{\|(\Kbb\delta_u)_j\|}.
\end{equation}

Likewise, the approximation can be made using zero-inactive points. In this case $(\Kbb\delta_u)_j \neq 0$, $(q_k)_j=0$ and $(\alpha_k)_j=0$. From this sequence we can derive $(\delta_q)_j=0$, $(\delta_\alpha)_j=0$ and $(\Kbb \delta_u)_j\in\R^2$, which is included in <ref>.

Moving forward, we can also approximate through strongly active points, i.e., $(\Kbb u_k)_j=0$ and $\|(q_k)_j\|<(\alpha_k)_j$. From this sequence we know $(\Kbb \delta_u)_j=0$ and $(\delta_\alpha)_j\ge 0$ and the dual variable will have the following form

\begin{equation}\label{eq:tan_cone_I0_3}
\|(\delta_q)_j\| = \left\|\lim_{t_k\to 0}\frac{(q_k)_j}{t_k}\right\| = \lim_{t_k\to 0} \frac{1}{t_k}\|(q_k)_j\| \le \lim_{t_k\to 0}\frac{1}{t_k}(\alpha_k)_j = (\delta_\alpha)_j,
\end{equation}

yielding $\|(\delta_q)_j\|\le(\delta_\alpha)_j$.

Finally, we consider an approximation through biactive points, meaning $(\Kbb u_k)_j=0$ and $\|(q_k)_j\|=(\alpha_k)_j$. We the obtain $(\Kbb \delta_u)_j =0$, $(\delta_\alpha)_j\ge 0$ and

\begin{equation*}
\|(\delta_q)_j\| = \left\|\lim_{t_k\to 0}\frac{(q_k)_j}{t_k}\right\| = \lim_{t_k\to 0} \frac{1}{t_k}\|(q_k)_j\| = \lim_{t_k\to 0}\frac{1}{t_k}(\alpha_k)_j = (\delta_\alpha)_j,
\end{equation*}

which is a particular case of <ref>. Consequently, $\|(\delta_q)_j\|=(\delta_\alpha)_j$.

The Fréchet normal cone to the graph of $Q$, described in <ref>, is given by

\begin{multline}\label{eq:frechet_normal_cone}
N^F_{\gph Q}(\alpha,u,\Kbb^\top q) =\\ \left\{(\vartheta,\Kbb^\top \mu,p):
\left\{
\begin{aligned}
& \mu_j + \alpha_jT_j(\Kbb p)_j = 0,&&\text{if }j\in\Ical,\\
& \vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}=0,&&\text{if }j\in\Ical,\\
& \vartheta_j=0,\;(\Kbb p)_j=0  &&\text{if }j\in\Acal_s,\\
& \vartheta_j+c\alpha_j = 0,(\Kbb p)_j = c q_j (c\ge 0),\scalar{\mu_j}{q_j} \le 0,&&\text{if }j\in\Bcal,\\
& \vartheta_j + \frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|} \le 0,\quad \mu_j = 0,&&\text{if }j\in\Ical_0,\\
& \vartheta_j\le 0,\;(\Kbb p)_j = 0,\;\mu_j=0,&&\text{if }j\in\Tcal.
\end{aligned}
\right.
\right\}
\end{multline}

Using the definition of the Fréchet normal cone for this problem

\begin{equation*}
N^F_{\gph Q}(\alpha,u,\Kbb^\top q) = \{(\vartheta,\Kbb^\top \mu,p)\in\R^n\times\R^m\times\R^m:\scalar{(\vartheta,\Kbb^\top \mu,p)}{(\delta_\alpha,\delta_u,\Kbb^\top \delta_q)}\le 0\}.
\end{equation*}

Indeed, we can rewrite this inequality as

\begin{equation*}
\sum_{j=1}^n ((\delta_\alpha)_j\vartheta_j + \scalar{(\Kbb \delta_u)_j}{\mu_j}+\scalar{(\delta_q)_j}{(\Kbb p)_j})\le 0.
\end{equation*}

Using this representation, along with the tangent cone characterization from <ref>, we analyze the different cases according to their index set.

Case 1: $\bm{j\in\Ical.}$
Using the characterization of the elements in the tangent cone, we have

\begin{equation*}
(\delta_\alpha)_j\vartheta_j+\scalar{(\Kbb \delta_u)_j}{\mu_j} + \bigscalar{(\delta_\alpha)_j\frac{(\Kbb u)_j}{\|(\Kbb u)_j\|}}{(\Kbb p)_j} + \scalar{\alpha_j T_j(\Kbb \delta_u)_j}{(\Kbb p)_j}\le 0.
\end{equation*}

Rearranging the terms and using the symmetry of $T_j$,

\begin{equation*}
(\delta_\alpha)_j \left(\vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}\right) + \scalar{(\Kbb \delta_u)_j}{\mu_j+\alpha_jT_j(\Kbb p)_j}\le 0.
\end{equation*}

Since $(\Kbb \delta_u)_j\in\R^2$ and $(\delta_\alpha)_j\in\R$, it necessarily must hold

\begin{equation*}
\vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}=0,\qquad \mu_j+\alpha_jT_j(\Kbb p)_j = 0.
\end{equation*}

Case 2: $\bm{j\in\Acal_s.}$
In this index set we know that $(\Kbb \delta_u)_j=0$, $(\delta_\alpha)_j\in\R$ and $(\delta_q)_j\in\R^2$. Consequently the product reads

\begin{equation*}
(\delta_\alpha)_j\vartheta_j + \scalar{(\delta_q)_j}{(\Kbb p)_j}\le 0,
\end{equation*}

and we obtain that $(\Kbb p)_j=0$ and $\vartheta_j=0$.

Case 3: $\bm{j\in\Bcal.}$
In this index set there are two conditions on the normal directions. For the first one we take $(\Kbb \delta_u)_j=0$ and the cone inequality reads

\begin{equation}\label{eq:frechet_bcal_6}
(\delta_\alpha)_j\vartheta_j + \scalar{(\delta_q)_j}{(\Kbb p)_j}\le 0,\quad\forall (\delta_\alpha)_j,(\delta_q)_j\text{ s.t. }\scalar{(\delta_q)_j}{q_j}\le \alpha_j(\delta_\alpha)_j.
\end{equation}

Taking in particular, $(\delta_\alpha)_j=0$ we get

\begin{equation*}
\scalar{(\delta_q)_j}{(\Kbb p)_j}\le 0,\forall (\delta_q)_j\text{ s.t. }\scalar{(\delta_q)_j}{q_j}\le 0.
\end{equation*}

Therefore, $(\Kbb p)_j = c q_j$ with $c\ge 0$. Using this result in <ref> for the particular case $(\delta_\alpha)_j = \frac{1}{\alpha_j}\scalar{(\delta_q)_j}{q_j}$, we obtain

\begin{equation*}
0\ge(\delta_\alpha)_j\vartheta_j + \scalar{(\delta_q)_j}{c q_j} =(\delta_\alpha)_j(\vartheta_j + c\alpha_j),
\end{equation*}

from where $\vartheta_j +c\alpha_j=0$ holds. The resulting cone reads

\begin{equation}\label{eq:frechet_bcal_1}
\vartheta_j+c\alpha_j=0,\quad (\Kbb p)_j = c q_j,\quad c\ge 0, \quad \mu_j\in\R^2.
\end{equation}

For the second case we take $(\Kbb\delta_u) = c q_j\;(c\ge 0)$ in <ref>,

\begin{equation}\label{eq:frechet_bcal_7}
(\delta_\alpha)_j\vartheta_j + \scalar{c q_j}{\mu_j}+\scalar{(\Kbb p)_j}{(\delta_q)_j}\le 0,\forall (\delta_\alpha)_j\in\R,\;(\delta_q)_j\text{ s.t. }\scalar{(\delta_q)_j}{q_j}= \alpha_j(\delta_\alpha)_j,
\end{equation}

Again, considering $(\delta_\alpha)_j=0$ and $(\Kbb\delta_u)_j=0$, we get

\begin{equation*}
\scalar{(\delta_q)_j}{(\Kbb p)_j}\le 0,\forall (\delta_q)_j\text{ s.t. }\scalar{(\delta_q)_j}{q_j}= 0.
\end{equation*}

Consequently, $(\Kbb p)_j = c q_j$ with $c\in\R$. Using this result in <ref>, while keeping $(\delta_\alpha)_j=0$, yields $\tilde{c}\scalar{q_j}{\mu_j}\le 0$. Thanks to the positiveness of $\tilde{c}$ we then get $\scalar{q_j}{\mu_j}\le 0$. Now, applying all previous results in <ref> we get

\begin{equation*}
0\ge(\delta_\alpha)_j\vartheta_j +\scalar{c q_j}{(\delta_q)_j} = (\delta_\alpha)_j\vartheta_j +c\alpha_j(\delta_\alpha)_j,\forall (\delta_\alpha)_j\in\R, c\in\R,
\end{equation*}

yielding $\vartheta_j\in\R$. Therefore, the resulting cone for the second case reads

\begin{equation}\label{eq:frechet_bcal_2}
\vartheta_j \in\R,\quad (\Kbb p)_j = c q_j,\quad c\in\R.
\end{equation}

Finally, considering both cases we obtain

\begin{equation*}
\vartheta_j+c\alpha_j = 0\wedge\scalar{\mu_j}{q_j} \le 0 \wedge (\Kbb p)_j = c q_j \wedge c\ge 0.
\end{equation*}

Case 4: $\bm{j\in\Ical_0.}$
By using the characterization of the tangent cone in this index set, we have

\begin{equation*}
(\delta_\alpha)_j\left(\vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}\right) + \scalar{(\Kbb \delta_u)_j}{\mu_j} \le 0.
\end{equation*}

This relationship must hold for all $(\Kbb \delta_u)_j\in\R^2$ and $(\delta_\alpha)_j\ge 0$, which implies

\begin{equation*}
\vartheta_j + \frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}\le 0,\quad \scalar{(\Kbb \delta_u)_j}{\mu_j} \le 0.
\end{equation*}

Since $(\Kbb \delta_u)_j\in\R^2$ we get $\mu_j=0$.

Case 5: $\bm{j\in\Tcal.}$
For the first case in this index set we know the elements of the tangent cone satisfy

\begin{equation*}
(\delta_q)_j = (\delta_\alpha)_j\frac{(\Kbb \delta_u)_j}{\|(\Kbb \delta_u)_j\|},\;(\delta_\alpha)_j\ge 0,\;(\Kbb\delta_u)_j\in\R^2\backslash\{0\}.
\end{equation*}

Replacing these terms into the normal cone inequality,

\begin{equation*}
(\delta_\alpha)_j\left(\vartheta_j+\frac{\scalar{(\Kbb p)_j}{(\Kbb\delta_u)_j}}{\|(\Kbb\delta_u)_j\|}\right) + \scalar{\mu_j}{(\Kbb\delta_u)_j}\le 0,\;\forall (\delta_\alpha)_j\ge 0,(\Kbb\delta_u)_j\in\R^2\backslash\{0\}.
\end{equation*}

In particular, for $(\delta_\alpha)_j = 0$, we get that $\scalar{\mu_j}{(\Kbb\delta_u)_j}\le 0$, for all $(\Kbb \delta_u)_j\in\R^2\backslash\{0\}$, which implies that $\mu_j=0$. Moreover, thanks to the positiveness of $(\delta_\alpha)_j\ge 0$, we get

\begin{equation}\label{eq:frechet_bcal_4}
\vartheta_j + \frac{\scalar{(\Kbb p)_j}{(\Kbb\delta_u)_j}}{\|(\Kbb\delta_u)_j\|}\le 0,\forall (\Kbb\delta_u)_j\in\R^2\backslash\{0\}.
\end{equation}

Testing this inequality with $\pm(\Kbb\delta_u)_j$, we get $\vartheta_j\le\scalar{(\Kbb p)_j}{(\Kbb\delta_u)_j}\le -\vartheta_j$, which implies that the product $\scalar{(\Kbb p)_j}{(\Kbb\delta_u)_j}=0$ for all $(\Kbb\delta_u)_j\in\R^2\backslash\{0\}$, consequently $(\Kbb p)_j=0$. Using this result in <ref> yields also $\vartheta_j\le 0$.

Now, regarding the second condition, we know $\|(\delta_q)_j\| - (\delta_\alpha)_j\le 0, ~(\Kbb \delta_u)_j =0$ and $(\delta_\alpha)_j\ge 0.$
Since $(\Kbb \delta_u)_j=0$, it follows $\mu_j\in\R^2$. Using the normal cone inequality,

\begin{equation*}
0 \ge (\delta_\alpha)_j\vartheta_j + \scalar{(\delta_q)_j}{(\Kbb p)_j}\ge (\delta_\alpha)_j\vartheta_j - \|(\delta_q)_j\|\|(\Kbb p)_j\|\ge (\delta_\alpha)_j\vartheta_j - (\delta_\alpha)_j\|(\Kbb p)_j\|,
\end{equation*}

from where we derive $(\delta_\alpha)_j(\vartheta_j - \|(\Kbb p)_j\|)\le 0$. Along with the positivity of $(\delta_\alpha)_j$, this implies $\vartheta_j \le \|(\Kbb p)_j\|$. Finally, considering both conditions, the result is obtained.

The Mordukhovich normal cone to the graph of $Q$, described in <ref>, is given by

\begin{multline*}
N_{\gph Q}^{M}(\alpha,u,\mathbb K^\top q) =\\ \left\{ (\vartheta,\mathbb K^\top \mu,p):
\left \{
\begin{aligned}
& \mu_j + \alpha_jT_j(\Kbb p)_j = 0,&&\text{if }j\in\Ical,\\
& \vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}=0,&&\text{if }j\in\Ical,\\
& \vartheta_j=0,(\Kbb p)_j=0, &&\text{if }j\in\Acal_s,\\
& \left.
\begin{aligned}
& \vartheta_j=0,(\Kbb p)_j=0,\;\vee\\
& (\mathbb K p)_j=c q_j (c\in\R), ~\langle \mu_j,q_j \rangle =0,\;\vee \\
& \vartheta_j+c\alpha_j=0,(\mathbb K p)_j = c q_j (c\ge 0),\langle \mu_j,q_j \rangle\le 0.
\end{aligned} \right\}  &&\text{if }j\in\Bcal,\\
& \vartheta_j + \frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|} \le 0,\quad \mu_j = 0,&&\text{if }j\in\Ical_0,\\
& \left.
\begin{aligned}
& \vartheta_j=0,(\Kbb p)_j = 0,\;\vee \\
& \vartheta_j\le\|(\Kbb p)_j\|,\mu_j=0,\;\vee \\
& |\vartheta_j|\le \|(\Kbb p)_j\|,\scalar{\mu_j}{(\Kbb p)_j}\le 0
\end{aligned} \right\}  &&\text{if }j\in\Tcal.
%& \vartheta_j=0,(\Kbb p)_j = 0,&&\text{if }j\in\Bcal_0.
\end{aligned} \right.
\right\}
\end{multline*}

Let us recall the definition of the Mordukhovich normal cone for our problem

\begin{multline*}
N^M_{\gph Q}(\alpha,u,\Kbb^\top q) = \{(\vartheta,\Kbb^\top \mu,p):(\vartheta_k,\Kbb^\top \mu_k,p_k)\in N^F_{\gph Q}(\alpha_k,u_k,\Kbb^\top q_k):\\(\vartheta_k,\Kbb^\top \mu_k,p_k)\to (\vartheta,\Kbb^\top \mu,p),(\alpha_k,u_k,\Kbb^\top q_k)\to (\alpha,u,\Kbb^\top q)\}.
\end{multline*}

Considering limiting sequences to the inactive, strongly active and zero-inactive sets, the same directions as for the Fréchet normal cone are obtained. The differences lie in the biactive and triactive sets, where several approximations may be considered.

Case 1: $\bm{j\in\Bcal.}$
By taking approximation sequences in the inactive set, from <ref> we know

\begin{equation}\label{eq:m_cone_b0_6}
0 = (\mu_k)_j + (\alpha_k)_j\frac{(\Kbb p_k)_j}{\|(\Kbb u_k)_j\|} - (\alpha_k)_j\frac{(\Kbb u_k)\scalar{((\Kbb u_k)_j)}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|^3}.
\end{equation}

Multiplying <ref> with $(\Kbb p_k)_j$ yields

\begin{equation*}
\scalar{(\mu_k)_j}{(\Kbb p_k)_j} = (\alpha_k)_j\frac{\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}^2}{\|(\Kbb u_k)_j\|^3} - (\alpha_k)_j\frac{\scalar{(\Kbb p_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|}
\end{equation*}

Again, multiplying by $(\alpha_k)_j\|(\Kbb u_k)_j\|$ on both sides and recalling that $(q_k)_j = (\alpha_k)_j\frac{(\Kbb u_k)_j}{\|(\Kbb u_k)_j\|}$, we get

\begin{equation*}
(\alpha_k)_j\|(\Kbb u_k)_j\|\scalar{(\mu_k)_j}{(\Kbb p_k)_j} = \scalar{(q_k)_j}{(\Kbb p_k)_j}^2 - (\alpha_k)_j^2\|(\Kbb p_k)_j\|^2.
\end{equation*}

Taking the limit as $k\to\infty$ and recalling $(\alpha_k)_j=\|(q_k)_j\|$ in this index set, we obtain

\begin{equation*}
\scalar{q_j}{(\Kbb p)_j}^2 = \|q_j\|^2\|(\Kbb p)_j\|^2,
\end{equation*}

which implies that $(\Kbb p)_j = cq_j(c\in\R)$. Now, multiplying <ref> with $(q_k)_j$ we get the following product

\begin{align*}
\scalar{(\mu_k)_j}{(q_k)_j} &= (\alpha_k)_j\frac{\scalar{(q_k)_j}{(\Kbb u_k)_j}\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|^3} - (\alpha_k)_j\frac{\scalar{(q_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|},\\
&= (\alpha_k)_j^2\frac{\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|^2} - (\alpha_k)_j^2\frac{\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|^2}=0.\numberthis \label{eq:m_cone_b0_7}
\end{align*}

Taking the limit we get that $\scalar{\mu_j}{q_j}=0$.

Regarding $\vartheta_j$ we have $(\alpha_k)_j(\vartheta_k)_j + \scalar{(q_k)_j}{(\Kbb p)_j} = 0.$
Taking the limit as $k\to \infty$ we obtain

\begin{equation*}
0 = \alpha_j\vartheta_j + \scalar{q_j}{(\Kbb p)_j}= \alpha_j\vartheta_j + c\|q_j\|^2,(c\in\R),
\end{equation*}

which implies that $\vartheta_j=-c\alpha_j\in\R$.

Finally, when taking the approximation through the strongly active and biactive sets, the cone directions coincide with the Fréchet normal ones.

Case 2: $\bm{j\in\Tcal.}$
This index set can be approximated by sequences belonging either to the inactive, biactive, strongly active or zero-inactive sets.
Considering strongly active sequences, $(\Kbb p_k)_j=0$ and $(\vartheta_k)_j=0$. Taking the limit as $k\to\infty$ we get $(\Kbb p)_j=0$ and $\vartheta_j=0$ as well.

Likewise, when taking biactive sequences we get $(\vartheta_k)_j+c(\alpha_k)_j=0$, $(\Kbb p_k)_j=c(q_k)_j(c\ge 0)$ and $\scalar{(\mu_k)_j}{(q_k)_j}\le 0$. Again, taking the limit as $k\to\infty$ we get, since $q_j=0$ and $\alpha_j=0$, that $\mu_j\in\R^2$, $(\Kbb p)_j=0$ and $\vartheta_j=0$.

Furthermore, taking sequences in the zero-inactive set we have $(\mu_k)_j=0$, which implies that $\mu_j=0$. Using the Cauchy-Schwarz inequality we get

\begin{equation*}
0\ge (\vartheta_k)_j + \frac{\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|} \ge (\vartheta_k)_j -\|(\Kbb p_k)_j\|.
\end{equation*}

Taking the limit as $k\to\infty$ yields $\vartheta_j\le\|(\Kbb p)_j\|$.

When taking inactive sequences, we know that

\begin{align}
(\mu_k)_j + (\alpha_k)_j\left(\frac{I}{\|(\Kbb u_k)_j\|}-\frac{(\Kbb u_k)_j(\Kbb u_k)_j^\top}{\|(\Kbb u_k)_j\|^3}\right)(\Kbb p_k)_j &= 0,\label{eq:m_cone_b0_1}\\
(\vartheta_k)_j + \frac{\scalar{(\Kbb u_k)_j}{(\Kbb p_k)_j}}{\|(\Kbb u_k)_j\|}&= 0.\label{eq:m_cone_b0_2}
\end{align}

Applying Cauchy-Schwarz in <ref>, yields $|(\vartheta_k)_j|\le\|(\Kbb p_k)_j\|.$

Now, multiplying <ref> with $(\Kbb p_k)_j$ we obtain

\begin{equation}\label{eq:m_cone_b0_4}
\scalar{(\mu_k)_j}{(\Kbb p_k)_j} \le (\alpha_k)_j\left(\frac{\|(\Kbb p_k)_j\|^2\|(\Kbb u_k)_j\|^2}{\|(\Kbb u_k)_j\|^3}-\frac{\|(\Kbb p_k)_j\|^2}{\|(\Kbb u_k)_j\|}\right)=0.
\end{equation}

Taking the limit as $k\to\infty$ we get $\scalar{\mu_j}{(\Kbb p)_j}\le 0$, finishing the proof.

Let $J:\R^m\to\R$ be continuously differentiable, $\phi:\R^m\to\R$ twice continuously differentiable and strongly convex, and $(\alpha^*,u^*,q^*)$ be a local solution to <ref>. Then there exist KKT multipliers $(\vartheta,\Kbb^\top \mu,p)$ such that

\begin{align}
\phi'(u^*)+\Kbb^\top q^* = 0,\\
\scalar{q^*_j}{(\Kbb u^*)_j}-\alpha^*_j\|(\Kbb u^*)_j\| =0, &&&\forall j=1,\dots,n,\\
\|q^*_j\| \le \alpha_j^*, &&&\forall j=1,\dots,n,\\
\phi''(u^*)^\top p -\Kbb^\top \mu - \nabla J(u^{*})=0,\\
\vartheta + \rho =0,\\
0 \leq \alpha^* \perp \rho \geq 0,\\
(\vartheta,\Kbb^\top \mu,p)\in N_{\gph Q}^M(\alpha^*,u^*,\Kbb^\top q^*)
% \mu_j + \alpha_jT_j(\Kbb p)_j = 0,&&&\text{if }j\in\Ical,\\
% \vartheta_j+\frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|}=0,&&&\text{if }j\in\Ical,\\
% \vartheta_j =0,(\Kbb p)_j=0,&&&\text{if }j\in\Acal_s,\\
% \left.
% \begin{aligned}
%   \vartheta_j=0,(\Kbb p)_j=0,\;\vee\\
%   (\Kbb p)_j=c q_j (c\in\R), ~\langle \mu_j,q_j \rangle =0,\;\vee \\
%   \vartheta_j+c\alpha_j=0,(\Kbb p)_j = c q_j (c\ge 0),\langle \mu_j,q_j \rangle\le 0.
% \end{aligned} \right\}  &&&\text{if }j\in\Bcal,\\
% \vartheta_j + \frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|} \ge 0,\; \mu_j = 0,&&&\text{if }j\in\Ical_0,\\
% \left.
% \begin{aligned}
%   \vartheta_j=0,(\Kbb p)_j = 0,\;\vee \\
%   \vartheta_j\le\|(\Kbb p)_j\|,\mu_j=0,\;\vee \\
%   |\vartheta_j|\le \|(\Kbb p)_j\|,\scalar{\mu_j}{(\Kbb p)_j}\le 0
% \end{aligned} \right\} &&&\text{if }j\in\Tcal.
\end{align}

Referring to <ref> let us take $F_1(\alpha,u)=\phi'(u)\in\R^m$ and $F_2(\alpha,u) = (\alpha,u)\in\R^n_+\times\R^m$. Existence of the KKT multipliers is guaranteed if the following constraint qualification condition holds for $(\vartheta,\Kbb^\top \mu,p)\in N_{\gph Q}^M(\alpha^*,u^*,\Kbb^\top q^*)$

\begin{equation}\label{eq:outrata_optimality_system_proof_1}
\begin{bmatrix}
I & \bm{0} & \bm{0} \\
\bm{0} & I & -\phi''(u^*)^\top
\end{bmatrix} \begin{bmatrix}
\vartheta\\
\Kbb^\top \mu\\
\end{bmatrix}\in -N^M_{\R^n_+}(\alpha^*)\times\{0\} \text{ implies }\vartheta = 0,\;\Kbb^\top\mu = 0,\;p = 0.
\end{equation}

Recalling <ref> and using the expression of the Mordukhovich normal cone $N_{R^n_+}^M(\alpha^*)=N_{R^n_+}(\alpha^*)=\{v\in\R^2:\scalar{v}{\alpha^*} = 0,\;v\le 0\}$, condition <ref> can also be written as

\begin{align}
\Kbb^\top \mu - \phi''(u^*)^\top p &= 0,\label{eq:constraint_qualification_1}\\
\scalar{\alpha^*}{\vartheta} &= 0,\label{eq:constraint_qualification_2}\\
\vartheta &\ge 0.\label{eq:constraint_qualification_3}
\end{align}

Let us take $(\vartheta,\Kbb^\top\mu,p)\in N^M_{\gph Q}(\alpha^*,u^*,\Kbb^\top q^*)$ and let us multiply <ref> by $p$ on the left. Recalling $(\Kbb p)_j=0$ in $\Acal_s$ and $\mu_j=0$ in $\Ical_0$, we have for each remaining index set

\begin{align*}
\scalar{p}{\phi''(u^*)^\top p}&=\scalar{p}{\Kbb^\top\mu} = \sum_{j\in\Ical}\scalar{\mu_j}{(\Kbb p)_j} + \sum_{j\in\Bcal}\scalar{\mu_j}{(\Kbb p)_j} + \sum_{j\in\Bcal_0}\scalar{\mu_j}{(\Kbb p)_j},\\
&= \sum_{j\in\Ical}-\alpha_j\scalar{(\Kbb p)_j}{T_j(\Kbb p)_j} + \sum_{j\in\Bcal}c\underbrace{\scalar{\mu_j}{q_j}}_{\le 0} + \sum_{j\in\Bcal_0}\underbrace{\scalar{\mu_j}{(\Kbb p)_j}}_{\le 0} \le 0,
\end{align*}

where we used the characterization of the Mordukhovich normal cone. Furthermore, using the strong convexity of the function $\phi$ we have $\scalar{p}{\phi''(u^*)^\top p}\ge 0$. Both inequalities imply $p=0$ and, according to <ref>, it also yields $\Kbb^\top\mu = 0$.
Moreover, if we consider the index set $\Ical\cup\Acal_s\cup\Bcal$, we know in all these sets $\alpha^*_j>0$, and therefore, to satisfy equation <ref> it must hold $\vartheta_j=0$. Since $p=0$ in $\Tcal$ we know $\vartheta_j = 0$ or $\vartheta_j \le \|(\Kbb p)_j\|$ for this index set, in both cases it leads to $\vartheta_j=0$. In $\Ical_0$, we have $\vartheta_j \le - \frac{\scalar{(\Kbb u)_j}{(\Kbb p)_j}}{\|(\Kbb u)_j\|} = 0$ and <ref> yields $\vartheta_j=0$. Therefore, $\vartheta_j = 0$ for all $j$. Consequently, the existence of multipliers is guaranteed and there exists a vector $\rho\in N^M_{R^n_+}(\alpha^*)$ and KKT multipliers $(\vartheta,\Kbb^\top\mu,p)\in N_{\gph Q}^M(\alpha^*,u^*,\Kbb^\top q^*)$ such that

\begin{align}
0 &= \nabla_u J(u^*) + (\nabla_u F_2(\alpha,u))^\top \begin{bmatrix}
\vartheta\\
\Kbb^\top\mu
\end{bmatrix} - (\nabla_u F_1(\alpha,u))^\top p,\\
0 &= (\nabla_\alpha F_2(\alpha,u))^\top \begin{bmatrix}
\vartheta\\
\Kbb^\top\mu
\end{bmatrix}- (\nabla_\alpha F_1(\alpha,u))^\top p + \rho.
\end{align}

§ TRUST REGION ALGORITHM

In this section we will describe the numerical algorithm used for finding optimal parameters of <ref>. Thanks to the Bouligand subdifferential characterization given in <ref>, it is possible to compute a linearized representative of the directional derivative via the linear system <ref>. Using this information we can make use of a descent-like algorithm for its numerical solution. Indeed, by using the uniqueness properties of the solution operator, we can write a reduced optimization problem

where $u(\alpha)$ is the image reconstruction corresponding to a particular value of $\alpha$. With this reduced problem, we can make use of the stationarity condition for the bilevel problem described in <ref> and the directional derivative characterization <ref>. By using the definition of the directional derivative for the reduced optimization problem, we get

\begin{equation}\label{eq:generalized_gradient_rep}
\scalar{j'(\alpha)}{h} = \scalar{\nabla J(u)}{S'(\alpha;h)} = \scalar{\nabla J(u)}{\tilde{\eta}}.
\end{equation}

where $\tilde{\eta}$ is a solution of system <ref> for a particular partition of the biactive set $\Bcal = \Bcal_1\cup\Bcal_2$. Let us now define a generalized adjoint $p\in\R^m$ as the solution of the following system

\begin{align*}
\scalar{\phi''(u)^\top p}{v} + \sum_{j\in\Ical}\scalar{\mu_j}{(\Kbb v)_j} - \scalar{\nabla J(u)}{v} &= 0,&\forall v\in V,\\
\mu_j - \alpha_jT_j(\Kbb p)_j &= 0,&\forall j\in\Ical,
\end{align*}

where $V$ is defined as in <Ref>. Moreover, using the results in <Ref>, we know that $\tilde{\eta}\in V$ is a linear representative of the directional derivative. Consequently, <ref> reads

\begin{equation*}
\scalar{j'(\alpha)}{h} = \scalar{\nabla J(u)}{\tilde{\eta}} = \scalar{\phi''(u)^\top p}{\tilde{\eta}} + \sum_{j\in\Ical}\scalar{\alpha_j T_j(\Kbb p)_j}{(\Kbb\tilde{\eta})_j},
\end{equation*}

Rearranging the terms we get

\begin{equation*}
\scalar{j'(\alpha)}{h} = \scalar{p}{\phi''(u)\tilde{\eta}} + \sum_{j\in\Ical}\scalar{(\Kbb p)_j}{\alpha_j T_j(\Kbb\tilde{\eta})_j} = \scalar{p}{\phi''(u)\tilde{\eta}} + \sum_{j\in\Ical}\scalar{(\Kbb p)_j}{\tilde{\lambda}_j},
\end{equation*}

Finally, using <ref> we get

\begin{equation}
\scalar{j'(\alpha)}{h}= -\sum_{j\in\Ical\cup\Bcal_2}\frac{h_j}{\alpha_j}\scalar{q_j}{(\Kbb p)_j}\label{eq:b_subgrad_equation}
\end{equation}

In this work, we will rely on a nonsmooth trust region method in the spirit of <cit.>. In general, a trust-region algorithm works by defining for an iteration $k$ a radius $\Delta_k$ and replaces the function with a model function $m_k(\alpha_k)$ within a trust-region. This region in this paper will be represented by a $l_\infty$ ball, and the radius $\Delta_k$ is updated according to a quality measure based on the predicted decrease of the model function pred and the actual cost function reduction ared.

The model function we will be using has the following form

\begin{equation*}
m_k(\alpha_k) = j(\alpha_k) + g_k^\top d_k + \frac{1}{2}d_k^\top B_k d_k,
\end{equation*}

where $B_k$ is a second order matrix built using BFGS updates, and $g_k$ is obtained either by using the Bouligand subdifferential element obtained using <ref> or a gradient obtained from the regularized model, $g_{\gamma,k}$, if the radius falls below a threshold value, see <Ref>. The idea of differentiating two phases for the trust region subproblem approximation was taken from <cit.>, where this strategy is presented for optimizing nonsmooth, non-convex and locally-Lipschitz functions. By defining an appropriate model function, dependent on the radius value, convergence to a Clarke (C)-stationary points may be verified <cit.>. The complete steps are provided in <ref>.

Non-smooth Trust-Region Algorithm
[1]
Choose initial parameter $\alpha_0$, radius $\Delta_0$, $0<\eta_1\le\eta_2<1$, $0<\gamma_1\le\gamma_2<1$ and $tol>0$
Choose initial second order matrix $B_0$ and a threshold radius $\Delta_t$
$\Delta_k > tol$
$m_k(\alpha_k) = j(\alpha_k) + g_k^\top d_k + \frac{1}{2}d_k^\top B_kd_k$
$m_k(\alpha_k) = j(\alpha_k) + g_{\gamma,k}^\top d_k + \frac{1}{2}d_k^\top B_kd_k$
Compute a step $s_k$ that “sufficiently” reduces the model $m_k$ such that $\alpha_k+s_k\in B_{\Delta_k}$
Update second order matrix $B_k$ using limited memory BFGS.
Compute the quality measure $\rho_k$
$\alpha_{k+1} = \begin{cases}
\alpha_k &\text{if }\rho_k\le\eta_1,\\
\alpha_k+s_k &\text{otherwise}.
\end{cases}$
$\Delta_{k+1} = \begin{cases}
[\Delta_k,\infty)&\text{if }\rho_k\ge\eta_2,\\
[\gamma_2\Delta_k,\Delta_k]&\text{if }\rho_k\in(\eta_1,\eta_2),\\
[\gamma_1\Delta_k,\gamma_2\Delta_k]&\text{if }\rho_k<\eta_1.
\end{cases}$
$k \leftarrow k+1$

§.§ Trust Region Subproblem

Regarding step <ref> in <ref>, we need to approximate the solution of the following trust region sub-problem

s∈^nj(α)+g^⊤s + 1/2s^⊤Bs

which corresponds to a classical trust-region sub-problem with additional positivity constraints (see, e.g., <cit.>). The main idea is to reformulate the problem taking advantage of the $l_\infty$ norm used for the ball at the point $x_k$

s∈^nj(α)+g^⊤s + 1/2s^⊤Bs
max(-α_j,-Δ) ≤s_j ≤Δ, ∀j=1,…,n.

For performance purposes it is desirable to solve this problem approximately in such a way that we can guarantee a descent on the cost function. With that goal in mind we will make use of a dogleg strategy that takes into account a Newton step $s_N$ and a Cauchy step $s_C$. In the context of this constrained problem, let us take $\tilde{B}_\Delta = B_\Delta\cap\R^n_+$ and distinguish the following three cases

* $s_N\in \tilde{B}_\Delta$,

* $s_C\in \tilde{B}_\Delta$ and $s_N\notin \tilde{B}_\Delta$,

* $s_C\notin \tilde{B}_\Delta$ and $s_N\notin \tilde{B}_\Delta$.

For case 1  we take the Newton step; for case 2 a dogleg strategy for box constraints is used; for case 3 we make use of a scaled Cauchy direction as described in <ref>.

$s_N\in \tilde{B}_\Delta$

$s_C\in \tilde{B}_\Delta$, $s_N\notin \tilde{B}_\Delta$

$s_C\notin \tilde{B}_\Delta$, $s_N\notin \tilde{B}_\Delta$

Projected Cauchy Step:
The figure depicts a two dimensional example for the three different possible cases when approximating the trust region subproblem using a dogleg strategy with $l_\infty$ norm and positivity constraints.

Dogleg Step for Box Constraints
[1]
Calculate Newton step by solving the linear system $B_k s_N=-g_k$.
$s_N \in \tilde{B}_\Delta$
Calculate $s_C = -\frac{\|g_k\|^2}{g_k^\top*B_k*g_k}g_k$
$s_C \in \tilde{B}_\Delta$
Calculate the intersection of $s_N-s_C$ with $\tilde{B}_\Delta$ and obtain $s_{DL}$
Find $t$ such that $t*s_C/\|s_C\|$ remains in $\tilde{B}_\Delta$.

§.§ Local Regularization

As a safeguard, for the cases where the trust-region radius fall below a threshold value $\Delta_t$ we will rely on a switching strategy between the nonsmooth model presented earlier and a regularized model. Let us consider the following optimization problem

u̅(α) = _u ϕ(u) + ∑_j=1^n α_j(u)_j_γ

where $\|\cdot\|_\gamma$ is a local regularization of the euclidean norm given by

\begin{equation*}
\|z\|_\gamma = \begin{cases}
\|z\| - \frac{1}{2\gamma} &\text{if }\|z\| \ge \frac{1}{\gamma},\\
\frac{\gamma}{2}\|z\|^2 &\text{if } \|z\| < \frac{1}{\gamma}.
\end{cases}
\end{equation*}

Since the problem is now differentiable, we can define the following optimality condition for the lower level problem is obtained

\begin{equation*}
\scalar{\phi'(u)}{v} + \sum_{j=1}^n\alpha_j\scalar{h_\gamma ((\Kbb u)_j)}{(\Kbb v)_j} = 0,\;\forall v\in\R^m.
\end{equation*}

Here, $h_\gamma$ correspond to the first derivative of the regularized euclidean norm. Moreover, regarding the bilevel problem, it is possible to formulate the KKT conditions for this problem. Let $(\alpha^*,u^*)\in\R^n_+\times\R^m$ be a stationary point for the regularized bilevel problem <ref>. Then there exist Lagrange multipliers $(p,\sigma)\in\R^m\times\R^n$ such that the following optimality system holds true

\begin{align}
\scalar{\phi'(u)}{v} + \sum_{j=1}^n\alpha_j^*\scalar{h_\gamma ((\Kbb u^*)_j)}{(\Kbb v)_j} &= 0,&\;\forall v\in\R^m,\\
\scalar{\phi''(u^*)^\top p}{v} + \sum_{j=1}^n\alpha_j\scalar{h'^*_\gamma ((\Kbb u^*)_j)(\Kbb p)_j}{(\Kbb v)_j} &= -\scalar{J'(u)}{v},&\;\forall v\in\R^m,\label{eq:regularized_adjoint}\\
\scalar{h_\gamma((\Kbb u^*)_j)}{(\Kbb p)_j} &= \sigma_j,&\;\forall j=1,\dots,n,\\
0\le \sigma \perp \alpha^* \ge 0.
\end{align}

Moreover, with help of the adjoint equation <ref> it is possible to derive a gradient formula for the reduced cost function $j(\alpha) = J(u_\gamma(\alpha))$ as

\begin{equation*}
(j'(\alpha))_j = \scalar{h_\gamma((\Kbb u)_j)}{(\Kbb p)_j}.
\end{equation*}

According to <ref>, if the value for the radius runs below a predetermined threshold, the algorithm will make use of this gradient for its calculations.

§ EXPERIMENTS

In this section we report on the performance of <ref> presented in <Ref> to find optimal parameters. With this goal in mind we prepared two training image datasets.

§.§ Single Training Pair

The first dataset we will explore is a single 128 by 128 pixel image pair dataset based on the cameraman image and a corrupted version, obtained by adding gaussian noise with zero mean and standard deviation $\sigma=0.05$. <Ref> shows this training image pair along with the optimal parameter obtained using the trust-region algorithm for both a scalar and a two-dimensional patch parameter. The improvement when using a two-dimensional patch parameter, according to the SSIM value of the image reconstructions, in <Ref>. Moreover, the cost function when using a scalar and a two dimensional parameter is shown, respectively. These plots show the non-convexity of the cost function, and the cost corresponding to the optimal parameter.

(image) at (0,0)

[green,very thick] (-1.35,0.0) rectangle (1.35,1.35) node[below left,black,fill=green]$\alpha_1$;
[green,very thick] (-1.35,-1.35) rectangle (1.35,0.0) node[below left,black,fill=green]$\alpha_2$;

$\alpha^*=[0.0233; 0.0126]$

Cameraman Dataset

Optimal reconstructions using a scalar regularization parameter and a 2 dimensional regularization parameter.

Cost Function - Cameraman Dataset

Values for the $l_2$ squared cost function using a scalar regularization parameter and a two dimensional regularization parameter using the Cameraman single image dataset.

In order to explore the behavior of the algorithm with more degrees of freedom we will make use of a so called patch-based parameter. This parameter will be a piecewise constant parameter. We split the image size into a grid of different size according to the number of patches, i.e., the patch size will be smaller as we require more patches within the image domain. For instance a 4x4 patch will calculate an optimal parameter of size 16.

Regarding the behavior of the algorithm, <Ref> shows the performance on the cameraman training dataset when using different number of patches, along with the number of iterations, cost, quality measures and residue. Indeed, the improvement on the reconstructed images when using more patches can be verified from the data. Furthermore, <Ref> shows the optimal reconstructions along with the optimal parameter for different number of patches. In this experiment it can be seen how the regularization parameter can adjust different regularization levels according to the training pair.

2cItem   2c   2cReconstruction
Patch     Iterations     $\|\alpha_{k+1}-\alpha_k\|$     COST     SSIM     PSNR
Scalar     19     7.629e-7     34.1206     0.7712     26.3794

2x1     20     4.883e-5     33.6915      0.8007     25.6222

2x2     19     4.639e-5     33.6625      0.8034     25.8982

4x4     22     1.461e-4     33.0175     0.8189     25.8303

8x8     28     5.614e-5     32.2040     0.8283     24.4066

16x16     36     9.119e-4     31.6976     0.8455     25.8264

Trust Region Algorithm behavior on the Cameraman dataset.

Optimal Patch Parameter - Cameraman Dataset

Values for the $l_2$ squared cost function using a scalar regularization parameter and a two dimensional regularization parameter using the Cameraman dataset.

§.§ Multiple Training Pairs

For the second experiment, we used ten image pairs containing images of faces to generate a training dataset and ten different image pairs to generate a testing dataset; both datasets were based on the CelebA dataset <cit.>. These images are of size 128 by 128 pixels and in both datasets, the degenerated pairs were generated by adding gaussian noise with zero-mean and standard deviation $\sigma = 0.1$. A subset of the training dataset is depicted in <Ref>.

Faces Dataset

A subset of the CelebA dataset corrupted with gaussian noise.

In <ref> we plot the cost function corresponding to a scalar parameter and two-dimensional patch parameter along with the cost function corresponding to the optimal value calculated by the algorithm.
It is worth mentioning that when considering a patch-dependent parameter, as it was the case with the single image dataset, a scale dependent structure appears (see <ref>).

Cost Function - Faces Dataset

Values for the $l_2$ squared cost function using a scalar regularization parameter and a two dimensional regularization parameter using the Faces dataset.

Optimal Patch Parameter - Faces Dataset

Values for the optimal parameters calculated for different parameter patch sizes.

For the training dataset proposed, the quality reconstruction results for different number of patches is shown in <Ref>. This table shows the mean values for the SSIM and PSNR quality metric for the reconstructed images from the training dataset. Again, an improvement on the reconstruction quality can be seen as the degrees of freedom for the regularization parameter increases.

2cItem   2c   2cReconstruction
Patch     Iterations     $\|\alpha_{k+1}-\alpha_k\|$     COST     MSSIM     MPSNR
Scalar     14     1.562e-3     152.3726     0.8141     27.0689

2x1     17     1.049e-3     152.3209     0.8142     27.0589

2x2     21     1.762e-4     152.3120      0.8146     27.0739

4x4     21     7.421e-4     152.0911     0.8154     27.0766

8x8     24     7.050e-4     151.2011     0.8173     27.0787

16x16     25     1.410e-3     149.9673     0.8193     27.1047

32x32     35     3.285e-4     147.6573     0.8223     27.2141

Trust Region Algorithm behavior on the Faces dataset.

Finally, we can estimate the denoiser performance in images from the testing dataset. This experiment will show the overfitting phenomena that may occur when dealing with large number of patches, including the case of scale-dependent parameters ($\alpha\in\R^n$). Indeed, it can be seen in the testing dataset an increment on the mean SSIM (MSSIM) for the reconstructed images from the testing dataset up to a 8x8 patch size. Any higher number of patches results in quality degradation. This is indeed the expected behavior when dealing with overfitting problems.

img num     noisy     scalar     2x2     4x4     8x8     16x16     32x32
1     0.5247     0.7951     0.7948     0.7937     0.7956     0.7940     0.7930

2     0.4588     0.6614     0.6618      0.6616     0.6625     0.6625     0.6592

3     0.4267     0.7719     0.7721     0.7723     0.7747     0.7743     0.7710

4     0.3836     0.7237     0.7229     0.7214     0.7235     0.7233     0.7202

5     0.4580     0.7812     0.7799      0.7783     0.7798     0.7799     0.7768

6     0.4263     0.7400     0.7403     0.7403     0.7420     0.7428     0.7423

7     0.4547     0.6321     0.6322     0.6308     0.6305     0.6298     0.6262

8     0.4117     0.7381     0.7386      0.7405     0.7416     0.7413     0.7398

9     0.4655     0.7430     0.7409     0.7399     0.7405     0.7388     0.7354

10     0.5081     0.8210     0.8210     0.8204     0.8211     0.8197    0.8181
MSSIM         $\bm{0.7408}$     $\bm{0.7405}$     $\bm{0.7400}$     $\bm{0.7412}$     $\bm{0.7406}$     $\bm{0.7382}$

Faces Dataset SSIM Quality Measures in the validation dataset.

§.§ Learning Optimal Total Variation Discretization

The selection of an adequate discretization of the total variation seminorm in the context of image reconstruction problems is still an open problem <cit.>. In <cit.>, the authors propose a methodology for finding optimal discretization where instead of using hand-crafted discretization schemes for the total variation seminorm, a learning strategy is proposed.

The bilevel framework presented in this work, can also be used to learn optimal gradient discretization, by using different discretization schemes and their corresponding regularization parameters. We will make use of a training dataset to estimate the optimal regularization parameters for the contributions of each discretization scheme into the final solution. Let us define the following variational denoising model

\begin{equation}\label{eq:sumregs_lower_level}
\min_{u\in\R^n} \Ecal(u) := \frac{1}{2}\|u-f\|^2+\sum_{j=1}^n(\alpha_1)_j\|(\Kbb_1 u)_j\|+\sum_{j=1}^n(\alpha_2)_j\|(\Kbb_2 u)_j\|+\sum_{j=1}^n(\alpha_3)_j\|(\Kbb_3 u)_j\|,
\end{equation}

where $\Kbb_1,\Kbb_2$ and $\Kbb_3$ are the forward, backward and centered finite differences discretization of the gradient operator respectively. The goal is to determine optimal parameters $(\alpha_1,\alpha_2,\alpha_3)^\top \in\R^{3\times n}_+$ that lead to an optimal discretization of the total variation operator. Indeed, we consider the following bilevel learning strategy

u̅ = _u∈^m(u),

By extending the model presented previously, we can make use of a similar analysis using the Bouligand candidate presented in <ref> by defining the following adjoint state

\begin{align*}
\scalar{p}{v} + \sum_{i=1}^3\sum_{j\in\Ical^i}\scalar{\mu_j^i}{(\Kbb^i v)_j} - \scalar{\nabla J(u)}{v} &= 0,\forall v\in \mathcal{V}\\
\mu^i_j - \alpha_j^i T_j^i(\Kbb^i p)_j &= 0,\;\forall j\in\Ical^i,\;i=1,2,3,\\
\end{align*}

where $\mathcal{V}:=\{v\in\R^m:(\Kbb^i v)_j=0,\;\forall \Acal_s^i\cup\Bcal_1^i,\;(\Kbb^i v)_j\in span(q_j^i),\;\forall j\in\Bcal_2^i,\;i=1,2,3\}$, and the following gradient form

\begin{equation*}
\scalar{j'(\alpha^i)}{h^i} = -\sum_{j\in\Ical^i\cup\Bcal_2^i}\frac{h_j^i}{\alpha_j^i}\scalar{q_j^i}{(\Kbb^i p)_j}\;\text{for }i=1,2,3.
\end{equation*}

The same procedure applies for the regularized problem. Using the KKT conditions for the differentiable problem we obtain

\begin{equation*}
\scalar{p}{v} + \sum_{i=1}^3\sum_{j=1}^n \alpha_j^i\scalar{h'^*_\gamma ((\Kbb^i u)_j)(\Kbb^i p)_j}{(\Kbb^i v)_j}  = -\scalar{\nabla J(u)}{v},\;\forall v\in\R^m,\\
\end{equation*}

and its corresponding gradient characterization

\begin{equation*}
(j'(\alpha^i))_j = \scalar{h_\gamma((\Kbb^i u)_j)}{(\Kbb^i p)_j},\;\forall i=1,2,3.
\end{equation*}

The optimal regularization parameters for 4x4 and 8x8 patch-parameters for each of the three regularizers considered is presented in <ref> for the cameraman dataset and <ref> for the faces dataset. Again, we can see an improvement on the image quality by using more patches in the training datasets. Finally, in <ref> we can see a comparison between different patch sizes when using these three regularizers. Again, it is worth noticing that in the validation dataset the best reconstruction quality corresponds to an 8x8 patch-parameter size, and any further enlargement of the patch dimension drives a lower reconstruction quality.

Optimal Patch Parameters - Cameraman Dataset

Values for the optimal parameter calculated for the Cameraman dataset for different patch sizes.

Optimal Patch Parameters - Faces Dataset

Values for the optimal parameters calculated for different parameter patch sizes.

img num     noisy     scalar     2x2     4x4     8x8     16x16

1     0.5247     0.7922     0.7952     0.7951     0.7977     0.7965

2     0.4588     0.6568     0.6676     0.6663     0.6649     0.6653

3     0.4267     0.7745     0.7742     0.7729     0.7748     0.7741

4     0.3836     0.7237     0.7239     0.7226     0.7244     0.7206

5     0.4580     0.7814     0.7796     0.7779     0.7806     0.7809

6     0.4263     0.7394     0.7406     0.7414     0.7432     0.7441

7     0.4547     0.6262     0.6372     0.6357     0.6338     0.6333

8     0.4117     0.7377     0.7380     0.7413     0.7429     0.7421

9     0.4655     0.7438     0.7481     0.7445     0.7432     0.7420

10     0.5081     0.8185     0.8225     0.8217     0.8219     0.8211

MSSIM         $\bm{0.7395}$     $\bm{0.7420}$     $\bm{0.7419}$     $\bm{0.7428}$     $\bm{0.7424}$

Faces Dataset SSIM Quality Measures - Optimal Gradient Discretization in the validation dataset

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

titleError estimates for finite differences approximations of the total variation

booktitleIFIP Conference on System Modeling and Optimization
titleDynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints

Walter de Gruyter GmbH

booktitleVariational Methods
titleBilevel approaches for learning of variational imaging models

journaltitleJournal of Mathematical Imaging and Vision
titleCrouzeix–Raviart approximation of the total variation on simplicial meshes

titleLearning Consistent Discretizations of the Total Variation

journaltitleSIAM Journal on Optimization
titleA Nonsmooth Trust-Region Method for Locally Lipschitz Functions with Application to Optimization Problems Constrained by Variational Inequalities

journaltitleJournal of Mathematical Imaging and Vision
titleBilevel Parameter Learning for Nonlocal Image Denoising Models

titleNumerical PDE-Constrained Optimization

journaltitleJournal of Mathematical Imaging and Vision
titleBilevel parameter learning for higher-order total variation regularisation models

journaltitleJournal of Mathematical Analysis and Applications
titleThe structure of optimal parameters for image restoration problems

journaltitleInverse Problems & Imaging
titleImage Denoising: Learning the noise model via nonsmooth pde-constrained optimization.

journaltitleHandbook of Mathematical Models and Algorithms in Computer Vision and Imaging
titleBilevel Optimization Methods in Imaging

journaltitleJournal of Mathematical Imaging and Vision
titleAutomated regularization parameter selection in multi-scale total variation models for image restoration

titleConvex analysis and variational problems

journaltitleJournal of Applied Mechanics
titleNumerical methods for nonlinear variational problems

journaltitleJournal of Mathematical Analysis and Applications
titleAnalytical aspects of spatially adapted total variation regularisation

booktitleHandbook of Numerical Analysis
titleGenerating structured nonsmooth priors and associated primal-dual methods

Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät

titleBilevel optimization for calibrating point spread functions in blind deconvolution

Springer Nature

titleIntroduction to the theory of nonlinear optimization

booktitleComputer Vision Winter Workshop (CVWW)
titleContinuous hyper-parameter learning for support vector machines

journaltitleSIAM Journal on Imaging Sciences
titleA bilevel optimization approach for parameter learning in variational models

booktitleProceedings of International Conference on Computer Vision (ICCV)
titleDeep Learning Face Attributes in the Wild

Cambridge University Press

titleMathematical programs with equilibrium constraints

journaltitleJournal of Mathematical Imaging and Vision
titleTechniques for gradient-based bilevel optimization with non-smooth lower level problems

journaltitleSIAM Journal on Control and Optimization
titleA generalized mathematical program with equilibrium constraints

booktitleGerman Conference on Pattern Recognition
titleA deep variational model for image segmentation

Springer Science & Business Media

titleVariational analysis

Springer Science & Business Media

titleIntroduction to piecewise differentiable equations

IOP Publishing

journaltitleInverse problems
titleEdge-preserving and scale-dependent properties of total variation regularization

booktitleDiusion in Image Processing, UCLA Math Department CAM Report
titleSpatially and scale adaptive total variation based regularization and anisotropic diffusion in image processing

IOP Publishing

journaltitleInverse Problems
titleLearning optimal spatially-dependent regularization parameters in total variation image denoising

booktitleWSEAS International Conference on Applied Mathematics
titleA rectangular trust region dogleg approach for unconstrained and bound constrained nonlinear optimization

journaltitleOptimization online, July
titleASTRAL: An active set linfinity-trust-region algorithm for box constrained optimization

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

titleError estimates for finite differences approximations of the total variation

booktitleIFIP Conference on System Modeling and Optimization
titleDynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints

Walter de Gruyter GmbH

booktitleVariational Methods
titleBilevel approaches for learning of variational imaging models

journaltitleJournal of Mathematical Imaging and Vision
titleCrouzeix–Raviart approximation of the total variation on simplicial meshes

titleLearning Consistent Discretizations of the Total Variation

journaltitleSIAM Journal on Optimization
titleA Nonsmooth Trust-Region Method for Locally Lipschitz Functions with Application to Optimization Problems Constrained by Variational Inequalities

journaltitleJournal of Mathematical Imaging and Vision
titleBilevel Parameter Learning for Nonlocal Image Denoising Models

titleNumerical PDE-Constrained Optimization

journaltitleJournal of Mathematical Imaging and Vision
titleBilevel parameter learning for higher-order total variation regularisation models

journaltitleJournal of Mathematical Analysis and Applications
titleThe structure of optimal parameters for image restoration problems

journaltitleInverse Problems & Imaging
titleImage Denoising: Learning the noise model via nonsmooth pde-constrained optimization.

journaltitleHandbook of Mathematical Models and Algorithms in Computer Vision and Imaging
titleBilevel Optimization Methods in Imaging

journaltitleJournal of Mathematical Imaging and Vision
titleAutomated regularization parameter selection in multi-scale total variation models for image restoration

titleConvex analysis and variational problems

journaltitleJournal of Applied Mechanics
titleNumerical methods for nonlinear variational problems

journaltitleJournal of Mathematical Analysis and Applications
titleAnalytical aspects of spatially adapted total variation regularisation

booktitleHandbook of Numerical Analysis
titleGenerating structured nonsmooth priors and associated primal-dual methods

Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät

titleBilevel optimization for calibrating point spread functions in blind deconvolution

Springer Nature

titleIntroduction to the theory of nonlinear optimization

booktitleComputer Vision Winter Workshop (CVWW)
titleContinuous hyper-parameter learning for support vector machines

journaltitleSIAM Journal on Imaging Sciences
titleA bilevel optimization approach for parameter learning in variational models

booktitleProceedings of International Conference on Computer Vision (ICCV)
titleDeep Learning Face Attributes in the Wild

Cambridge University Press

titleMathematical programs with equilibrium constraints

journaltitleJournal of Mathematical Imaging and Vision
titleTechniques for gradient-based bilevel optimization with non-smooth lower level problems

journaltitleSIAM Journal on Control and Optimization
titleA generalized mathematical program with equilibrium constraints

booktitleGerman Conference on Pattern Recognition
titleA deep variational model for image segmentation

Springer Science & Business Media

titleVariational analysis

Springer Science & Business Media

titleIntroduction to piecewise differentiable equations

IOP Publishing

journaltitleInverse problems
titleEdge-preserving and scale-dependent properties of total variation regularization

booktitleDiusion in Image Processing, UCLA Math Department CAM Report
titleSpatially and scale adaptive total variation based regularization and anisotropic diffusion in image processing

IOP Publishing

journaltitleInverse Problems
titleLearning optimal spatially-dependent regularization parameters in total variation image denoising

booktitleWSEAS International Conference on Applied Mathematics
titleA rectangular trust region dogleg approach for unconstrained and bound constrained nonlinear optimization

journaltitleOptimization online, July
titleASTRAL: An active set linfinity-trust-region algorithm for box constrained optimization