Iterative Potts minimization for the recovery of signals with discontinuities from indirect measurements -- the multivariate case

TL;DR

This paper introduces an iterative method for minimizing the multivariate Potts functional to recover signals with discontinuities from noisy, indirect measurements, with proven convergence and supporting numerical experiments.

Contribution

A novel iterative minimization strategy for the multivariate Potts functional addressing indirect, noisy data, with theoretical convergence analysis.

Findings

Convergence of the proposed method is established.

Numerical experiments demonstrate effectiveness in signal recovery.

Method outperforms existing approaches in handling noise and indirect measurements.

Abstract

Signals and images with discontinuities appear in many problems in such diverse areas as biology, medicine, mechanics, and electrical engineering. The concrete data are often discrete, indirect and noisy measurements of some quantities describing the signal under consideration. A frequent task is to find the segments of the signal or image which corresponds to finding the discontinuities or jumps in the data. Methods based on minimizing the piecewise constant Mumford-Shah functional -- whose discretized version is known as Potts functional -- are advantageous in this scenario, in particular, in connection with segmentation. However, due to their non-convexity, minimization of such functionals is challenging. In this paper we propose a new iterative minimization strategy for the multivariate Potts functional dealing with indirect, noisy measurements. We provide a convergence analysis and underpin our findings with numerical experiments.

Figures

Equations

$${\mathrm{argmin}}_u\gamma \Vert \nabla u{\Vert }_0+{\Vert Au-f\Vert }_2^2.$$

$$P_{\gamma }(u)=\Vert Au-f{\Vert }_2^2+\gamma \sum _{s=1}^S{\omega }_s{\Vert {\nabla }_{a_s}u\Vert }_0,$$

$$P_{\gamma }(u_1,\dots ,u_S)\to \min ,\text{ s.t. }{u}_1=\dots =u_S,$$

$$P_{\gamma }(u_1,\dots ,u_S)=\sum _{s=1}^S\frac{1}{S}{\Vert A{u}_s-f\Vert }_2^2+\gamma \sum _{s=1}^S{\omega }_s{\Vert {\nabla }_{a_s}{u}_s\Vert }_0.$$

$$P_{\gamma ,\rho }(u_1,\dots ,u_S)=\sum _{s=1}^S\frac{1}{S}{\Vert A{u}_s-f\Vert }_2^2+\gamma \sum _{s=1}^S{\omega }_s{\Vert {\nabla }_{a_s}{u}_s\Vert }_0+\rho \sum _{1\le s<s^{\prime }\le S}{c}_{s,s^{\prime }}\Vert u_s-u_{s^{\prime }}{\Vert }_2^2.$$

$${\rho }_{s,s^{\prime }}=\rho c_{s,s^{\prime }}.$$

$${\rho }_{s,s^{\prime }}=\rho \text{ for all }s,s^{\prime },\text{ or}{\rho }_{s,s^{\prime }}=\rho {\delta }_{((s+1)\mathrm{mod}S),s^{\prime }},$$

$$B\ =\ \begin{pmatrix}S^{-1/2}A&0&&\cdots&&0\\ 0&S^{-1/2}A&&\cdots&&0\\ \vdots&&&\ddots&&\vdots\\ 0&0&&\cdots&S^{-1/2}A&0\\ 0&0&&\cdots&0&S^{-1/2}A\vspace{.8em}\\ \rho_{1,2}^{1/2}I&-\rho_{1,2}^{1/2}I&0&\ldots&0&0\\ \rho_{1,3}^{1/2}I&0&-\rho_{1,3}^{1/2}I&\ldots&0&0\\ &\vdots&&&\vdots&\\ \rho_{1,S}^{1/2}I&0&0&\ldots&0&-\rho_{1,S}^{1/2}I\\ 0&\rho_{2,3}^{1/2}I&-\rho_{2,3}^{1/2}I&\ldots&0&0\\ &\vdots&&&\vdots&\\ 0&\rho_{2,S}^{1/2}I&0&\ldots&0&-\rho_{2,S}^{1/2}I\\ &&&\vdots&&\\ &&&\vdots&&\\ 0&0&0&\ldots&\rho_{S-1,S}^{1/2}I&-\rho_{S-1,S}^{1/2}I\\ \end{pmatrix},\qquad g\ =\ \begin{pmatrix}S^{-1/2}f\\ S^{-1/2}f\\ \vdots\\ S^{-1/2}f\\ S^{-1/2}f\vspace{.8em}\\ 0\\ 0\\ \vdots\\ 0\\ 0\\ \vdots\\ 0\\ \vdots\\ \vdots\\ 0\\ \end{pmatrix}.\$$

$$D(u_1,\dots ,u_S)=\left(\begin{array}{c}{\nabla }_{a_1}{u}_1\\ ⋮\\ {\nabla }_{a_S}{u}_S\end{array}\right)$$

$$\Vert D(u_1,\dots ,u_S){\Vert }_{0,\omega }=\sum _{s=1}^S{\omega }_s{\Vert {\nabla }_{a_s}{u}_s\Vert }_0.$$

$$P_{\gamma,\rho}(u_{1},\ldots,u_{S})=\left\|B(u_{1},\ldots,u_{S})^{T}-g\right\|_{2}^{2}+\gamma\ \Big{\|}\ D(u_{1},\ldots,u_{S})\ \Big{\|}_{0,\omega}.$$

$$P_{\gamma ,\rho }^{\mathrm{surr}}$$

$$-\frac{1}{L_{\rho }^2}{\Vert B(u_1,\dots ,u_S{)}^T-B(v_1,\dots ,v_S{)}^T\Vert }_2^2+{\Vert (u_1,\dots ,u_S{)}^T-(v_1,\dots ,v_S{)}^T\Vert }_2^2.$$

$$L_{\rho }^2>\Vert A{\Vert }_2^2/S+2\underset{s\in \{1,\dots ,S\}}{\max }\sum _{s^{\prime }:s^{\prime }\ne s}^S{\rho }_{s,s^{\prime }}.$$

$$P_{\gamma ,\rho }^{\mathrm{surr}}(u_1,\dots ,u_S,v_1,\dots ,v_S)=$$

$$\displaystyle+\frac{\gamma}{L_{\rho}^{2}}\Big{\|}D(u_{1},\ldots,u_{S})\Big{\|}_{0,\omega}+R(v_{1},\ldots,v_{S}),$$

$$P_{\gamma ,\rho }^{\mathrm{surr}}$$

$$=$$

$$\left(u_1^{(n+1)},\dots ,u_S^{(n+1)}\right)\in \underset{u_1,\dots ,u_S}{\mathrm{argmin}}\sum _{s=1}^S\left[{\Vert u_s-h_s^{(n)}\Vert }_2^2+\frac{\gamma {\omega }_s}{L_{\rho }^2}{\Vert {\nabla }_{a_s}{u}_s\Vert }_0\right]$$

$$h_s^{(n)}=u_s^{(n)}+\frac{1}{S{L}_{\rho }^2}{A}^{\ast }f-\frac{1}{S{L}_{\rho }^2}{A}^{\ast }A{u}_s^{(n)}-\sum _{s^{\prime }:s^{\prime }\ne s}\frac{{\rho }_{s,s^{\prime }}}{L_{\rho }^2}(u_s^{(n)}-u_{s^{\prime }}^{(n)}),\text{ for all }s\in \{1,\dots ,S\}.$$

$$\Vert u_s-u_{s^{\prime }}{\Vert }_2^2=\sum _{i,j}|(u_s{)}_{ij}-(u_{s^{\prime }}{)}_{ij}{|}^2<\frac{{\epsilon }^2}{c_{s,s^{\prime }}},$$

$$h_s^{(n)}=u_s^{(n)}+\frac{1}{S{L}_{\rho }^2}{A}^{\ast }f-\frac{1}{S{L}_{\rho }^2}{A}^{\ast }A{u}_s^{(n)}-\sum _{s^{\prime }:s^{\prime }\ne s}\frac{{\rho }_{s,s^{\prime }}}{L_{\rho }^2}(u_s^{(n)}-u_{s^{\prime }}^{(n)}),s=1,\dots ,S,$$

$$\left(u_1^{(n+1)},\dots ,u_S^{(n+1)}\right)\in \underset{u_1,\dots ,u_S}{\mathrm{argmin}}\sum _{s=1}^S\left[{\Vert u_s-h_s^{(n)}\Vert }_2^2+\frac{\gamma {\omega }_s}{L_{\rho }^2}{\Vert {\nabla }_{a_s}{u}_s\Vert }_0\right].$$

$${\cup }_{i=1}^{N_{\mathcal{P}}}{\mathcal{P}}_i=\Omega ,{\mathcal{P}}_i\cap {\mathcal{P}}_j=\varnothing \text{for all i,j=1,…,NP.}$$

$$\mathcal{I}(\mathcal{P})={\mathcal{I}}_{\mathcal{P}},\mathcal{P}(\mathcal{I})={\mathcal{P}}_{\mathcal{I}},$$

