A convex variational model for learning convolutional image atoms from incomplete data

TL;DR

This paper introduces a convex variational model for learning convolutional image atoms from incomplete or corrupted data, enabling simultaneous image reconstruction and atom learning with proven stability and numerical efficiency.

Contribution

It presents a novel convex variational approach for joint image reconstruction and atom learning, including a semi-convex variant for improved numerical performance.

Findings

Proven well-posedness and stability in inverse problems.

Numerical computation of globally optimal solutions.

Effective handling of incomplete, noisy, and blurry data.

Abstract

A variational model for learning convolutional image atoms from corrupted and/or incomplete data is introduced and analyzed both in function space and numerically. Building on lifting and relaxation strategies, the proposed approach is convex and allows for simultaneous image reconstruction and atom-learning in a general, inverse problems context. Further, motivated by an improved numerical performance, also a semi-convex variant is included in the analysis and the experiments of the paper. For both settings, fundamental analytical properties allowing in particular to ensure well-posedness and stability results for inverse problems are proven in a continuous setting. Exploiting convexity, globally optimal solutions are further computed numerically for applications with incomplete, noisy and blurry data and numerical results are shown.

Figures

Equations

$$\underset{(c_i{)}_i,(p_i{)}_i}{\min }\sum _{i=1}^k\Vert c_i{\Vert }_1\text{s.t. }\{\begin{array}{l}u={\sum }_{i=1}^k{c}_i\ast p_i,\\ \Vert p_i{\Vert }_2\le 1\text{ for }i=1,\dots ,k.\end{array}$$

$$\underset{C}{\min }\Vert C{\Vert }_{1,2}\text{s.t. }u=KC,$$

$$N_{\nu }(u)=\underset{C}{\min }\nu \Vert C{\Vert }_{1,2}+(1-\nu )\Vert C{\Vert }_{\ast }\text{ such that }u=KC,MC=0.$$

$$\underset{u}{\min }\frac{\lambda }{2}\Vert Au-u_0{\Vert }_2^2+N_{\nu }(u),$$

$$\underset{c,p}{\min }\frac{\lambda }{2}\Vert cp-D{\Vert }_2^2+\Vert c{\Vert }_1+R(p),$$

$$\underset{C:\text{rank}(C)\le k}{\min }\frac{\lambda }{2}\Vert C-D{\Vert }_2^2+\Vert C{\Vert }_{1,2},$$

$$\underset{C}{\min }\frac{\lambda }{2}\Vert C-D{\Vert }_2^2+\nu \Vert C{\Vert }_{1,2}+(1-\nu )\Vert C{\Vert }_{\ast },$$

$$\underset{C}{\min }\frac{\lambda }{2}\Vert C-K^{\ast }u{\Vert }_2^2+\nu \Vert C{\Vert }_{1,2}+(1-\nu )\Vert C{\Vert }_{\ast },$$

$$\underset{C}{\min }\frac{\lambda }{2}\Vert KC-u{\Vert }_2^2+\nu \Vert C{\Vert }_{1,2}+(1-\nu )\Vert C{\Vert }_{\ast }.$$

$$\underset{u,v}{\min }\frac{\lambda }{2}\Vert Au-u_0{\Vert }_2^2+{\mu }_{1}R(u-v)+{\mu }_2{N}_{\nu }(v),$$

$${\Omega }_{\Sigma }:=\{x\in {\mathbb{R}}^d|\text{there exists }y\in \Sigma \text{ s.t. }x-y\in \Omega \},$$

$$K(c,p)(x):={\int }_{{\Omega }_{\Sigma }}p(x-y)\mathrm{d}c(y),$$

$$K_{c,p}(h):={\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\tilde{h}(z+y)\tilde{p}(z)\mathrm{d}z\mathrm{d}\tilde{c}(y),$$

$$\left|K_{c,p}(h)\right|\le {\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}|\tilde{h}(z+y)||\tilde{p}(z)|\mathrm{d}z\mathrm{d}|\tilde{c}|(y)\le \Vert \tilde{h}{\Vert }_{L^{q^{\prime }}({\mathbb{R}}^d)}\Vert \tilde{p}{\Vert }_{L^q({\mathbb{R}}^d)}\Vert \tilde{c}{\Vert }_{\mathcal{M}({\mathbb{R}}^d)}\le D\Vert p{\Vert }_s\Vert c{\Vert }_{\mathcal{M}}\Vert h{\Vert }_{q^{\prime }}.$$

$$K_{c,p}(h)={\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\tilde{h}(x)\tilde{p}(x-y)\mathrm{d}x\mathrm{d}\tilde{c}(y)={\int }_{\Omega }h(x)\left({\int }_{{\Omega }_{\Sigma }}p(x-y)\mathrm{d}c(y)\right)\mathrm{d}x.$$

$$K(c,p):=[K_{c,p}]$$

$$N_{\text{cl,s}}(u)=\underset{\begin{array}{c}(c_i{)}_{i=1}^k\subset \mathcal{M}({\Omega }_{\Sigma })\\ (p_i{)}_{i=1}^k\subset L^s(\Sigma )\end{array}}{\inf }\sum _{i=1}^k\Vert c_i{\Vert }_{\mathcal{M}}\text{s.t. }\{\begin{array}{rl}& \Vert p_i{\Vert }_s\le 1,M{p}_i=0i=1,\dots ,k,\\ & u=\sum _{i=1}^{k}K(c_i,p_i)\text{ in }\Omega ,\end{array}$$

$$N_{\text{cl,s}}(u)=\underset{C\in X_s}{\inf }\Vert C{\Vert }_{\pi ,k,M}\text{s.t. }\begin{array}{rl}& u=\hat{K}C\text{ in }\Omega ,\end{array}$$

$$\Vert C{\Vert }_{\pi ,k,M}:=\inf \left\{\sum _{i=1}^k\Vert c_i{\Vert }_{\mathcal{M}}\Vert p_i{\Vert }_s|C=\sum _{i=1}^k{c}_i\otimes p_i\text{ with }M{p}_i=0\text{ for }i=1,\dots ,k\right\}.$$

$$C\mapsto \Vert C{\Vert }_{\pi }+{\mathcal{I}}_{\ker (\hat{M})}(C),$$

$$\Vert C{\Vert }_{\pi }=\inf \left\{\sum _{i=1}^{\infty }\Vert c_i{\Vert }_{\mathcal{M}}\Vert p_i{\Vert }_s|C=\sum _{i=1}^{\infty }{c}_i\otimes p_i\right\}.$$

$$\Vert \cdot {\Vert }_{\pi }+{\mathcal{I}}_{\ker (\hat{M})}\le \Vert \cdot {\Vert }_{\pi ,k,M}\le \Vert \cdot {\Vert }_{\pi ,1,M},$$

$$\Vert \cdot {\Vert }_{\pi }+{\mathcal{I}}_{\ker (\hat{M})}\le {\left(\Vert \cdot {\Vert }_{\pi ,k,M}\right)}^{\ast \ast }\le {\left(\Vert \cdot {\Vert }_{\pi ,1,M}\right)}^{\ast \ast }.$$

$$0=\hat{M}C=\sum _{i=1}^{\infty }{c}_i\otimes M{p}_i=\sum _{i=1}^{\infty }{c}_i\otimes M(p_i-P{p}_i).$$

$$0=\hat{B}(\hat{M}C)=\sum _{i=1}^{\infty }\hat{B}(c_i\otimes M(p_i-P{p}_i))$$

$$=\varphi \left(\sum _{i=1}^{\infty }\psi (c_i)M(p_i-P{p}_i)\right)$$

$$=\varphi \left(M\left(\sum _{i=1}^{\infty }\psi (c_i)(p_i-P{p}_i)\right)\right)$$

$${\left(\Vert \cdot {\Vert }_{\pi ,M}\right)}^{\ast }(B)=\underset{C\in X_s}{\sup }⟨B,C⟩-\Vert C{\Vert }_{\pi ,M}=\underset{n}{\lim }⟨B,C_n⟩-\Vert C_n{\Vert }_{\pi ,M}$$

$$C_n=\sum _{i=1}^{\infty }{c}_i^n\otimes p_i^n\text{and}\sum _{i=1}^{\infty }\Vert c_i^n{\Vert }_{\mathcal{M}}\Vert p_i^n{\Vert }_s\le \Vert C_n{\Vert }_{\pi ,M}+1/n.$$

$${\left(\Vert \cdot {\Vert }_{\pi ,M}\right)}^{\ast }(B)$$

$$=\underset{n}{\lim }\underset{m}{\lim }⟨B,\sum _{i=1}^m{c}_i^n\otimes p_i^n⟩-\sum _{i=1}^m\Vert c_i^n{\Vert }_{\mathcal{M}}\Vert p_i^n{\Vert }_s+1/n$$

$$\le \underset{n}{\lim }\underset{m}{\sup }\underset{\begin{array}{c}(c_i{)}_{i=1}^m,(p_i{)}_{i=1}^m\\ M{p}_i=0\end{array}}{\sup }⟨B,\sum _{i=1}^m{c}_i\otimes p_i⟩-\sum _{i=1}^m\Vert c_i{\Vert }_{\mathcal{M}}\Vert p_i{\Vert }_s+1/n$$

Full text

Mathematical subject classification: 94A08 49M29 65F22 49K30

Keywords: Variational methods, learning approaches, inverse problems, functional lifting, convex relaxation, convolutional lasso, machine learning, texture reconstruction.

1 Introduction

An important task in image processing is to achieve an appropriate regularization or smoothing of images or image related data. In particular, this is indispensable for most application-driven problems in the field, such as denoising, inpainting, reconstruction, segmentation, registration or classification. Also beyond imaging, for general problem settings in the field of inverse problems, an appropriate regularization of unknowns plays a central role as it allows for a stable inversion procedure.

Variational methods and partial-differential-equation-based methods can now be regarded as classical regularization approaches of mathematical image processing (see for instance [42, 3, 40, 30]). An advantage of such methods is the existence of a well-established mathematical theory and, in particular for variational methods, a direct applicability to general inverse problems with provable stability and recovery guarantees [24, 25]. While in particular piecewise smooth images are typically well-described by such methods, their performance for oscillatory or texture-like structures, however, is often limited to pre-described patterns (see for instance [21]).

Data-adaptive methods such as patch-based methods (see for instance [28, 17, 11, 18]) on the other hand are able to exploit redundant structures in images independent of an a-priory description and are, at least for some specific tasks, often superior to variational- and PDE-based methods. In particular machine-learning-based methods have advanced the state-of-the art significantly in many typical imaging applications in the past years. A disadvantage of such methods, however, is a current lack of mathematical understanding in particular compared to existing results for variational methods in the context of inverse problems. Data-adaptive methods typically build on learning some image atoms in some way or the other and neither the mapping of training data to image atoms nor the application of learned atoms to reconstruction is known to enjoy the same stability and regularization properties as with classical methods. In particular, for typical deep neuronal networks, the training step corresponds to the minimization of a non-convex energy, but neither is this minimization carried out until convergence nor are the properties of local minimizers well understood mathematically.

To goal of this work is to provide a first step towards bridging the gap between data-driven methods and variational methods. Building on a tensorial-lifting approach, we introduce a convex variational method for learning image atoms from noisy and/or incomplete data in an inverse problems context. We further extend this model by a semi-convex variant that improves the performance in some applications. For both settings, we are able to prove well-posedness results in function space and, for the convex version, to compute globally optimal solutions numerically. In particular, classical stability and convergence results for inverse problems such as the ones of [24, 25] are applicable to our model, providing a stable recovery of both learned atoms and images from given, incomplete data.

Our approach is motivated by a sparse, convolutional representation of images via a few image atoms and conceptually allows for a joint learning of image atoms and image reconstruction in a single step. Nevertheless, it can also be regarded purely as an a approach for learning image atoms from potentially incomplete data in a training step, after which the learned atoms can be further incorporated in a second step, e.g., for reconstruction or classification. It should also be noted that, while we show some examples where our approach achieves a good practical performance for image reconstruction compared to existing methods, the main purpose of this paper is to provide a mathematical understanding rather than an algorithm that achieves the best performance in practice.

Regarding existing literature in the context of data-adaptive variational learning approaches in imaging, we note that there are many recent approaches that aim to infer either parameter or filters for variational methods from given training data, see e.g. [27, 12, 23]. A continuation of such techniques more towards the architecture of neuronal networks are so called variational networks [26, 2] where not only model parameters but also components of the solution algorithm such as stepsizes or proximal mappings are learned. We also refer to [29] for a recent work on combining variational methods and neuronal networks. While for some of those methods also a function space theory is available, the learning step is still non-convex and one can in general only expect to obtain locally optimal solutions.

1.1 Outline of the paper

In Section 2 we present the main ideas for our approach in a formal setting. This is done from two perspectives, once from the perspective of a convolutional lasso approach and once from the perspective of patch-based methods. In Section 3 we then carry out an analysis of the proposed model in function space, were we motivate our approach via convex relaxation and derive well-posedness results. Section 4 then presents the model in a discrete setting and the numerical solution strategy and Section 5 provides numerical results and a comparison to existing methods. At last, an appendix provides a brief overview on some results for tensor spaces that are used in Section 3. We note that, while the analysis of Section 3 is an important part of our work, the paper is structured in a way such that readers only interested in the conceptual idea and the realization of our approach can skip Section 3 and focus on Sections 2 and 4.

2 A convex approach to image atoms

In this section we present the proposed approach to image-atom-learning and texture reconstruction, where we focus on explaining the main ideas rather than precise definitions of the involved terms. For the latter, we refer to Section 3 for the continuous model and Section 4 for the discrete setting.

Our starting point is the convolutional lasso problem [44, 15], which aims to decompose a given image $u$ as a sparse linear combination of basic atoms $(p_{i})_{i=1}^{k}$ with coefficient images $(c_{i})_{i=1}^{k}$ by inverting a sum of convolutions as follows

[TABLE]

It is important to note that, by choosing the $(c_{i})_{i}$ to be composed of delta peaks, this allows to place the atoms $(p_{i})_{i}$ and any position in the image. In [44], this model has been used in the context of convolutional neural networks for generating image atoms and other image related tasks. Subsequently, many works have dealt with the algorithmic solution of the resulting optimization problem, where the main difficulty lies in the non-convexity of the atom learning step, and we refer to [22] for a recent review.

Our goal is to obtain a convex relaxation of this model that can be used for both, learning image atoms from potentially noisy data as well as image reconstruction tasks such as inpainting, deblurring or denoising. To this aim, we lift the model to the tensor product space of coefficient images and image atoms, i.e., the space of all tensors $C=\sum_{i}c_{i}\otimes p_{i}$ with $c_{i}\otimes p_{i}$ being a rank-1 tensor such that $(c_{i}\otimes p_{i})(x,y)=c_{i}(x)p_{i}(y)$. We refer to Figure 1 for a visualization of this lifting in a one-dimensional setting, where both coefficients and image atoms are vectors and $c_{i}\otimes p_{i}$ corresponds to a rank-one matrix. Notably, in this tensor product space, the convolution $c_{i}*p_{i}$ can be written as linear operator $K$ such that $KC(x)=\sum_{i}K(c_{i}\otimes p_{i})(x)=\sum_{i}\int p_{i}(x-y)c_{i}(y)$. Exploiting redundancies in the penalization of $(\|c_{i}\|_{1})_{i}$ and the constraint $\|p_{i}\|_{2}\leq 1$, $i=1\ldots,k$ and re-writing the above setting in the lifted tensor space, as discussed in Section 3, we obtain the following minimization problem as convex relaxation of the convolutional lasso approach

[TABLE]

where $\|\cdot\|_{1,2}$ takes the $1$-norm and $2$-norm of $C$ in coefficient and atom direction, respectively. Now while a main feature of the original model was that the number of image atoms was fixed, this is no longer the case in the convex relaxation and would correspond to constraining the rank of the lifted variable $C$ (defined as the minimal number of simple tensors needed to decompose $C$) to be below a fixed number. As convex surrogate, we add an additional penalization of the nuclear norm of $C$ in the above objective functional (here we refer to the nuclear norm of $C$ in the tensor product space which, in the discretization of our setting, coincides with the classical nuclear norm of a matrix-reshaping of $C$). Allowing also for additional linear constraints on $C$ via a linear operator $M$, we arrive at the following convex norm that measures the decomposability of a given image $u$ into a sparse combination of atoms as

[TABLE]

Interestingly, this provides a convex model for learning image atoms, which for simple images admitting a sparse representation seems quite effective. In addition, this can in principle also be used as a prior for image reconstruction tasks in the context of inverse problems via solving for example

[TABLE]

with $u_{0}$ given some corrupted data, $A$ a forward operator and $\lambda>0$ a parameter.

Both the original motivation for our model as well as its convex variant have many similarities with dictionary learning and patch-based methods. The next section strives to clarify similarities and difference and provides a rather interesting, different perspective on our model.

2.1 A dictionary-learning/patch-based methods perspective

In classical dictionary-learning-based approaches, the aim is to represent a resorted matrix of image patches as a sparse combination of dictionary atoms. That is, with $u\in\mathbb{R}^{NM}$ a vectorized version of an image and $D=(D_{1},\ldots,D_{l})^{T}\in\mathbb{R}^{l\times nm}$ a patch-matrix containing $l$ vectorized (typically overlapping) images patches of size $nm$, the goal is to obtain a decomposition $D=cd$, where $c\in\mathbb{R}^{l\times k}$ is a coefficient matrix and $p\in\mathbb{R}^{k\times nm}$ is a matrix of $k$ dictionary atoms such that $c_{i,j}$ is the coefficient for the atom $p_{j,\cdot}$ in the representation of the patch $D_{i}$. In order to achieve a decomposition in this form, using only a sparse representation of dictionary atoms, a classical approach is to solve

[TABLE]

where $R$ potentially puts additional constraints or cost on the dictionary atoms, e.g. $R(p)=0$ if $\|p_{j,\cdot}\|_{2}\leq 1$ for all $j$ and $R(p)=\infty$ else.

A difficulty with such an approach is again the bilinear and hence non-convex nature of the optimization problem, leading to potentially many non-optimal stationary points and making the approach sensitive to initialization.

As a remedy, one strategy is to consider a convex variant (see for instance [4]). That is, re-writing the above minimization problem (and again using the ambiguity in the product $cd$ to eliminate the $L^{2}$ constraint) we arrive at the problem

[TABLE]

where $\|C\|_{1,2}=\sum_{i}\|C_{i,\cdot}\|_{2}$. A possible convexification is then given as

[TABLE]

where $\|\cdot\|_{*}$ is the nuclear norm of the matrix $C$.

A disadvantage of such an approach is that the selection of patches is a-priory fixed and that the lifted matrix $C$ has to approximate each patch. In the situation of overlapping patches, this means that individual rows of $C$ have to represent different shifted version of the same patch several times, which inherently contradicts the low-rank assumption.

It is now interesting to see our proposed approach in relation to this dictionary learning methods and the above-described disadvantage: Denote again by $K$ the lifted version of the convolution operator, which in the discrete setting takes a lifted patch-matrix as input and provides an image composted of overlapping patches as output. It is than easy to see that $K^{*}$, the adjoint of $K$, is in fact a patch selection operator and it holds that $KK^{*}=I$. Now using $K^{*}$, the approach in (1) can be re-written as

[TABLE]

where we remember that $u$ is the original image. Considering the problem of finding an optimal patch-based representation of an image as the problem if inverting $K$, we can see that the previous approach in fact first applies a right-inverse of $K$ and then decomposes the data. Taking this perspective, however, it seems much more natural to consider instead an adjoint formulation as

[TABLE]

Indeed, this means that we do not fix the patch-decomposition of the image a-priory but rather allow the method itself to optimally select the position and size of patches. In particular, given a particular patch at an arbitrary location, this patch can be represented by using only one line of $C$ and the other lines (corresponding to shifted versions) can be left empty. Figure 2 shows the resulting improvement by solving both of the above optimization problems for a particular test image, where the parameters are set such that the data error of both methods, i.e., $\|KC-U\|_{2}^{2}$ is approximately the same. As can be seen, solving (2), which we call patch denoising, does not yield meaningful dictionary atoms as the dictionary elements need to represent different, shifted version of the single patch that makes up the image. In contrast to that, solving (3), which we call patch reconstruction, allows to identify the underlying patch of the image and the corresponding patch-matrix is indeed row-sparse.

2.2 The variational model

Now while the proposed model can, in principle, describe any kind of images, in particular its convex relaxation seems best suited for situations where the image can be described by only a few, repeating atoms, as would be for instance the case with texture images. In particular, since we do not include rotations in the model, there are many simple situations, such as $u$ being the characteristic function of a disk, which would in fact require an infinite number of atoms. To overcome this, it seems beneficial to include an additional term which is rotationally invariant and takes care of piecewise smooth parts of the image. Denoting $R$ to be any regularization functional for piecewise smooth data and taking the infimal convolution of this functional with our atom-based norm, we then arrive at the convex model

[TABLE]

for learning image data and image atoms kernels from potentially noisy or incomplete measurements.

A particular example of this model can be given when choosing $R=\operatorname{TV}$, the total variation function [37]. In this setting, a natural choice for the operator $M$ in the definition of $N_{\nu}$ is to take the point-wise mean of the lifted variable in atom direction, which corresponds to constraining the learned atoms to have zero mean and enforces some kind of orthogonality between the cartoon and the texture part in the spirit of [31]. In our numerical experiments, in order to obtain an even richer model for the cartoon part, we use the second order Total Generalized Variation function ($\operatorname{TGV}_{\alpha}^{2})$ [7, 5] as cartoon prior and, in the spirit of a dual $\operatorname{TGV}_{\alpha}^{2}$ norm, use $M$ to constrain the [math]th and $1$st moments of the atoms to be zero.

We also remark that, as shown in the analysis part of Section 3, while an $\ell^{1}/\ell^{2}$-type norm on the lifted variables indeed arises as convex relaxation of the convolutional lasso approach, the addition of the nuclear norm is to some extend arbitrary and in fact in the context of compressed sensing it is known that a summation of two norms is is suboptimal for a joint penalization of sparsity and rank [32] (we refer to Remark 5 below for an extended discussion). Indeed, our numerical experiments also indicate that the performance of our method is to some extend limited by a sub-optimal relaxation of a joint sparsity and rank penalization. To account for that, we also tested with semi-convex potential functions (instead of the identity) for a penalization of the singular values in the nuclear norm. Since this provided a significant improvement in some situations, we also include this more general setting in our analysis and the numerical results.

3 The model in a continuous setting

The goal of this section is to define and analyze the model introduced in Section 2 in a continuous setting. To this aim, we regard images as functions in the Lebesgue spaces $L^{q}(\Omega)$ with a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, $d\in\mathbb{N}$ and $1<q\leq 2$. Image atoms are regarded as functions in $L^{s}(\Sigma)$, with $\Sigma\subset\mathbb{R}^{d}$ a second (smaller) bounded Lipschitz domain (either a circle or a square around the origin) and $s\in[q,\infty]$ an exponent that is a-priori allowed to take any value in $[q,\infty]$, but will be further restricted below. We also refer to the appendix for further notation and results, in particular in the context of tensor product spaces, that will be used in this section.

As described in Section 2 above, the main idea is to synthesize an image via the convolution of a small number of atoms with corresponding coefficient images, where we think of the latter as a sum delta peaks that define the locations where atoms are placed. For this reason, and also due to compactness properties, the coefficient images are modeled as Radon measures in the space $\mathcal{M}(\Omega_{\Sigma})$, the dual of $C_{0}(\Omega_{\Sigma})$, where we denote

[TABLE]

i.e., the extension of $\Omega$ by $\Sigma$. The motivation for using this extension of $\Omega$ is to allow atoms also to be placed arbitrarily close to be boundary (see Figure 1). We will further use the notation $r^{\prime}=r/(r-1)$ for an exponent $r\in(1,\infty)$ and denote duality pairings between $L^{r}$ and $L^{r^{\prime}}$ and between $\mathcal{M}(\Omega)$ and $C_{0}(\Omega)$ by $(\cdot,\cdot)$, while other duality pairings (e.g. between tensor spaces) are denoted by $\langle\cdot,\cdot\rangle$. By $\|\cdot\|_{r},\|\cdot\|_{\mathcal{M}}$ we denote standard $L^{r}$ and Radon norms whenever the domain of definition is clear from the context, otherwise we write $\|\cdot\|_{L^{r}(\Omega_{\Sigma})},\|\cdot\|_{\mathcal{M}(\Omega_{\Sigma})}$ etc.

3.1 The convolutional lasso prior

As first step, we deal with the convolution operator that synthesizes an image from a pair of a coefficient image and an image atom in function space. Formally, we aim to define $K:\mathcal{M}(\Omega_{\Sigma})\times L^{s}(\Sigma)\rightarrow L^{q}(\Omega)$ as

[TABLE]

where we extend $p$ by zero outside of $\Sigma$. An issue with this definition is that, in general, $p$ is only defined Lebesgue almost everywhere and so we have to give a rigorous meaning to the integration of $p$ with respect to an arbitrary Radon measure. To this aim, we define the convolution operator via duality (see [38]). For $c\in\mathcal{M}(\Omega_{\Sigma})$, $p\in L^{s}(\Sigma)$ we define by $K_{c,p}$ the functional on $C(\overline{\Omega})$ as dense subset of $L^{q^{\prime}}(\Omega)$ as

[TABLE]

where $\tilde{g}$ always denotes the zero extension of the function or measure $g$ outside their domain of definition. Now we can estimate with $D>0$

[TABLE]

Hence, by density we can uniquely extend $K_{c,p}$ to a functional in $L^{q^{\prime}}(\Omega)^{*}\simeq L^{q}(\Omega)$ and we denote by $[K_{c,p}]$ the associated function in $L^{q}(\Omega)$. Now in case $p$ is integrable w.r.t. $c$ and $x\mapsto\int_{\Omega_{\Sigma}}p(x-y)\>\mathrm{d}c(y)\in L^{q}(\Omega)$, we get by a change of variables and Fubini’s theorem that for any $h\in C(\overline{\Omega})$

[TABLE]

Hence we get that in this case, $[K_{c,p}](x)=\int_{\Omega_{\Sigma}}p(x-y)\>\mathrm{d}c(y)$ and defining $K:\mathcal{M}(\Omega_{\Sigma})\times L^{s}(\Omega)\rightarrow L^{q}(\Omega)$ as

[TABLE]

we get that $K(c,p)$ coincides with the convolution of $c$ and $p$ whenever the latter is well defined. Note that $K$ is bilinear and, as the previous estimate shows, there exists $D>0$ such that $\|K(c,p)\|_{q}\leq D\|c\|_{\mathcal{M}}\|p\|_{s}$. Hence $K\in\mathcal{B}(\mathcal{M}(\Omega_{\Sigma})\times L^{s}(\Sigma),L^{q}(\Omega))$, the space of bounded bilinear operators (see the appendix).

Using the bilinear operator $K$ and denoting by $k\in\mathbb{N}$ a fixed number for kernels, we now define the convolutional lasso prior for an exponent $s\in[q,\infty]$ and for $u\in L^{q}(\Omega)$ as

[TABLE]

and set $N_{\text{cl,s}}(u)=\infty$ if the constraint set above is empty. Here, we include an operator $M\in\mathcal{L}(L^{s}(\Sigma),\mathbb{R}^{m})$ in our model that optionally allows to enforce additional constraints on the atoms. A simple example of $M$ that we have in mind is an averaging operator, i.e., $Mp:=|\Sigma|^{-1}\int_{\Sigma}p(x)\>\mathrm{d}x$, hence the constraint that $Mp=0$ corresponds to a zero-mean constraint.

3.2 A convex relaxation

Our goal is now to obtain a convex relaxation of the convolutional lasso prior. To this aim, we introduce by $\hat{K}$ and $\hat{M}:=I\otimes M$ the lifting of the bilinear operator $K$ and the linear operators $I$ and $M$, with $I\in\mathcal{L}(\mathcal{M}(\Omega_{\Sigma}),\mathcal{M}(\Omega_{\Sigma}))$ being the identity, to the projective tensor product space $X_{s}:=\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}L^{s}(\Sigma)$ (see the appendix). In this space, we consider a reformulation as

[TABLE]

where

[TABLE]

Note that this reformulation is indeed equivalent. Next we aim to derive the convex relaxation of $N_{\text{cl,s}}$ in this tensor product space. First we consider a relaxation of the functional $\|\cdot\|_{\pi,k,M}$. To this aim, we need an additional assumption on the constraint set $\ker(M)$, which is satisfied for instance if $s=2$ or for $M=0$, in particular will be fulfilled by the concrete setting we use later on.

Lemma 1**.**

Assume that there exists a continuous, linear, norm-one projection onto $\ker(M)$. Then, the convex, lower semi-continuous relaxation of $\|\cdot\|_{\pi,k,M}:X_{s}\rightarrow\overline{\mathbb{R}}$ is given as

[TABLE]

where $\mathcal{I}_{\ker(\hat{M})}(C)=0$ if $\hat{M}C=0$ and $\mathcal{I}_{\ker(\hat{M})}(C)=\infty$ else, and $\|\cdot\|_{\pi}$ is the projective norm on $X_{s}$ given as

[TABLE]

Proof.

At first note that, for a general function $g:X_{s}\rightarrow\overline{\mathbb{R}}$, its convex, lower-semi continuous relaxation can be computed as the biconjugate $g^{**}:X_{s}\rightarrow\overline{\mathbb{R}}$, where $g^{*}(x^{*})=\sup_{x\in X_{s}}\langle x^{*},x\rangle-g(x)$ and $g^{**}(x)=\sup_{x^{*}\in X^{*}_{s}}\langle x^{*},x\rangle-g^{*}(x^{*})$. Keeping this in mind, we note that

[TABLE]

and consequently

[TABLE]

Hence the assertion follows if we show that $\left(\|\cdot\|_{\pi,1,M}\right)^{**}\leq\|\cdot\|_{\pi}+\mathcal{I}_{\ker(\hat{M})}.$ To this aim, we first show that $\|\cdot\|_{\pi,M}\leq\|\cdot\|_{\pi}+\mathcal{I}_{\ker(\hat{M})}$, where we set $\|C\|_{\pi,M}=\|C\|_{\pi,\infty,M}$. Let $C\in X_{s}$ be such that $\hat{M}C=0$ and take $(c_{i})_{i}$, $(p_{i})_{i}$ be such that $\|C\|_{\pi}\geq\sum_{i=1}^{\infty}\|c_{i}\|_{\mathcal{M}}\|p_{i}\|_{s}-\epsilon$ for some $\epsilon>0$ and $C=\sum_{i=1}^{\infty}c_{i}\otimes p_{i}$. Then, with $P$ the projection to $\ker(M)$ as in the assumption, we get that

[TABLE]

Now remember that, according to [39, Theorem 2.9], we have $(\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}\mathbb{R}^{m})^{*}=\mathcal{B}(\mathcal{M}(\Omega_{\Sigma})\times\mathbb{R}^{m})$ with the norm $\|B\|_{\mathcal{B}}:=\sup\{|B(x,y)|\,\left|\right.\,\|x\|_{\mathcal{M}}\leq 1,\,\|y\|_{2}\leq 1\}$. Taking arbitrary $\psi\in\mathcal{M}(\Omega)^{*}$, $\phi\in(\mathbb{R}^{m})^{*}$, we get that $B:(c,p)\mapsto\psi(c)\phi(p)\in\mathcal{B}(\mathcal{M}(\Omega_{\Sigma})\times\mathbb{R}^{m})$ hence

[TABLE]

and since $\phi$ was arbitrary it follows that $M\left(\sum_{i=1}^{\infty}\psi(c_{i})(p_{i}-Pp_{i})\right)=0$. Finally, by closedness of $\operatorname{Rg}(I-P)$ we get that $\sum_{i=1}^{\infty}\psi(c_{i})(p_{i}-Pp_{i})=0$ and, since $\mathcal{M}(\Omega_{\Sigma})$ has the approximation property (see [19, Section VIII.3]), from [39, Proposition 4.6], it follows that $\sum_{i=1}^{\infty}c_{i}\otimes(p_{i}-Pp_{i})=0$, hence $C=\sum_{i=1}^{\infty}c_{i}\otimes Pp_{i}$ and by assumption $\sum_{i=1}^{\infty}\|c_{i}\|_{\mathcal{M}}\|p_{i}\|_{s}\geq\sum_{i=1}^{\infty}\|c_{i}\|_{\mathcal{M}}\|Pp_{i}\|_{s}$. Consequently, $\|\cdot\|_{\pi}+\mathcal{I}_{\ker(\hat{M})}\geq\|\cdot\|_{\pi,M}-\epsilon$ and since $\epsilon$ was arbitrary, the claimed inequality follows.

Now we show that $\left(\|\cdot\|_{\pi,M}\right)^{*}\leq\left(\|\cdot\|_{\pi,1,M}\right)^{*}$, from which the claimed assertion follows by the previous estimate and taking the convex conjugate on both sides. To this aim, take $(C_{n})_{n}\subset X_{s}$ such that

[TABLE]

and take $(c_{i}^{n})_{i}$, $(p_{i}^{n})_{i}$ such that $Mp^{n}_{i}=0$ for all $n,i$ and

[TABLE]

We then get

[TABLE]

Now it can be easily seen that the last expression equals [math] in case $|B(c,p)|\leq\|c\|_{\mathcal{M}}\|p\|_{s}$ for all $c,p$ with $Mp=0$. In the other case, we can pick $\hat{c},\hat{p}$ with $M\hat{p}=0$ and $r>1$ such that $B(\hat{c},\hat{p})>r\|\hat{c}\|_{\mathcal{M}}\|\hat{p}\|_{s}$ and get for any $\lambda>0$ that

[TABLE]

Hence, the last line of the above equation is either 0 or infinity and equals

[TABLE]

This result suggests that the convex, lower semi-continuous relaxation of (5) will be obtained by replacing $\|\cdot\|_{\pi,k,M}$ with the projective tensor norm $\|\cdot\|_{\pi}$ on $X_{s}$ and the constraint $\hat{M}C=0$. Our approach to show this will in particular require us ensure lower semi-continuity of this candidate for the relaxation, which in turn requires us to ensure a compactness property of the sublevelsets of the energy appearing in (5) and closedness of the constraints. To this aim, we consider a weak* topology on $X_{s}$ and rely on a duality result for tensor product spaces (see the appendix), which states that, under some conditions, the projective tensor product $X_{s}=\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}L^{s}(\Sigma)$ can be identified with the dual of the so called injective tensor product $C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{s^{\prime}}(\Sigma)$. Different from one what one would expect from the individual spaces, however, this can only be ensured for the case $s<\infty$ which excludes the space $L^{\infty}(\Sigma)$ for the image atoms. This restriction will also be required later on in order to show well-posedness of a resulting regularization approach for inverse problems, hence we will henceforth always consider the case that $s\in[q,\infty)$ and use the identification $(C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{s^{\prime}}(\Sigma))^{*}\hat{=}\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}L^{s}(\Sigma)$ (see the appendix).

