$\epsilon$-Subgradient Algorithms for Bilevel Convex Optimization

TL;DR

This paper develops and analyzes a new class of explicit $psilon$-subgradient algorithms for bilevel convex optimization, demonstrating their effectiveness in large-scale practical problems like tomographic image reconstruction.

Contribution

It introduces a novel $psilon$-subgradient method with proven convergence for bilevel convex problems, applicable to large-scale applications.

Findings

Effective in solving large-scale problems

Convergence properties established

Applicable to tomographic image reconstruction

Abstract

This paper introduces and studies the convergence properties of a new class of explicit $\epsilon$-subgradient methods for the task of minimizing a convex function over the set of minimizers of another convex minimization problem. The general algorithm specializes to some important cases, such as first-order methods applied to a varying objective function, which have computationally cheap iterations. We present numerical experimentation regarding certain applications where the theoretical framework encompasses efficient algorithmic techniques, enabling the use of the resulting methods to solve very large practical problems arising in tomographic image reconstruction.

Equations

$$\begin{array}{rl}\min & f_1(\mathbf{x})\\ s.t.& \mathbf{x}\in \underset{\mathbf{y}\in X_0}{\mathrm{argmin}}{f}_0(\mathbf{y})\text{,}\end{array}$$

$$-\frac{{\mathbf{x}}_{k+1}-{\mathbf{x}}_k}{{\lambda }_k}\in {\partial }_{{ϵ}_k}(f_0+{\eta }_k{f}_1)({\mathbf{x}}_{k+1})\text{,}$$

$${\partial }_{ϵ}f(\mathbf{x}):=\{\mathbf{v}:f(\mathbf{y})\ge f(\mathbf{x})+{\mathbf{v}}^T(\mathbf{y}-\mathbf{x})-ϵ,\forall \mathbf{y}\in {\mathbb{R}}^n\},$$

$${\mathbf{x}}_{k+1}\in {ϵ}_k\text{mbox}-\underset{\mathbf{x}\in X_0}{\mathrm{argmin}}\left\{\frac{1}{2{\lambda }_k}\Vert \mathbf{x}-{\mathbf{x}}_k{\Vert }^2+f_0(\mathbf{x})+{\eta }_k{f}_1(\mathbf{x})\right\}\text{.}$$

$$\dot{\boldsymbol{x}}(t)+\nabla f_{0}\bigl{(}\boldsymbol{x}(t)\bigr{)}+\eta(t)\nabla f_{1}\bigl{(}\boldsymbol{x}(t)\bigr{)}=0.$$

$$-\frac{{\mathbf{x}}_{k+1}-{\mathbf{x}}_k}{{\lambda }_k}\in {\partial }_{{ϵ}_k}(f_0+{\eta }_k{f}_1)({\mathbf{x}}_k)\text{.}$$

$$\begin{array}{rl}-\frac{{\mathbf{x}}_{k+1/2}-{\mathbf{x}}_k}{{\lambda }_k}& \in {\partial }_{{ϵ}_k^0}{f}_0({\mathbf{x}}_k)\\ -\frac{{\mathbf{x}}_{k+1}-{\mathbf{x}}_{k+1/2}}{{\lambda }_k{\eta }_k}& \in {\partial }_{{ϵ}_k^1}{f}_1({\mathbf{x}}_{k+1/2}).\end{array}$$

$$\boldsymbol{x}_{k+1}:=\boldsymbol{x}_{k}-\lambda_{k}\bigl{(}\tilde{\nabla}f(\boldsymbol{x}_{k})+\alpha(\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1})\bigr{)}.$$

$$\sum _{k=0}^{\infty }{\lambda }_k=\infty \text{and}{\lambda }_k\to 0^{+}\text{.}$$

$$0<\underline{\lambda }\le {\lambda }_k\le \overline{\lambda },$$

$$\sum _{k=0}^{\infty }{\eta }_k=\infty \text{and}{\eta }_k\to 0^{+},$$

$$\sum _{k=0}^{\infty }{\lambda }_k{\eta }_k=\infty \text{and}{\eta }_k\to 0^{+}\text{.}$$

$$\begin{array}{rl}{\mathbf{x}}_{k+1/3}& :={\mathcal{O}}_{f_0}({\lambda }_k,{\mathbf{x}}_k)\text{;}\\ {\mathbf{x}}_{k+2/3}& :={\mathcal{O}}_{f_1}({\mu }_k,{\mathbf{x}}_{k+1/3})\text{;}\\ {\mathbf{x}}_{k+1}& :={\mathcal{P}}_{X_0}({\mathbf{x}}_{k+2/3})\text{,}\end{array}$$

$${\mathcal{P}}_X(\mathbf{x}):=\underset{\mathbf{y}\in X}{\mathrm{argmin}}\Vert \mathbf{x}-\mathbf{y}\Vert .$$

$$\sum _{k=0}^{\infty }{\mu }_k=\infty \text{and}\frac{{\mu }_k}{{\lambda }_k}\to 0^{+}\text{.}$$

$$\|\mathcal{O}_{f_{i}}(\lambda,\boldsymbol{x}_{k+i/3})-\boldsymbol{y}\|^{2}\leq\|\boldsymbol{x}_{k+i/3}-\boldsymbol{y}\|^{2}-\beta\lambda\bigl{(}f_{i}\bigl{(}\mathcal{O}_{f_{i}}(\lambda,\boldsymbol{x}_{k+i/3})\bigr{)}-f_{i}(\boldsymbol{y})\bigr{)}\\ {}+\lambda\rho_{i}(\lambda,k)\text{,}$$

$$\Vert \mathbf{x}-{\mathcal{O}}_{f_i}(\lambda ,\mathbf{x}){\Vert }_2\le \lambda \gamma \text{.}$$

$${\mathcal{S}}_f(\lambda ,\mathbf{x}):=\mathbf{x}-\lambda {\tilde{\nabla }}_{ϵ}f(\mathbf{x})\text{,}$$

$$\|\mathcal{S}_{f}(\lambda,\boldsymbol{x})-\boldsymbol{y}\|_{2}^{2}\leq\|\boldsymbol{x}-\boldsymbol{y}\|_{2}^{2}-2\lambda\bigl{(}f(\boldsymbol{x})-f(\boldsymbol{y})\bigr{)}+\lambda\bigl{(}\lambda\|\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})\|_{2}^{2}+2\epsilon\bigr{)}\text{.}$$

$$f(\mathbf{x})\ge f(\mathbf{z})+\tilde{\nabla }f(\mathbf{z}{)}^T(\mathbf{x}-\mathbf{z})=f(\mathbf{z})+\lambda \tilde{\nabla }f(\mathbf{z}{)}^T{\tilde{\nabla }}_{ϵ}f(\mathbf{x}),$$

$$-2\lambda f(\boldsymbol{x})\leq-2\lambda\Bigl{(}f\bigl{(}\mathcal{S}_{f}(\lambda,\boldsymbol{x})\bigr{)}+\lambda\tilde{\nabla}f\bigl{(}\mathcal{S}_{f}(\lambda,\boldsymbol{x})\bigr{)}^{T}\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})\Bigr{)}$$

$$\|\mathcal{S}_{f}(\lambda,\boldsymbol{x})-\boldsymbol{y}\|_{2}^{2}\leq\|\boldsymbol{x}-\boldsymbol{y}\|_{2}^{2}-2\lambda\bigl{(}f(\mathcal{S}_{f}(\lambda,\boldsymbol{x}))-f(\boldsymbol{y})\bigr{)}\\ {}+\lambda\bigl{(}\lambda\|\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})\|_{2}^{2}+2\epsilon-2\lambda\tilde{\nabla}f\bigl{(}\mathcal{S}_{f}(\lambda,\boldsymbol{x})\bigr{)}^{T}\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})\bigr{)}\text{,}$$

$$\|\mathcal{O}_{f}(\lambda,\boldsymbol{x})-\boldsymbol{y}\|_{2}^{2}\leq\|\boldsymbol{x}-\boldsymbol{y}\|_{2}^{2}-2\lambda\bigl{(}f(\boldsymbol{x})-f(\boldsymbol{y})\bigr{)}+\lambda\varrho(\lambda)\text{,}$$

$$\Vert {\mathcal{O}}_f(\lambda ,\mathbf{x})-\mathbf{x}\Vert \le \lambda \gamma ,$$

$$\|\mathcal{O}_{f}(\lambda,\boldsymbol{x})-\boldsymbol{y}\|_{2}^{2}\leq\|\boldsymbol{x}-\boldsymbol{y}\|_{2}^{2}-2\lambda\bigl{(}f(O_{f}(\lambda,\boldsymbol{x}))-f(\boldsymbol{y})\bigr{)}\\ {}+\lambda\Bigl{(}\varrho(\lambda)+2\lambda\gamma\bigl{\|}\tilde{\nabla}f\bigl{(}\mathcal{O}_{f}(\lambda,\boldsymbol{x})\bigr{)}\bigr{\|}\Bigr{)}\text{,}$$

$$\Vert P_{X_0}(\mathbf{x})-\mathbf{y}\Vert \le \Vert \mathbf{x}-\mathbf{y}\Vert \text{.}$$

$$\underset{k\to \infty }{\lim }{d}_{X_1}({\mathbf{x}}_k)=0\text{.}$$

$$\begin{split}\|\boldsymbol{x}_{k+1}-\boldsymbol{y}\|^{2}&{}\leq\|\boldsymbol{x}_{k+2/3}-\boldsymbol{y}\|^{2}\\ &{}\leq\|\boldsymbol{x}_{k+1/3}-\boldsymbol{y}\|^{2}-\beta\mu_{k}\bigl{(}f_{1}(\boldsymbol{x}_{k+2/3})-f_{1}(\boldsymbol{y})\bigr{)}+\mu_{k}\rho_{1}(\mu_{k},k)\\ &{}\leq\|\boldsymbol{x}_{k}-\boldsymbol{y}\|^{2}-\beta\lambda_{k}\bigl{(}f_{0}(\boldsymbol{x}_{k+1/3})-f_{0}^{*}\bigr{)}+\lambda_{k}\rho_{0}(\lambda_{k},k)\\ &\qquad\qquad\qquad\qquad{}-\beta\mu_{k}\bigl{(}f_{1}(\boldsymbol{x}_{k+2/3})-f_{1}(\boldsymbol{y})\bigr{)}+\mu_{k}\rho_{1}(\mu_{k},k)\text{.}\end{split}$$