$${\pi }_{\mathcal{P}}(u){|}_{{\mathcal{P}}_i}=\frac{{\sum }_{x\in {\mathcal{P}}_i}{\sum }_{s=1}^S{u}_s(x)}{S\#{\mathcal{P}}_i},$$

$$t>2{\sigma }_1^{-1/2}{S}^{-1/2}\Vert A\Vert \Vert f\Vert ,$$

$$t>\frac{2}{S}\Vert A\Vert \Vert f\Vert ,\text{ and }t>2(2-2\cos (2\pi /S){)}^{-1/2}{S}^{-1/2}\Vert A\Vert \Vert f\Vert ,\text{ respectively.}$$

$$\Vert u_s^{(k,n)}-u_{s^{\prime }}^{(k,n)}\Vert >\frac{t}{\rho \sqrt{c_{s,s^{\prime }}}},\text{ or }\Vert u_s^{(k,n)}-u_s^{(k,n-1)}\Vert >\frac{\delta }{L_{\rho }}$$

$$h_s^{(k,n)}=u_s^{(k,n)}+\frac{1}{S{L}_{\rho }^2}{A}^{\ast }f-\frac{1}{S{L}_{\rho }^2}{A}^{\ast }A{u}_s^{(k,n)}-\sum _{s^{\prime }:s^{\prime }\ne s}\frac{{\rho }_{s,s^{\prime }}}{L_{\rho }^2}(u_s^{(k,n)}-u_{s^{\prime }}^{(k,n)}),s=1,\dots ,S,$$

$$\left(u_1^{(k,n+1)},\dots ,u_S^{(k,n+1)}\right)\in \underset{u_1,\dots ,u_S}{\mathrm{argmin}}\sum _{s=1}^S\left[{\Vert u_s-h_s^{(k,n)}\Vert }_2^2+\frac{\gamma {\omega }_s}{L_{\rho }^2}{\Vert {\nabla }_{a_s}{u}_s\Vert }_0\right],$$

$$u_s^{(k+1)}=u_s^{(k+1,0)}=u_s^{(k,n)},$$

$$\begin{array}{r}\underset{u_s:\Omega \to \mathbb{R}}{\mathrm{argmin}}\Vert u_s-f{\Vert }_2^2+{\gamma }_s^{\prime }\Vert {\nabla }_{a_s}u{\Vert }_0\end{array}$$

Full text

Keywords: Piecewise-constant Mumford-Shah model, Potts model, majorization-minimization methods, image segmentation, joint reconstruction and segmentation, ill-posed inverse problems, Radon transform, deconvolution.

AMS subject classifications. 94A08, 68U10, 65D18, 65K10, 90C26, 90C39.

1 Introduction

Problems involving reconstruction tasks for functions with discontinuities appear in various biological and medical applications. Examples are the steps in the rotation of the bacterial flagella motor [79, 78, 70], the cross-hybridization of DNA [77, 30, 44], x-ray tomography [74], electron tomography [49] and SPECT [51, 93]. An engineering example is crack detection in brittle material in mechanics [3]. Further examples may for instance be found in the papers [58, 59, 25, 22, 34] and the references therein. In general, signals with discontinuities appear in many applied problems. A central task is to restore the jumps, edges, change points or segments of the signals or images from the observed data. These observed data are usually indirectly measured. Furthermore, they consist of measurements on a discretized grid and are typically corrupted by noise.

In many scenarios, nonconvex nonsmooth variational methods are a suitable choice for the partitioning task, i.e., the task of finding the jumps/edges/change points; see for example [58, 70, 13]. In particular, methods based on piecewise constant Mumford-Shah functionals [62, 63] have been used in various different applications. The piecewise-constant Mumford-Shah model also appears in statistics and image processing where it is often called Potts model [36, 13, 14, 15, 72, 91]; this is a tribute to Renfrey B. Potts and his work in statistical mechanics [73]. The variational formulation of the piecewise-constant Mumford-Shah/Potts model (with an indirect measurement term) is given by

[TABLE]

Here, $A$ is a linear operator modeling the measurement process, e.g., the Radon transform in computed tomography (CT), or the point-spread function of the microscope in microscopy. Further, $f$ is an element of the data space, e.g., a sinogram or part of it in CT, or the blurred microscopy image in microscopy. The mathematically precise definition of the jump term $\|\nabla u\|_{0}$ in the general situation is rather technical. However, if $u$ is piecewise-constant and the discontinuity set of $u$ is sufficiently regular, say, a union of $C^{1}$ curves, then $\|\nabla u\|_{0}$ is just the total arc length of this union. In general, the gradient $\nabla u$ is given in the distributional sense and the boundary length is expressed in terms of the $(d-1)$-dimensional Hausdorff measure. When $u$ is not piecewise-constant, the jump penalty is infinite [75]. The second term measures the fidelity of a solution $u$ to the data $f.$ The parameter $\gamma>0$ controls the balance between data fidelity and jump penalty. (A wider class of Mumford-Shah models can be obtained by replacing the squared $L^{2}$ distance by more general data terms such as other norm-based expressions or divergences.)

The piecewise-constant Mumford-Shah/Potts model can be interpreted in two ways. On the one hand, if the imaged object is (approximately) piecewise-constant, then the solution is an (approximate) reconstruction of the imaged object. On the other hand, since a piecewise-constant solution directly induces a partitioning of the image domain, it can be seen as joint reconstruction and segmentation. Executing reconstruction and segmentation jointly typically leads to better results than performing the two steps successively [51, 74, 75, 85]. We note that in order to deal with the discrete data, the energy functional is typically discretized; see Section 2.1. Some references concerning Mumford-Shah functionals are [2, 18, 75, 45, 8, 33, 67] and also the references therein; see also the book [1]. The piecewise constant Mumford-Shah functionals are among the maybe most well-known representatives of the class of free-discontinuity problems introduced by De Giorgi [29].

The analysis of the nonsmooth and nonconvex problem (1) is rather involved. We discuss some analytic aspects. We first note that without additional assumptions the existence of minimizers of (1) is not guaranteed in a continuous domain setting [33, 32, 75, 84]. To ensure the existence of minimizers, additional penalty terms such as an $L^{p}$ ($1<p<\infty$) term of the form $\|u\|_{p}^{p}$ [74, 75] or pointwise boundedness constraints [45] have been considered. We note that the existence of minimizers is guaranteed in the discrete domain setup for typical discretizations [33, 84]. Another important topic is to verify that the Potts model is a regularization method in the sense of inverse problems. The first work dealing with this task is [75]: The authors assume that the solution space consists of non-degenerate piecewise-constant functions with at most $k$ (arbitrary, but fixed) different values which are additionally bounded. Under relatively mild assumptions on the operator $A$, they show stability. Further, by giving a suitable parameter choice rule, they show that the method is a regularizer in the sense of inverse problems. Related references are [50, 45] with the latter including (non-piecewise-constant) Mumford-Shah functionals. We note that Mumford-Shah approaches (including the piecewise constant Mumford-Shah variant) also regularize the boundaries of the discontinuity set of the underlying signal [45].

Solving the Potts problem is algorithmically challenging. For $A=\mathrm{id},$ it is NP-hard for multivariate domains [87, 13], and, for general linear operators $A,$ it is even NP-hard for univariate signals [84]. Thus, finding a global minimizer within reasonable time seems to be unrealistic in general. Nevertheless, due to its importance, many approximative strategies for multivariate Potts problems with $A=\mathrm{id}$ have been proposed. (We note that the case $A=\mathrm{id}$ is important as well since it captures the partitioning problem in image processing.) For the Potts problem with general $A$ there are still some but not that many existing approaches, in particular in the multivariate situation. For a more detailed discussion, we refer to the paragraph on algorithms for piecewise constant Mumford-Shah problems below. A further discussion of methods for reconstructing piecewise constant signals may be found in [59]. In [90], we have considered the univariate Potts problem for a general operator $A$ and have proposed a majorization-minimization strategy which we called iterative Potts minimization in analogy to iterative thresholding schemes. In this work, we will develop iterative Potts minimization schemes for the more demanding multivariate situation which is important for multivariate applications as appearing in imaging problems.

Existing algorithmic approaches to the piecewise constant Mumford-Shah problem and related problems.

We start to consider the Potts problem for general operator $A.$ In [5], Bar et al. consider an Ambrosio-Tortorelli-type approximation. Kim et al. use a level-set based active contour method for deconvolution in [48]. Ramlau and Ring [74] employ a related level-set approach for the joint reconstruction and segmentation of x-ray tomographic images; further applications are electron tomography [49] and SPECT [51]. The authors of the present paper have proposed a strategy based on the alternating methods of multipliers in [84] for the univariate case and in [85] for the multivariate case.

Fornasier and Ward [33] rewrite Mumford-Shah problems as a pointwise penalized problem and derive generalized iterative thresholding algorithms for the rewritten problems in the univariate situation. Further, they show that their method converges to a local minimizer in the univariate case. Their approach principally carries over to the piecewise constant Mumford-Shah functional as explained in [84, 90] and then results in a $\ell^{0}$ sparsity problem. In the univariate situation, this NP-hard optimization problem is unconstrained and may be addressed by iterative hard thresholding algorithms for $\ell^{0}$ penalizations, analyzed by Blumensath and Davies in [9, 10]. (Note that related algorithms based on iterative soft thresholding for $\ell^{1}$ penalized problems have been considered by Daubechies, Defrise, and De Mol in [28].) Artina et al. [3] in particular consider the multivariate discrete Mumford-Shah model using the pointwise penalization approach of [33]. In the multivariate setting, this results in a corresponding nonconvex and nonsmooth problem with linear constraints. The authors successively minimize local quadratic and strictly convex perturbations (depending on the previous iterate) of a (fixed) smoothed version of the objective by augmented Lagrangian iterations which themselves can be accomplished by iterative thresholding via a Lipschitz continuous thresholding function. They show that the accumulation points of the sequences produced by their algorithm are constraint critical points of the smoothed problem. In the multivariate situation, a similar approach for rewriting the Potts problem results in an $\ell^{0}$ sparsity problem with additional equality constraints. Algorithmic approaches for such $\ell^{0}$ sparsity problem with equality constraints are the penalty decomposition methods of [60, 61, 96]. The connection with iterative hard thresholding is that the inner loop of the employed two stage process usually is of iterative hard thresholding type. The difference of the hard thresholding based methods to our approach in this paper is that we do not have to deal with constraints and the full matrix $A$ but with the nonseparable regularizing term $\|\nabla u\|_{0}$ instead of its separable counterpart $\|u\|_{0}.$ Hence we cannot use hard thresholding.

Another frequently appearing method in the context of restoration of piecewise constant images is total variation minimization [76]. There the jump penalty $\|\nabla u\|_{0}$ is replaced by the total variation $\|\nabla u\|_{1}.$ The arising minimization problem is convex and therefore numerically tractable with convex optimization techniques [21, 26]. Candès, Wakin, and Boyd [17] use iteratively reweighted total variation minimization for piecewise constant recovery problems. Results of compressed sensing type related to the Potts problem have been derived by Needell and Ward [65, 66]: under certain conditions, minimizers of the Potts function agree with total variation minimizers. However, in the presence of noise, total variation minimizers might significantly differ from minimizers of the Potts problem. But the minimizers of the Potts problem are the results frequently desired in practice. Further, algorithms based on convex relaxations of the Potts problem (1) have gained a lot of interest in recent years; see, e.g., [56, 4, 20, 16, 37, 86].

We next discuss approaches for the multivariate Potts problem for the situation $A=\mathrm{id}$ which is particularly interesting in image processing and for which there are some further approaches. The first class of approaches is the approach via graph cuts. Here, the range space of $u$ is a priori restricted to a relatively small number of values. The problem remains NP-hard, but it then allows for an approach by sequentially solving binary partitioning problems via minimal graph cut algorithms [13, 52, 12]. We here point out that this approach also can deal with (possibly non-convex) data fidelity terms more general than the squared $L^{2}$ data term employed in (1) (in the case $A=\mathrm{id}$). Another approach is to limit the number $k$ of different values which $u$ may take without discretizing the range space a priori. For $k=2,$ active contours were used by Chan and Vese [24] to minimize the corresponding binary Potts model. They use a level set function to represent the partitions which evolves according to the Euler-Lagrange equations of the Potts model. A globally convergent strategy for the binary segmentation problem is presented in [23]. The active contour method for $k=2$ was extended to larger $k$ in [88]. Note that, for $k>2$ the problem is NP hard. We refer to [27] for an overview on level set segmentation. In [40, 41, 42], Hirschmüller proposes a non-iterative strategy for the Potts problem which is based on cost aggregation. It has lower computational cost, but comes with lower quality reconstructions compared with graph cuts. Due to the small number of potential values of $u,$ these methods mainly appear in connection with image segmentation. Methods for restoring piecewise constant images without restricting the range space are proposed in Nikolova et al. [69, 68]. They use non-convex regularizers which are algorithmically approached using a graduated non-convexity approach. We note that the Potts problem (1) does not fall into the class of problems considered in [69, 68]. Last but not least, Xu et al. [94] proposed a piecewise constant model reminiscent of the Potts model that is approached by a half-quadratic splitting using a pixelwise iterative thresholding type technique. It was later extended to a method for blind image deconvolution [95].

Contributions.

The contributions of this paper are threefold: (i) We propose a new iterative minimization strategy for multivariate piecewise constant Mumford-Shah/Potts objective functions as well as a (still NP-hard) quadratic penalty relaxation. (ii) We provide a convergence analysis of the proposed schemes. (iii) We show the applicability of our schemes in several experiments.

Concerning (i), we propose two schemes which are based on majorization-minimization or forward-backward splitting methods of Douglas-Rachford type [57]. The one scheme addresses the Potts problem directly, whereas the other scheme treats a quadratic penalty relaxation. The solutions of the relaxed problem themselves are not feasible for the Potts problem but near to a feasible solution of the Potts problem where nearness can be quantified. In particular, when a given tolerance in applications is acceptable the relaxed scheme is applicable. In contrast to the approaches in [33, 9] and [60, 61] for sparsity problems which lead to thresholding algorithms, our approach leads to non-separable yet computationally tractable problems in the backward step.

Concerning (ii), we first analyze the proposed quadratic penalty relaxation scheme. In particular, we show convergence towards a local minimizer. Due to the NP hardness of the quadratic penalty relaxation, the convergence result is in the range of what can be expected best. Concerning the scheme for the non-relaxed Potts problem we also perform a convergence analysis. In particular, we obtain results on the convergence towards local minimizers on subsequences. The quality of the convergence results is comparable with the ones in [60, 61]. We note that compared with [60, 61] we face the additional challenge to deal with the non-separability of the backward step. (We note that in practice we observe convergence of the whole sequence, not on a subsequence.)

Concerning (iii) we consider problems with full and partial data. We begin to apply our algorithms to deconvolution problems. In particular, we consider deblurring and denoising Gaussian blur images and motion blur images, respectively. We further consider noisy and undersampled Radon data, together with the task of joint reconstruction, denoising and segmentation. Finally, we use our method in the situation of pure image partitioning (without blur) which is a widely considered problem in computer vision.

Organization of the paper.

In Section 2 we derive the proposed algorithmic schemes. In Section 3 we provide a convergence analysis for the proposed schemes. In Section 4 we apply the algorithms derived in the present paper to concrete reconstruction problems. In Section 5 we draw conclusions.

2 Majorization-minimization algorithms for multivariate Potts problems

2.1 Discretization

We use the following finite difference type discretiziation of the multivariate Potts problem (1) given by

[TABLE]

where the $a_{s}\in\mathbb{Z}^{2}$ come from a finite set of directions and the symbol $\nabla_{a_{s}}u\,(i,j)$ denotes the directional difference $u_{(i,j)+a_{s}}-u_{i,j}$ with respect to the direction $a_{s}$ at the pixel $(i,j)$. The symbol $\|\nabla_{a_{s}}u\|_{0}$ denotes the number of nonzero entries of $\nabla_{a_{s}}u.$ The simplest set of directions consists of the unit vectors $a_{1}=(0,1),$ $a_{2}=(1,0)$ along with unit weights. Unfortunately, when refining the grid, this discretization converges to a limit that measures the boundary in terms of the $\ell^{1}$ analogue of the Hausdorff measure [18]. The practical consequences are unwanted block artifacts in the reconstruction (geometric staircasing). More isotropic results are obtained by adding the diagonals $a_{3}=(1,1),a_{4}=(1,-1)$ to the directions $a_{1}$ and $a_{2};$ a near isotropic discretization can be achieved by extending this system by the knight moves $a_{5}=(1,2),a_{6}=(2,1),a_{7}=(1,-2),a_{8}=(2,-1).$ (The name is inspired by the possible moves of a knight in chess.) Weights $\omega_{s}$ for the system $\{a_{1},a_{2},a_{3},a_{4}\}$ of coordinate directions and diagonal directions can be chosen as $\omega_{s}=\sqrt{2}-1$ for the coordinate part $s=1,2$ and $\omega_{s}=1-\tfrac{\sqrt{2}}{2}$ for diagonal part $s=3,4$. When additionally adding knight-move directions, weights $\omega_{s}$ for the system $\{a_{1},\ldots,a_{8}\}$ can be chosen as $\omega_{s}=\sqrt{5}-2$ for the coordinate part $s=1,2,$ $\omega_{s}=\sqrt{5}-\frac{3}{2}\sqrt{2}$ for diagonal part $s=3,4$, and $\omega_{s}=\frac{1}{2}(1+\sqrt{2}-\sqrt{5})$ for diagonal part $s=5,\ldots,8.$ There are several ways to derive weights $\omega_{s}$ for the neighborhood systems: the method of [19] is based on an optimization approach, the method of [11] is based on the Cauchy-Crofton formula, and the approach of [85] is based on equating the euclidean lengths of straight lines and the lengths of their digital counterparts. We note that for the system $\{a_{1},a_{2},a_{3},a_{4}\}$ of coordinate directions and diagonal directions the weights of [19] and in [85] coincide; the weights displayed for the knight-move case above are the ones derived by the scheme in [85]. For further details we refer to these references.

We record that the considered problem (2) has a minimizer.

Theorem 1**.**

The discrete multivariate Potts problem (2) has a minimizer.

The validity of Theorem 1 can be seen by following the lines of the proof of [43, Theorem 2.1] where an analogous statement is shown for the (non-piecewise constant) Mumford-Shah problem.

Vector-valued images.

We briefly discuss the extension of (2) to vector-valued images and multi-channel data, e.g., (blurred) RGB color images. To this end, we assume multi-channel data $f=(f_{1},\ldots,f_{C})$ consisting of $C$ channels and images $u=(u_{1},\ldots,u_{C})$. In this situation, the role of the first summand on the right-hand side of (2) is taken by the channel-wise sum $\sum_{c=1}^{C}\|Au_{c}-f_{c}\|_{2}^{2}$. The symbol $\nabla_{a_{s}}u(i,j)$ now denotes the vector of directional differences with entries $u_{(i,j)+a_{s},c}-u_{i,j,c}$, $c=1,\ldots,C$ and the entirety of these vectors form the rows of $\nabla u_{a_{s}}$. Consequently, $\|\nabla_{a_{s}}u\|_{0}$ denotes the number of nonzero rows of $\nabla_{a_{s}}u$. As a result, introducing a jump between two pixels in all channels has the same costs as opening a jump in a single channel only. This enforces the jumps to be aligned across the channels which is in contrast to a channel-wise application of the single-channel Potts model (2).

2.2 Derivation of the proposed algorithmic schemes

We start out with the discretization (2) of the multivariate Potts problem. We introduce $S$ versions $u_{1},\ldots,u_{S}$ of the target $u$ and link them via equality constraints in the following consensus form to obtain the problem

[TABLE]

where the function $P_{\gamma}(u_{1},\ldots,u_{S})$ of the $S$ variables $u_{1},\ldots,u_{S}$ is given by

[TABLE]

Note that solving (3) is equivalent to solving the discrete Potts problem (2). Further, note that we have overloaded the symbol $P_{\gamma}$ which, for one argument $u,$ denotes the Potts function of (2) and for $S$ arguments $u_{1},\ldots,u_{S}$ denotes the energy function of (4); we have the relation $P_{\gamma}(u,\ldots,u)=P_{\gamma}(u)$.

A majorization-minimization approach to the quadratic penalty relaxation of the Potts problem.

The quadratic penalty relaxation of (4) is given by

[TABLE]

Here, the soft constraints which replace the equalities $u_{1}=\ldots=u_{S}$ are realized via the squared Euclidean norms $\sum_{1\leq s<s^{\prime}\leq S}c_{s,s^{\prime}}\ \|u_{s}-u_{s^{\prime}}\|_{2}^{2},$ where the nonnegative numbers $c_{s,s^{\prime}}$ denote weights (which may be set to zero if no direct coupling between the particular $u_{s},u_{s^{\prime}}$ is desired.) The symbol $\rho$ denotes a positive penalty parameter promoting the soft constraint, i.e., increasing $\rho$ enforces the $u_{i}$ to be closer to each other w.r.t. the Euclidean distance. We note that we later analytically quantify the size of $\rho$ which is necessary to obtain an a priori prescribed tolerance in the $u_{i};$ see (18) below. Frequently, we use the short-hand notation

[TABLE]

Typical choices of the $\rho_{s,s^{\prime}}$ are

[TABLE]

i.e., the constant choice ($c_{s,s^{\prime}}=1$), as well as the coupling between consecutive variables with constant parameter ($\delta_{s,t}=1$ if and only if $s=t,$ and $\delta_{s,t}=0$ otherwise.) We note that in these situations only one additional positive parameter $\rho$ appears, and that this parameter is tied to the tolerance one is willing to accept as a distance of the $u_{i};$ see Algorithm 1.

For the majorization-minimization approach, we derive a surrogate functional [28] of the function $P_{\gamma,\rho}(u_{1},\ldots,u_{S})$ of (5). For this purpose, we introduce the block matrix $B$ and the vector $g$ given by

[TABLE]

Here $I$ denotes the identity matrix and [math] the zero matrix; The matrix $B$ has $S$ block columns and $S+S(S-1)/2$ block rows. Further, we introduce the difference operator $D$ given by

[TABLE]

which applies the difference w.r.t. the $i$th direction to the $i$th component of $u.$ We employ the weights $\omega_{1},$ $\ldots,\omega_{S}$ to define the quantity $\|D(u_{1},\ldots,u_{S})\|_{0,\omega}$ which counts the weighted number of jumps by

[TABLE]

With all this comprehensive notation at hand, we may rewrite the function of (5) as

[TABLE]

Using the representation (11), the surrogate functional in the sense of [28] of $P_{\gamma,\rho}$ is given by

[TABLE]

Here $L_{\rho}\geq 1$ denotes a constant which is chosen larger than the spectral norm $\|B\|$ of $B$ (i.e., the operator norm w.r.t. the $\ell^{2}$ norm.) This scaling is made to ensure that $B/L_{\rho}$ is contractive. In terms of $A$ and the penalties $\rho_{s,s^{\prime}},$ we require that

[TABLE]

For the particular choice $\rho_{s,s^{\prime}}=\rho$ as on the left-hand side of (7) we can choose $L_{\rho}^{2}$ smaller, i.e., $L_{\rho}^{2}>\|A\|_{2}^{2}/S+S\rho.$ For only coupling neighboring $u_{s}$ with the same constant $\rho$, i.e., the right-hand coupling of (7), we have $L_{\rho}^{2}>\|A\|_{2}^{2}/S+\alpha\rho,$ where $\alpha=4,$ if $S$ is even, and $\alpha=2-2\cos\left(\frac{\pi(S-1)}{S}\right)$ if $S$ is odd. These choices ensure that $B/L_{\rho}$ is contractive by Lemma 8. Basics on surrogate functionals as we need them for this paper are gathered in Section 3.4. Further details on surrogate functionals can be found in [28, 9, 10].

Using elementary properties of the inner product shows that

[TABLE]

where $R(v_{1},\ldots,v_{S})$ is a rest term which is irrelevant when minimizing $P_{\rm surr}$ w.r.t. $u_{1},\ldots,u_{S}$ for fixed $v_{1},\ldots,v_{S}.$ Writing this down in terms of the original system matrix $A$ and the data $f$ yields

[TABLE]

For the quadratic penalty relaxation of the Potts problem, i.e., for minimizing the problem (5) we propose to use the surrogate iteration, i.e. $u^{(n+1)}_{1},\ldots,u^{(n+1)}_{S}$ $\in\operatorname*{argmin}_{u_{1},\ldots,u_{S}}P_{\gamma,\rho}^{\rm surr}(u_{1},\ldots,$ $u_{S},u^{(n)}_{1},\ldots,u^{(n)}_{S}).$ Applied to (15), this surrogate iteration reads

[TABLE]

where $h^{(n)}_{s}$ is given by

[TABLE]

Note that in Section 2.3 below, we derive an efficient algorithm which computes an exact minimizer of (16). Now assume that we are willing to accept a deviation between the $u_{s}$ which is small, i.e.,

[TABLE]

for $\varepsilon>0$ and for indices $s,s^{\prime}$ with $c_{s,s^{\prime}}\neq 0.$ The following algorithm computes a result fulfilling (18).

Algorithm 1**.**

We consider the quadratic penalty relaxed Potts problem (5) and tolerance $\varepsilon$ for the targets $u_{s}$ we are willing to accept. We propose the following algorithm for the relaxed Potts problem (5) (which yields a result with targets $u_{s}$ deviating from each other by at most $\varepsilon/\sqrt{c_{s,s^{\prime}}}$).

• •

Set $\rho$ according to (34), set $L_{\rho}$ according to (13) (or, in the special cases of (7), as below (34) and (13).)

• Initialize $u^{(n)}_{s}$ as discussed in the corresponding paragraph below, (e.g., $u^{(n)}_{s}=0$ for all $s.$)
• •

Iterate until convergence:

[TABLE]

We will see in Theorem 4 that this algorithm converges to a local minimizer of the quadratic penalty relaxation (5) and that the $u_{s}$ are $\varepsilon$-close, i.e., (18) is fulfilled.

The relation between the Potts problem and its quadratic penalty relaxation and

obtaining a feasible solution for the Potts problem (4) from the output of Algorithm 1.

As pointed out above, we show in Theorem 4 that Algorithm 1 produces a local minimizer of the quadratic penalty relaxation (5) of the Potts problem (4) and that the corresponding variables of a resulting solution are close up to an a priori prescribed tolerance. This may in practice be already enough. However, strictly speaking a local minimizer of the quadratic penalty relaxation (5) is not feasible for the Potts problem (4).

We will now explain a projection procedure to derive a feasible solution for the Potts problem (4) from a local minimizer of (5) with nearby variables $u_{s}$ (as produced by Algorithm 1.) Related theoretical results are stated as Theorem 5. In particular, we will see that in case the image operator $A$ is lower bounded, the projection procedure applied to the output of Algorithm 1 yields a feasible point which is close to a local minimizer of the original Potts problem (4).

In order to explain the averaging procedure, we need some notions on partitionings. Recall that a partitioning $\mathcal{P}$ consists of a (finite number of) segments $\mathcal{P}_{i}$ which are pairwise disjoint sets of pixel coordinates whose union equals the image domain $\Omega,$ i.e.,

[TABLE]

Here, we assume that each segment $\mathcal{P}_{i}$ is connected w.r.t. the neighborhood system $a_{1},\ldots,a_{S}$ in the sense that there is a path connecting any two elements in $\mathcal{P}_{i}$ with steps in $a_{1},\ldots,a_{S}.$

We will need the following proposed notion of a directional partitioning. A directional partition w.r.t. a set of $S$ directions $a_{1},\ldots,a_{S}$ consists of a set $\mathcal{I}$ of (discrete) intervals $I$, where each interval $I$ is associated with exactly one of the directions $a_{1},\ldots,a_{S};$ here, an interval $I$ associated with the direction $a_{s}$ has to be of the form $I=\{(i,j)+ka_{s}:k=0,\ldots,K-1\},$ where $K\in\mathbb{N}$ and $I$ belongs to the discrete domain. (For each direction $a_{s}$, the corresponding intervals form an ordinary partition.) We note that Algorithm 1 which produces output $u=(u_{1},\ldots,u_{S}):\Omega\to\mathbb{R}^{s}$ induces a directional partitioning as follows. We observe that each variable $u_{s}$ is associated with a direction $a_{s}.$ For any $s\in\{1,\ldots,S\},$ we let each (maximal) interval of constance of $u_{s}$ be an interval in $\mathcal{I}$ associated with $a_{s}.$

Each partitioning induces a directional partitioning $\mathcal{I}$ by letting the intervals $I$ of $\mathcal{I}$ be the stripes with direction $a_{s}$ obtained from segment $\mathcal{P}_{i}$ for each direction $s=1,\ldots,S$ and each segment $\mathcal{P}_{i},i=1,\ldots,N_{\mathcal{P}}.$ Furthermore, each directional partitioning $\mathcal{I}$ induces a partitioning by the following merging process.

Definition 2**.**

We say that pixels $x,y$ are related, in symbols, $x\sim y$, if there is a path $x_{0}=x,\ldots,x_{N}=y$ connecting $x,y$ in the sense that for any consecutive members $x_{i},x_{i+1},$ $i=1,\ldots,N-1,$ of the path there is an interval $I$ of the directional partitioning $\mathcal{I}$ containing both $x_{i},x_{i+1}.$

The relation $x\sim y$ obviously defines an equivalence relation and the corresponding equivalence classes $\mathcal{P}_{i}$ yield a partitioning on $\Omega.$ We use the symbols

[TABLE]

to denote the mappings assigning a partitioning a directional partitioning and vice versa, respectively.

As a final preparation we consider a function $u=(u_{1},\ldots,u_{S}):\Omega\to\mathbb{R}^{s}$ as produced by Algorithm 1 and a partitioning $\mathcal{P}$ of $\Omega$ and define the following projection to a function $\pi_{\mathcal{P}}(u):\Omega\to\mathbb{R}$ by

[TABLE]

where the symbol $\#\mathcal{P}_{i}$ denotes the number of elements in the segment $\mathcal{P}_{i}.$ Hence the projection $\pi$ defined via (22) averages w.r.t. all components of $u$ and all members of the segment $\mathcal{P}_{i}$ and so produces a piecewise constant function w.r.t. the partitioning $\mathcal{P}.$

Using these notions we propose the following projection procedure.

Procedure 1** (Projection Procedure).**

We consider output $u=(u_{1},\ldots,u_{S}):\Omega\to\mathbb{R}^{s}$ of Algorithm 1 together with its induced directional partitioning $\mathcal{I}.$

Compute the partitioning $\mathcal{P}(\mathcal{I})=\mathcal{P}_{\mathcal{I}}$ induced by the directional partitioning $\mathcal{I}$ as explained above (21).

• 2.

Project $u=(u_{1},\ldots,u_{S}):\Omega\to\mathbb{R}^{s}$ to $\pi_{\mathcal{P}_{\mathcal{I}}}(u)$ using (22) for the partitioning $\mathcal{P}(\mathcal{I})=\mathcal{P}_{\mathcal{I}},$ and return $\pi_{\mathcal{P}_{\mathcal{I}}}(u)$ as output.

We notice that when having a partitioning $\mathcal{P}_{\mathcal{I}}$ solving the normal equation in the space of functions constant on $\mathcal{P}_{\mathcal{I}}$ would be an alternative to the above second step which, however, might be more expensive.

A penalty method for the Potts problem based on

a majorization-minimization approach for its quadratic penalty relaxation.

Intuitively, increasing the parameters $\rho$ during the iterations should tie the $u_{s}$ closer together such that the constraint of (3) should be ultimately fulfilled which results in an approach for the initial Potts problem (2). Recall that $\rho_{s,s^{\prime}}=\rho\ c_{s,s^{\prime}},$ was defined by (6), where the $c_{s,s^{\prime}}$ are nonnegative numbers weighting the constraints. We here increase $\rho$ while leaving the $c_{s,s^{\prime}}$ fixed during this process.

Algorithm 2**.**

We consider the Potts problem (3) in $S$ variables (which is equivalent to (2) as explained above). We propose the following algorithm for the Potts problem (3).

• Let $\rho^{(k)}$ be a strictly increasing sequence (e.g., $\rho^{(k)}=\tau^{k}\rho^{(0)},$ with $\rho_{0},\tau>1$) and $\delta_{k}\to 0$ be a strictly decreasing sequence converging to zero (e.g., $\delta_{k}=\delta_{0}/\tau^{k}.$) Further, let

[TABLE]

where $\sigma_{1}$ is the smallest non-zero eigenvalue of $C^{T}C$ with $C$ given by (49). For the particular choice of coupling given by the left-hand and right hand side of (7) we let

[TABLE]

• •

Initialize $u^{(0)}_{s}:=u^{(0,0)}_{s}$ as discussed in the corresponding paragraph below, (e.g., $u^{(0)}_{s}=0$ for all $s.$)

• Set  $\rho=\rho^{(0)},\ \rho_{s,s^{\prime}}=\rho^{(0)}c_{s,s^{\prime}},\ \delta=\delta_{0},\ k,n=0;$ set $L_{\rho}$ according to (13) (or, in the special cases of (7), as explained below (13))
• A.

While

[TABLE]

do

[TABLE]

and set $n=n+1.$

• B.

Set

[TABLE]

set $k=k+1,n=0,$ and let $\rho=\rho^{(k)},\rho_{s,s^{\prime}}=\rho^{(k)}\ c_{s,s^{\prime}},\delta=\delta_{k};$ set $L_{\rho}$ according to (13) (or, in the special cases of (7), as below (13)) and goto A.

This approach is inspired by [60] which considers quadratic penalty methods in the sparsity context. There, the authors are searching for a solution with only a few nonzero entries. The corresponding prior is separable. In contrast to this work, the present work considers a non-separable prior.

Initialization.

Although the initialization of Algorithm 1 and of Algorithm 2 is not relevant for its convergence properties (cf. Section 3), the choice of the initialization influences the final result. (Please note that this also might happen for convex but not strictly convex problems.) We discuss different initialization strategies. The simplest choice is the all-zero initialization $(u_{1}^{(0)},...,u_{s}^{(0)})=(0,...,0).$ Likewise, one can select the right hand side of the normal equations of the underlying least squares problem, that is $A^{T}f$. A third reasonable choice is the solution of the normal equation itself or an approximation of it. Using an approximation might in particular be reasonable to get a regularized approximation of the normal equation. A possible strategy to obtain such a regularized initialization is to apply a fixed number of Landweber iterations [54] or of the conjugate gradient method to the underlying least square problem. (In our experiments, we initialized Algorithm 1 with the result of 1000 Landweber iterations and Algorithm 2 with $A^{T}f$.)

2.3 A non-iterative algorithm for minimizing the Potts subproblem (16)

Both proposed algorithms require solving the Potts subproblem (16) in the backward step, see (19),(26). We first observe that (16) can be solved for each of the $u_{s}$ separately. The corresponding $s$ minimization problems are of the prototypical form

[TABLE]

with given data $f$, the jump penalty $\gamma^{\prime}_{s}=\tfrac{\gamma\omega_{s}}{L_{\rho}^{2}}>0$ and the direction $a_{s}\in\mathbb{Z}^{2}$. As a next step, we see that (28) decomposes into univariate Potts problems for data along the paths in $f$ induced by $a_{s}$, e.g., for $a_{s}=e_{1}$ those paths correspond to the rows of $f$ and we obtain a minimizer $u_{s}^{\ast}$ of (28) by determining each of its rows individually. The univariate Potts problem amounts to minimizing

[TABLE]

where the data $g$ is given by the restriction of $f$ to the pixels in $\Omega$ of the form $v+a_{s}z,$ for $z\in\mathbb{Z}$, i.e., $g(z)=f(v+a_{s}z)$.

Here the offset $v$ is fixed when solving each univariate problem, but varied afterwards to get all lines in the image with direction $a_{s}.$ The target to optimize is denoted by $x\in{\mathbb{R}^{n}}$ and, in the resulting univariate situation, $\|\nabla x\|_{0}=|\{i:x_{i}\neq x_{i+1}\}|$ denotes the number of jumps of $x$.

It is well-known that the univariate direct problem (29) has a unique minimizer. Further these particular problems can be solved exactly by dynamic programming [35, 92, 18, 62, 63] which we briefly describe in the following. For further details we refer to [35, 82]. Assume we have computed minimizers $x^{l}$ of (29) for partial data $(g_{1},...,g_{l})$ for each $l=1,...,r$, $r<n$. Then the minimum value of (29) for $(g_{1},...,g_{r+1})$ can be found by

[TABLE]

where we let $x^{0}$ be the empty vector, $P^{\mathrm{id,1d}}_{\gamma}(x^{0})=-\gamma$ and $\mathcal{E}^{l:r+1}$ be the quadratic deviation of $(g_{l},...,g_{r+1})$ from its mean. By denoting the minimizing argument in (30) by $l^{\ast}$ the minimizer $x^{r+1}$ is given by

[TABLE]

where $\mu_{[l^{\ast},r]}$ is the mean value of $(g_{l^{\ast}},...,g_{r})$. Thus, we obtain a minimizer for full data $g$ by successively computing $x^{l}$ for each $l=1,...,n$. By precomputing the first and second moments of data $g$ and storing only jump locations the described method can be implemented in $\mathcal{O}(n^{2})$, [35]. Another way to achieve $\mathcal{O}(n^{2})$ is based on the QR decomposition of the design matrix by means of Givens rotations, see [82]. Furthermore, the search space can be pruned to speed up computations [47, 83].

We briefly describe the extensions of the above scheme necessary to approach (29) for vector valued-data $g\in\mathbb{R}^{n\times C}$ (e.g., the row of a color image). In this situation, the symbol $\mathcal{E}^{l:r+1}$ in (30) denotes the sum of the quadratic deviations of $(g_{l},\ldots,g_{r+1})$ from its channel-wise means. Further, $\mu_{[l^{\ast},r]}\in\mathbb{R}^{C}$ in (31) is the vector of channel-wise means of the data $(g_{l^{\ast}},\ldots,g_{r})$. On the computational side, the first and second moments of each channel have to be precomputed separately. It is worth mentioning that the theoretical computational costs of the described method grows only linearly in the number of channels [83]. Thus, the proposed algorithm can be efficiently applied to vector-valued images with a high-dimensional codomain.

3 Analysis

3.1 Analytic results

In the course of the derivation of the proposed algorithms above, we consider the quadratic penalty relaxation (5) of the multivariate Potts problem. Although it is more straight-forward to access algorithmically via our approach, we first note that this problem is still NP hard (as is the original problem).

Theorem 3**.**

Finding a (global) minimizer of the quadratic penalty relaxation (5) of the multivariate Potts problem is an NP hard problem.

The proof is given in Section 3.3 below. In Section 2.2 we have proposed Algorithm 1 to approach the quadratic penalty relaxation of the multivariate Potts problem. We show that the proposed algorithm converges to a local minimizer and that a feasible point of the original multivariate Potts problem is nearby.

Theorem 4**.**

We consider the iterative Potts minimization Algorithm 1 for the quadratic penalty relaxation (5) of the multivariate Potts problem:

1. i.

Algorithm 1 computes a local minimizer of the quadratic penalty relaxation (5) of the multivariate Potts problem for any starting point. The convergence rate is linear. 2. ii.

We have the following relation between local minimizers ${\mathcal{L}}$, global minimizers ${\mathcal{G}}$ and the fixed points $\mathrm{Fix}(\mathbb{I})$ of the iteration of Algorithm 1,

[TABLE] 3. iii.

Assume a tolerance $\varepsilon$ we are willing to accept for the distance between the $u_{s},$ i.e.,

[TABLE]

Running Algorithm 1 with the choice of the parameter $\rho$ by

[TABLE]

(where $\sigma_{1}$ is the smallest non-zero eigenvalue of $C^{T}C$ with $C$ given by (49); for the particular choice of the coupling given by (7), $\sigma_{1}=S$ and $\sigma_{1}=(2-2\cos(2\pi/S)),$ respectively) yields a local minimizer of the quadratic penalty relaxation (5) such that the $u_{s}$ are close up to $\varepsilon,$ i.e., (33) is fulfilled.

The proof is given in Section 3.5 below. A solution of Algorithm 1 is not a feasible point for the initial Potts problem (3). However, we see below that it produces a $\delta$-approximative solution $u^{\ast}$ in the sense that there is $\mu^{*}$ and a partitioning $\mathcal{P}^{\ast}$ such that

[TABLE]

where $L(\mu^{\ast})$ is given by (53) below. In this context note that the conditions for a local minimizer are given by $\sum_{s,s^{\prime}}c_{s,s^{\prime}}\|u^{\ast}_{s}-u^{\ast}_{s^{\prime}}\|^{2}_{2}=0$ and the Lagrange multiplier condition $L(\mu^{\ast})=0.$ So (35) intuitively means that both the constraint and the Lagrange multiplier condition are approximately fulfilled for the partitioning induced by $u^{\ast}$.

Further, given a solution of Algorithm 1 we find a feasible point for the Potts problem (3) (or, equivalently,(2)) which is nearby as detailed in the following theorem.

Theorem 5**.**

We consider the iterative Potts minimization Algorithm 1 for the quadratic penalty relaxation (5) in connection with the (non-relaxed) Potts problem (3).

1. i.

Algorithm 1 produces an approximative solution in the sense of (35) of the Potts problem (3). 2. ii.

The projection procedure (Procedure 1) proposed in Section 2.2 applied to the solution $u^{\prime}=(u^{\prime}_{1},\ldots,u^{\prime}_{S})$ of Algorithm 1 produces a feasible image $\hat{u}$ (together with a valid partitioning) for the Potts problem (3) which is close to $u^{\prime}$ in the sense that

[TABLE]

where $\varepsilon=\max_{s,s^{\prime}}\|u^{\prime}_{s}-u^{\prime}_{s^{\prime}}\|$ quantifies the deviation between the $u_{s}.$ Here $C_{1}=\#\Omega/4,$ where the symbol $\#\Omega$ denotes the number of elements in $\Omega.$ If the imaging operator $A$ is lower bounded, i.e., there is a constant $c>0$ such that $\|Au\|\geq c\|u\|$, a local minimizer $u^{\ast}$ of the Potts problem (3) is nearby, i.e.,

[TABLE]

where

[TABLE]

The proof of Theorem 5 can be found at the end of Section 3.4, where most relevant statements are already shown in Section 3.3. Theorem 5 theoretically underpins the fact that, on the application side, we may use Algorithm 1 for the Potts problem (3) (accepting some arbitrary small tolerance we may fix in advance).

In addition, in Section 2.2, we have proposed Algorithm 2 to approach the Potts problem (3). We first show that Algorithm 2 is well-defined.

Theorem 6**.**

Algorithm 2 is well-defined in the sense that the inner iteration governed by (25) terminates, i.e., for any $k\in\mathbb{N},$ there is $n\in\mathbb{N}$ such that the termination criterium given by (25) holds.

The proof of Theorem 6 is given in Section 3.6. Concerning the convergence properties of Algorithm 2 we obtain the following results.

Theorem 7**.**

We consider the iterative Potts minimization algorithm (Algorithm 2) for the Potts problem (3).

• •

Any cluster point of the sequence $u^{(k)}$ is a local minimizer of the Potts problem (3) (which implicitly implies that the components of each limit $u^{\ast}$ are equal, i.e., $u_{s}^{\ast}=u_{s^{\prime}}^{\ast}$ for all $s,s^{\prime}.$)

• •

If $A$ is lower bounded, the sequence $u^{(k)}$ produced by Algorithm 2 has a cluster point and the produced cluster points are local minimizers of the Potts problem (3).

The proof of Theorem 7 can be found in Section 3.6.

3.2 Estimates on operator norms and Lagrange multipliers

Lemma 8**.**

The spectral norm of the block matrix $B$ given by (8) fulfills

[TABLE]

For the particular choice of constant $\rho_{s,s^{\prime}}=\rho$ (independent of $s,s^{\prime}$) as on the left-hand side of (7) we have the improved estimate

[TABLE]

For only coupling neighboring $u_{s}$ with the same constant $\rho$, i.e., the right-hand coupling of (7), we have

[TABLE]

Proof.

We decompose the matrix $B$ according to $B=\begin{pmatrix}S^{-1/2}\tilde{A}\\ \tilde{P}\end{pmatrix}.$ Here $\tilde{A}$ denotes an $S\times S$- block diagonal matrix with each diagonal entry being equal to $A,$ where $A$ is the matrix representing the forward/imaging operator; see (8). The matrix $\tilde{P}$ is given as the lower ${S\choose 2}\times S$- block in (8) which represents the soft constraints.

Using this decomposition of $B,$ we may decompose the symmetric and positive (semidefinite) matrix $B^{T}B$ according to

[TABLE]

where $\tilde{A}^{T}\tilde{A}$ is an $S\times S$- block diagonal matrix with each diagonal entry being equal to $A^{T}A,$ and $\tilde{P}^{T}\tilde{P}$ is an $S\times S$- block diagonal matrix with block entries given by

[TABLE]

with $\rho_{l,k}:=\rho_{k,l}$ for $l>k.$ Using Gerschgorin’s Theorem (see for instance [81]), the eigenvalues of $\tilde{P}$ are contained in the union of the balls with center $x_{r}=\sum_{k=1,k\neq r}^{S}\rho_{r,k}$ and radius $x_{r}=\sum_{k=1,k\neq r}^{S}|-\rho_{r,k}|.$ These balls are all contained in the larger ball with center [math] and radius $2\cdot\max_{r}x_{r}.$ This implies the general estimate (39).

For seeing (40) we decompose an argument $u=(u_{1},\ldots,u_{S})$ according to $u=\bar{u}+u^{0}$ with an “average” part $\bar{u}=(\tfrac{1}{S}\sum_{i=1}^{S}u_{i},\ldots,\tfrac{1}{S}\sum_{i=1}^{S}u_{i})$ and $u^{0}:=u-\bar{u}$ such that $u^{0}$ has average $0,$ i.e., $\sum_{i=1}^{S}u^{0}_{i}=0,$ where [math] denotes the vector containing only zero entries here. In the situation of (40), the matrix $\tilde{P}^{T}\tilde{P}$ has the form $\tilde{P}^{T}\tilde{P}=\rho(S\cdot I-(1,\ldots,1)(1,\ldots,1)^{T})$ We have $\tilde{P}^{T}\tilde{P}\bar{u}=0.$ Further, $\tilde{P}^{T}\tilde{P}u^{0}=\rho Su^{0}.$ Hence, the largest modulus of an eigenvalue of $\tilde{P}^{T}\tilde{P}$ equals $\rho S$ which in turn shows the estimate (40).

For seeing (41), we notice that in case of (41), the matrix $\tilde{P}^{T}\tilde{P}$ has cyclic shift structure with three nonzero entries in each line. The discrete Fourier matrix w.r.t. the cyclic group of order $S$ diagonalizes $\tilde{P}^{T}\tilde{P}.$ The corresponding eigenvalues are given by $\lambda_{k}=\rho\left(2-2\cos\left(2\pi\frac{k}{S}\right)\right),$ where $k=0,\ldots,S-1.$ The largest modulus of an eigenvalue is thus given by $4\ \rho,$ if $S$ is even, and by $\rho\cdot\left(2-2\cos\left(\frac{\pi(S-1)}{S}\right)\right).$ ∎

Note that the problem of estimating the operator norm of $B$ in (39) involves computing the operator norm of $\tilde{P}$ given by (43). This problem is intimately related to computing the spectral norm of the Laplacian of a corresponding weighted graph (e.g., [80, 38]), in particular, we conclude from this link that the general estimate (39) is sharp in the sense that the factor of $2$ in front of the sum cannot be made smaller. This is because, for a general graph, the spectral radius of the (normalized) Laplacian has spectral norm smaller than two and this factor of two is sharp; cf. [80, 38].

We recall that we have introduced the concept of a directional partitioning $\mathcal{I}$ and discussed its relation with the concept of a partitioning near (21) above. For a function $f:\Omega\to\mathbb{R}^{S}$ (representing its $S$ component functions $f_{1},\ldots,f_{S}:\Omega\to\mathbb{R}$) defined on a grid $\Omega,$ we consider the orthogonal projection $P_{\mathcal{I}}$ associated with a directional partition $\mathcal{I}$ by first sorting the intervals $\mathcal{I}$ into $\mathcal{I}_{1},\ldots,\mathcal{I}_{S}$ according to their associated directions $a_{s},$ $s=1,\ldots,S,$ and then letting

[TABLE]

i.e., the function $P_{\mathcal{I}_{s}}f_{s}$ on the interval $I$ is given as the arithmetic mean of $f_{i}$ on the interval $I$ for all intervals $I\in\mathcal{I}_{s},$ and for all $s=1,\ldots,S.$ Here, the symbol $\#I$ denotes the number of elements in $I.$ We note that $P_{\mathcal{I}}$ defines an orthogonal projection on the corresponding $\ell^{2}$ space of discrete functions $f:\Omega\to\mathbb{R}^{S}$ with the norm $\|f\|^{2}=\sum_{s,i}|(f_{s})_{i}|^{2}$ where $i$ iterates through all the indices of $f_{s}.$

We consider a partitioning $\mathcal{P}$ of $\Omega,$ its induced directional partitioning $\mathcal{I}^{\mathcal{P}}$ w.r.t. a set of $S$ directions $a_{1},\ldots,a_{S},$ and the subspace

[TABLE]

of functions which are constant on the intervals of the induced directional partitioning $\mathcal{I}^{\mathcal{P}}$ (which equal the image of the orthogonal projection $P_{\mathcal{I}^{\mathcal{P}}}.$)

Functions $g:\Omega\to\mathbb{R}$ which are piecewise constant w.r.t. a partitioning $\mathcal{P},$ i.e., they are constant on each segment $\mathcal{P}_{i}$ are in one-to-one correspondence with the linear subspace $\mathcal{B}^{\mathcal{P}}$ of $\mathcal{A}^{\mathcal{P}}$ given by

[TABLE]

as shown by the following lemma.

Lemma 9**.**

There is a one-to-one correspondence between the linear space of piecewise constant mappings w.r.t. the partitioning $\mathcal{P},$ and the subspace $\mathcal{B}^{\mathcal{P}}$ of $\mathcal{A}^{\mathcal{P}}$ via the mapping $\iota:g\mapsto(g,\ldots,g).$

Proof.

Let $g$ be a piecewise constant mapping w.r.t. the partitioning $\mathcal{P},$ then $(g,\ldots,g)$ is constant on each interval $I$ of the induced directional partitioning $\mathcal{I}^{\mathcal{P}},$ and $(g,\ldots,g)\in\mathcal{B}^{\mathcal{P}}.$ This shows that $\iota$ is well-defined in the sense that its range is contained in $\mathcal{B}^{\mathcal{P}}.$ Obviously, $\iota$ is an injective linear mapping so that it remains to show that any $f\in\mathcal{B}^{\mathcal{P}}$ is the image under $\iota$ of some $g:\Omega\to\mathbb{R}$ which is piecewise constant w.r.t. the partitioning $\mathcal{P}.$ To this end, let $f\in\mathcal{B}^{\mathcal{P}}.$ By definition, $f$ has the form $f=(g,\ldots,g)$ for some $g:\Omega\to\mathbb{R}.$ Now, towards a contradiction, assume there is a segment $\mathcal{P}_{i}$ and points $x,y\in\mathcal{P}_{i}$ with $g(x)\neq g(y).$ Since there is a path $x_{0}=x,\ldots,x_{N}=y$ connecting $x,y$ in $\mathcal{P}_{i}$ with steps in $a_{1},\ldots,a_{S},$ we have that for any $i$ there is an interval $I$ in the induced partitioning $\mathcal{I}^{\mathcal{P}}$ containing $x_{i}$ together with $x_{i+1}.$ Since $g$ is constant on each $I$ in $\mathcal{I}^{\mathcal{P}}$ we get $g(x_{i})=g(x_{i+1})$ for all $i$ which implies $g(x)=g(y).$ This contradicts our assumption and shows the lemma. ∎

Using the identification given by Lemma 9, we define, for a given partitioning $\mathcal{P},$ the projection $Q_{\mathcal{P}}$ onto $\mathcal{B}^{\mathcal{P}}$ by

[TABLE]

i.e., we average w.r.t. the segment and to all component functions as given by (22). Since the components of $Q_{\mathcal{P}}f$ are all identical, we will not distinguish $Q_{\mathcal{P}}$ and $\pi_{\mathcal{P}}$ in the following. This means that we also use the symbol $Q_{\mathcal{P}}f$ to denote the scalar-valued function which is piecewise constant on the partitioning $\mathcal{P}.$

On $\mathcal{A}^{\mathcal{P}},$ we consider the problem

[TABLE]

i.e., given the directional partitioning we are searching for a solution which belongs to $\mathcal{B}^{\mathcal{P}}.$ Here, $C$ denotes the matrix

[TABLE]

where the $c_{s,s^{\prime}}$ are as in (5); if $c_{s,s^{\prime}}=0,$ the corresponding line is removed from the constraint matrix $C.$ For the special choices of (7), we have

[TABLE]

which reflects the constraints $u_{1}=\ldots=u_{S}.$ We recall that $\mu_{\mathcal{P}}$ is a Lagrange multiplier of the problem in (48) if

[TABLE]

We note that for quadratic problems such as in (48) Lagrange multipliers always exist [7]. We have that

[TABLE]

or, in other form,

[TABLE]

where $\tilde{A}$ is the block diagonal matrix with constant entry $A$ on each diagonal component, $\tilde{f}$ is a block vector of corresponding dimensions with entry $f$ in each component, and $u_{\mathcal{P}}^{\ast}$ is a minimizer of the constraint problem in $\mathcal{B}^{\mathcal{P}}$. We note that the last equality $C^{T}\mu_{\mathcal{P}}=P_{\mathcal{I}^{\mathcal{P}}}C^{T}\mu_{\mathcal{P}}$ in (52) holds since the left-hand side of (52) is contained in the image of $P_{\mathcal{I}^{\mathcal{P}}}.$

Lemma 10**.**

We consider a partitioning $\mathcal{P}$ of the discrete domain $\Omega,$ and the corresponding problem (48). There is a Lagrange multiplier $\mu_{\mathcal{P}}$ for (48) with

[TABLE]

Here, $\sigma_{1}$ is the smallest non-zero eigenvalue of $C^{T}C$ with $C$ given by (49). For the particular choice of $C$ given by the left-hand side of (50) we have

[TABLE]

and, for the particular choice of $C$ given by the right-hand side of (50) we have

[TABLE]

(e.g., for $S=4,$ an eight neigborhood, $2-2\cos(2\pi/S)=\sigma_{1}=2.$) In particular, the right-hand side and the constants in all these estimates are independent of the particular partitioning $\mathcal{P}$.

Proof.

For any minimizer $u_{\mathcal{P}}^{\ast}$ of the constraint problem in $\mathcal{B}^{\mathcal{P}},$ we have that

[TABLE]

where we recall that $\tilde{A}$ is the block diagonal matrix with constant entry $A,$ and $\tilde{f}$ is a block vector with entry $f$ in each component. The first inequality is a consequence of the fact that $P_{\mathcal{I}^{\mathcal{P}}}$ is an orthogonal projection. The second inequality may be seen by evaluating the term for the constant zero function (which always belongs to $\mathcal{B}^{\mathcal{P}}$ ) as a candidate and by noting that $\|A^{T}\|=\|A\|.$

Using (52), we have $\|C^{T}\mu_{\mathcal{P}}\|\leq\tfrac{2}{\sqrt{S}}\|A\|\|f\|.$ Choosing $\mu_{\mathcal{P}}$ in the complement of the zero space of $C^{T},$ we get

[TABLE]

We observe that finding the infimum in (58) corresponds to finding the square root of the smallest nonzero eigenvalue of $C^{T}C.$ This is because (i) the nonzero eigenvalues of $C^{T}C$ equal the nonzero eigenvalues of $CC^{T},$ i.e.,

[TABLE]

where $\sigma_{1}$ is the smallest non-zero eigenvalue of $C^{T}C$. Further, (ii) for $x\in\left(\operatorname{ker}\left(C^{T}\right)\right)^{\perp},$ $\|C^{T}x\|^{2}=\langle x,CC^{T}x\rangle\geq\min\left\{\sigma:\sigma\in\mathrm{spectrum}(CC^{T})\setminus\{0\}\right\}\|x\|^{2}.$ Hence, using (59) in (58) we get that $\|C^{T}\mu_{\mathcal{P}}\|\geq\sqrt{\sigma_{1}}\|\mu_{\mathcal{P}}\|,$ and together with (52) and (57), we obtain

[TABLE]

which shows (54).

Now we consider the particular choice of $C$ given by the left-hand side of (50). Similar to the derivation in (43), we have that $C^{T}C=S\cdot I-(1,\ldots,1)(1,\ldots,1)^{T}).$ Further, the constants constitute the kernel of $C^{T}C$ and any vector $u$ in its orthogonal complement is mapped to $Su.$ Hence, $\sigma_{1}=S$ which shows (55).