As first step towards the final relaxation result, and also as crucial ingredient for well-posedness results below, we show weak* continuity of the operator $\hat{K}$ on the space $X_{s}$.

Lemma 2**.**

Let $s\in[q,\infty)$. Then the operator $\hat{K}:X_{s}\rightarrow L^{q}(\Omega)$ is continuous w.r.t. weak convergence in $X_{s}$ and weak convergence in $L^{q^{\prime}}(\Omega)$. Also, for any $\phi\in C_{0}(\Omega)\subset L^{q^{\prime}}(\Omega)$ it follows that $\hat{K}^{*}\phi\in C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{s^{\prime}}(\Sigma)$ and, via the identification $C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{s^{\prime}}(\Sigma)\hat{=}C_{0}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$ (see the appendix), can be given as $K^{*}\phi(t)=[x\mapsto\phi(t+x)]$.*

Proof.

First we note that for any $\psi\in C_{c}(\Omega)$, the function $\hat{\psi}$ defined as $\hat{\psi}(t)=[x\mapsto\psi(t+x)]$ (where we extend $\psi$ by 0 to $\mathbb{R}^{d}$) is contained in $C_{c}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$. Indeed, continuity follows since by uniform continuity for any $\epsilon>0$ there exists a $\delta>0$ such that for any $r\in\mathbb{R}^{d}$ with $|r|\leq\delta$ and $t\in\Omega_{\Sigma}$ with $t+r\in\Omega_{\Sigma}$