$$\|\boldsymbol{x}_{k+1}-\boldsymbol{y}\|^{2}\leq\|\boldsymbol{x}_{k}-\boldsymbol{y}\|^{2}-\beta\lambda_{k}\bigl{(}f_{0}(\boldsymbol{x}_{k+1/3})-f_{0}^{*}\bigr{)}\\ {}+\lambda_{k}\rho_{0}(\lambda_{k},k)+\mu_{k}\bigl{(}\beta N+\rho_{1}(\mu_{k},k)\bigr{)}\text{.}$$

$$f_0({\mathbf{x}}_{k+1/3})\ge f({\mathbf{x}}_k)-M\Vert {\mathbf{x}}_{k+1/3}-{\mathbf{x}}_k\Vert \text{,}$$

Full text

We present numerical experimentation regarding certain applications where the theoretical framework encompasses efficient algorithmic techniques, enabling the use of the resulting methods to solve very large practical problems arising in tomographic image reconstruction.

1 Introduction

Our aim in the present paper is to solve a bilevel or hierarchical optimization problem of the form

[TABLE]

where $f_{i}:\mathbb{R}^{n}\to\mathbb{R}$ ($i\in\{0,1\}$) are convex functions and $X_{0}$ is a nonempty closed convex set.

Bilevel problems like (1) have already been considered in the literature. For example, for the case $X_{0}=\mathbb{R}^{n}$, Cabot Cabot (2005) suggests the use of the following algorithm:

[TABLE]

where $\partial_{\epsilon}f(\boldsymbol{x})$ is the $\epsilon$-subdifferential of $f$ at $\boldsymbol{x}$:

[TABLE]

$\eta_{k}\to 0^{+}$ and $\lambda_{k}$ is a nonnegative stepsize. Such iterations are reminiscent of approximate proximal methods, in the sense that (2) is equivalent to (we follow the notation of Cabot (2005)):

[TABLE]

While method (2) is powerful and conceptually simple, its application may be complicated by the implicit formulation. Assuming differentiability, iteration (2) can also be interpreted as a discretization of the continuous dynamical system

[TABLE]

Another way of discretizing this system is to do it explicitly, that is, to use iterations similar to

[TABLE]

or in two steps:

[TABLE]

Among the consequences of the results we will present in this paper, there are sufficient conditions on the sequences $\{\lambda_{k}\}$, $\{\epsilon_{k}\}$ and $\{\eta_{k}\}$ for the convergence of iterations (5) to the solution of problem (1). In fact, convergence of iterations (4) could also be proven using our abstract results, but would require more restrictive conditions on $\epsilon_{k}$ and we will keep this topic off the present paper.

Despite the fact that algorithms (2) and (5) are formally very similar, they differ significantly in both practical and theoretical aspects. At the practical side, implementation of an explicit iteration like (5) requires little more than evaluation of suitable $\epsilon$-subgradients. On the other hand, while strong convexity makes this minimization perhaps be more computationally amenable than approximately optimizing $f_{0}+\eta_{k}f_{1}$, computing (3) is still a nontrivial task. Furthermore, even if it allows for some tolerance $\epsilon_{k}$ in the optimization subproblem, current theory requires $\sum_{k=0}^{\infty}\epsilon_{k}<\infty$ to ensure convergence of (2), which means that this tolerance will decrease quickly.

Arguably, an algorithm such as (5) would be among the easiest to implement methodology for solving (1) in the general case. For example, Solodov Solodov (2007, 2008) has also provided algorithms for such bilevel problems: in Solodov (2007), where only the differentiable case is considered, the proposed algorithm has the form (4) with $\epsilon_{k}\equiv 0$ and $\lambda_{k}$ selected through a line search with a sufficient decrease criterion based on the value of $f_{0}+\eta_{k}f_{1}$. While this procedure is still simple to implement on the differentiable case, descent directions are harder to be found in the presence of nondifferentiability and in Solodov (2008) a bundle technique is used to this end. This approach requires sophisticated quadratic solvers, which complicates implementations. Another, partial and perturbed, descent method was developed by Helou and De Pierro Helou Neto and De Pierro (2011), where only sufficient decrease of $f_{0}$, assumed to be smooth, is enforced, alleviating the need of a descent direction for $f_{1}$. But this technique may still require multiple evaluations of $f_{0}$, while an iterative process like (5), differently, does not require any $f_{i}$ value because it does not rely on descent criteria. Furthermore, it should be remarked that the fact that (5) allows for inexactness in the computation has positive impact in algorithmic performance, which we will illustrate through experimental work. An approach currently available in the literature which is similar to a special case of the techniques that we can analyze within our framework can be found in Bredies and Zhariy (2013). In this work the results appear to be restricted to monotone methods for a quadratic residual function, but a practical stopping criterion based on the discrepancy principle is given. Another recent work dealing with this kind of problem is Beck and Sabach (2014) where a first order algorithm is proposed and convergence analysis including rates is provided. While we do not provide convergence rates, the theory we present requires less hypothesis on the objective functions and seems to give rise to a wider range of practical algorithms.

1.1 Contributions and Outline of the Paper

The main contribution of the present paper is as follows. The thread lead by Cabot and followed by Solodov is based on the idea of applying classical convex minimization algorithms to the ever-changing objective function $f_{0}+\eta_{k}f_{1}$. It started with the “tight” near-minimization from Cabot (2005) and evolved to the less stringent sufficient decrease policy of Solodov (2007, 2008). We here pave this way one step further by showing that the same principle is applicable to the more anarchic nonmonotone $\epsilon$-subgradient techniques.

Several unconstrained optimization algorithms have iterations that can be described as $\epsilon$-subgradient steps, among which we can mention the incremental subgradient methods Nedić and Bertsekas (2001), the aggregated incremental gradient of Blatt et al. (2007) when applied to the nondiferentiable case Helou Neto and De Pierro (2011), the recent incremental proximal method Bertsekas (2011), and Polyak’s heavy ball method (let $\tilde{\nabla}f(\boldsymbol{x})\in\partial f(\boldsymbol{x})$):

[TABLE]

The theory we develop will, therefore, cope with all the just mentioned cases simultaneously. That is, we show that application of any of these algorithms to the varying objective function $f_{0}+\eta_{k}f_{1}$ will converge to the solution of the bilevel problem (1) under assumptions on the stepsize $\lambda_{k}$ which are not much different from those required in the one level case.

We present numerical experimentation showing the effectiveness of the technique when applied to high-resolution micro-tomographic image reconstruction from simulated and from real synchrotron radiation illumination projection data. In this case the amount of data and the number of variables is very large, requiring efficient algorithms with computationally cheap iterations. In this context, important practical contributions are the introduction of certain perturbed, fista-inspired Beck and Teboulle (2009) algorithms, resulting in effective methods for problems like (1) with Lipschitz-differentiable $f_{0}$. Furthermore, when solving an instance with a non-differentiable $f_{0}$ composed as a sum of many convex functions, incremental techniques are very efficient in the first iterations, also resulting in good algorithmic performance. Both the theoretical analysis and the application of these practical algorithms to the bilevel problem (1) are new.

Our methods can be seen as perturbations of classical algorithms, in the spirit of the superiorization approach Garduño and Herman (2014). However, we show more powerful convergence results because we impose some structure on what would otherwise be called a superiorization sequence. We believe that this is a major contribution of the paper because opens the possibility of pursuing bilevel results alongside with superiorization techniques.

2 Theoretical Analysis

2.1 Stepsize Considerations

Recall from the theory of $\epsilon$-subgradient methods Correa and Lemaréchal (1993) for the one level case (that is, problem (1) with $f_{1}\equiv 0$) that convergence of iterations (4) to a solution, under mild extra assumptions of subgradient boundedness, can be ensured with slowly diminishing stepsizes satisfying:

[TABLE]

The non-summability hypothesis seems necessary. Vanishing stepsizes, however, may have the negative effect of slowing down asymptotic convergence of the algorithm. Therefore, owing to its computationally cheap iteration, methods like (4) are usually thought to be most competitive when the problem size is very large, or when highly accurate solutions are not required. However, in some important particular (with smooth primary objective function $f_{0}$) cases supported by the theory developed here, the stepsize $\lambda_{k}$ does not necessarily have to vanish and in such applications we obtain reasonably fast algorithms. Both stepsize regimes (decreasing and non-decreasing) are evaluated in the experimental work we present and we shall see that incremental techniques are efficient too in the cases where its good characteristics apply, even if with vanishing stepsizes.

For the classical (one level) convex optimization problem, proximal methods Rockafellar (1976) require $\lambda_{k}\geq\underline{\lambda}$ for some $\underline{\lambda}>0$. For the bilevel case, the same stepsize requirement, with an extra upper boundedness assumption, i.e.,

[TABLE]

would ensure convergence of (2) to the optimizer of (1). For such results to hold Cabot (2005), the extra assumption of a slow control:

[TABLE]

was made in order to ensure that the influence of $f_{1}$ throughout the iterations was strong enough while still becoming arbitrarily small. The alternative form below is more appropriate to us this time, as it will generalize immediately to our algorithms:

[TABLE]

Notice that if (6) holds, then (8) is equivalent to (7). However, because the net contribution to each iteration from $\partial f_{1}(\boldsymbol{x}_{k})$ in algorithm (4) is actually $O(\lambda_{k}\eta_{k})$, (8) generalizes to the case $\lambda_{k}\to 0^{+}$, while (7) does not.