Finally, we consider the particular choice of $C$ given by the right-hand side of (50). As already explained in the proof of Lemma 8 the discrete Fourier transform shows that the corresponding eigenvalues are given by $\lambda_{k}=\rho\left(2-2\cos\left(2\pi\frac{k}{S}\right)\right),$ where $k=0,\ldots,S-1.$ The smallest nonzero eigenvalue is thus given by $2-2\cos(2\pi/S).$ This shows (56) which completes the proof of the lemma. ∎

3.3 The quadratic penalty relaxation of the Potts problem

and its relation to the Potts problem

In this subsection we reveal some relations between the Potts problem and its quadratic penalty relaxation; in particular, we show Theorem 3 and parts of Theorem 5. We start out to show that the quadratic penalty relaxation of the Potts problem is NP hard which was formulated as Theorem 3.

Proof of Theorem 3.

We consider the quadratic penalty relaxation (5) of the multivariate Potts problem in its equivalent form (11) which reads

[TABLE]

with $B$ and $g$ given by (8) and $D$ given by (9). We serialize $u:(u_{1},\ldots,u_{S}):\Omega\to\mathbb{R}^{S}$ into a function $\hat{u}:X\to\mathbb{R}$ with $X\subset\mathbb{Z}$ being a discrete interval of size $S\#\Omega$ as follows: For $u_{s},$ we consider the discrete lines in the image with direction $a_{s}$ and interpret $u$ on these lines as a vector; then we concatenate these vectors starting with the one corresponding to the leftmost upper line to obtain a vector of length $\#\Omega;$ for each $s,$ we obtain such a vector and we again concatenate these vectors starting with index $s=1,2,\ldots$ to obtain the resulting object which we denote by $\hat{u}.$ Using this serialization we may arrange $B,g$ and $D$ accordingly to obtain the univariate Potts problem