[TABLE]

Also, taking $K\subset\subset\Omega$ to be the support of $\psi$ we get, with $K_{\Sigma}$ the extension of $K$ by $\Sigma$, for any $t\in\Omega_{\Sigma}\setminus K_{\Sigma}$ that $\psi(t+x)=0$ for any $x\in\Sigma$ and hence $\hat{\psi}=0$ in $L^{s^{\prime}}(\Sigma)$ and $\hat{\psi}\in C_{c}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$.

Now for $\phi\in C_{0}(\Omega_{\Sigma})$, taking $(\phi_{n})_{n}\in C_{c}(\Omega_{\Sigma})$ to be a sequence converging to $\phi$, we get that

[TABLE]

Thus $\hat{\phi}$ can be approximated by a sequence of compactly supported functions and hence $\hat{\phi}\in C_{0}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$. Fixing now $u=c\otimes p\in X_{s}$ we note that for any $\psi\in C_{0}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$, the function $t\rightarrow\int_{\Sigma}\psi(t)(x)p(x)\>\mathrm{d}x$ is continuous, hence we can define the linear functional

[TABLE]

and get that $F_{u}$ is continuous on $C_{0}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$. Then, since $\hat{\phi}\in C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{s^{\prime}}(\Sigma)$ it can be approximated by a sequence of simple tensors $(\sum_{i=1}^{m_{n}}x_{i}^{n}\otimes y_{i}^{n})_{n}$ in the injective norm, which coincides with the norm in $C_{0}(\Omega_{\Sigma},L^{s^{\prime}}(\Sigma))$ and, using Lemma 21 in the appendix, we get