2.2 Formal Algorithm Description

The stepping stone of our analysis will be an abstract three-step algorithm given as follows:

[TABLE]

where the operators $\mathcal{O}_{f_{i}}$ for $i\in\{0,1\}$ have specific conceptual roles and must satisfy certain corresponding properties, which we will discuss right next, and, for a nonempty convex and closed set $X$, $\mathcal{P}_{X}$ is the projector:

[TABLE]

Sequences $\{\lambda_{k}\}$ and $\{\mu_{k}\}$ are stepsize sequences. The first one plays, in this abstract setting, the same role it plays in iterative scheme (4), while sequence $\{\mu_{k}\}$ should be identified with $\{\lambda_{k}\eta_{k}\}$. Therefore, application of $\mu_{k}=\lambda_{k}\eta_{k}$ in (8) leads immediately to

[TABLE]

Let us then describe the properties required for the optimality operators $\mathcal{O}_{f_{i}}$. The imposed characteristics are easy to meet, as we later illustrate.

Property 1

There is $\beta>0$ such that for any $\lambda\geq 0$ and for all $\boldsymbol{x}_{k+i/3}$, $\boldsymbol{y}\in\mathbb{R}^{n}$, and $i\in\{0,1\}$:

[TABLE]

where $\rho_{i}(\lambda,k)$ represents an error term, with properties to be describe later.

Below the description of the next property, we give an example of a class of operators which satisfy this condition. Furthermore, Subsections 2.4, 2.5, and 3.4 bring four other instances that will be used in the experimental part of the paper: the projected gradient, the incremental subgradient, the proximal map and the iterated subgradient step. In fact, the key utility of this abstract definition is to be able to encompass several useful classical optimization steps while still ensuring sufficient qualities in order to provide convergence results. For this to be true, the error term will have to be controlled in a specific way, but, for every case we have found, the error term magnitude is bounded by a constant times the stepsize and this way we can always obtain convergent algorithms by selecting proper stepsize sequences.

Property 2

There exists $\gamma>0$ such that

[TABLE]

Property 1 guarantees that, going from some fixed $\boldsymbol{x}$, the operator $\mathcal{O}_{f_{i}}$ will approach a point $\boldsymbol{y}$ with improved $f_{i}$ value if only the result of the operation does not have a better $f_{i}$ value than $\boldsymbol{y}$ and the stepsize $\lambda$ is small enough. Property 2 is no more than a boundedness assumption on the operators which makes sure that the stepsize controls the magnitude of the movement.

These can be derived from somewhat standard hypothesis for $\epsilon$-subgradient algorithms (see, e.g., Correa and Lemaréchal (1993)) and, as such, a plethora of concrete realizations of such operators $\mathcal{O}_{f}$ is possible, the most obvious being $\epsilon$-subgradient steps:

[TABLE]

where $\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})\in\partial_{\epsilon}f(\boldsymbol{x})$. In this case we have:

[TABLE]

Denote $\boldsymbol{z}=\boldsymbol{x}-\lambda\tilde{\nabla}_{\epsilon}f(\boldsymbol{x})$. Then, convexity leads to

[TABLE]

where $\tilde{\nabla}f(\boldsymbol{z})\in\partial f(\boldsymbol{z})$. Multiplying the above inequality by $-2\lambda$ we get

[TABLE]

This with (10) gives

[TABLE]

so that we can satisfy Properties 1 and 2 for $\mathcal{O}_{f}=\mathcal{S}_{f}$ if we further assume $\epsilon$-subgradient boundedness (and consequently subgradient boundedness).

A straightforward generalization of the argument leading from (10) to (11), omitted for brevity, results in the following statement, which will be useful later:

Proposition 1

Assume an operator $\mathcal{O}_{f}:\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ satisfies, for $\lambda>0$

[TABLE]

where $\varrho(\lambda)$ is an error term, and

[TABLE]

for some $\gamma>0$. Then, we have:

[TABLE]

where \tilde{\nabla}f\bigl{(}\mathcal{O}_{f}(\lambda,\boldsymbol{x})\bigr{)}\in\partial f\bigl{(}\mathcal{O}_{f}(\lambda,\boldsymbol{x})\bigr{)}.

Our analysis will focus on algorithms more general than (4), allowing simple constraint sets $X_{0}$ to be handled. We recall that the projection onto a nonempty convex closed set $X_{0}$ satisfies:

Property 3

For all $\boldsymbol{x}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\in X_{0}$, we have

[TABLE]

2.3 Convergence Results

We introduce some simplifying notations:

• •

$f_{i}^{*}$ is the optimal value of $f_{i}$ over $X_{i}$ for $i\in\{0,1\}$ where;

• •

$X_{0}$ is given and $X_{i+1}:=\{\boldsymbol{x}\in X_{i}:f_{i}(\boldsymbol{x})=f_{i}^{*}\}$ for $i\in\{0,1\}$;

• •

$[x]_{+}:=\max\{0,x\}$;

• •

$d_{X}(\boldsymbol{x}):=\|\boldsymbol{x}-\mathcal{P}_{X}(\boldsymbol{x})\|$.

Our first result shows convergence of the iterates to the set of minimizers of $f_{0}$ over $X_{0}$. We next prove convergence to the set of minimizers of $f_{1}$ over $X_{1}$. Both of these preliminary results contain certain technical and some apparently strong hypothesis. We subsequently weaken and clarify such ad hoc requirements in order to obtain our main results.

Proposition 2

Assume that $X_{1}\neq\emptyset$, $X_{1}$ is bounded (or $\{\boldsymbol{x}_{k}\}$ is bounded) $\sum_{i=0}^{\infty}\lambda_{k}=\infty$, $\mathcal{O}_{f_{0}}$ and $\mathcal{O}_{f_{1}}$ satisfy Property 1, $\mathcal{O}_{f_{1}}$ satisfy also Property 2, $f_{1}(\boldsymbol{x}_{k+2/3})\geq\underline{f}>-\infty$, $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$, $\rho_{0}(\lambda_{k},k)\to 0$, $\rho_{1}(\mu_{k},k)\leq\overline{\rho}_{1}<\infty$, $\mu_{k}\to 0$, and $\mu_{k}/\lambda_{k}\to 0$. Suppose also that there exists $M$ such that $\forall k\in\mathbb{N}$ there is $\boldsymbol{v}_{k}\in\partial f_{0}(\boldsymbol{x}_{k})$ for which $\|\boldsymbol{v}_{k}\|<M$, then we have

[TABLE]

Proof

First notice that $\mu_{k}\to 0$ and Property 2 imply $\|\boldsymbol{x}_{k+1/3}-\boldsymbol{x}_{k+2/3}\|\to 0$. Then we take into consideration the non-expansiveness of the projection and of Property 1 of $\mathcal{O}_{f_{0}}$ and $\mathcal{O}_{f_{1}}$, there holds, for $\boldsymbol{y}\in X_{1}$:

[TABLE]

Denote $N=f_{1}(\boldsymbol{y})-\underline{f}$, thereby simplifying the above expression to:

[TABLE]

Then, the boundedness of $\partial f_{0}(\boldsymbol{x}_{k})$ leads to

[TABLE]

which together with (14) gives

[TABLE]

We shall denote, for $\delta\geq 0$:

[TABLE]

Notice that if $X_{1}$ is bounded (or $\{\boldsymbol{x}_{k}\}$ is bounded), then $X_{1}^{\delta}$ is bounded. Therefore, the following quantity is well defined:

[TABLE]

Furthermore, we have

[TABLE]

Let $\delta$ be any positive real number and consider, with $\rho_{1}(\mu_{k},k)\leq\overline{\rho}_{1}$, $\mu_{k}/\lambda_{k}\to 0$, $\|\boldsymbol{x}_{k+1/3}-\boldsymbol{x}_{k}\|\to 0$, and $\rho_{0}(\lambda_{k},k)\to 0$ in mind, that $k_{0}$ is large enough such that $k\geq k_{0}$ implies

[TABLE]

Then, two situations may occur:

$d_{X_{1}}(\boldsymbol{x}_{k})\geq\Delta_{1}(\delta)$; 2. 2.

$d_{X_{1}}(\boldsymbol{x}_{k})<\Delta_{1}(\delta)$.

Let us first suppose that Case 1 holds, that is $f_{0}(\boldsymbol{x}_{k})-f_{0}^{*}\geq\delta$. Then, for $k\geq k_{0}$, from (15) and (16) we get:

[TABLE]

In particular,

[TABLE]

Therefore, since $\sum_{k=0}^{\infty}\lambda_{k}=\infty$, there must exist an arbitrarily large $k_{1}\geq k_{0}$ such that $d_{X_{1}}(\boldsymbol{x}_{k})<\Delta_{1}(\delta)$.

Now, let us notice that, because of (12)

[TABLE]

Given the hypothesis, we may assume that $k_{0}$ is large enough such that, in addition to (16), we have also

[TABLE]

Therefore, for $k>k_{1}$, there holds:

[TABLE]

Since $\delta>0$ was arbitrary and $\lim_{\delta\to 0}\Delta_{1}(\delta)=0$, the claim is proven.∎

Proposition 3

Assume $X_{2}\neq\emptyset$, $X_{2}$ is bounded (or $\{\boldsymbol{x}_{k}\}$ is bounded), that $\mu_{k}\to 0$, $\sum_{i=0}^{\infty}\mu_{k}=\infty$, $\mathcal{O}_{f_{0}}$ and $\mathcal{O}_{f_{1}}$ satisfy Property 1, $\mathcal{O}_{f_{1}}$ also satisfies Property 2, $\lambda_{k}[f_{0}^{*}-f_{0}(\boldsymbol{x}_{k+1/3})]_{+}/\mu_{k}\to 0$, $d_{X_{1}}(\boldsymbol{x}_{k})\to 0$, $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$, $\lambda_{k}\rho_{0}(\lambda_{k},k)/\mu_{k}\to 0$ and $\rho_{1}(\mu_{k},k)\to 0$. Suppose also that there exists an $M$ such that $\forall k\in\mathbb{N}$ there are $\boldsymbol{v}_{k}\in\partial f_{0}(\boldsymbol{x}_{k})$ and \boldsymbol{w}_{k}\in\partial f_{1}\bigl{(}\operatorname{\mathcal{P}}_{X_{0}}(\boldsymbol{x}_{k})\bigr{)} for which $\|\boldsymbol{v}_{k}\|<M$ and $\|\boldsymbol{w}_{k}\|\leq M$, then we have