[TABLE]

is a weight vector, $\omega\nabla\hat{u}$ denotes pointwise multiplication, and $\hat{B},\hat{g}$ are the matrix and the vector corresponding to $B,g$ w.r.t. the serialization. The weight vector may be zero which in particular happens at the line breaks, i.e., those indices where two vectors have been concatenated in the above procedure. More precisely, constant data induce a directional segmentation on $\Omega$ and the image of the directional segmentation under the above serialization procedure induces a partitioning of the unvariate domain $D;$ precisely between these segments, the weight vectors equals zero. Now, for each segment $[d_{1},\ldots,d_{r}]$ in $D,$ we transform the basis $\delta_{d_{1}},\ldots,\delta_{d_{r}}$ to the basis $\delta_{d_{2}}-\delta_{d_{1}},\ldots,\delta_{d_{r}}-\delta_{d_{r-1}},\tfrac{1}{r}\sum_{l=1}^{r}\delta_{d_{l}}$ obtained by neighboring differences and the average. As a result (which is in detail elaborated in [84]), we obtain a problem of the form

[TABLE]

which is a sparsity problem and which is known to be NP hard; see, for instance, [84]. This shows the assertion. ∎

We next characterize the local minimizers of the relaxed Potts problem (5) and of the Potts problem (2).