[TABLE]

Now by density of simple tensors in the projective tensor product, it follows that $K^{*}\phi=\hat{\phi}$. In order to show the continuity assertion, take $(u_{n})_{n}$ weak * converging to some $u\in X_{s}$. Then by the previous assertion we get for any $\phi\in C_{c}(\Omega_{\Sigma})$ that

[TABLE]

hence $(Ku_{n})_{n}$ weakly converges to $Ku$ on a dense subset of $L^{s^{\prime}}(\Omega)$ which, together with boundedness of $(Ku_{n})_{n}$, implies weak convergence. ∎

We will also need weak*-to-weak* continuity of $\hat{M}$, which is shown in the following lemma in a slightly more general situation than needed.

Lemma 3**.**

Take $s\in[q,\infty)$ and assume that $M\in\mathcal{L}(L^{s}(\Sigma),Z)$ with $Z$ a reflexive space and define $\hat{M}:=I\otimes_{\pi}M\in\mathcal{L}(X_{s},\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}Z)$, where $I$ is the identity on $\mathcal{M}(\Omega_{\Sigma})$. Then $\hat{M}$ is continuous w.r.t. weak star convergence in both spaces.

Proof.

Take $(u_{n})_{n}\in X$ weak* converging to some $u\in X$ and write $u_{n}=\lim_{k}\sum_{i=1}^{k}x_{i}^{n}\otimes y_{i}^{n}$. We note that, since $Z$ is reflexive, it satisfies in particular the Radon Nikodým property (see the appendix) and hence $C(\Omega_{\Sigma})\otimes_{\mathfrak{i}}Z^{*}$ can be regarded as predual of $\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}Z$ and we test with $\phi\otimes\psi\in C(\Omega_{\Sigma})\otimes_{\mathfrak{i}}Z^{*}$. Then

[TABLE]