[TABLE]

Proof

Notice for later reference that just like in Proposition 2, the hypotheses imply that $\|\boldsymbol{x}_{k+1/3}-\boldsymbol{x}_{k+2/3}\|\to 0$. Now, if we use (13) with $\boldsymbol{y}\in X_{2}\subset X_{1}$, we get:

[TABLE]

Now let $\tilde{\nabla}f_{1}(\boldsymbol{x}_{k})\in\partial f_{1}(\boldsymbol{x}_{k})$ and then notice that convexity of $f_{1}$, Cauchy-Schwarz inequality and the boundedness assumption on $\partial f_{1}(\boldsymbol{x}_{k})$ lead to:

[TABLE]

Then, using (19) in (18) it is possible to obtain:

[TABLE]

Similarly to the $\Delta_{1}$ notation introduced above, we will denote, for $\delta\geq 0$:

[TABLE]

Notice that if $X_{2}$ is bounded (or $\{\boldsymbol{x}_{k}\}$ is bounded) and $d_{X_{1}}(\boldsymbol{x}_{k})$ is also bounded, then $X_{2}^{\delta}$ is bounded. Therefore, the following quantity is well defined:

[TABLE]

Furthermore, we have

[TABLE]

Given the hypothesis, for any fixed $\delta>0$, there is $k_{0}$ such that $k\geq k_{0}$ implies that

[TABLE]

We from now on assume $k>k_{0}$ and split in two different possibilities:

$f_{1}(\boldsymbol{x}_{k})>f_{1}^{*}+\delta$; 2. 2.

$f_{1}(\boldsymbol{x}_{k})\leq f_{1}^{*}+\delta$.

We start by analyzing Case 1. Using (21) in (20) we get, for $\boldsymbol{y}\in X_{2}$:

[TABLE]

In particular:

[TABLE]

Because of $\sum_{k=0}^{\infty}\mu_{k}=\infty$, this inequality means that there is an arbitrarily large $k_{1}$ such that $f_{1}(\boldsymbol{x}_{k_{1}})\leq f_{1}^{*}+\delta$.

Let us then focus on Case (2). We first notice that the assumed boundedness of \partial f_{1}\bigl{(}\operatorname{\mathcal{P}}_{X_{1}}(\boldsymbol{x}_{k})\bigr{)} leads to

[TABLE]

Therefore, $f_{1}(\boldsymbol{x}_{k})\leq f_{1}^{*}+\delta$ implies

[TABLE]

Thus, (17) now reads

[TABLE]

Then, because we have assumed $d_{X_{1}}(\boldsymbol{x}_{k})\to 0$, $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$, and $\|\boldsymbol{x}_{k+1/3}-\boldsymbol{x}_{k+2/3}\|\to 0$, we can recall $\lim_{\eta\to 0}\Delta_{2}(\delta+\eta)=\Delta_{2}(\delta)$, so that the argumentation above leads to the conclusion that

[TABLE]

Finally, because $\delta>0$ was arbitrary and $\lim_{\delta\to 0}\Delta_{2}(\delta)=0$, we have just proven the claimed result.∎

We now present two different algorithms and prove their convergence based on the above general results. Next section contains numerical experimentation regarding some of these methods in four different bilevel models arising in high-resolution micro-tomographic image reconstruction from synchrotron illumination.

2.4 An Algorithm for Lipschitz-Differentiable Primary Objective Functions

In this Subsection we suppose $f_{0}$ in the bilevel optimization problem (1) is differentiable with uniformly bounded and Lipschitz continuous gradient. For this kind of problem, we will consider Algorithm 1, which we name Fast Iterative Bilevel Algorithm (fiba). fiba first performs a projected gradient descent step, with a stepsize that does not change unless the magnitude of this operation is larger than a control sequence. This computation is then followed by the application of an optimality operator of the kind described by Properties 1 and 2. Such optimality operator is actually applied to a perturbation of the point obtained by the projected gradient descent, in a fashion similar to the Fast Iterative Soft-Thresholding Algorithm (fista) Beck and Teboulle (2009), but with the magnitude of the perturbation bounded by $\mu_{k}\zeta_{k}$, where $\{\zeta_{k}\}$ is a positive vanishing sequence.

In order to analyze convergence of Algorithm 1 through our previous results, we first look at the simple projected gradient descent

[TABLE]

as an instance of the optimality operators considered above. Let $L_{f}$ denote the Lipischtz constant of $\nabla f$. Then, if $\lambda\leq 1/L_{f}$, it is possible to show (see, e.g., Beck and Teboulle (2009) and references therein) that

[TABLE]

Let now $\iota_{X_{0}}$ be the indicator function:

[TABLE]

Then, inequality (22) leads to

[TABLE]

Therefore, (Beck and Teboulle, 2009, Lemma 2.3) can be used with $L=1/\lambda$, $g=\iota_{X_{0}}$, $\mathbf{y}=\boldsymbol{x}$ and $\mathbf{x}=\boldsymbol{y}$ in order to get, for $\boldsymbol{y}\in X_{0}$:

[TABLE]

Thus, for $\lambda\leq 1/L_{f}$, we can see that $\mathcal{G}_{f}(\lambda,\boldsymbol{x})$ satisfies Property 1 with $\beta=2$ and $\rho(\lambda,k)\equiv 0$. Notice that the operation described at line 7 of Algorithm 1 is actually:

[TABLE]

Consequently, according to (24), Algorithm 1 is an instance of (9) with an optimality operator $\mathcal{O}_{f_{0}}$ which satisfies Property 1 with $\rho_{0}(\lambda,k)\equiv 0$, whenever $\lambda\leq 1/L_{f}$ or if, which is weaker, (22) holds with $\boldsymbol{x}=\boldsymbol{x}_{k}$ and $\boldsymbol{y}=\boldsymbol{x}_{k+1/3}$. This characteristic of the error term implies that it is possible to have a precise enough operator without requiring $\lambda_{i_{k}}\to 0$. However, $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$ would still require $\lambda_{i_{k}}\to 0$ if the projection on line 7 of Algorithm 1 were not performed and this is the reason why there are two projections in this method.

Now, let us recall that $\boldsymbol{x}_{k}\in X_{0}$, so that by the definition of $\mathcal{G}_{f}$ and (12)

[TABLE]

Therefore, if the sequence $\{\nabla f_{0}(\boldsymbol{x}_{k})\}$ is bounded, and if $\lambda_{k}\to 0$ and $\zeta_{k}\to 0$, then the procedure in lines 8–12 of Algorithm 1 implies that $\|\boldsymbol{x}_{k+1/3}-\boldsymbol{x}_{k}\|\to 0$. Observe also that if $\sum_{k}\lambda_{k}=\infty$, then $\sum_{k}\lambda_{i_{k}}=\infty$ too.

Now, we consider the fact that the optimization operator for the secondary function $f_{1}$ is used in a perturbed point $\boldsymbol{y}_{k+1/3}$, instead of at $\boldsymbol{x}_{k+1/3}$. Our goal is to verify that the relevant properties of $\mathcal{O}_{f_{1}}$ are maintained. The first observation is that, given the way that $\boldsymbol{y}_{k+1/3}$ is defined, we have

[TABLE]

We then define a new operator $\tilde{\mathcal{O}}_{f_{1}}$, based on $\mathcal{O}_{f_{1}}$, as follows:

[TABLE]

Notice that $\tilde{\mathcal{O}}_{f_{1}}$ role is to hide the perturbation from the analysis. Also, this kind of operator is the reason why we use an iteration-dependent error term in Property 1, as we will see just below. Let us then assume that Property 1 holds for $O_{f_{1}}$, therefore:

[TABLE]

Now, by taking (26) into consideration, a straightforward computation leads to

[TABLE]

Then, using the above bound in (27) we have:

[TABLE]

That is, $\tilde{\mathcal{O}}_{f_{1}}$ satisfies Property 1 with $\rho_{1}$ replaced by

[TABLE]

where we notice that the set of points $\boldsymbol{y}$ where Property 1 is applied in the convergence proofs is bounded if $\{\boldsymbol{x}_{k}\}$ is bounded.

Given the above considerations, we are ready to provide the convergence results for Algorithm 1.

Theorem 2.1

Assume $f_{0}$ is differentiable with Lipschitz-continuous gradient and has Lipschitz constant $L_{0}$. Suppose too that $f_{0}$ has a bounded gradient and that $f_{1}$ has a bounded subgradient and $\{f_{1}(\boldsymbol{x}_{k+1/3})\}$ bounded from below. Assume $\{\lambda_{k}\}$, $\{\mu_{k}\}$ and $\{\zeta_{k}\}$ are non-negative vanishing scalar sequences such that $\sum_{k=0}^{\infty}\lambda_{k}=\infty$, $\lambda_{k}\leq 1/L_{0}$ (or each $\lambda_{i_{k}}$ satisfies (22)), $\sum_{k=0}^{\infty}\mu_{k}=\infty$, and $\mu_{k}/\lambda_{k}\to 0$. Then, if $X_{2}\neq\emptyset$, $\{\boldsymbol{x}_{k}\}$ is bounded, and $\mathcal{O}_{f_{1}}$ satisfies Properties 1 and 2 with $\rho_{1}(\mu_{k},k)\to 0$, we have

[TABLE]

Proof

First let us notice that Algorithm 1 can be written as

[TABLE]