Lemma 11**.**

A local minimizer $u=(u_{1},\ldots,u_{S})$ of the quadratic penalty relaxation (5) is characterized as follows: let $\mathcal{I}$ be the directional partitioning induced by the minimizer $u,$ and $\mathcal{P}=\mathcal{P}_{\mathcal{I}}$ be the induced partitioning, then $u$ is a minimizer of the problem

[TABLE]

Conversely, if $u$ minimizes (62) on $\mathcal{A}^{\mathcal{P}},$ then $u$ is a minimizer of the relaxed Potts problem (5).

Proof.

Let $u=(u_{1},\ldots,u_{S})$ be a local minimizer of the quadratic penalty relaxation (5). Hence, there is a neigborhood $\mathcal{U}$ of $u$ such that, for any $v\in\mathcal{U},$ $P_{\gamma,\rho}(v)\geq P_{\gamma,\rho}(u).$ Now if $v\in\mathcal{A}^{\mathcal{P}}$ and $\|v-u\|$ is small, then $\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}u_{s}\,\right\|_{0}$ $=\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}v_{s}\,\right\|_{0}$ which implies that

[TABLE]

This shows that $u$ minimizes (62). Conversely, we assume that $u$ minimizes (62). If the directional partitioning $\mathcal{I}^{\prime}$ induced by $u$ is coarser than $\mathcal{I}$ consider the coarser directional partitioning $\mathcal{I}^{\prime}$ instead of $\mathcal{I}.$ Let the maximum norm of $h=(h_{1}\ldots,h_{S})$ be smaller than the height of the smallest jump of $u,$ then, for $u+h,$

[TABLE]

If inequality holds in (64), the continuity of $F_{\rho}$ implies that $F_{\rho}(u+h)\geq F_{\rho}(u)-\varepsilon$ for small enough $h$ and arbitrary $\varepsilon.$ Hence,

[TABLE]

if we choose $\varepsilon$ small enough. If equality holds in (64) we have that $u+h\in\mathcal{A}^{\mathcal{P}}$ which implies $F_{\rho}(u)\leq F_{\rho}(u+h)$ since $u$ is a minimizer of $F_{\rho}$ on $\mathcal{A}^{\mathcal{P}}.$ This in turn implies $P_{\gamma,\rho}(u)\leq P_{\gamma,\rho}(u+h)$ by the assumed equality in (64). Together, in any case, $P_{\gamma,\rho}(u)\leq P_{\gamma,\rho}(u+h)$ for any small perturbation $h.$ This shows that $u$ is a local minimizer of $P_{\gamma,\rho}$ which completes the proof. ∎

Lemma 12**.**

We consider a function $u^{\ast}:\Omega\to\mathbb{R}$ and its induced partitioning $\mathcal{P}.$ Then $u$ is a local minimizer of the Potts problem (2), if and only if $(u^{\ast},\ldots,u^{\ast})$ minimizes (48) w.r.t. $\mathcal{P}.$

Proof.

Since the proof of this statement is very similar to the proof of Lemma 11 we keep it rather short and refer to the proof of Lemma 11 if more explanation is necessary. Let $u$ be a minimizer of (2) which is equivalent to $\bar{u}=(u,\ldots,u)$ being a minimizer of (4). There is a neigborhood $\mathcal{U}$ of $\bar{u}$ such that, for any $\bar{v}=(v,\ldots,v)\in\mathcal{U},$ $P_{\gamma}(v)\geq P_{\gamma}(u).$ For $\bar{v}\in\mathcal{B}^{\mathcal{P}}$ with small $\|\bar{v}-\bar{u}\|,$ we have $\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}u\,\right\|_{0}$ $=\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}v\,\right\|_{0}.$ Hence, by the definition of $P_{\gamma}$ in (4) $\|Au-f\|_{2}^{2}\leq\|Av-f\|_{2}^{2}$ which shows that $(u^{\ast},\ldots,u^{\ast})$ minimizes (48).