where the convergence follows since $M^{*}\psi\in L^{s^{\prime}}(\Sigma)$, the predual of $L^{s}(\Sigma)$, and hence $\phi\otimes M^{*}\psi\in C(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{q^{\prime}}(\Sigma).$ ∎

Now we can obtain the convex, lower semi-continuous relaxation of $N_{\text{cl,s}}$.

Lemma 4**.**

With the assumptions of Lemma 1 and $s\in[q,\infty)$, the convex, l.s.c. relaxation of $N_{\text{cl,s}}$ is given as

[TABLE]

Proof.

Again we first compute the convex conjugate:

[TABLE]

Similarly, we see that $N^{*}(v)=\left(\|\cdot\|_{\pi}+\mathcal{I}_{\ker(\hat{M})}\right)^{*}(\hat{K}^{*}v)$, where

[TABLE]

Now in the proof of Lemma 1, we have in particularly shown that $\left(\|\cdot\|_{\pi}+\mathcal{I}_{\ker(\hat{M})}\right)^{*}=\|\cdot\|_{\pi,k,M}^{*}$, hence, if we show that $N$ is convex and lower semi-continuous, the assertion follows from $N(u)=N^{**}(u)=N_{\text{cl,s}}^{**}(u)$. To this aim, take a sequence $(u_{n})_{n}$ in $L^{q}(\Omega)$ converging weakly to some $u$ for which, without loss of generality, we assume that

[TABLE]

Now with $(C_{n})_{n}$ such that $\|C_{n}\|_{\pi}\leq N(u_{n})+n^{-1}$, $\hat{M}C_{n}=0$ and $u_{n}=\hat{K}C_{n}$ we get that $(\|C_{n}\|_{\pi})_{n}$ is bounded. Since $X_{s}$ admits a separable predual (see the appendix), this implies that $(C_{n})_{n}$ admits a subsequence $(C_{n_{i}})_{i}$ weak* converging to some $C$. By weak* continuity of $\hat{K}$ and $\hat{M}$ we get that $u=\hat{K}C$ and $\hat{M}C=0$, respectively, and by weak* lower semi-continuity of $\|\cdot\|_{\pi}$ it follows that

[TABLE]

which concludes the proof. ∎

This relaxation results suggest to use $N(\cdot)$ as in Equation (6) as convex texture prior in the continuous setting. There is, however, an issue with that, namely that such a functional cannot be expected to penalize the number of used atoms at all. Indeed, taking some $C=\sum_{i=1}^{l}c_{i}\otimes p_{i}$ and assume that $\|C\|_{\pi}=\sum_{i=1}^{l}\|c_{i}\|_{\mathcal{M}}\|p_{i}\|_{s}$. Now note that we can split any summand $c_{i_{0}}\otimes p_{i_{0}}$ as follows: Write $c_{i_{0}}=c^{1}_{i_{0}}+c^{2}_{i_{0}}$ with disjoint support such that $\|c_{i_{0}}\|_{\mathcal{M}}=\|c^{1}_{i_{0}}\|_{\mathcal{M}}+\|c^{2}_{i_{0}}\|_{\mathcal{M}}$. Then we can re-write

[TABLE]

which gives a different representation of $C$ by increasing the number of atoms without changing the cost of the projective norm. Hence, in order to maintain the original motivation of the approach to enforce a limited number of atoms, we need to add an additional penalty on $C$ for the lifted texture prior.

3.3 Adding a rank penalization

Considering the discrete setting and the representation of the tensor $C$ as a matrix, the number of used atoms corresponds to the rank of the matrix, for which it is well known that the nuclear norm constitutes a convex relaxation [20]. This construction can in principle also be transferred to general tensor products of Banach spaces via the identification (see Proposition 22 in the appendix)

[TABLE]

and the norm

[TABLE]

It is important to realize, however, that the nuclear norm of operators depends on the underlying spaces and in fact coincides with the projective norm in the tensor product space (see Proposition 22). Hence adding the nuclear norm in $X_{s}$ does not change anything and, more generally, whenever one of the underlying spaces is equipped with and $L^{1}$-type norm we cannot expect a rank-penalizing effect (consider the example of the previous section).

On the other hand, going back to the nuclear norm of a matrix in the discrete setting, we see that it relies on orthogonality and an inner product structure and that the underlying norm is the Euclidean inner product norm. Hence an appropriate generalization of a rank penalizing nuclear norm needs to be built on a Hilbert space setting. Indeed, it is easy to see that any operator between Banach spaces with a finite nuclear norm is compact, and in particular for any $T\in\mathcal{L}(H_{1},H_{2})$ with finite nuclear norm and $H_{1},H_{1}$ Hilbert spaces, there are orthonormal systems $(x_{i})_{i}$, $(y_{i})_{i}$ and uniquely defined singular values $(\sigma_{i})_{i}$ such that

[TABLE]

Motivated by this, we aim to define an $L^{2}$-based nuclear norm as extended real valued function on $X_{s}$ as convex surrogate of a rank penalization. To this aim, we consider from now on the case $s=2$. Remember that the tensor product $X\otimes Y$ of two spaces $X,Y$ is defined as the vector space spanned by linear mappings $x\otimes y$ on the space of bilinear forms on $X\times Y$, which are given as $x\otimes y(B)=B(x,y)$. Now since $L^{2}(\Omega_{\Sigma})$ can be regarded as subspace of $\mathcal{M}(\Omega_{\Sigma})$, also $L^{2}(\Omega_{\Sigma})\otimes L^{2}(\Sigma)$ can be regarded as subspace of $\mathcal{M}(\Omega_{\Sigma})\otimes L^{2}(\Sigma).$ Further, defining for $C\in L^{2}(\Omega_{\Sigma})\otimes L^{2}(\Sigma)$,

[TABLE]

we get that $\|\cdot\|_{\pi}\leq B\|\cdot\|_{\pi,L^{2}\otimes L^{2}}$ for a constant $B>0$, hence also the completion $L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$ can be regarded as subspace of $\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$. Further, $L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$ can be identified with the space of nuclear operators $\mathcal{N}(L^{2}(\Omega_{\Sigma}),L^{2}(\Sigma))$ as above such that

[TABLE]

Using this, and introducing a potential function $\phi:[0,\infty)\rightarrow[0,\infty)$, we define for $C\in X_{2}$,

[TABLE]

We will mostly focus on the case $\phi(x)=x$, in which $\|\cdot\|_{\text{nuc},\phi}$ coincides with an extension of the nuclear norm and can be interpreted as convex relaxation of the rank. However, since we observed a significant improvement in some cases in practice by choosing $\phi$ to be a semi-convex potential function, i.e., a function such that $\phi+\tau|\cdot|^{2}$ is convex for $\tau$ sufficiently small, we include the more general situation in the theory.

Remark 5** (Sparsity and low-rank).**

It is important to note that $\|C\|_{\text{nuc},\phi}<\infty$ restricts $C$ to be contained in the smoother space $L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$ and in particular does not allow for simple tensors $\sum_{i=1}^{k}c_{i}\otimes p_{i}$ with the $c_{i}$’s being composed of delta peaks. Thus we observe some inconsistency of a rank penalization via the nuclear norm and a pointwise sparsity penalty, which is only visible in the continuous setting trough the regularity of functions. Nevertheless, such an inconsistency has already been observed in the finite dimensional setting in the context of compressed sensing for low-rank AND sparse matrices, manifested via a poor performance of the sum of a nuclear norm and $\ell^{1}$ norm for exact recovery (see [32]). As a result, there exists many literature on improved, convex priors for the recovery of low-rank and sparse matrices, see for instance [36, 16, 35]. While such improved priors can be expected to be highly beneficial for our setting, the question does not seem to be solved in such a way that that can be readily applied in our setting.

One direct way to circumvent this inconsistency would be to include an additional smoothing operator for $C$ as follows: Take $S\in\mathcal{L}(\mathcal{M}(\Omega_{\Sigma}),\mathcal{M}(\Omega_{\Sigma}))$ such that $\text{range}(S)\subset L^{2}(\Omega_{\Sigma})$ to be a weak to weak* continuous linear operator and define the operator $\hat{S}:X_{2}\rightarrow X_{2}$ as $\hat{S}:=S\otimes I_{L^{2}}$, were $I_{L^{2}}$ denotes the identity in $L^{2}(\Sigma)$. Then one could alternatively also use $\|SC\|_{\text{nuc}}$ as alternative for penalizing the rank of $C$ while still allowing $C$ to be a general measure. Indeed, in the discrete setting, by choosing $S$ also to be injective, we even obtain the equality $\text{rank}(SC)=\text{rank}(C)$ (where we interpret $C$ and $SC$ as matrices). In practice, however, we did not observe an improvement by including such a smoothing and thus do not include $\hat{S}$ in our model.*

3.4 Well-posedness and a cartoon-texture model

Including $\|C\|_{\text{nuc},\phi}$ for $C\in X_{2}$ as additional penalty in our model, we ultimately arrive at the following variational texture prior in the tensor product space $X_{2}:=\mathcal{M}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$, which is convex whenever $\phi$ is convex, in particular for $\phi(x)=|x|$.

[TABLE]

where $\nu\in(0,1)$ is a parameter balancing the sparsity and the rank penalty.

In order to employ $N_{\nu}$ as a regularization term in an inverse problems setting, we need to obtain some lower-semi continuity and coercivity properties. As first step, the following lemma, which is partially inspired by techniques used in [8, Lemma 3.2], shows that, under some weak conditions on $\phi$, $\|\cdot\|_{\text{nuc},\phi}$ defines a weak* lower semi-continuous function on $X_{2}$.

Lemma 6**.**

Assume that $\phi:[0,\infty)\rightarrow[0,\infty)$ is lower semi-continuous, non-decreasing, that

• •

$\phi(x)\rightarrow\infty$* for $x\rightarrow\infty$ and that*

• •

there exist $\epsilon,\eta>0$ such that $\phi(x)\geq\eta x$ for $0\leq x<\epsilon$.

Then the functional $\|\cdot\|_{\text{nuc},\phi}:X_{2}\rightarrow\overline{R}$ defined as in (7) is lower semi-continuous w.r.t. weak convergence in $X_{2}$.*

Proof.

Take $(C_{n})_{n}\subset X_{2}$ weak* converging to some $C\in X_{2}$ for which, w.l.o.g., we assume that

[TABLE]

We only need to consider the case that $(\|C_{n}\|_{\text{nuc},\phi})_{n}$ is bounded, otherwise the assertion follows trivially. Hence we can write $T_{C_{n}}(x)=\sum_{i=1}^{\infty}\sigma_{i}^{n}(x_{i}^{n},x)y_{i}^{n}$ such that $\|C_{n}\|_{\text{nuc},\phi}=\sum_{i=1}^{\infty}\phi(\sigma_{i}^{n})$. Now we aim to bound $(\|C_{n}\|_{\pi,L^{2}\otimes L^{2}})_{n}$ in terms of $(\|C_{n}\|_{\text{nuc},\phi})_{n}$. To this aim, first note that the assumptions in $\phi$ imply that for any $\epsilon^{\prime}>0$ there is $\eta^{\prime}>0$ such that $\phi(x)\geq\eta^{\prime}x$ for all $x<\epsilon^{\prime}$. Also, $\phi(\sigma_{i}^{n})\leq\|C_{n}\|_{\text{nuc},\phi}$ for any $i,n$ and via a direct contradiction argument it follows that there exists $\hat{\epsilon}>0$ such that $\sigma_{i}^{n}<\hat{\epsilon}$ for all $i,n$. Picking $\hat{\eta}$ such that $\phi(x)\geq\hat{\eta}x$ for all $x<\hat{\epsilon}$ we obtain

[TABLE]

hence $(C_{n})_{n}$ is also bounded as sequence in $L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$ and admits a (non-relabeled) subsequence weak* converging to some $\hat{C}\in L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$, with $L^{2}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{2}(\Sigma)$ being the predual space. By the inclusion $C_{0}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{2}(\Sigma)\subset L^{2}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{2}(\Sigma)$ and uniqueness of the weak* limit, we finally get $\hat{C}=C\in L^{2}(\Omega_{\Sigma})\otimes_{\mathfrak{i}}L^{2}(\Sigma)$ and can write $T_{C}x=\sum_{i=1}^{\infty}\sigma_{i}(x_{i},x)y_{i}$ and $\|C\|_{\text{nuc},\phi}=\sum_{i=1}^{\infty}\sigma_{i}$. By lower semi-continuity of $\|\cdot\|_{\text{nuc}}$ this would suffice to conclude in the case $\phi(x)=x$. For the more general case, we need to show a point-wise lim-inf property of the singular values. To this aim, note that by the Courant-Fischer min-max principle (see for instance [10, Problem 37]) for any compact operator $T\in\mathcal{L}(H_{1},H_{2})$ with $H_{1},H_{2}$ Hilbert spaces and $\lambda_{k}$ the $k$-th singular value of $T$ sorted in descending order we have

[TABLE]

Now consider $k\in\mathbb{N}$ fixed. For any subspace $V$ with $\dim(V)=k$, the minimum in the equation above is achieved, hence we can denote $x_{V}$ to be a minimizer and define $F_{V}(T):=\|Tx_{V}\|_{H_{2}}$ such that $\lambda_{k}=\sup_{\dim(V)=k}F_{V}(T)$. Since weak* convergence of a sequence $(T_{n})$ to $T$ in $L^{2}(\Omega_{\Sigma})\otimes_{\pi}L^{2}(\Sigma)$ implies in particular $T_{n}(x)\rightharpoonup T(x)$ for all $x$, by lower semi-continuity of the norm $\|\cdot\|_{H_{2}}$ it follows that $F_{V}$ is lower semi-continuous with respect to weak* convergence. Hence this is also true for the function $T\mapsto\lambda_{k}(T)$ is being the pointwise supremum of a family of lower semi-continuous functional. Consequently, for the sequence $(T_{C_{n}})_{n}$ it follows that $\sigma_{k}\leq\liminf_{n}\sigma^{n}_{k}$. Finally, by monotonicity and lower semi-continuity of $\phi$ and Fatou’s lemma we conclude

[TABLE]

∎

The lemma below now establishes the main properties of $N_{\nu}$ that in particular allow to employ it as regularization term in an inverse problems setting.

Lemma 7**.**

The infimum in the definition of (8) is attained and $N_{\nu}:L^{q}(\Omega)\rightarrow\overline{\mathbb{R}}$ is convex and lower semi-continuous. Further, any sequence $(v_{n})_{n}$ such that $N_{\nu}(v_{n})$ is bounded admits a subsequence converging weakly in $L^{q}(\Omega)$.

Proof.

The proof is quite standard, but we provide it for the readers convenience. Take $(v_{n})_{n}$ to be a sequence such that $(N_{\nu}(v_{n}))_{n}$ is bounded. Then we can pick a sequence $(C_{n})_{n}$ in $X_{2}$ such that $\hat{M}C_{n}=0$, $v_{n}=\hat{K}C_{n}$ and

[TABLE]

This implies that $(C_{n})_{n}$ admits a subsequence $(C_{n_{i}})_{i}$ weak* converging to some $C\in X_{2}$. Now by continuity of $\hat{M}$ and $\hat{K}$ we have that $\hat{M}C=0$ and that $(v_{n_{i}})_{i}=(\hat{K}C_{n_{i}})_{i}$ is bounded. Hence also $(v_{n_{i}})_{i}$ admits a (non-relabeled) subsequence converging weakly to some $v=\hat{K}C$. This already shows the last assertion. In order to show lower semi-continuity, assume that $(v_{n})_{n}$ converges to some $v$ and, without loss of generality, that

[TABLE]

Now this is a particular case of the argumentation above, hence we can deduce with $(C_{n})_{n}$ as above that

[TABLE]

which implies lower semi-continuity. Finally, specializing even more to the case that $(v_{n})_{n}$ is the constant sequence $(v)_{n}$, also the claimed existence follows. ∎

In order to model a large class of natural images and to keep the number of atoms needed in the above texture prior low, we combine it with a second part that models cartoon-like images. Doing so, we arrive at the following model

[TABLE]

where we assume $R$ to be a functional that models cartoon images, $D(\cdot,f_{0}):Y\rightarrow\overline{\mathbb{R}}$ is a given data discrepancy, $A\in\mathcal{L}(L^{q}(\Omega),Y)$ a forward model and we define the parameter balancing function

[TABLE]

Now we get the following general existence result.

Proposition 8**.**

Assume that $R:L^{q}(\Omega)\rightarrow\overline{R}$ is convex, lower semi-continuous and that there exists a finite dimensional subspace $U\subset L^{q}(\Omega)$ such that for any $u\in L^{q}(\Omega)$, $v\in U^{\perp}$, $w\in U$,

[TABLE]

with $B>0$ and $U^{\perp}$ denoting the complement of $U$ in $L^{q}(\Omega)$. Further assume that $A\in\mathcal{L}(L^{q}(\Omega),Y)$, $D(\cdot,f_{0})$ is convex, lower semi-continuous and coercive on the finite dimensional space $A(U)$ in the sense that for any two sequences $(u_{n}^{1})_{n}$, $(u_{n}^{2})_{n}$ such that $(u_{n}^{1})_{n}\subset U$, $(u_{n}^{2})_{n}$ is bounded and $(D(A(u_{n}^{1}+u_{n}^{2}),f_{0}))_{n}$ is bounded, also $(\|Au^{1}_{n}\|_{q})_{n}$ is bounded. Then there exists a solution to (P).

Remark 9**.**

Note that for instance in case $D$ satisfies a triangle inequality, the sequence $(u^{2}_{n})$ in the coercivity assumption is not needed, i.e., can be chosen to be zero.

Proof.

The proof is rather standard and we provide only a short sketch. Take $((u_{n},v_{n}))_{n}$ a minimizing sequence for (P). From Lemma 7 we get that $(v_{n})_{n}$ admits a (non-relabeled) weakly convergent subsequence. Now we split $u_{n}=u_{n}^{1}+u_{n}^{2}\in U+U^{\perp}$ and $v_{n}=v_{n}^{1}+v_{n}^{2}\in U+U^{\perp}$ and by assumption get that $\|u^{2}_{n}-v^{2}_{n}\|_{q}$ is bounded. But since $(\|v_{n}\|_{q})_{n}$ is bounded, so is $(\|v^{2}_{n}\|_{q})_{n}$ and consequently also $(\|u^{2}_{n}\|_{q})_{n}$. Now we split again $u_{n}^{1}=u_{n}^{1,1}+u_{n}^{1,2}\in\ker(A)\cap U+(\ker(A)\cap U)^{\perp}$ and note that also $(u_{n}^{1,2}+u_{n}^{2},v_{n})$ is a minimizing sequence for (P). Hence it remains to show that $(u_{n}^{1,2})_{n}$ is bounded in order to get a bounded minimizing sequence. To this aim, we note that $(u_{n}^{1,2})_{n}\subset(\ker(A)\cap U)^{\perp}\cap U$ and that $A$ is injective on this finite dimensional space. Hence $\|u_{n}^{1,2}\|_{q}\leq B\|Au_{n}^{1,2}\|_{q}$ for some $B>0$ and by the coercivity assumption on the data term we finally get that $(\|u_{n}^{1,2}\|_{q})_{n}$ is bounded. Hence also $(u_{n}^{1,2}+u_{n}^{2})_{n}$ admits a weakly convergent subsequence in $L^{q}(\Omega)$ and by continuity of $A$ as well as lower semi-continuity of all involved functionals existence of a solution follows. ∎

Remark 10** (Choice of regularization).**

A particular choice of regularization for $R$ in (P) that we consider in this paper is $R=\operatorname{TGV}_{\alpha}^{2}$, with $\operatorname{TGV}_{\alpha}^{2}$ the second order total generalized variation functional [7], $q\leq d/(d-1)$ and

[TABLE]

Since in this case [9, 5]

[TABLE]

with $B>0$ and for all $u\in\mathcal{P}_{1}(\Omega)^{\perp}$, the complement of the first order polynomials, and $\operatorname{TGV}_{\alpha}^{2}$ is invariant on first order polynomials, the result of Proposition 8 applies.

Remark 11** (Norm-type data terms).**

We also note that the result of Proposition 8 in particular applies to $D(w,f_{0}):=\frac{1}{r}\|w-f_{0}\|_{r}^{r}$ for any $r\in[1,\infty)$ or $D(w,f_{0}):=\|w-f\|_{\infty}$, where we extend the norms by infinity to $L^{q}(\Omega)$ whenever necessary. Indeed, lower semi-continuity of these norms is immediate for both $r\leq q$ and $r>q$ and since the coercivity is only required on a finite dimensional space, it also holds by equivalence of norms.

Remark 12** (Inpainting).**

At last we also remark that the assumptions of Proposition 8 also hold for an inpainting data term defined as

[TABLE]

whenever $\omega$ has nonempty interior. Indeed, lower semi-continuity follows from the fact that $L^{q}$ convergent sequences admit pointwise convergent subsequences and the coercivity follows from finite dimensionality of $U$ and the fact that $\omega$ has non-empty interior.

Remark 13** (Regularization in a general setting).**

We also note that Lemma 7 provides the basis for employing either $N_{\nu}$ directly or its infimal-convolution with a suitable cartoon prior as in Proposition 8 for the regularization of general (potentially non-linear) inverse problems and with multiple data fidelities, see for instance [24, 25] for general results in that direction.

4 The model in a discrete setting

This section deals with the discretization of the proposed model and its numerical solution. For the sake of brevity, we provide only the main steps and refer to the publicly available source code [13] for all details.

We define $U=\mathbb{R}^{N\times M}$ to be the space of discrete grayscale images, $W=\mathbb{R}^{(N+n-1)\times(M+n-1)}$ to be the space of coefficient images and $Z=\mathbb{R}^{n\times n}$ to be the space of image atoms for which we assume $n<\min\{N,M\}$ and, for simplicity, only consider a square domain for the atoms. The tensor product of a coefficient image $c\in W$ and a atom $p\in Z$ is given as $(c\otimes p)_{i,j,r,s}=c_{i,j}p_{r,s}$ and the lifted tensor space is given as the four-dimensional space $X=\mathbb{R}^{(N+n-1)\times(M+n-1)\times n\times n}$.

Texture norm. The forward operator $K$ being the lifting of the convolution $c\ast p$ and mapping lifted matrices to the vectorized image space is then given as

[TABLE]

and we refer Figure 1 for a visualization in the one dimensional case. Note that by extending the first two dimensions of the tensor space to $N+n-1$, $M+n-1$ we allow to place an atom at any position where it still effects the image, also partially outside the image boundary.

Also we note that, in order to reduce dimensionality and accelerate the computation, we introduce a stride parameter $\eta\in\mathbb{N}$ in practice which introduces a stride on the possible atom positions. That is, the lifted tensor space and forward operator is reduced in such a way that the grid of possible atom positions in the image is essentially $\{(\eta i,\eta j)\,\left|\right.\,i,j\in\mathbb{N},(\eta i,\eta j)\in\{1,\ldots,N\}\times\{1,\ldots,M\}\}$. This reduces the dimension of the tensor space by a factor $\eta^{-2}$, while for $\eta>1$ it naturally does not allow for arbitrary atom positions anymore and for $\eta=n$ it corresponds to only allowing non-overlapping atoms. In order to allow for atoms being placed next to each other it is important to choose $\eta$ to be a divisor of the atom-domain size $n$ and we used $n=15$ and $\eta=3$ in all experiments of the paper. In order to avoid extensive indexing and case distinctions, however, we only consider the case $\eta=1$ here and refer to the source code [13] for the general case.

As a straightforward computation shows that, in the discrete lifted tensor space, the projective norm corresponding to discrete $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$ norms for the coefficient images and atoms, respectively, is given as a mixed 1-2 norm as

[TABLE]

The nuclear norm for a potential $\phi$ on the other hand reduces to the evaluation of $\phi$ on the singular values of a matrix-reshaping of the lifted tensors and is given as

[TABLE]

where $[A_{(NM,nn)}]$ denotes a reshaping of the tensor $A$ to a matrix of dimensions $NM\times nn$. For the potential function $\phi$ we consider two choices: Mostly we are interested in $\phi(x)=x$ which yields a convex texture model and enforces sparsity of the singular values. A second choice we consider is $\phi:[0,\infty)\rightarrow[0,\infty)$ given as

x$$\phi(x)

$\phi(x)=\begin{cases}x-\epsilon\delta x^{2}&x\in[0,\frac{1}{2\epsilon}]\\ (1-\delta)x+\frac{\delta}{4\epsilon}&\text{else,}\end{cases}$

(10)

where $\delta<1$, $\delta\approx 1$ and $\epsilon>0$. It is easy to see that $\phi$ fulfills the assumptions of Lemma 6 and that $\phi$ is semi-convex, i.e., $\phi+\rho|\cdot|^{2}$ is convex for $\rho>\delta\epsilon$. While the results of Section 3 hold for this setting even without the semi-convexity assumption, we can in general not expect to obtain an algorithm that provably delivers a globally optimal solution in the semi-convex (or generally non-convex) case. The reason for using a semi-convex potential rather than a arbitrary non-convex one is twofold: First, for a suitably small stepsize $\tau$ the proximal mapping

[TABLE]

is well defined and hence proximal-point type algorithms are applicable at least conceptually. Second, since we employ $\phi$ on the singular values of the lifted matrices $A$, it will be important for numerical feasibility of the algorithm that the corresponding proximal mapping on $A$ can be reduced to a proximal mapping on the singular values. While this is not obvious for a general choice of $\phi$, it is true (see Lemma 14 blow) for semi-convex $\phi$ with suitable parameter choices.

Cartoon prior As cartoon-prior we employ the second-order total generalized variation functional which we define for fixed parameters $(\alpha_{0},\alpha_{1})=(\sqrt{2},1)$ and a discrete image $u\in U$ as

[TABLE]

Here $\nabla$ and $\mathcal{E}$ denote discretized gradient and symmetrized Jacobian operators, respectively and we refer to [6] and the source code [13] for details on a discretization of $\operatorname{TGV}_{\alpha}^{2}$. To ensure a certain orthogonality of the cartoon and texture part, we further define the operator $M$ that incorporates atom-constraints, to evaluate the 0th and 1st moments of the atoms, which in the lifted setting yields

[TABLE]

The discrete version of (P) is then given as

[TABLE]

where the parameter balancing functions $s_{1},s_{2}$ are given as in (9) and the model depends on three parameters $\lambda,\mu,\nu$, with $\lambda$ defining the trade-off between data and regularization, $\mu$ defining the trade-off between the cartoon and the texture part and $\nu$ defining the trade-off between sparsity and low-rank of the tensor $C$.

Numerical solution. For the numerical solution of (DP) we employ the primal-dual algorithm of [14]. Since the concrete form of the algorithm depends on whether the proximal mapping of the data term $u\mapsto D(Au,f_{0})$ is explicit or not, we replace the data term $D(Au,f_{0})$ by

[TABLE]

where we assume the proximal mappings of $v\mapsto D_{i}(v,f_{0})$ to be explicit and, depending on the concrete applications, set either $D_{1}$ or $D_{2}$ to be the constant zero function.

Denoting by $g^{*}(v):=\sup_{w}(v,w)-g(w)$ the convex conjugate of a function $g$, with $(\cdot,\cdot)$ being the standard inner product of the sum of all pointwise-products of entries of $v$ and $w$, we reformulate (DP) to a saddle-point problem as

[TABLE]

Here, the dual variables $(p,q,d,r,m)\in(U^{2},U^{3},A(U),U,U^{3})$ are in the image space of the corresponding operators, $\mathcal{I}_{S}(z)=0$ if $z\in S$ and $\mathcal{I}_{S}(z)=\infty$ else, $\{\|\cdot\|_{\infty}\leq\delta\}:=\{z\,\left|\right.\,\|z\|_{\infty}\leq\delta\}$ with $\|z\|_{\infty}=\|(z_{1},\ldots,z_{l})\|_{\infty}=\sup_{i,j}\sqrt{\sum_{s=1}^{l}(z_{i,j}^{s})^{2}}$ a point-wise infinity norm on $z\in U^{l}$. The operator $B$ and the functional $G$ are given as

[TABLE]

and $F^{*}(y)=F^{*}(p,q,d,r,m)$ summarizes all the conjugate functionals as above. Applying the algorithm of [14] to this reformulation yields the numerical scheme as in Algorithm 1.

Note that there, we set either $D_{1}(\cdot,f_{0})\equiv 0$ such that the dual variable $d$ is constant 0 and line $8$ of the algorithm can be skipped, or we set $D_{2}(\cdot,f_{0})\equiv 0$ such that the proximal mapping in line $12$ reduces to the identity. The concrete choice of $D_{1},D_{2}$ and the proximal mappings will be given in the corresponding experimental sections. All other proximal mappings can be computed explicitly and reasonably fast: The mappings $\operatorname{proj}_{\alpha_{1}}$ and $\operatorname{proj}_{\alpha_{1}}$ can be computed as pointwise projections to the $L^{\infty}$-ball (see for instance [6]) and the mapping $\operatorname{prox}_{\sigma,(s_{2}(\mu)\nu\|\cdot\|_{1,2})^{*}}$ is a similar projection given as

[TABLE]

Most of the computational effort lies in the computation of $\operatorname{prox}_{\tau,s_{2}(\mu)(1-\nu)\|\cdot\|_{\text{nuc},\phi}}$, which, as the following lemma shows, can be computed via an SVD and a proximal mapping on the singular values.

Lemma 14**.**

Let $\phi:[0,\infty)\rightarrow[0,\infty)$ be a differentiable and increasing function and $\tau,\rho>0$ be such that $x\mapsto\frac{x^{2}}{2\tau}+\rho\phi(x)$ is convex on $[0,\infty)$. Then the proximal mapping of $\rho\|\cdot\|_{\text{nuc},\phi}$ for parameter $\tau$ is given as

[TABLE]

where $[A_{(NM,nn)}]=U\Sigma V^{T}$ is the SVD of $[A_{(NM,nn)}]$ and for $x_{0}\geq 0$

[TABLE]

In particular, in case $\phi(x)=x$ we have

[TABLE]

and in case

[TABLE]

we have that that $x\mapsto\frac{x^{2}}{2\tau}+\rho\phi(x)$ is convex whenever $\tau\leq\frac{1}{2\epsilon\delta\rho}$ and in this case

[TABLE]

Proof.

At first note that it suffices to consider $\rho\|\cdot\|_{\text{nuc},\phi}$ as a function on matrices and show the assertion without the reshaping operation. For any matrix $A$, we denote by $A=U_{A}\Sigma_{A}V^{T}_{A}$ the SVD of $A$ and $\Sigma_{A}=\text{diag}((\sigma^{A}_{i})_{i})$ contains the singular values sorted in non-increasing order, where $\Sigma_{A}$ is uniquely determined by $A$ and $U_{A},V_{A}$ are chosen to be suitable orthonormal matrices.

We first show that $G(A):=\frac{\|A\|^{2}_{2}}{2\tau}+\rho\|A\|_{\text{nuc},\phi}$ is convex. For $\lambda\in[0,1]$, $A,B$ matrices we get by sub-additivity of the singular values (see for instance [41]) that

[TABLE]

Now with $H(A):=\frac{\|A-A_{0}\|^{2}_{2}}{2\tau}+\rho\|A\|_{\text{nuc},\phi}$ we get that $H(A)=G(A)-\frac{1}{2\tau}(2(A,A_{0})+\|A_{0}\|_{2}^{2})$, thus also $H$ is convex. Hence first order optimality conditions are necessary and sufficient and we get (using the derivative of the singular values as in [33]) with $DH$ the derivative of $H$ that $A=\operatorname{prox}_{\tau,\rho\|\cdot\|_{\text{nuc},\phi}}(A_{0})$ is equivalent to

[TABLE]

and consequently to

[TABLE]

which is equivalent to

[TABLE]

as claimed. The other results follow by direct computation. ∎

Note also that, in Algorithm 1, $KC_{+}$ returns the part of the images that is represented by the atoms (the “texture part”) and $\text{RSV}(C_{+})$ stand for right-singular values of $[(C_{+})_{(NM,nn)}]$ and returns the image atoms. For the sake of simplicity, we use a rather high, fixed number of iterations in all experiment but note that, alternatively, a duality-gap based stopping criterion (see for instance [6]) could be used.

5 Numerical results

In this section we present numerical results obtained with the proposed method as well as its variants and compare to existing methods. We will mostly focus on the setting of (DP), where $\phi(x)=|x|$ and we use different data terms $D$. Hence, the regularization term is convex and consists of $\operatorname{TGV}_{\alpha}^{2}$ for the cartoon part and a weighted sum of a nuclear norm and $\ell^{1,2}$ norm for the texture part. Besides this choice of regularization (called CT-cvx), we will compare to pure $\operatorname{TGV}_{\alpha}^{2}$ regularization (called TGV), the setting of (DP) with the semi-convex potential $\phi$ as in (10) (call CT-scvx) and the setting of (DP) with $\operatorname{TGV}$ replaced by $\mathcal{I}_{\{0\}}$, i.e., only the texture norm is used for regularization, and $\phi(x)=|x|$ (called TXT). Further, in the last subsection, we also compare to other methods as specified there. For CT-cvx and CT-scvx we use the algorithm described in the previous section (where convergence can only be ensured for CT-cvx) and for the other variants we use a adaption of the algorithm to the respective special case.

We fix the size of the atom domain to $15\times 15$ pixel and the stride to $3$ pixel (see Section 4) for all experiments, and use four different test images (see Figure 3): The first two are synthetic images of size $120\times 120$, containing four different blocks of size $60\times 60$, whose size is a multiple of the chosen atom-domain size. The third and fourth image have size $128\times 128$ (not being a multiple of the atom-domain size) and the third image contains four sections of real images of size $64\times 64$ each (again not a multiple of the atom-domain size). All but the first image contain a mixture of texture and cartoon parts. The first four subsections consider only convex variants of our method ($\phi(x)=|x|$) and the last one considers the improvement obtained with a non-convex potential $\phi$ and also compares to other approaches.

Regarding the choice of parameters for all methods, we generally aimed to reduce the number of varying parameters for each method as much as possible such that for each method and type of experiment, at most two parameters need to be optimized. Whenever we incorporate the second order TGV functional for the cartoon part, we fix the parameters $(\alpha_{0},\alpha_{1})$ to $(\sqrt{2},1)$. The method CT-cvx then essentially depends on the three parameters $\lambda,\mu,\nu$. We experienced that the choice of $\nu$ is rather independent of the data and type of experiments, hence we leave it fixed for all experiments with incomplete or corrupted data, leaving our method with two parameters to be adapted: $\lambda$ defining the tradeoff between data and regularization and $\mu$ defining the tradeoff between cartoon and texture regularization. For the semi-convex potential we choose $\nu$ as with the convex one, fix $\delta=0.99$ and use two different choice of $\epsilon$, depending on the type of experiment, hence again leaving two parameters to be adapted. A summary of the parameter choice for all methods is provided in Table 11 below.

We also note that, whenever we tested a range of different parameters for any method presented below, we show the visually best results in the figure. Those are generally not the ones delivering the best result in terms of peak-signal-to-noise ratio, and for the sake of completeness we also provide in Table 10 the best PSNR result obtained with each method and each experiment over the range of tested parameters.

5.1 Image-atom learning and texture separation

As first experiment we test the method CT-cvx for learning image atoms and texture separation directly on the ground truth images. To this aim, we use

[TABLE]

The results can be found in Figure 4, where for the pure texture image we used only the texture norm (i.e. the method TXT) without the TGV part for regularization.

It can be observed that the proposed method achieves a good decomposition of cartoon and texture and also is able to learn the most important image structure effectively. While there are some repetitions of shifted structures in the atoms, the different structures are rather well-separated and the first nine atoms corresponding to the nine largest singular values still contain the most important features of the texture parts.

5.2 Inpainting and leaning from incomplete data

This section deals with the task of inpainting a partially available image and learning image atoms from this incomplete data. For reference, we also provide results with pure $\operatorname{TGV}_{\alpha}^{2}$ regularization (the method TGV). The data fidelity in this case is

[TABLE]

with $\mathcal{M}$ the index set of known pixels. Again we use only the texture norm for the first image (the method TXT) and the cartoon-texture functional for the others.

The results can be found in Figure 5. For the first and third image, $20\%$ of the pixels where given while for the other two, $30\%$ were given. It can be seen that our method is generally still able to identify the underlying pattern of the texture part and to reconstruct it reasonably well. Also the learned atoms are reasonable and are in accordance with the ones learned from the full data as in the previous section. In contrast to that, pure $\operatorname{TGV}$ regularization (which assumes piecewise smoothness) has no chance to reconstruct the texture patterns. For the cartoon part, both methods are comparable. It can also be observed that the target-like structure in the bottom right of the second image is not reconstructed well and also not well identified with the atoms (only the 8th one contains parts of this structure). The reason might be that due to the size of the repeating structure there is not enough redundant information available to reconstruct it from the missing data. Concerning the optimal PSNR values of Table 10, we can observe a rather strong improvement with CT-cvx compared to TGV.

5.3 Learning and separation under noise

In this section we test our method for image-atom-learning and de-noising with data corrupted by Gaussian noise (with standard deviation 0.5 and 0.1 times the image range for the Texture and the other images, respectively). Again we compare to TGV regularization in this section (but also to other methods in Section 5.5 below) and use the texture norm for the first image (the method TXT). The data fidelity in this case is

[TABLE]

It can be observed that also under the presence of rather strong noise, our method is able to learn some of the main features of the image within the learned atoms. Also the quality of the reconstructed image is improved compared to TGV, in particular for the right-hand side of the Mix image, where the top left structure is only visible in the result obtained with CT-cvx. On the other hand, the circle of the Patches image obtained with CT-cvx contains some artifacts of the texture part. Regarding the optimal PSNR values of Table 10, the improvement with CT-cvx compared to TGV is still rather significant.

5.4 Deconvolution

This section deals with the learning of image features and image reconstruction an an inverse problem setting, where the forward operator is given as a convolution with a Gaussian kernel (standard deviation 0.2, kernel size $9\times 9$ pixels) and the data is degraded by Gaussian noise with standard deviation 0.05 times the image range. The data fidelity in this case is

[TABLE]

with $A$ being the convolution operator.

We show results for the Mix and the Barbara image and compare to TGV. It can be seen that the improvement is comparable to the denoising case. In particular, the method is still able to learn reasonable atoms from the given, blurry data and in particular for the texture parts the improvement is quite significant. Regarding the optimal PSNR values, the improvement is roughly around 1 to 1.5 decibel.

5.5 Comparison

This section compares the method CT-cvx to its semi-convex variant CT-scvx and to other methods. At first, we consider the learning of atoms from incomplete data and image inpainting in Figure 8. It can be seen there that for the Patches image, the semi-convex variant achieves an almost perfect results: It is able to learn exactly the three atoms that compose the texture part of the image and to inpaint the image very well. For the Barabara image, where more atoms are necessary to synthesize the texture part, the two methods yield similar results and also the atoms are similar. These results are also reflected in the PSNR values of Table 10, where CT-scvx is more that 7 decibel better for the Patches image and achieve only a slight improvement for Barbara.

Next we consider the semi-convex variant CT-scvx for denoising the Patches and Barbara images of Figure 6. In this setting, also other methods are applicable and we compare to a costume implementation of a variant of the convolutional lasso algorithm (called CL) and to BM3D denoising [17] (called BM3D). For the former, we strive to solve the non-convex optimization problem

[TABLE]

where $(c_{i})_{i}$ are coefficient images, $p_{i}$ are atoms and $k$ is the number of used atoms. Note that we use the same boundary extension, atom-domain-size and stride variable than in the methods CT-cvx, CT-scvx, and that $\operatorname{TV}_{\rho}$ denotes a discrete TV functional with a slight smoothing of the $L^{1}$ norm to make it differentiable (see the source code [13] for details). For the solution we use an adaption of the algorithm of [34]. For BM3D we use the implementation obtained from [1].

Remark 15**.**

We note that, while we provide the comparison to BM3D in order to have a reference on achievable denoising quality, we do not aim to propose an improved denoising method that is comparable to BM3D. In contrast to BM3D, our method constitutes a variational (convex) approach, that is generally applicable for inverse problems and for which we were able to provide a detailed analysis in function space such that in particular stability and convergence results for vanishing noise can be proven. Furthermore, beyond mere image reconstruction, we regard the ability of simultaneous image-atom-learning and cartoon-texture decomposition as an important feature of our approach.

Results for the Patches and Barbara image can be found in Figure 9, where for CL we allowed for three atoms for the Patches images an tested 3,5 and 7 atoms for the Barbara image, showing the best result that was obtained with 7 atoms. It can be seen that, as with the inpainting results, CT-scvx achieves a very strong improvement compared to CT-cvx for the Patches image (obtaining the atoms almost perfectly) and only a slight improvement for the Barbara image. The CL method performs similar but slightly worse than CT-scvx. While for the Patches image also the three main features are identified correctly, they are not centered which leads to artifacts in the reconstruction and might be explained by the method being stuck in a local minimum. For the Patches image, the result of BM3D are comparable but slightly smoother than the ones of CT-scvx. In particular, the target-like structure in the bottom left is not very well reconstructed with BM3D but suffers from less remaining noise. For the Barbara image, BM3D delivers the best and result, but a slight overs-smoothing is visible. Regarding the PSNR values of Table 10, BM3D performs best and CT-scvx second best (better that CL), where in accordance with the visual results the difference of BM3D and CT-scvx is not as high as with Barbara.

6 Discussion

Using lifting techniques, we have introduced a (potentially convex) variational approach for learning image atoms from corrupted and/or incomplete data. An important part of our work is the analysis of the proposed model, which shows well-posedness results for the proposed model in function space for a general inverse problem setting. The numerical part shows that indeed our model can effectively learn image atoms from different types of data. While this works well also in a convex setting, moving to a semi-convex setting (which is also captured by our theory) yields a further, significant improvement. While the proposed method can also be regarded solely as image reconstruction method, we believe its main feature is in fact the ability to learn image atoms from incomplete data in a mathematically well understood framework.

Interesting future research questions are for instance the exploration of our method for classification problems or the exploration of similar lifting techniques for a mathematical understanding of deep neural networks.

7 Acknowledgements

MH acknowledges support by the Austrian Science Fund (FWF) (Grant J 4112). TP is supported by the European Research Council under the Horizon 2020 program, ERC starting grant agreement 640156.

Appendix A Appendix: Tensor spaces

We recall here some basic results on tensor products of Banach spaces that will be relevant for our work. Most of these results are obtained from [39, 19], to which we refer to for further information and a more complete introduction to the topic.

Throughout this section, let always $(X,\|\cdot\|_{X}),(Y,\|\cdot\|_{Y}),(Z,\|\cdot\|_{Z})$ be Banach spaces. By $X^{*}$ we denote the analytic dual of $X$, i.e., the space of bounded linear functionals from $X$ to $\mathbb{R}$. By $\mathcal{L}(X,Y)$ and $\mathcal{B}(X\times Y,Z)$ we denote the spaces of bounded linear and bilinear mappings, respectively, where the norm for the latter is given by $\|B\|_{\mathcal{B}}=\sup\{\|B(x,y)\|_{Z}\,\left|\right.\,\|x\|_{X}\leq 1,\,\|y\|_{Y}\leq 1\}$. In case the image space is the reals, we write $\mathcal{L}(X)$ and $\mathcal{B}(X\times Y)$.

Algebraic tensor product. The tensor product $x\otimes y$ of two elements $x\in X$, $y\in Y$ can be defined as a linear mapping on the space of bilinear forms on $X\times Y$ via

[TABLE]

The algebraic tensor product $X\otimes Y$ is then defined as the subspace of the space of linear functionals on $B(X,Y)$ spanned by elements $x\otimes y$ with $x\in X$, $y\in Y$.

Tensor norms. We will use two different tensor norms, the projective and the injective tensor norm (also known as the largest and smallest reasonable cross norm, respectively). The projective tensor norm on $X\otimes Y$ is defined for $u\in X\otimes Y$ as

[TABLE]

Note that indeed $\|\cdot\|_{\pi}$ is a norm and $\|x\otimes y\|_{\pi}=\|x\|_{X}\|y\|_{Y}$ (see [39, Proposition 2.1]). We denote by $X\otimes_{\pi}Y$ the completion of the space $X\otimes Y$ equipped with this norm. The following result gives a useful representation of elements in $X\otimes_{\pi}Y$ and their projective norm.

Proposition 16**.**

For $u\in X\otimes_{\pi}Y$ and $\epsilon>0$ there exist bounded sequences $(x_{n})_{n}\subset X$, $(y_{n})_{n}\subset Y$ such that

[TABLE]

In particular,

[TABLE]

Now for the injective tensor norm, we note that elements of the tensor product $X\otimes Y$ can be viewed as bounded bilinear forms on $X^{*}\times Y^{*}$ by associating to a tensor $u=\sum_{i=1}^{n}x_{i}\otimes y_{i}$ the bilinear form $B_{u}(\phi,\psi)=\sum_{i=1}^{n}\phi(x_{i})\psi(y_{i})$, where this association is unique (see [39, Section 1.3]). Hence $X\otimes Y$ can be regarded as a subspace of $\mathcal{B}(X^{*}\times Y^{*})$ and the injective tensor norm is the norm induced by this space. Thus for $u=\sum_{i=1}^{n}x_{i}\otimes y_{i}$ the injective tensor norm $\|\cdot\|_{\mathfrak{i}}$ is given as

[TABLE]

and the injective tensor product $X\otimes_{\mathfrak{i}}Y$ is defined as the completion of $X\otimes Y$ with respect to this norm.

Tensor lifting. The next result (see [39, Theorem 2.9]) shows that there is a one-to-one correspondence between bounded bilinear mappings from $X\times Y$ to $Z$ and bounded linear mappings from $X\otimes_{\pi}Y$ to $Z$.

Proposition 17**.**

For $B\in\mathcal{B}(X\times Y,Z)$ there exists a unique linear mapping $\hat{B}:X\otimes_{\pi}Y\rightarrow Z$ such that $\hat{B}(x\otimes y)=B(x,y)$. Further $\hat{B}$ is bounded and the mapping $B\mapsto\hat{B}$ is an isometric isomorphism between the Banach spaces $\mathcal{B}(X\times Y,Z)$ and $\mathcal{L}(X\otimes_{\pi}Y,Z)$.

Using this isometry, for $B\in\mathcal{B}(X\times Y,Z)$ we will always denote by $\hat{B}$ the corresponding linear mapping on the tensor product.

The following result is provided in [39, Proposition 2.3] and deals with the extension of linear operators to the tensor product.

Proposition 18**.**

Let $S\in\mathcal{L}(X,W)$, $T\in\mathcal{L}(Y,Z)$. Then there exists a unique operator $S\otimes_{\pi}T:X\otimes_{\pi}Y\rightarrow W\otimes_{\pi}Z$ such that $S\otimes_{\pi}T(x\otimes y)=(Sx)\otimes(Ty)$. Furthermore, $\|S\otimes_{\pi}T\|=\|S\|\|T\|$.

Tensor space isometries. The following proposition deals with duality of the injective and the projective tensor products. To this aim, we need the notion of Radon Nikodým property and approximation property, which we will not define here but rather refer to [39, Sections 4 and 5] and [19]. For our purposes, it is important to note that both properties hold for $L^{r}$-spaces with $r\in(1,\infty)$, the Radon Nikodým property holds for reflexive spaces, but while we cannot expect the Radon Nikodým property to hold for $L^{\infty}$ and $\mathcal{M}$, the approximation property does.

Lemma 19**.**

Assume that either $X^{*}$ or $Y^{*}$ has the Radon Nikodým property and that either $X^{*}$ or $Y^{*}$ has the approximation property. Then

[TABLE]

and for simple tensors $u=\sum_{i=1}^{n}x_{i}\otimes y_{i}\in X\otimes_{\mathfrak{i}}Y$ and $u^{*}=\sum_{i=1}^{m}x^{*}_{i}\otimes y^{*}_{i}\in X^{*}\otimes_{\pi}Y^{*}$ the duality pairing is given as

[TABLE]

Proof.

The identification of the duals is shown in [39, Theorem 5.33]. For the duality paring, we first note that the action of an element $u^{*}\in X^{*}\otimes_{\pi}Y^{*}$ on $X\otimes_{\mathfrak{i}}Y$ is given as the action of the associated bilinear form $B_{u^{*}}$ [39, Section 3.4], which for simple tensors $u=\sum_{i=1}^{n}x_{i}\otimes y_{i}$ can be given as

[TABLE]

Now in case also $u^{*}$ is a simple tensor, i.e., $u^{*}=\sum_{i=1}^{m}x_{i}^{*}\otimes y_{i}^{*}$, the action of this bilinear form can be given more explicitly [39, Section 1.3], which yields

[TABLE]

∎

The duality between the injective and projective tensor product will be used for compactness assertions on subsets of the latter. To this aim, we note in the following lemma that separability of the individual space transfers to the tensor product. As consequence, in case $X^{*}$ and $Y^{*}$ satisfy the assumption of Lemma 19 above and both admit a separable predual, also $X^{*}\otimes_{\pi}Y^{*}$ admits a separable predual and hence bounded sets are weakly* compact.

Lemma 20**.**

Let $X,Y$ be separable spaces. Then both $X\otimes_{\mathfrak{i}}Y$ and $X\otimes_{\pi}Y$ are separable.

Proof.

Take $X^{\prime}$ and $Y^{\prime}$ to be dense countable subsets of $X$ and $Y$, respectively. First note that it suffices to show that any simple tensor $x\otimes y$ can be approximated arbitrarily close by $x^{\prime}\otimes y^{\prime}$ with $x^{\prime}\in X^{\prime}$, $y^{\prime}\in Y^{\prime}$. But this is true since (using [39, Propositions 2.1 and 3.1])

[TABLE]

where $\|\cdot\|$ denotes either the projective or the injective norm.

The following result, which can be obtained by direct modification of the result shown at the beginning of [39, Section 3.2], provides an equivalent representation of the injective tensor product in a particular case.

Lemma 21**.**

Denote by $C_{c}(\Omega_{\Sigma},X)$ the space of compactly supported continuous functions mapping from $\Omega_{\Sigma}$ to $X$ and denote by $C_{0}(\Omega_{\Sigma},X)$ its completion with respect to the norm $\|\phi\|_{\infty}:=\sup_{t\in\Omega_{\Sigma}}\|X\|_{X}$. Then we have that

[TABLE]

where the isometry is given as the completion of the isometric mapping $J:C_{0}(\Omega_{\Sigma})\otimes X\rightarrow C_{0}(\Omega_{\Sigma},X)$ defined for $u=\sum_{i=1}^{n}f_{i}\otimes x_{i}$ as

[TABLE]

Next we consider the identification of tensor products with linear operators which is provided in the following proposition [39, Corollary 4.8].

Proposition 22**.**

Define the mapping $J:X^{*}\otimes_{\pi}Y\rightarrow\mathcal{L}(X,Y)$ as

[TABLE]

Then, $J$ is well-defined and has unit norm. Defining $\mathcal{N}(X,Y)\subset\mathcal{L}(X,Y)$ as the range of $J$, equipped with the norm

[TABLE]

we get that $\mathcal{N}(X,Y)$ is a Banach space, called the space of nuclear operators. If further either $X^{*}$ or $Y$ has the approximation property, then $J$ is an isometric isomorphism, that is, we can identify

[TABLE]

It is easy to see that nuclear operators are compact that that we can equivalently write

[TABLE]

Also, in a Hilbert space setting (see [43] for details), it is a classical result that for any compact $T\in\mathcal{L}(H_{1},H_{2})$ with $(H_{1},(\cdot,\cdot)),(H_{2}(\cdot,\cdot))$ Hilbert spaces there exist orthonormal systems $(x_{i})_{i}$, $(y_{i})_{i}$ and uniquely defined singular values $(\sigma_{i})_{i}:=(\sigma_{i}(T))_{i}$ such that

[TABLE]

In addition, in case $T$ has finite nuclear norm, it follows that $\|T\|_{\text{nuc}}=\sum_{i=1}^{\infty}\sigma_{i}.$