where ${\tilde{\lambda}}_{k}:=\lambda_{i_{k}}$. Since the construction of the algorithm guarantees that $i_{k}\leq k$, we have $\sum_{k=0}^{\infty}{\tilde{\lambda}}_{k}=\infty$ and $\mu_{k}/{\tilde{\lambda}}_{k}\to 0$. Because ${\tilde{\lambda}_{k}}\leq 1/L_{0}$, (24) holds and therefore, $\mathcal{G}_{f_{0}}$ satisfies Property 1 with $\rho_{0}({\tilde{\lambda}}_{k},k)\equiv 0$. Furthermore, because (25) and the algorithm definition, as already argued, we have $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$. Also, if $\mathcal{O}_{f_{1}}$ satisfies Property 2 so does ${\tilde{\mathcal{O}}}_{f_{1}}$. Furthermore, as shown above, if $\mathcal{O}_{f_{1}}$ satisfies Property 1 so does ${\tilde{\mathcal{O}}}_{f_{1}}$, with the error term given by the factor multiplying $\mu_{k}$ in the second line of (28). Thus, the assumed boundedness of $\{\boldsymbol{x}_{k}\}$ and of $f_{1}$ ensure that Proposition 2 can be applied so that

[TABLE]

Now, because $\boldsymbol{x}_{k+1/3}\in X_{0}$, we then have $f_{0}(\boldsymbol{x}_{k+1/3})\geq f_{0}^{*}$. Furthermore, because of the boundedness assumptions and of $\mu_{k}\to 0$, it is possible to see that $\tilde{\rho}_{1}(\mu_{k},k)\to 0$. Therefore, Proposition 3 can be applied, which leads to the desired conclusion.∎

2.5 Incremental Algorithms for Non-Differentiable Problems

Here we specialize Algorithm (9) to the case where $f_{0}$ is the sum of several non-differentiable convex functions:

[TABLE]

In this situation we propose the use, for the primary optimization problem, of the incremental subgradient operator, denoted as $\mathcal{I}_{f}:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$, given by:

[TABLE]

Incremental operators are well known for its fast initial convergence rate and, accordingly, several variations of it have been thoroughly analyzed in the literature Solodov and Zavriev (1998); Solodov (1998); Nedić and Bertsekas (2001); Blatt et al. (2007); De Pierro and Yamagishi (2001); Bertsekas and Tsitsiklis (2000); Bertsekas (1997). We will use here the result (Nedić and Bertsekas, 2001, Lemma 2.1):

Lemma 1

Assume the subgradients of the convex functions $f_{0}^{i}$ are bounded in the following sense:

[TABLE]

Then, the incremental subgradient operator satisfies, for every $\lambda\in\mathbb{R}_{+}$, and $\boldsymbol{y},\boldsymbol{x}\in\mathbb{R}^{n}$:

[TABLE]

where $f_{0}:=\sum_{i=1}^{m}f_{0}^{i}$.

Now, notice that the boundedness condition on the subdifferentials leads, for every $\boldsymbol{x}\in\mathbb{R}^{n}$, to

[TABLE]

Thus, applying Proposition 1 we are lead to the following result:

Corollary 1

Assume the subgradients of the convex functions $f_{0}^{i}$ satisfy (29). Then, the incremental subgradient operator satisfies, for every $\lambda\in\mathbb{R}_{+}$, and $\boldsymbol{y},\boldsymbol{x}\in\mathbb{R}^{n}$:

[TABLE]

where, again, $f_{0}:=\sum_{i=1}^{m}f_{0}^{i}$.

The second algorithm we propose in this work will be called iiba, from Incremental Iterative Bilevel Algorithm, and is described in Algorithm 2 below. For this algorithm we have the following convergence result.

Theorem 2.2

Assume that $f_{0}$ is of the form $f_{0}:=\sum_{i=1}^{m}f_{0}^{i}$ and satisfies (29), that $f_{1}$ has a bounded subgradient, and that $\{f_{1}(\boldsymbol{x}_{k+1/3})\}$ is bounded from below. Assume $\{\lambda_{k}\}$ and $\{\mu_{k}\}$ are non-negative vanishing scalar sequences such that $\sum_{k=0}^{\infty}\lambda_{k}=\infty$, $\sum_{k=0}^{\infty}\mu_{k}=\infty$, $\mu_{k}/\lambda_{k}\to 0$ and $\lambda_{k}^{2}/\mu_{k}\to 0$. Then, suppose $X_{2}\neq\emptyset$ and $X_{2}$ is bounded (or $\{\boldsymbol{x}_{k}\}$ is bounded), and $\mathcal{O}_{f_{1}}$ satisfies Properties 1 and 2 with $\rho_{1}(\mu_{k},k)\to 0$. Then, the sequence $\{\boldsymbol{x}_{k}\}$ generated by Algorithm 2 satisfies

[TABLE]

Proof

Notice that Lemma 1 together with Proposition 1 and the subgradient boundedness assumption imply that $\mathcal{I}_{f_{0}}$ satisfies the desired Property 1. Also, $\lambda_{k}\to 0$ implies, together with the subgradient boundedness assumption, that $\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k+1/3}\|\to 0$ and that $\rho_{0}(\lambda_{k},k)\to 0$, the latter because of Lemma 1. Therefore, Proposition 2 implies that

[TABLE]

Now, notice the subgradient boundedness assumption and (31) imply that $[f_{0}(\boldsymbol{x}_{k})-f_{0}(\boldsymbol{x}_{k+1/3})]_{+}=O(\lambda_{k})$, and, therefore, since $\boldsymbol{x}_{k}\in X_{0}$ for $k>0$, we have $[f_{0}^{*}-f(\boldsymbol{x}_{k+1/3})]_{+}=O(\lambda_{k})$. Thus, $\lambda_{k}^{2}/\mu_{k}\to 0$ implies $\lambda_{k}[f_{0}^{*}-f(\boldsymbol{x}_{k+1/3})]_{+}/\mu_{k}\to 0$. Furthermore, since, by (30), $\rho_{0}(\lambda_{k},k)=O(\lambda_{k})$, $\lambda_{k}^{2}/\mu_{k}$ also implies $\lambda_{k}\rho_{0}(\lambda_{k},k)/\mu_{k}\to 0$. So, finally, Proposition 3 can be applied, which proves the result.∎

2.6 Stopping Criterion

Here we devise a stopping criterion for the proposed bilevel methods, based on inequalities (15) and (20) coupled to the following Lemma:

Lemma 2

Let $\{a_{k}\}$ and $\{b_{k}\}$ be non-negative sequences such that $a_{k}/b_{k}\to 0$ and $\sum_{k=0}^{\infty}b_{k}=\infty$. Then

[TABLE]

Proof

Choose any $\alpha>0$ and let $k_{0}$ be such that $a_{k}/b_{k}\leq\alpha$ for every $k\geq k_{0}$. Then

[TABLE]

But

[TABLE]

Therefore, $\lim_{\to\infty}\left|\sum_{k=0}^{n}a_{k}\middle/\sum_{k=0}^{n}b_{k}\right|\leq\alpha$ for any $\alpha>0$.∎

In order to explain the stopping criterion, we resort to the concept of best-so-far iteration. Let us denote as $\phi_{i}^{k_{0},k}$, for $i\in\{0,1\}$ and $k_{0}\leq k$ integers in $\{0,1,\dots\}$, the smallest value in the set $\{f_{i}(\boldsymbol{x}_{k_{0}}),f_{i}(\boldsymbol{x}_{1}),\dots,f_{i}(\boldsymbol{x}_{k})\}$. For the special case $k_{0}=0$, we simplify the notation by $\phi_{i}^{k}:=\phi_{i}^{0,k}$. Then, successive application of (15) together with $\phi_{i}^{k}\leq f_{i}(\boldsymbol{x}_{j})$ for all $j\leq k$ leads to

[TABLE]

Therefore, under the hypothesis of Proposition 2, Lemma 2 ensures that the right-hand side of the above inequality vanishes as the iterations proceed. Thus, the quantity $\sigma_{0}^{k}$ can be used as a measure of convergence to $X_{1}$ bounding the difference between the best $f_{0}$ function value to the optimal $f_{0}^{*}$, as long as it is possible to estimate the distance $\|\boldsymbol{x}_{0}-\boldsymbol{x}^{*}\|$, the behavior of the error terms and the subgradient bounding constants. Notice that, in principle, the constant $N$ require knowledge of the optimal value in this case, but it can be replaced by an upper bound for it, which should not be difficult to obtain in many cases.

We can use a very similar reasoning in order to estimate optimality of the secondary objective function too. Applying (20) repeatedly and recalling that $\phi_{1}^{k_{0},k}\leq f_{1}(\boldsymbol{x}_{i})$ for every $i\in\{k_{0},k_{0}+1,\dots,k\}$ we have

[TABLE]

Now, under the hypothesis of Proposition 3, Lemma 2 ensures that for any $k_{0}$, we have $\lim_{k\to\infty}\sigma_{0}^{k_{0},k}=0$.

We can now describe how to stop the algorithm at a non-negative integer $\kappa_{\epsilon_{0},\epsilon_{1}}$ such that we have

[TABLE]

for any pair of positive numbers $\epsilon_{0}$ and $\epsilon_{1}$. Let us consider the following procedure:

Iterate the algorithm until $\sigma_{0}^{k}\leq\epsilon_{0}$; 2. 2.

$k_{0}\leftarrow k$, $\kappa\leftarrow k$; 3. 3.

Iterate the algorithm until $\sigma_{1}^{\kappa,k}\leq\epsilon_{1}$; 4. 4.

Let $k_{1}\geq\kappa$ be such that $f_{1}(\boldsymbol{x}_{k_{1}})=\phi_{1}^{\kappa,k}$; 5. 5.

If $f_{0}(\boldsymbol{x}_{k_{1}})\leq\phi_{0}^{k_{0}}$: STOP; 6. 6.

$\kappa\leftarrow k$; Go to step 3.

Notice that once the procedure has stopped, then

[TABLE]

and

[TABLE]