Conversely, let $\bar{u}=(u,\ldots,u)$ be a minimizer of (48) with the partitioning $\mathcal{P}$ induced by $u$. For $\bar{h}=(h,\ldots,h)$ with absolute value smaller than the minimal height of a jump of $u,$ we have the estimate $\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}(u+h)\,\right\|_{0}$ $\geq\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}u\,\right\|_{0}.$ If inequality holds in this estimate, the continuity of $F_{\rho}$ implies that $\|A(u+h)-f\|_{2}^{2}\geq\|Au-f\|_{2}^{2}-\varepsilon$ for small enough $h$ and arbitrary $\varepsilon.$ Hence, $P_{\gamma}(\bar{u})\leq\|A(u+h)-f\|_{2}^{2}-\gamma\min_{s}\omega_{s}+\gamma\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}(u+h)\,\right\|_{0}+\varepsilon$ $\leq P_{\gamma,\rho}(\bar{u}+\bar{h})$ if $\varepsilon$ is small. If equality holds above, i.e, $\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}(u+h)\,\right\|_{0}$ $=\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}u\,\right\|_{0},$ then $\bar{u}+\bar{h}\in\mathcal{B}^{\mathcal{P}}$ which implies that $\|Au-f\|_{2}^{2}\leq\|A(u+h)-f\|_{2}^{2}$ since $\bar{u}$ is a minimizer of the corresponding function on $\mathcal{B}^{\mathcal{P}}.$ As a consequence $P_{\gamma}(\bar{u})\leq P_{\gamma}(\bar{u}+\bar{h})$ for any small perturbation $h.$ This shows that $u$ is a local minimizer of $P_{\gamma}$ which completes the proof. ∎

Proposition 13**.**

Any local minimizer of the quadratic penalty relaxation (5) is an approximate local minimizer in the sense of (35) of the Potts problem (3).

Proof.

By Lemma 11, a local minimizer $u=(u_{1},\ldots,u_{S})$ of the quadratic penalty relaxation (5) is a minimizer of the problem (62). Let us thus consider a local minimizer $u$ of (5) with induced partitioning $\mathcal{P}=\mathcal{P}_{\mathcal{I}}.$ Since $u$ minimizes (62), we have

[TABLE]

since the gradient projected to $\mathcal{A}^{\mathcal{P}}$ equals zero for any local minimizer of the restricted problem on the subspace $\mathcal{A}^{\mathcal{P}}.$ (The notation is chosen as in (53) above.) We define $\mu$ by $\mu=\rho CP_{\mathcal{I}}u$ and obtain

[TABLE]

by (66). It remains to show that $\|Cu\|$ becomes small. To this end, we observe that, by Lemma 14, for arbitrary $v=(v_{1},\ldots,v_{S})\in\mathcal{A}^{\mathcal{P}}$, $\|Cv\|=\|CP_{\mathcal{I}}v\|\leq\tfrac{1}{\rho}\|\mu^{\ast}\|+\sqrt{\tfrac{F_{\rho}(v)-\min_{x\in\mathcal{A}^{\mathcal{P}}}F_{\rho}(x)}{\rho}},$ where $\mu^{\ast}$ is an arbitrary Lagrange multiplier of (48). Plugging in the minimizer $u$ for $v$ yields $\|Cu\|<\tfrac{1}{\rho}\|\mu^{\ast}\|.$ Thus letting $\delta=\tfrac{1}{\rho}\|\mu^{\ast}\|,$ we have

[TABLE]

and $L(\mu)=0$ by (67) which by (35) shows the assertion and completes the proof. ∎

For the proof of Proposition 13 as well as in the following we need the next lemma. Similar statements are [53, Proposition 13] and [60, Lemma 2.5]. However, since there are differences concerning the precise estimate in these references, and the setup here is slightly different, we provide a brief proof here for the readers convenience.

Lemma 14**.**

Let $\mathcal{P}$ be a partitioning and $\mathcal{I}=\mathcal{I}_{\mathcal{P}}$ be the corresponding induced partitioning. For arbitrary $v=(v_{1},\ldots,v_{S})\in\mathcal{A}^{\mathcal{P}}$,

[TABLE]

where $\mu^{\ast}$ is an arbitrary Lagrange multiplier of (48).

Proof.

By [53, Corollary 2], we have for arbitrary $v=(v_{1},\ldots,v_{S})\in\mathcal{A}^{\mathcal{P}}$ that

[TABLE]

Then,

[TABLE]

For the first inequality we wrote down the definition of $F_{\rho}$ and restricted the set with respect to which the minimum is formed which results in a potentially larger function value. For the second inequality we notice that, for $(y,\ldots,y)\in\mathcal{B}^{\mathcal{P}}$, we have $C(y,\ldots,y)=0,$ and for the last inequality we employed (70). Now, writing $z^{2}-\frac{\|\mu^{\ast}\|}{\|\rho\|}z$ $=z^{2}-\frac{\|\mu^{\ast}\|}{\rho}z+\left(\tfrac{\|\mu^{\ast}\|}{2\rho}\right)^{2}-\left(\tfrac{\|\mu^{\ast}\|}{2\rho}\right)^{2}$ $=(z-\tfrac{\|\mu^{\ast}\|}{2\rho})^{2}-\left(\tfrac{\|\mu^{\ast}\|}{2\rho}\right)^{2}$ and plugging this into (3.3) with $z:=\|Cv\|$ yields

[TABLE]

and hence

[TABLE]

where the last inequality is a consequence of the fact that the unit ball w.r.t. the $\ell^{1}$ norm is contained in the the unit ball w.r.t. the $\ell^{2}$ norm. As a consequence, $\|Cv\|\leq\sqrt{\tfrac{F_{\rho}(v)-\min_{x\in\mathcal{A}^{\mathcal{P}}}F_{\rho}(x)}{\rho}}+\tfrac{\|\mu^{\ast}\|}{2\rho}+\tfrac{\|\mu^{\ast}\|}{2\rho}$ which completes the proof. ∎

Next, we see that for any local minimizer of the quadratic penalty relaxation (5), we can find a nearby feasible point using the projection procedure (Procedure 1) proposed in Section 2.2. Further, if the imaging operator $A$ is lower bounded, we find a nearby minimizer.

Proposition 15**.**

Procedure 1 applied to a local minimizer $u^{\prime}=(u^{\prime}_{1},\ldots,u^{\prime}_{S})$ of the quadratic penalty relaxation (5) produces a feasible image $\hat{u}$ (together with a valid partitioning) for the Potts problem (3) which is close to $u^{\prime}$ in the sense that

[TABLE]

where $\varepsilon=\max_{s,s^{\prime}}\|u^{\prime}_{s}-u^{\prime}_{s^{\prime}}\|$ quantifies the deviation between the $u_{s}.$ Here $C_{1}=\#\Omega/4,$ where the symbol $\#\Omega$ denotes the number of elements in $\Omega.$

If the imaging operator $A$ is lower bounded, i.e., there is a constant $c>0$ such that $\|Au\|\geq c\|u\|$, a local minimizer $u^{\ast}$ of the Potts problem (3) is nearby, i.e.,

[TABLE]

where

[TABLE]

Proof.

We denote the directional partitioning induced by $u^{\prime}$ by $\mathcal{I}$ and the corresponding induced partitioning by $\mathcal{P}=\mathcal{P}_{\mathcal{I}}.$ We note that Procedure 1 applied to $u^{\prime}$ precisely produces

[TABLE]

with the projection $Q_{\mathcal{P}}$ given by (47). We first note, that the average $(\bar{u})_{ij}=\tfrac{1}{S}\sum_{s=1}^{S}(u_{s}^{\prime})_{ij}$ fulfills $|(\bar{u})_{ij}-(u_{s}^{\prime})_{ij}|<\varepsilon.$ Further, the function value of $\hat{u}$ which is piecewise constant w.r.t. $\mathcal{P}$ is obtained by $\hat{u}|_{\mathcal{P}_{i}}=\sum_{x\in\mathcal{P}_{i}}\bar{u}(x)/\#\mathcal{P}_{i}.$ Hence, we may estimate

[TABLE]

where $L$ is the maximal length of a path connecting any two pixels as given by Definition 2. As a worst case estimate, we get $L\leq C_{1}$ where we define $C_{1}$ as one fourth of the number of elements in $\Omega,$ i.e., $C_{1}=\tfrac{\#\Omega}{4}.$ This shows (74).

For $F_{\rho}$ given by (62), we have

[TABLE]

with

[TABLE]

as given in (76). The first inequality holds since as a local minimizer of the quadratic penalty relaxation (5), $u^{\prime}$ is the global minimizer of $F_{\rho}$ on $\mathcal{A}^{\mathcal{P}}$ by Lemma 11 and since $(\hat{u},\ldots,\hat{u})\in\mathcal{A}^{\mathcal{P}}$ by construction. The next inequalities apply the triangle inequality and estimates on matrix norms. The last inequality is a consequence of the fact that $\sum\nolimits_{s=1}^{S}\tfrac{1}{S}\left\|Au_{s}^{\prime}-f\right\|_{2}\leq\|f\|_{2}.$ Otherwise, if $\|Au_{s}^{\prime}-f\|_{2}>\|f\|_{2},$ choosing $u_{s}^{\prime}=0$ would yield a lower function value which would contradict the minimality of $u^{\prime}.$

Now consider the partitioning $\mathcal{P}^{\prime}$ induced by $\hat{u},$ and the corresponding minimizer $u^{\ast},$ i.e.,

[TABLE]

where, for $(u,\ldots,u)\in\mathcal{B}^{\mathcal{P}^{\prime}},$ we have $F_{\rho}(u,\ldots,u)=\left\|Au-f\right\|_{2}^{2}.$ By Lemma 12, $u^{\ast}$ is a local minimizer of the Potts problem (2). On the other hand, by orthogonality in an inner product space, we have

[TABLE]

where $P_{A\left(\mathcal{B}^{\mathcal{P}^{\prime}}\right)}$ denotes the orthogonal projection onto the image of $\mathcal{B}^{\mathcal{P}^{\prime}}$ under the linear mapping $A.$ Thus,

[TABLE]

Inserting $u^{\ast}$ in the estimate (3.3), we get

[TABLE]

This allows us to further estimate

[TABLE]

If now the operator $A$ is lower bounded, then

[TABLE]

which completes the proof. ∎

3.4 Majorization-minimization for multivariate Potts problems

In this part we build the basis for the convergence analysis of Algorithm 1 and Algorithm 2.

We first recall some basics on surrogate functionals. We consider functionals $F(u)$ of the form $F(u)=\|Xu-z\|^{2}+\gamma J(u),$ where $X$ is a given (measurement) matrix with operator norm $\|X\|<1$ (with the operator norm formed w.r.t. the $\ell^{2}$ norm), $z$ is a given vector (of data), $J$ is an arbitrary (not necessarily convex) lower semicontinuous functional, and $\gamma>0$ is a parameter. In general, the surrogate functional $F^{\rm surr}(u,v)$ of $F(u)$ is given by

[TABLE]

Lemma 16**.**

Consider the functionals $F(u)=\|Xu-z\|^{2}+\gamma J(u)$ as above with $\|X\|<1.$ (For our purposes, $J$ is the regularizer $\|D(u)\|_{0,\omega}$ given by (10).) Then, we get for the associated surrogate functional $F^{\rm surr}$ given by (87) (with $J$ as regularizer), that

1. i.

the inequality

[TABLE]

holds for all $v;$ and $F^{\rm surr}(u,v)=F(u)$ holds if and only if $u=v;$ 2. ii.

the functional values $F(u^{k})$ of the sequence $u^{k}$ given by the surrogate iteration $u^{k+1}=\operatorname*{argmin}_{u}F^{\rm surr}(u,u^{k})$ are non-increasing, i.e.,

[TABLE] 3. iii.

the distance between consecutive members of the previous surrogate sequence $u^{k}$ converges to $0,$ i.e.,

[TABLE]

We note that –when minimizing $F$– the condition $\|X\|<1$ can always be achieved by rescaling, i.e., dividing the functional $F$ by a number which is larger than $\|X\|^{2}.$ Proofs of the general statements above on surrogate functionals (which do not rely on the specific structure of the problems considered here) may for instance be found in the above mentioned papers [28, 9, 33].

We now employ properties of the quadratic penalty relaxation $P_{\gamma,\rho}(u_{1},\ldots,u_{S})$ of the Potts energy given by (5). The strategy is similar to the authors’ approach for the univariate case in [90]. We first show that the minimizers of $P_{\gamma,\rho}(u_{1},\ldots,u_{S})$ (with $B=\mathrm{id}$ in (11)) which are precisely the solutions of (16) have a minimal directional jump height which only depends on the scale parameter $\gamma,$ the directional weights $\omega_{s}$ and the constant $L_{\rho}$ but not on the particular input data. Here, for the multivariate discrete function $u=(u_{1},\ldots,u_{S})$ (and the directional system $a_{s},$ $s=1,\ldots,S$) a directional jump is a jump in the $s$th component $u_{s}$ in direction $a_{s}$ for some $s.$ In particular, jumps of $u_{s}$ in directions $a_{s^{\prime}}$ with $s^{\prime}\neq s$ are not considered.

Lemma 17**.**

We consider the function $P_{\gamma,\rho}(u_{1},\ldots,u_{S})$ of (11) for the choice $B=\mathrm{id}$ and data $h=(h_{1},\ldots,h_{S})$. In other words, we consider the problem (16) for arbitrary data $h=(h_{1},\ldots,h_{S})$. Then there is a constant $c>0$ which is independent of the minimizer $u^{\ast}=(u_{1}^{\ast},\ldots,u_{S}^{\ast})$ of (16) and the data $h$ such that the minimal directional jump height $j_{\min}(u^{\ast})$ (w.r.t. the directional system $a_{s},$ $s=1,\ldots,S,$) of a minimzer $u^{\ast}$ fulfills

[TABLE]

The constant $c$ depends on $\gamma,$ the directional weights $\omega_{s}$ and the constant $L_{\rho}.$

Proof.

Writing $u=(u_{1},\ldots,u_{S})$ we restate (16) as the problem of minimizing

[TABLE]

where we use the notation $\|D(u_{1},\ldots,u_{S})\|_{0,\omega}=\sum_{s=1}^{S}\omega_{s}\left\|\nabla_{a_{s}}u_{s}\,\right\|_{0}$ introduced in (10). We let

[TABLE]

where $W$ denotes the maximal length of the signal $u$ per dimension (e.g., if $u$ denotes an $l\times b$ image, then $W=\max(l,b)$.) We now assume that $h_{\min}(u^{\ast})<c,$ which means that the minimizer $u^{\ast}$ has a directional jump of height smaller than $c.$ For such $u^{\ast},$ we construct an element $u^{\prime}$ with a smaller $P^{\mathrm{id}}_{\gamma/L_{\rho}^{2}}$ value which yields a contradiction since $u^{\ast}$ is a minimizer of $P^{\mathrm{id}}_{\gamma/L_{\rho}^{2}}$. To this end, we let $a_{s}$ be a direction such that the component $u^{\ast}_{s}$ of $u^{\ast}$ has a jump of height smaller than $c.$ We denote the (discrete) directional intervals in direction $a_{s}$ near the directional jump by $I_{1},I_{2}$ and the corresponding points near the directional jump of $u_{s}^{\ast}$ by $x_{1}$ and $x_{2}.$ We let $m_{1},m_{2}$ and $m$ be the mean of $h_{s}$ on $I_{1},I_{2}$ and $I_{1}\cup I_{2}$, respectively. We define

[TABLE]

By construction, $\|\nabla u^{\prime}\|_{0}=\|\nabla u^{\ast}\|_{0}-1.$ and thus

[TABLE]

Since $u^{\ast}$ is a minimizer of $P^{\mathrm{id}}_{\gamma/L_{\rho}^{2}},$ its $s$th component $u^{\ast}_{s}$ equals $m_{1}$ on $I_{1}$ and $m_{2}$ on $I_{2}.$ Further, as $u_{s^{\prime}}^{\prime}=u_{s^{\prime}}^{\ast}$ if $s^{\prime}\neq s$ and $u_{s}^{\ast}$ and $u_{s}^{\prime}$ only differ on $I_{1}\cup I_{2},$ we have that

[TABLE]

where $l_{1},l_{2}$ denote the length of $I_{1}$ and $I_{2},$ respectively. Employing (94) together with (3.4) we get

[TABLE]

The validity of the last inequality follows by (92). Together, $u^{\prime}$ has a smaller function value than $u^{\ast}$ which is a contradiction to $u^{\ast}$ being a minimizer which shows the assertion. ∎

Proposition 18**.**

The iteration (19) of Algorithm 1 converges to a local minimizer of the quadratic penalty relaxation $P_{\gamma,\rho}$ of the Potts objective function given by (5). The convergence rate is linear.

Proof.

We divide the proof into three parts. First, we show that the directional partitionings induced by the iterates $u^{(n)}$ become fixed after sufficiently many iterations. In a second part, we derive the convergence of Algorithm 1 and, in a third part, we show that the limit point is a local minimizer of $P_{\gamma,\rho}$.

(1) We first show that the directional partitioning $\mathcal{I}^{n}$ induced by the iterates $u^{(n)}$ gets fixed for large $n$. For every $n\in\mathbb{N},$ the iterate $u^{(n)}$ of Algorithm 1 is a minimizer of the function $P_{\gamma,\rho}$ of (11) for the choice $B=\mathrm{id}$ as it appears in (16). Here, the data $h=(h_{1},\ldots,h_{S})$ is given by (17). By Lemma 17 there is a constant $c>0$ which is independent of the particular $u^{(n)}=(u_{1}^{(n)},\ldots,u_{S}^{(n)})$ of (16) and the data $h$ such that the minimal directional jump height $j_{\min}(u^{(n)})$ fulfills

[TABLE]

We note that the parameter $\gamma,$ the directional weights $\omega_{s}$ and the constant $L_{\rho}$ which the constant $c$ depends on by Lemma 17 do not change during the iteration of Algorithm 1.

If two iterates $u^{(n)},u^{(n+1)}$ have different induced directional partitionings $\mathcal{I}^{n},\mathcal{I}^{n+1}$ their $\ell^{\infty}$ distance always fulfills $\|u^{(n)}-u^{(n+1)}\|_{\infty}>c/2$ since both $u^{(n)},u^{(n+1)}$ have minimal jump height of at least $c$ and different induced directional partitionings. This implies $\|u^{(n)}-u^{(n+1)}\|_{2}>c/2$ for the $\ell^{2}$ distance as well. This may only happen for small $n,$ since Lemma 16 by (89), we have $\|u^{(n)}-u^{(n+1)}\|_{2}\to 0$ as $n$ increases. Hence, there is an index $N$ such that, for all $n\geq N,$ the directional partitionings $\mathcal{I}^{n}$ are identical.

(2) We use the previous observation to show the convergence of Algorithm 1. We consider iterates $u^{(n)}$ with $n\geq N;$ they have the same induced directional partitionings which we denote by $\mathcal{I}^{\prime},$ and all jumps have minimal jump height $c$. Hence, for $n\geq N,$ the iteration of (16) can be written as

[TABLE]

with $P_{\mathcal{I}^{\prime}}$ being the orthogonal projection onto the $\ell^{2}$ space $\mathcal{A}^{\mathcal{P}}$ consisting of functions which are piecewise constant w.r.t. the directional partitioning $\mathcal{I}^{\prime},$ and where $h^{(n)}$ depends on $u^{(n)}$ via

[TABLE]

as given by (17). As introduced before, we use the symbols $\tilde{A}$ to denote the block diagonal matrix with constant entry $A$ on each diagonal component, and $\tilde{f}$ for the block vector of corresponding dimensions with entry $f$ in each component. With this notation we may write (97) as

[TABLE]

Since $u^{(n)}$ is piecewise constant w.r.t. the directional partitioning $\mathcal{I}^{\prime},$ we have $u^{(n)}=P_{\mathcal{I}^{\prime}}u^{(n)}.$ Using this fact and the fact that $P_{\mathcal{I}^{\prime}}$ is an orthogonal projection we obtain

[TABLE]

Since $C\tilde{A}^{T}\tilde{f}=0,$ the iteration (100) can be interpreted as Landweber iteration for the block matrix consisting of the upper block $(\tilde{A}P_{\mathcal{I}^{\prime}})/(\sqrt{S}L_{\rho})$ and the lower block $(\sqrt{\rho}CP_{\mathcal{I}^{\prime}})/(\sqrt{S}L_{\rho})$ and data $\tilde{f}/(\sqrt{S}L_{\rho})$ extended by $0.$ The Landweber iteration converges at a linear rate; cf., e.g., [31]. Thus, the iteration (97) convergences and, in turn, we get the convergence of Algorithm 1 at a linear rate to some limit $u^{\ast}$.

(3) We show that $u^{\ast}$ is a local minimizer. Since $u^{\ast}$ is the limit of the iterates $u^{(n)}$, the jumps of $u^{\ast}$ also have minimal height $c,$ the number of jumps are equal to those of the $u^{(n)}$ for all $n\geq N,$ and the induced directional partitioning $\mathcal{I}^{\ast}$ equals the partitioning $\mathcal{I}^{\prime}$ of the $u^{(n)}$ for $n\geq N.$ Since $u^{\ast}$ equals the limit of the above Landweber iteration, $u^{\ast}$ minimizes $F_{\rho}$ given by (62) on $\mathcal{A}^{\mathcal{P}_{\mathcal{I}^{\prime}}}.$ Then by Lemma 11 $u^{\ast}$ is a local minimizer of the relaxed Potts energy $P_{\gamma,\rho}$ which completes the proof. ∎

After having shown the convergence of Algorithm 1 to a local minimizer, we have now gathered all information to show Theorem 5.

Proof of Theorem 5.

Assertion i. was stated and shown as Proposition 13 in Section 3.3. By Proposition 18 Algorithm 1 produces a local minimizer. Then the assertion ii. is a consequence of Proposition 15. ∎

3.5 Estimating the distance between the objectives

The next lemma is a preparation for the proof of item (iii) of Theorem 4.

Lemma 19**.**

We consider Algorithm 1 for the quadratic penalty relaxation (5) of the multivariate Potts problem. For any output $u=(u_{1},\ldots,u_{S})$ of Algorithm 1 we have that

[TABLE]

Here, $\sigma_{1}$ denotes the smallest non-zero eigenvalue of $C^{T}C$ with $C$ given by (49).

Proof.

Since $u=(u_{1},\ldots,u_{S})$ is the output of Algorithm 1 it is a local minimizer of the relaxed Potts problem (5). In particular, there is a directional partitioning $\mathcal{I}$ with respect to which $u$ is piecewise constant. We denote the induced partitioning by $\mathcal{P}=\mathcal{P}_{\mathcal{I}}.$ By Lemma 14 we have

[TABLE]

where $\mu^{\ast}$ is an arbitrary Lagrange multiplier of (48). By Lemma 10 we have that $\|\mu^{\ast}\|\leq 2\sigma_{1}^{-1/2}S^{-1/2}\|A\|\|f\|,$ for any partitioning of the discrete domain $\Omega,$ and in particular for the partitioning $\mathcal{P}=\mathcal{P}_{\mathcal{I}}.$ This shows that

[TABLE]

Since $u$ is a local minimizer of the relaxed Potts problem (5), it is a minimizer of $F_{\rho}$ on $\mathcal{A}^{\mathcal{P}}$ by Lemma 11, and the second summand on the right hand side equals zero. This shows (101) and completes the proof. ∎

We have now gathered all information necessary to show Theorem 4.

Proof of Theorem 4.

Part (i) is shown by Proposition 18. Concerning (ii) we first show that any global minimizer of the relaxed Potts energy $P_{\gamma,\rho}$ given by (5) appears as a stationary point of Algorithm 1. To this end we start Algorithm 1 with a global minimizer $\bar{u}^{\ast}=(u_{1}^{\ast},\ldots,u_{S}^{\ast})$ as initialization. Then, we have for all $\bar{v}=(v_{1},\ldots,v_{S})$ with $\bar{v}\neq\bar{u}^{\ast},$

[TABLE]

Here, $B$ is given by (8). The estimate (103) means that $\bar{u}^{\ast}$ is the minimizer of the surrogate functional w.r.t. the first component, i.e., it is the minimizer of the mapping $\bar{v}\mapsto$ $P_{\gamma,\rho}^{\rm surr}(\bar{v},\bar{u}^{\ast}).$ Thus, the iterate $\bar{u}^{(1)}=(u^{(1)}_{1},\ldots,u^{(1)}_{S})$ of Algorithm 1 equals $\bar{u}^{\ast}$ when the iteration is started with $\bar{u}^{\ast}.$ Thus, the global minimizer $\bar{u}^{\ast}$ is a stationary point of Algorithm 1. It remains to show that each stationary point of Algorithm 1 is a local minimizer of the relaxed Potts energy $P_{\gamma,\rho}$. This has essentially already been done in the proof of Proposition 18: start the iteration given by (16) with a stationary point $u^{\prime};$ its limit equals $u^{\prime}$ and is thus a local minimizer by Proposition 18.

Concerning (iii), we use Lemma 19 to estimate

[TABLE]

The second inequality follows by our choice of $\rho$ in (34) as $\rho>2\varepsilon^{-1}\ \sigma_{1}^{-1/2}S^{-1/2}\|A\|\|f\|.$ This shows the validity of (iii) and completes the proof. ∎

3.6 Convergence Analysis of Algorithm 2

We start out showing that Algorithm 2 is well-defined in the sense that the inner iteration governed by (25) terminates. This result was formulated as Theorem 6.

Proof of Theorem 6.

We have to show that, for any $k\in\mathbb{N},$ there is $n\in\mathbb{N}$ such that

[TABLE]

To see the right-hand side of (105), we notice that, by Proposition 18, the iteration (19) converges to a local minimizer of the quadratic penalty relaxation $P_{\gamma,\rho}(u_{1},\ldots,u_{S})$ of the Potts energy. The inner loop of Algorithm 2 precisely computes the iteration (19) (for the penalty parameter $\rho$ which increases with $k.$) Thus, the distance between consecutive iterates $u^{(k,n)}_{s},u^{(k,n-1)}_{s}$ converges to zero as $n$ increases which implies the validity of the right-hand side of (105) for sufficiently large $n,$ and all $k\in\mathbb{N}.$

To see the left-hand inequality in (105), we notice that, by the considerations above, the inner loop of Algorithm 2 would converge to a minimizer $\bar{u}^{(k),\ast}=(u_{1}^{(k),\ast},\ldots,u_{S}^{(k),\ast})$ if it was not terminated by (105) for all $k\in\mathbb{N}.$ Since $\bar{u}^{(k),\ast}$ is a local minimizer of the relaxed Potts problem (5) for the parameter $\rho_{k}$, it is a minimizer of $F_{\rho_{k}}$ on $\mathcal{A}^{\mathcal{P}}$ (where $\mathcal{P}$ denotes the partitioning induced by $\bar{u}^{(k),\ast}$) by Lemma 11. Hence, for any $k\in\mathbb{N}$ and any $\xi>0$ there is $\bar{u}^{(k,n)}=(u^{(k,n)}_{1},\ldots,u^{(k,n)}_{S})$ such that $F_{\rho_{k}}(\bar{u}^{(k,n)})-F_{\rho_{k}}(\bar{u}^{(k),\ast})<\xi.$ We let $\tau=(t-2\sigma_{1}^{-1/2}S^{-1/2}\|A\|\ \|f\|)/\rho_{k},$ and choose $\xi=\rho_{k}\tau^{2}.$ Using this together with Lemma 14 we estimate

[TABLE]

where $\mu^{\ast}$ is an arbitrary Lagrange multiplier of (48). Here, the last inequality is true since by Lemma 10 we have that $\|\mu^{\ast}\|\leq 2\sigma_{1}^{-1/2}S^{-1/2}\|A\|\|f\|$ which implies that $\tau\leq(t-\mu^{\ast})/\rho_{k}.$ The estimate (3.6) shows the left-hand inequality in (105) and completes the proof. ∎

We have now gathered all information to prove Theorem 7 which deals with the convergence properties of Algorithm 2.

Proof of Theorem 7.

We start out to show that any accumulation point of the sequence $u^{(k)}$ produced by Algorithm 2 is a local minimizer of the Potts problem (3). Let $u^{\ast}$ be such an accumulation point and let $\mathcal{I}^{\ast}$ be the directional partitioning induced by $u^{\ast}.$ We may extract a subsequence $u^{(k_{l})}$ of the sequence $u^{(k)}$ such that $u^{(k_{l})}$ converges to $u^{\ast}$ as $l\to\infty,$ and such that the directional partitionings $\mathcal{I}^{k_{l}}$ induced by the $u^{(k_{l})}$ all equal the directional partitioning $\mathcal{I}^{\ast},$ i.e., $\mathcal{I}^{k_{l}}=\mathcal{I}^{\ast}$ for all $l\in\mathbb{N}.$ We let

[TABLE]

with the matrix $C$ given by (49), and estimate

[TABLE]

We recall that $\tilde{A}$ was the block diagonal matrix having the matrix $A$ as entry in each diagonal component and that $F_{\rho_{k_{l}}}$ was given by (62). We notice that the second before last inequality follows by the right hand side of (105). We further estimate

[TABLE]

which is a consequence of the left hand side of (105). Hence, the sequence $\mu^{k_{l}}$ is bounded and thus has a cluster point, say $\mu^{\ast},$ by the Bolzano-Weierstraß Theorem. By passing to a further subsequence (where we suppress the new indexation for better readability and still use the symbol $l$ for the index) we get that

[TABLE]

Now, on this subsequence, we have that $u^{(k_{l})}\to u^{\ast}$ and that $\mu^{k_{l}}\to\mu^{\ast}.$ Hence taking limits on both sides of (3.6) yields

[TABLE]

since $\delta_{k_{l}}\to 0$ as $l\to\infty.$ Further,

[TABLE]

This implies that the components of $u^{\ast}$ are equal, i.e., $u_{s}^{\ast}=u_{s^{\prime}}^{\ast}$ for all $s,s^{\prime}.$ In particular $u^{\ast}$ is a feasible point for the Potts problem (3). Or, letting $\mathcal{P}^{\ast}$ be the partitioning induced by $u^{\ast},$ we have that $u^{\ast}\in\mathcal{B}^{\mathcal{P}}.$ Then, (110) tells us that $u^{\ast}$ minimizes (48) which by Lemma 12 tells us that $u^{\ast}$ is a local minimizer of (3), or synonymously, that any component of $u^{\ast}$ (which are all equal) minimizes the Potts problem (2). This shows the first assertion of Theorem 7.

We continue by showing the second assertion of Theorem 7, i.e., if $A$ is lower bounded, then the sequence $u^{(k)}$ produced by Algorithm 2 has a cluster point. Then, by the above considerations, each cluster point is a local minimizer which shows the assertion. To this end, we show that, if $A$ is lower bounded, the sequence $u^{(k)}$ produced by Algorithm 2 is bounded which by the Heine-Borel property of finite dimensional Euclidean space implies that it has a cluster point. So we assume that $A$ is lower bounded, and consider the sequence $u^{(k)}=(u_{1}^{(k)},\ldots,u_{S}^{(k)})$ produced by Algorithm 2. As in the proof of Theorem 6 we see that, for any $k\in\mathbb{N}$, there is a local minimizer $u^{(k),\ast}=(u_{1}^{(k),\ast},\ldots,u_{S}^{(k),\ast})$ of (5) such that

[TABLE]

where $C_{2}$ is a constant independent of $k.$ By Lemma 11, $u^{(k),\ast}$ is a minimizer of $F_{\rho}$ on $\mathcal{A}^{\mathcal{P}}$ (where $\mathcal{P}$ denotes the partitioning induced by $u^{(k),\ast}$.) Hence,

[TABLE]

by choosing the function having the zero function as entry in each component as a candidate. This implies

[TABLE]

Then, since $A$ is lower bounded, there is a constant $c>0$ such that

[TABLE]

where we used (113) for the last inequality. Combining this estimate with (112) yields

[TABLE]

Since we have chosen $\delta_{k}$ as a sequence converging to zero, (115) shows that the sequence $u^{(k)}$ is bounded which implies that it has cluster points. This completes the proof. ∎

4 Numerical Results

In this section, we show the applicability of our methods to different imaging tasks. We start out by providing the necessary implementation details. Then we compare the results of the quadratic relaxation (5) (Algorithm 1) to the ones of the Potts problem (2) (Algorithm 2). Next, we apply Algorithm 2 to blurred image data and to image reconstruction from incomplete Radon data. Finally, we consider the image partitioning problem according to the classical Potts model.

Implementation details.

We implemented Algorithm 1 and Algorithm 2 for the coupling schemes in (7) and the set of compass and diagonal directions $(1,0),(0,1),(1,1),(1,-1)$ with weights $\omega_{1,2}=\sqrt{2}-1$ and $\omega_{3,4}=1-\frac{\sqrt{2}}{2}$.

Concerning Algorithm 1 we observed both visually and quantitatively appealing results if we use relaxed step-sizes $L_{\rho}^{\lambda}=L_{\rho}[\lambda+(1-(n+1)^{-1/2})(1-\lambda)]$ for an empirically chosen parameter $0<\lambda\leq 1$, where $L_{\rho}$ denotes the estimate in Lemma 8. The iterations were stopped when the nearness condition (18) was fulfilled and the iterates did not change anymore, i.e., when $\|u_{1}^{(n)}-u_{1}^{(n-1)}\|/(\|u_{1}^{(n)}\|+\|u_{1}^{(n-1)}\|)$ and $\|u_{2}^{(n)}-u_{2}^{(n-1)}\|/(\|u_{2}^{(n)}\|+\|u_{2}^{(n-1)}\|)$ were smaller than $10^{-6}$. The result of Algorithm 1 was transformed into a feasible solution of (3) by applying the projection procedure described in Section 2.2 (Procedure 1). As initialization we applied 1000 Landweber iterations with step-size $1/\|A\|^{2}$ to the least squares problem induced by the linear operator $A$ and data $f$.

Concerning Algorithm 2, we set $\rho^{(0)}=10^{-3}$ in all experiments which we incremented by the factor $\tau=1.05$ in each outer iteration. The $\delta$-sequence was chosen as $\delta^{(k)}=\frac{1}{\eta\rho^{(k)}}$ for $\eta=0.95$ when coupling all variables and $\eta=0.98$ when coupling consecutive variables. Similarly to Algorithm 1, step A of Algorithm 2 was performed using the relaxed step sizes $L_{\rho}^{\lambda}=L_{\rho}[\lambda+(1-(n+1)^{-1/2})(1-\lambda)]$ for an application-dependent parameter $0<\lambda\leq 1$ and for the estimate $L_{\rho}$ in Lemma 8. We stopped the iterations when the relative discrepancy of the first two splitting variables $\|u^{(k)}_{1}-u^{(k)}_{2}\|/(\|u^{(k)}_{1}\|+\|u^{(k)}_{2}\|)$ was smaller than $10^{-6}$. We initialized Algorithm 2 with $A^{T}f$.

Comparison of Algorithm 1 and Algorithm 2.

We compare Algorithm 1 and Algorithm 2 for blurred image data, that is, the linear operator $A$ in (1) amounts to the convolution by a kernel $K$. In the present experiment, we chose a Gaussian kernel with standard deviation $\sigma=3$ and of size $6\sigma+1$. Here, we coupled all splitting variables and chose the step-size parameter $\lambda=0.4$ for Algorithm 1 and $\lambda=0.35$ for Algorithm 2, respectively. In Figure 1 we applied both methods to a blurred natural image. While both algorithms yield reasonable partitionings, Algorithm 2 provides smoother edges than Algorithm 1. Further, Algorithm 1 produces some smaller segments (at the treetops).

Application to blurred data.

For the following experiments, we focus on Algorithm 2. In case of motion blur we set the step-size parameter to $\lambda=0.25$, while for Gaussian blur we set $\lambda=0.35$ as in Figure 1. We compare our method with the Ambrosio-Tortorelli approximation [2] of the classical Mumford-Shah model (which itself tends to the piecewise constant Mumford-Shah model for increasing variation penalty) given by

[TABLE]

The variable $v$ serves as an edge indicator and $\varepsilon>0$ is an edge smoothing parameter that is chosen empirically. The parameter $\gamma>0$ controls the weight of the edge length penalty and the parameter $\alpha>0$ penalizes the variation. In this respect, a higher value of $\alpha$ promotes solutions which are closer to being piecewise constant. In the limit $\alpha\to\infty$ minimizers of (116) are piecewise constant. Our implementation follows the scheme presented in [6]. The functional $A_{\varepsilon}$ is alternately minimized w.r.t. $u$ and $v$. To this end, we iteratively solve the Euler-Lagrange equations

[TABLE]

where $\widetilde{K}(x)=K(-x)$. The first line is solved w.r.t. $v$ using a MINRES solver and the second line is solved using the method of conjugate gradients [6]. The iterations were stopped when the relative change of both variables was small, i.e., if both $\|u^{k+1}-u^{k}\|/(\|u^{k}\|+10^{-6})<10^{-3}$ and $\|v^{k+1}-v^{k}\|/(\|v^{k}\|+10^{-6})<10^{-3}$.

Figure 4 shows the restoration of a traffic sign from simulated horizontal motion blur. For the Ambrosio-Tortorelli approximation we set $\alpha=10^{5}$ to promote a piecewise constant solution. We observe that both the Ambrosio-Tortorelli approximation and the proposed method restore the data to a human readable form. However, the Ambrosio-Tortorelli result shows clutter and blur artifacts. Our method provides sharp edges and it produces less artifacts.

In Figure 4 we partition a natural image blurred by a Gaussian kernel and corrupted by Gaussian noise. We observed that the Ambrosio-Tortorelli result was heavily corrupted by artifacts for $\alpha=10^{5}$. This might be attributed to the underlying linear systems in scheme (117) which become severely ill-conditioned for large choices of the variation penalty $\alpha$. Therefore, we chose the moderate variation penalty $\alpha=5$ which does only provide an approximately piecewise constant (rather piecewise smooth) result. The result does not fully separate the background from the fish in terms of edges. On the other hand, the result of the proposed method sharply differentiates between background and fish. Further, it highlights various segments of the fish.

Reconstruction from Radon data.

We here consider reconstruction from Radon data which appears for instance in computed tomography. We recall that the Radon transform reads

[TABLE]

where $s\in\mathbb{R},$ $\theta\in S^{1}$ and $\theta^{\perp}\in S^{1}$ is (counterclockwise) perpendicular to $\theta$; see [64]. For our experiments, we use a discretization of the Radon transform created using the AIR tools software package [39]. Regarding our method, we employed coupling of consecutive splitting variables and the step-size parameter was set to $\lambda=0.11$. To quantify the reconstruction quality we use the mean structural similarity index (MSSIM) [89] which is bounded from above by 1, where higher values indicate better results.

We compare the proposed method to filtered back projection (FBP) which is standard in practice [71]. The FBP is computed using its Matlab implementation with the standard Ram-Lak filter. Furthermore, we compare with total variation (TV) regularization [76] in the Lagrange form $\|Ru-f\|_{2}^{2}+\mu\|\nabla u\|_{1}$ with parameter $\mu>0$. Its implementation follows the Chambolle-Pock algorithm [21]. The corresponding parameter $\mu$ was tuned w.r.t. the MSSIM index.

In Figure 7 we show the reconstruction results for the Shepp-Logan phantom from undersampled (25 angles) and noisy Radon data. Standard FBP produces strong streak artifacts which are typical for angular undersampling, and the reconstruction suffers from noise. The TV regularization and the proposed method both provide considerably improved reconstruction results. The proposed method achieves a higher MSSIM value than the TV reconstruction, and it provides a reconstruction which is less grainy than the TV result.

Image partitioning.

Finally, we consider the classical Potts problem corresponding to $A=\mathrm{id}$ in (1). While the focus of the present paper is on a general imaging operator $A$, we next observe that it also works rather well for $A=\mathrm{id}$. We used the full coupling scheme and set the step-size parameter to $\lambda=0.55$.

To put our result in context we added the results of two other methods for $A=\mathrm{id}$: the $L_{0}$ gradient smoothing method of Xu et. al [94] and the state-of-the-art $\alpha$-expansion graph cut algorithm based on max-flow/min-cut of the library GCOptimization 3.0 of Veksler and Delong [13, 12, 52]. The method of [94] has a parameter $\kappa>1$ to control the convergence speed and a smoothing weight $\nu$. In our experiments, we set $\kappa=1.01$ and $\nu=0.1$. For the graph cuts the same neighborhood weights and jump penalty were used as for the proposed method. The discrete labels are computed via k-means.

In Figure 8, we show the results for a natural image corrupted by Gaussian noise. The Ambrosio-Tortorelli result suffers from clutter and remains noisy. The result of $L_{0}$ gradient smoothing over-segments the textured window area while it smooths out details of the cross. The state of the art graph cuts method and the proposed method both provide satisfying results which are visually comparable. Further, they yield solutions with comparable Potts energy values. For instance, on the IVC dataset [55] which consists of $10$ natural color images of size $512\times 512,$ for the model parameters, $\gamma=0.25$ and $\gamma=1,$ the mean values of the proposed approach are $7107.8$ and $13053.2$ compared with the respective mean energies of the graph cut approach $7093.2$ and $13008.7$ which differ by about half a percent. (For the results in Figure 8, the energy value of the proposed approach is 25067.7  compared with 25119.5  for the graph cuts approach.) Here, for the graph cut approach, we took the mean value of the input image on each computed segment before computing the Potts objective function. To sum up, while the proposed method can handle general linear operators $A$, the quality of the results for $A=\rm id$ is comparable with the state-of-the-art graph cut algorithm for $A=\rm id$.

5 Conclusion

In this paper, we have proposed a new iterative minimization strategy for multivariate piecewise constant Mumford-Shah/Potts energies as well as their quadratic penalty relaxations. Our schemes are based on majorization-minimization or forward-backward splitting methods of Douglas-Rachford type [57]. In contrast to the approaches in [33, 9, 60, 61] for sparsity problems which lead to thresholding algorithms, our approach leads to non-separable yet computationally tractable problems in the backward step.

As a second part, we have provided a convergence analysis for the proposed algorithms. For the proposed quadratic penalty relaxation scheme, we have shown convergence towards a local minimizer. Due to the NP hardness of the quadratic penalty relaxation, the convergence result is in the range of what can be expected best. Concerning the scheme for the non-relaxed Potts problem we have also performed a convergence analysis. In particular, we have obtained results on the convergence towards local minimizers on subsequences. The quality of these results is comparable with the results of [60, 61] where, compared with these papers, we had to deal with the non-separability of the backward step as an additional challenge.

Finally, we have shown the applicability of our schemes in several experiments. We have applied our algorithms to deconvolution problems including the problem of deblurring and denoising motion blur images. We have further dealt with noisy and undersampled Radon data for the task of joint reconstruction, denoising and segmentation. Finally, we have applied our approach to the situation of pure image partitioning (without blur) which is a widely considered problem in computer vision.

Acknowledgements

L. Kiefer and A. Weinmann were supported by the German Research Foundation (DFG) Grant WE5886/4-1. Additionally, A. Weinmann was supported by DFG Grant WE5886/3-1. M. Storath was supported by DFG Grant STO1126/2-1.