The procedure indeed stops because since Proposition 2 ensures convergence in norm to $X_{1}$, mild subgradient boundedness assumptions will guarantee also that $f_{0}(\boldsymbol{x}_{k})\to f_{0}^{*}$. Therefore we know that for large enough $k$, there holds $f_{0}(\boldsymbol{x}_{k})\leq\phi_{0}^{k_{0}}$. That is, the condition in step 5 will eventually be satisfied as the iterations proceed.

3 Application Problem Presentation

We have performed experiments involving tomographic image reconstruction. In tomography, the idealized problem is to reconstruct a function $\mu:\mathbb{R}^{2}\to\mathbb{R}$ given the values of its integrals along straight lines, that is, given its Radon transform denoted as $\mathcal{R}[\mu]$ and defined by the following equality:

[TABLE]

Figure 1 (adapted from Helou et al. (2014)) brings a graphical representation of this definition.

Data was obtained by transmitting synchrotron radiation through samples of eggs taken from a fish of the species Prochilodus lineatus collected at the Madeira River‘s bed, immersed in distilled water inside a capillary test tube. Data acquisition was performed at the Brazilian National Synchrotron Light Source (lnls)111http://lnls.cnpem.br/.

In a transmission tomography setup Kak and Slaney (1988); Natterer (1986); Herman (1980) like the one we have used, the value of the line integral is estimated through the emitted to detected intensity ratio according to Beer–Lambert law:

[TABLE]

where $I_{e}$ is the emitted intensity, $I_{d}$ is the detected intensity, $\mu$ gives the linear attenuation factor of the imaged object at each point in space, and $L$ is the straight line connecting detector to emitter. While in this case the reconstruction problem is essentially bi-dimensional, a simultaneous acquisition of a radiography of $2048$ parallel slices of the object to be imaged is made at each angle, which enables volumetric reconstruction, if desired, by the stacking of several bi-dimensional reconstructions. After each plain x-ray imaging, the sample is rotated and new samples of the Radon transform are estimated in the same manner at a new angle. Figure 2 depicts the process of assembling the $2048\times 200$ data array, which will be used for a slice reconstruction, from the $200$ images of size $2048\times 2048$.

In our application the imaged subject is sensitive in a way such that a overly long exposure time under a low energy x-ray beam may overheat or otherwise physically damage the sample. Therefore, because the exposure time for good radiographies under a monochromatic beam at LNLS’ facilities was experimentally found to be at least $20$ seconds, the Radon Transform was sampled at only $200$ evenly spaced angles covering the interval $[-\pi,0]$, a relatively small number if we are willing to reconstruct full resolution $2048\times 2048$ images from this data. In this case, it is likely that problem (34) will have many solutions and we need to select one of these, therefore the need of a bilevel model arises.

3.1 Primary Objective Functions

Assuming the original image $g:\mathbb{R}^{2}\to\mathbb{R}_{+}$ lies in a finite dimensional vector space generated by some basis $\{g^{1},g^{2},\dots,g^{n}\}$ and recognizing that the number of measurements is always finite in practice, one can reduce the problem of tomographic reconstruction to a linear system of equations:

[TABLE]

where the elements $r_{ij}$ of the matrix $R$ are given by

[TABLE]

and the elements $b_{i}$ of the vector $\boldsymbol{b}$ are the corresponding experimental data, that is, $b_{i}$ is an approximate sample of $\mathcal{R}[g](\theta_{i},t_{i})$, where $g=\sum_{i=1}^{n}x_{i}g^{i}$ is the desired image. Because actual microtomographic data from synchrotron illumination will always contain errors (either from micrometric misalignment of the experimental setup, small dust particles in the emitter-detector path, or from statistical fluctuation of photon emission and attenuation), the above linear system of equations may not have a solution.

3.2 Differentiable Primary Objective Function

System (33) can be replaced by a constrained least squares problem, in order to alleviate the likely lack of consistency:

[TABLE]

Therefore, in this case we will be lead to a bilevel problem of the form (1) with $X_{0}=\mathbb{R}_{+}^{n}$.

Another option for a feasibility set would be the bounded box $\{\boldsymbol{x}\in\mathbb{R}^{n}:0\leq x_{i}\leq u\}$ for some $u>0$. This may be a sensible idea if the maximum attenuation factor is known beforehand and we could thus include this information into the model. Because we do not make this kind of boundedness imposition, by using the least squares function we are at the risk of violating certain (sub)gradient boundedness assumption. We therefore use another continuously differentiable convex function $f_{0}$, which has an uniformly bounded gradient and is such that

[TABLE]

In order to build such function, take any $\tilde{\boldsymbol{x}}\in X_{0}$ (such as, in our case, $\tilde{\boldsymbol{x}}=\boldsymbol{0}$) and set

[TABLE]

then define

[TABLE]

where $r_{i}^{T}$ is the $i$-th row of $R$ and the function $h:\mathbb{R}\to\mathbb{R}_{+}$ is given as

[TABLE]

Notice that if $q(\boldsymbol{x})\leq q(\tilde{\boldsymbol{x}})=1/2\Delta^{2}$, then $f_{0}(\boldsymbol{x})=q(\boldsymbol{x})$ and, furthermore, if $f_{0}(\boldsymbol{x})\neq q(\boldsymbol{x})$, then $f_{0}(\boldsymbol{x})\geq 1/2\Delta^{2}=q(\tilde{\boldsymbol{x}})$. Therefore, (35) holds.

We need to show that the uniformly bounded continuous differentiability and convexity claims hold, but these facts follow from convexity and continuous differentiability of $h$. Its derivative can be computed to be

[TABLE]

with $\operatorname{sign}:\mathbb{R}\to 2^{\mathbb{R}}$ defined as

[TABLE]

Thus, the derivative of $h$ is uniformly bounded by $|h^{\prime}(x)|\leq 2\Delta$, is continuous as long as we use $\Delta>0$, and, since $h^{\prime}(x)$ is nondecreasing, $h$ is convex. Notice that $h$ is the well known Huber function Huber (1964).

In order to simplify the notation, we introduce the function $\boldsymbol{h}^{\prime}:\mathbb{R}^{m}\to\mathbb{R}^{m}$ given as:

[TABLE]

Under this notation, a straightforward calculation leads to

[TABLE]

The computationally expensive parts of this expression are products of the form $R\boldsymbol{x}$ and $R^{T}\boldsymbol{b}$. In our case, because of image resolution and dataset size, as explained in the previous Subsection, matrix $R$ has dimensions of over $4\cdot 10^{5}$ lines by $4\cdot 10^{6}$ columns.

In order to obtain a computationally viable algorithm, these matrix-vector products must be amenable to fast computation, and, in fact, there is a very effective approach to it based in a relation between the one-dimensional Fourier Transform of the data and one “slice” of the two-dimensional Fourier Transform of the original image. Recall that we denote our desired, unknown, image as $g:\mathbb{R}^{2}\to\mathbb{R}_{+}$ and define the projection of $f$ at an angle $t$ as

[TABLE]

Now, using the hat notation for the Fourier transform of an integrable function $f:\mathbb{R}^{n}\to\mathbb{C}$, as defined as follows, where $\jmath:=\sqrt{-1}$:

[TABLE]

there holds the so called Fourier-slice theorem Kak and Slaney (1988); Natterer (1986):

[TABLE]

This result implies that the action of $R$ can be evaluated first in the Fourier space, an operation which can be computed efficiently by Non-uniform Fast Fourier Transforms (nfft) Fourmont (2003) and later translated back to the original feature space by regular Fast Fourier Transforms (fft). We have used the pynfft binding for the nfft3 library Keiner et al. (2009) in order to compute samples of $\hat{p}_{\theta}$ following (3.2), and then used regular inverse fft routines from numpy for the final computation from the Fourier back to the feature space. Fast evaluation of the transpose operation is immediately available from the same set of libraries using the fft from numpy and the transpose of the nfft available from nfft3.

3.3 Non-Differentiable Primary Objective Function

Another option to circumvent the non-consistency of (33) under non-negativity constraints is to use the least 1-norm approach, which is useful because synchrotron illuminated tomographic data contains mostly very mild noise, but also has sparsely distributed highly perturbed points for example caused by small dust particles and other detector failures such as those caused by defective scintillator crystals and other mechanical, thermal or electronic causes Miqueles et al. (2014).

Therefore, our second option for primary optimization problem was given by

[TABLE]

Where, now, the primary objective function naturally has an everywhere bounded subdifferential and, as is well known, the $\|\cdot\|_{1}$ is more forgiving to sparse noise than is the euclidean norm Donoho and Logan (1992); Santosa and Symes (1986); Taylor et al. (1979); Claerbout and Muir (1973).

3.4 Secondary Objective Functions

For the purpose of selecting one among the solutions of (34) or (36), we will consider two particular cases of the bilevel program (1).

3.4.1 Haar 1-Norm

Here the secondary objective function $f_{1}$ is given by

[TABLE]

with $H\in\mathbb{R}^{n\times n}$ orthonormal, that is, satisfies $H^{T}H=I$, where $I$ is the identity matrix and superscript $T$ indicates transposition. Transformation $H$ is usually a sparsefying transform, in our case the Haar transform, and one is looking for a sparse (in some cases the sparsest Candès et al. (2006a, b); Donoho (2006a, b, c)), in the $H$-transformed space, optimizer of $f_{0}$ over $X_{0}$.

Notice that for this particular function we have

[TABLE]

where the set-valued function $\boldsymbol{\operatorname{sign}}:\mathbb{R}^{n}\to 2^{\mathbb{R}^{n}}$ is given by

[TABLE]

Let us now consider the soft-thresholding operator, given componentwise as

[TABLE]

We then define an intermediary optimality operator as

[TABLE]

A straightforward computation convinces us that

[TABLE]

where

[TABLE]

This means that $\overline{\boldsymbol{x}}=\mathcal{N}_{f_{\text{Haar}}}(\mu,\boldsymbol{x})$, and we have denoted $\tilde{\nabla}f_{1}(\overline{\boldsymbol{x}})\in\partial f_{\text{Haar}}(\overline{\boldsymbol{x}})$. Now we proceed:

[TABLE]

and we then can observe from (37) that $\partial f_{1}$ is bounded in a way such that $\|\tilde{\nabla}f_{1}(\boldsymbol{x})\|\leq\sqrt{n}$ for every $\boldsymbol{x}$, and from (40) that $\|\boldsymbol{x}-\overline{\boldsymbol{x}}\|\leq\mu\sqrt{n}$ for every $\boldsymbol{x}$. Using this facts followed by the subgradient inequality we have

[TABLE]

From where we see that Property 1 is satisfied with $\rho_{1}(\mu,k)=3\mu n$. Furthermore, since for every $\boldsymbol{x}$ we have $\|\tilde{\nabla}f_{\text{Haar}}(\boldsymbol{x})\|\leq\sqrt{n}$, formulation (39) of $\mathcal{O}_{f_{1}}$ implies that this operator satisfies Property 2 with $\gamma=\sqrt{n}$.

It is worth remarking that the soft-thresholding is an example of a proximal operator,

[TABLE]

In fact, it is possible to carry out an argumentation similar to the one that lead to inequality (41) above for more general proximal operators, because they too can be interpreted as examples of implicit subgradient algorithms or, alternatively, as appropriately controlled $\epsilon$-subgradient methods Correa and Lemaréchal (1993).

3.4.2 Total Variation

Another function we have used to select among the primary problem optimizers is the Total Variation Rudin et al. (1992):

[TABLE]

where we have assumed that $n$ is a perfect square, which is true for our reconstructed images, the lexicographical notation:

[TABLE]

and also the following boundary conditions:

[TABLE]

The $f_{\text{{tv}}}$ subdifferential is bounded and, therefore, application of a subgradient step would satisfy the required operator properties. However given the low computational time of the $f_{\text{{tv}}}$ subgradient step, it is reasonable to apply it several times instead of a single pass. In this case we can see that the resulting operation would still satisfy the required properties, as follows. Let us denote by $\hat{\mathcal{O}}_{f}^{J}$ the following operator:

[TABLE]

That is, the $J$-fold repetition of $\tilde{\mathcal{O}}_{f}$, with diminishing stepsizes. For this kind of operator there holds the following:

Lemma 3

Suppose that $f$ has a bounded subgradient and that $\tilde{\mathcal{O}}_{f}$ satisfies Properties 1 and 2. Then $\hat{\mathcal{O}}_{f}^{J}$ also satisfies Properties 1 and 2.

Proof

Repeated use of Property 2 and the triangle inequality leads to

[TABLE]

More generally, telescoping from $j$ to $l$ we have

[TABLE]

Now, repeated use of Property 1 gives

[TABLE]

Then, by subgradient boundedness and from (43), we have:

[TABLE]

Therefore, $\hat{\mathcal{O}}_{f}^{J}$ also satisfies Property 2, with error term given by

[TABLE]

Notice that if the original error term $\rho(\lambda,k)$ for $\tilde{\mathcal{O}}_{f}$ is $O(\lambda)$, then so is the one of the iterated operator $\hat{\mathcal{O}}_{f}^{J}$. This remark will be useful when considering convergence of the incremental algorithm to be presented in Subsection 2.5.

4 Numerical Experimentation

4.1 Smooth Primary Objective Function

In the present Subsection, we report the results obtained using Algorithm 1 in order to approximately solve the optimization problem

[TABLE]

with both $f_{1}(\boldsymbol{x})=\|H\boldsymbol{x}\|_{1}$ and $f_{1}(\boldsymbol{x})=f_{\text{{tv}}}(\boldsymbol{x})$. Below we describe the parameter selection and explain the rationale supporting each choice.

The bilevel algorithms were compared against the Fast Iterative Soft-Thresholding Algorithm (fista) Beck and Teboulle (2009) for the function

[TABLE]

where $\iota$ is the indicator function (23). We have used the starting image described in the next paragraph. We iterated fista with a constant stepsize $\lambda=3$, for the reasons explained in the paragraph following the next.

Starting image

The initial guess $\boldsymbol{x}_{0}$ used in the experiments presented in this subsection was the zero image.

Stepsize sequences

The sequence $\{\lambda_{k}\}$ was set to be

[TABLE]

where $\lambda$ was chosen to be just small enough so that the squared residual of the first iterations of the one-level algorithm is decreasing. In the examples, $\lambda=3$ worked well. This choice for $\lambda$ was noticed to provide the fastest primary objective function decrease during the iterations when using stepsize sequences of the form (44) above. Sequence $\{\eta_{k}\}$ was given by

[TABLE]

which was chosen to attain very high values in order not to negatively influence convergence speed by limiting the values of $\lambda_{i_{k}}$ or $\xi_{k}$ (both of which remained constant throughout the iterations of Algorithm 1, without detected need for decreasing). Finally, the stepsize sequence $\{\mu_{k}\}$ was given by

[TABLE]

Now, the initial stepsize $\mu$ has been selected in order to make the first pair of subiterations to have a specific relative strength. The precise procedure was the following: we first computed $\boldsymbol{x}_{1/3}$ as usual, following Algorithm 1 and using the already chosen $\lambda_{0}$. Then, a tentative subiteration $\tilde{\boldsymbol{x}}_{2/3}$ is computed using a tentative stepsize $\tilde{\mu}_{0}=1$. Finally, the value to be used as starting stepsize is computed by

[TABLE]

So, the step given by the first subiteration for the primary problem is about $10^{2}$ times the step given by the subiteration for the secondary optimization problem. This value provided a good compromise between primary problem convergence and secondary function value during the iterations.

Secondary Operators

When the problem being solved had $f_{1}=f_{\text{Haar}}$, we used fiba with $\mathcal{O}_{f_{1}}=\mathcal{N}_{f_{\text{Haar}}}$ according to (38). If the problem had $f_{1}=f_{\text{{tv}}}$, then we have used $\mathcal{O}_{f_{1}}=\hat{\mathcal{O}}_{f_{\text{{tv}}}}^{10}$ as defined in (42) with $\tilde{\mathcal{O}}_{f_{\text{{tv}}}}(\lambda,\boldsymbol{x})=\boldsymbol{x}-\lambda\tilde{\nabla}f_{\text{{tv}}}(\boldsymbol{x})$ where $\tilde{\nabla}f_{\text{{tv}}}(\boldsymbol{x})\in\partial f_{\text{{tv}}}(\boldsymbol{x})$.

Algorithm Convergence

Because all entries of the matrix $R$ are nonnegative and because every pixel of the reconstructed image is crossed by several rays during acquisition, the model has the property that $X_{1}$ is bounded. Therefore, given the subgradient boundedness of all involved objective functions and the fact that both secondary objective functions $f_{\text{Haar}}$ and $f_{\text{{tv}}}$ are bounded from below, Theorem 2.1 could be applied to show that Algorithm 1 converges in the cases covered in the present subsection, except for two issues. The first problem is that $\lambda_{0}$ may not be smaller than $1/L_{0}$, but we have observed during algorithm evaluation that the alternative sufficient decrease criterion (22) was satisfied in every iteration, which suffices, instead, for convergence. Furthermore, boundedness of the iterates was also observed, ensuring convergent behavior without forceful truncation.

Numerical Results

In the plots that follow, we refer to fiba when applied to the model with $f_{1}=f_{\text{Haar}}$ as fiba-h and when applied to the model with $f_{1}=f_{\text{{tv}}}$ as fiba-tv. Figure 3 shows that, as expected, the bilevel approach influences the iterates of the method, thereby resulting in lower secondary objective function value throughout the iterations. Figure 4, on the other hand and also unsurprisingly, shows that the convergence speed for the primary optimization problem is reduced as the secondary optimization step influences the iterations. Figure 5 displays some resulting images, all with the same prescribed primary objective function value, from the algorithms and Figure 6 brings a plot of the profiles across the same line of each image.

We can see the benefits of the enforced smoothing, in particular the good quality of the Total Variation image. The advantage of the bilevel technique here is that the regularization level is chosen by a stopping criterion instead of a parameter in an optimization problem. This latter situation implies that more computational effort would be required since reconstruction for each tentative parameter value would require the equivalent to a large fraction of the computation effort of performing some iterations of the bilevel techniques. We further would like to remark that although heuristic, the accelerating scheme does substantially speed up the method. While we do not show these results in the present paper, when using a non-accelerated technique (i.e, $\eta_{k}=0$), the convergence speed is reduced to a rate similar to the ista Beck and Teboulle (2009) method.

4.2 Non-Differentiable Primary Objective Function

Now we consider using Algorithm 2 applied to optimization problem

[TABLE]

with $f_{1}(\boldsymbol{x})=f_{\text{{tv}}}(\boldsymbol{x})$. We do not experiment using $f_{1}(\boldsymbol{x})=\|H\boldsymbol{x}\|_{1}$ here, we leave it to the simulated experiments presented in the next subsection.

Reconstructions were performed using algorithms where the data was divided in $s$ subsets with $s\in\{1,2,4,8,16,32\}$. Each subset was itself composed by data comprising several projection measurements in sequence, i.e., each subset was a vertical stripe of the image222Coincidently, the vertical colored lines delimit the subsets for the case $s=4$. shown in the bottom of Figure 2. At each iteration, the sequence of subset processing was selected by a pseudo-random shuffling. This makes the algorithm non-deterministic, but still covered by the theory because every subset was used once in every iteration and each of the corresponding subdifferentials is uniformly bounded. Furthermore, we have observed a consistent behavior among runs and we had not observed a run where sequential data processing led to better convergence than the random ordering.

The various incremental bilevel algorithms were pairwise compared against their pure projected incremental subgradient counterparts, i.e., a variation of Algorithm 2 with $\mu_{k}\equiv 0$. All of the methods where started with an uniform image as described for the differentiable model. We will denote the bilevel algorithms by iiba-$s$, where $s$ is the number of subsets, and the projected incremental method by inc-$s$.

Starting image

The initial guess $\boldsymbol{x}_{0}$ used in the experiments of the present section, for all algorithms tested, was a constant image such that $\sum_{i=1}^{m}(R\boldsymbol{x}_{0})_{i}=\sum_{i=1}^{m}b_{i}$. It is easy to compute the correct constant value $\alpha$ from $\alpha=\sum_{i=1}^{m}b_{i}/\sum_{i=1}^{m}(R\boldsymbol{1})_{i}$ where $\boldsymbol{1}$ is the vector, of appropriate dimension, with every coordinate equal to $1$. This choice makes sure that the Radon consistency condition $\sum_{i=1}^{m}(R\boldsymbol{x})_{i}=\sum_{i=1}^{m}b_{i}$ is satisfied for the first iteration, potentially avoiding large oscillations in the first steps of the algorithm.

Stepsize sequences

The sequence $\{\lambda_{k}\}$ for the incremental algorithms with $s$ subsets was set to be

[TABLE]

where $\lambda$ was of the form

[TABLE]

and $\tilde{\nabla}f_{0}(\boldsymbol{x}_{0})\in\partial f(\boldsymbol{x}_{0})$. The pair $(\alpha_{s},\epsilon_{s})$ was selected through a simple search procedure as follows. Let us denote by $\boldsymbol{x}_{(s,\alpha,\epsilon)}$ the first iteration to be completed past $4$ seconds of computation time by the pure projected incremental algorithm with $s$ subsets (that is, by inc-$s$) and using parameters $(\alpha,\epsilon)$. Then $(\alpha_{s},\epsilon_{s})$ was given by

[TABLE]

The parameters were only optimized for the non-bilevel case, and the same values were used for the corresponding (with relation to the number of subsets) bilevel algorithm. The second step sequence $\{\mu_{k}\}$ was prescribed as

[TABLE]

with

[TABLE]

where the rationale is the same than in (45), but with target relative importance between the first subiterations of $10$.

Secondary Operators

In the present experiments we have used $\mathcal{O}_{f_{1}}=\hat{\mathcal{O}}_{f_{\text{{tv}}}}^{5}$ as defined in (42) with $\tilde{\mathcal{O}}_{f_{\text{{tv}}}}(\lambda,\boldsymbol{x})=\boldsymbol{x}-\lambda\tilde{\nabla}f_{\text{{tv}}}(\boldsymbol{x})$ where $\tilde{\nabla}f_{\text{{tv}}}(\boldsymbol{x})\in\partial f_{\text{{tv}}}(\boldsymbol{x})$.

Algorithm Convergence

Notice that the stepsizes are of the form

[TABLE]

where $\epsilon\in[0.5,0.9]$ and $\lambda$ and $\mu$ are nonnegative. It is routine to check that for this range of $\epsilon$ there holds:

[TABLE]

Also, for the same reasons as in the differentiable primary problem case, the model has the property that $X_{1}$ is bounded. Given the subgradient boundedness of all involved objective functions and the fact that secondary objective function $f_{\text{{tv}}}$ is bounded from below, Theorem 2.2 can be applied to prove algorithm convergence.

Numerical Results

By denoting $R_{i}$, $i\in\{1,2,\dots,s\}$, the matrix with the rows corresponding to the $i$-th subset, we notice that the computationally demanding parts of the algorithm are products of the form

[TABLE]

because the partial subgradients are given by

[TABLE]

We were not able to make the nfft library as efficient for such partial matrix-vector products, which imposed a large overhead in the partial iterations. We have instead used a ray-tracing algorithm Han et al. (1999); Siddon (1985) implemented to run in gpus (Graphics Processing Units) under single precision floating point arithmetics.

A special scheme was used so that access to slow gpu memory is minimized by performing the computations in sub-images loaded to/from the gpu’s shared memory (a $64$kb fast L1-cache-like memory) in coalesced reads/writes, and summing up the partial results. In this setting, the sequential computation of the $s$ different partial subgradients takes longer than the computation of the subgradient itself, because there are more memory copy to/from shared memory. Yet another per-iteration overhead of the incremental methods are the multiple subiteration updates. Therefore one iteration of the incremental method with $s$ subsets still takes slightly longer than with $s-1$ subsets.

Even with mandatory overheads, a careful implementation was able to make the iteration-wise speed up provided by the incremental approach advantageous time-wise, as can be seen in Figure 7. An important feature in this plot is that this speed up is retained by the bilevel algorithms in a similar fashion to the non-bilevel incremental method. That is, if we take into account that considering the secondary objective function in the optimization does, expectedly, slow down the non-incremental method from the viewpoint of primary objective function decrease in comparison to the corresponding non-bilevel algorithm, it can be seen that incremental bilevel techniques too present a speed up in this convergence rate as the number of subsets grow. Computations were performed on a gtx 480 gpu and timing figures were obtained considering iterate updates only, disregarding both data input/output and objective value computation.

A particular point in the experimental results can be seen in Figure 8. In the application of tomographic reconstruction, incrementalism seems to induce more roughness and we therefore see that inc-$s$ achieves lower Total Variation for the same value of 1-norm of the residual than inc-$2s$. Looking at the iiba-$s$ curves in the same plot, we notice that the choice of secondary objective function, which as we have seen is conflicting with the incrementality idea in the one level case, also plays a similar role in the bilevel case. Consequently, at least for our algorithmic parameters selection, there is a incrementality level (or equivalently, primary objective function decrease speed) versus secondary objective function decrease trade-off. Even so, iiba-$32$ provides substantially lower Total Variation for a given 1-norm of the residual compared to inc-$s$ for every $s$ while still achieving faster experimental primary objective function decrease rate than inc-$s$ with $s\leq 4$.

For other bilevel models, at first glance, there seems not to be any reason for an increase in incrementality to lead to worse secondary to primary objective function ratios. However, such antagonism between fast algorithms (with relation to data adherence) and desirable solution properties may appear naturally in models for ill-posed inverse problems like the one we consider here, because in this case overly fit solutions to noisy data are instable and the secondary objective function is usually an attempt at instability prevention.

4.3 Simulated Data

The present set of experiments intends to highlight the practical differences and advantages of the bilevel approach

[TABLE]

over regularized techniques of the form

[TABLE]

In order to be able to quantify reconstructed image quality, we use simulated data in these experiments. The ideal image will be denoted by $\boldsymbol{x}^{\dagger}$ and was a $512\times 512$ pixels discretized and scaled version of the Shepp-Logan phantom that can be seen at the left of Figure 1. The scaling was such that the relative error after Poisson data simulation was around $10\%$. Tomographic data was computed at $64$ angular samples evenly spaced in $[0,\pi)$, each by its turn sampled in $512$ points in $[-1,1]$.

We have used fiba for solving (46) with just the same parameters (including the secondary operator $\mathcal{O}_{f_{1}}=\mathcal{N}_{f_{\text{Haar}}}$) of those used in Subsection 2.4 except that convergence required $\lambda=2$ and a reasonable starting point for the secondary stepsize sequence was $\mu=10^{2}$. Problem (47) was solved by the fista algorithm with a constant stepsize $\lambda=2$. We have tried $\gamma\in\{10^{2},10,1.5,1,0\}$. We have run both algorithms for $400$ iterations in this experiment and, as in Subsection 2.4, the algorithms were started with the zero image.

Figure 10 shows how the image quality, as measured by the relative error $\|\boldsymbol{x}_{k}-\boldsymbol{x}^{\dagger}\|/\|\boldsymbol{x}^{\dagger}\|$, evolves over the iterations for the tested methods. In Figure 9 we see the best image obtained by each of the methods throughout the iterations. We notice that the bilevel image is competitive with the fista for a certain range of values of $\gamma$, whereas if $\gamma$ is not within the reasonable range, reconstruction by solving (47) quickly degrades. In fact, had we used larger values of the starting secondary stepsize good images would still be obtained, maybe requiring more iterations, while larger than ideal values for $\gamma$ in the non-bilevel approach produce a wholly unusable sequence of iterates. On the other hand, smallish secondary stepsize sequences in the bilevel method have practically the same effect than using a small $\gamma$ in the non-bilevel approach, except that the notion of “small” includes a wider range of values in the bilevel approach because of the diminishing nature of the sequence $\{\mu_{k}\}$. Therefore, if we consider the starting secondary stepsize $\mu$ a parameter of the bilevel technique, it is considerably easier to choose than the parameter $\gamma$ of the traditional regularization approach. On the other hand, if a good and efficient procedure for selecting $\gamma$ is available and the regularization function is appropriate, solving (47) has the potential of delivering better reconstructions.

5 Conclusions

The present paper introduced an abstract class of explicit numerical methods for instances of bilevel non-differentiable convex optimization problems and proposed two first-order concrete representatives of these new algorithms with reduced per iteration computational cost. The proposed methods have computationally simple iterations. The reported numerical experimentation showed that, when used in tomographic reconstruction problems where the problem size is huge, the low computational complexity of the iterations and the advanced modeling allow to the technique to generate high quality images at a moderate computational cost from sparsely sampled and noisy data. Algorithmic flexibility was highlighted by two conceptually distinct implementations: on one side, an implementation in the high-level Python language, which is free and portable to most computing architectures and environments. On the other hand, the simplicity of the method was also suitable for a low-level hardware-specific implementation fully running on a gpu with no external software library dependency.

Acknowledgements

We would like to thank the LNLS for providing the beam time for the tomographic acquisition, obtained under proposal number 17338. We are also grateful to Prof. Marcelo dos Anjos (Federal University of Amazonas) for kindly providing the fish egg samples used in the presented experimentation and Dr. Eduardo X. Miqueles for invaluable help in data acquisition and for the discussions about tomographic reconstruction. We are also indebted to the anonymous referees who gave numerous suggestions that lead to the improvement of the original manuscript.