Multiscale Higher Order TV Operators for L1 Regularization

TL;DR

This paper introduces a multiscale higher order total variation (MHOTV) method for signal and image denoising that improves upon traditional $ ext{L}_1$ regularization by reducing artifacts and better promoting sparsity, leveraging connections to wavelets.

Contribution

The authors develop a novel MHOTV approach related to multiscale wavelets, with efficient computational tools, demonstrating improvements over existing wavelet and HOTV methods.

Findings

MHOTV reduces artifacts compared to traditional $ ext{L}_1$ methods.

MHOTV shows better sparsity promotion and denoising performance.

Numerical results confirm the advantages of MHOTV over wavelets and HOTV.

Abstract

In the realm of signal and image denoising and reconstruction, $\ell_1$ regularization techniques have generated a great deal of attention with a multitude of variants. A key component for their success is that under certain assumptions, the solution of minimum $\ell_1$ norm is a good approximation to the solution of minimum $\ell_0$ norm. In this work, we demonstrate that this approximation can result in artifacts that are inconsistent with desired sparsity promoting $\ell_0$ properties, resulting in subpar results in {some} instances. With this as our motivation, we develop a multiscale higher order total variation (MHOTV) approach, which we show is related to the use of multiscale Daubechies wavelets. We also develop the tools necessary for MHOTV computations to be performed efficiently, via operator decomposition and alternatively converting the problem into Fourier space. The relationship of higher order regularization methods with wavelets, which we believe has generally gone unrecognized, is shown to hold in several numerical results, although notable improvements are seen with our approach over both wavelets and classical HOTV.

Tables

Table 1. Table 1: Average relative reconstruction error over 100 simulations, as a function of the order of the method and number of levels in the multiscale approaches. The minimums for each SNR are emphasized in bold.

SNR		Order 1	Order 2	Order 3	Order 1.5	Order 2.5
		mhotv  Daub	mhotv  Daub	mhotv  Daub	mhotv	mhotv
	1 level	.1624, .1624	.2039, .1961	.2464, .2306	.1819	.2328
2	2 levels	.1612, .1617	.1742, .1852	.2135, .2223	.1782	.2183
	3 levels	.1699, .1615	.1513, .1778	.1745, .2149	.1776	.1975
	4 levels	.2001, .1647	.1584, .1745	.1764, .2104	.2031	.2102
	1 level	.0864, .0864	.0971, .0914	.1293, .1090	.1025	.1287
5	2 levels	.0858, .0857	.0761, .0838	.0946, .1004	.0987	.1172
	3 levels	.0926, .0864	.0668, .0805	.0766, .0982	.1016	.1133
	4 levels	.1100, .0894	.0742, .0801	.0828, .0981	.1186	.1276
	1 level	.0543, .0542	.0509, .0480	.0694, .0572	.0690	.0841
10	2 levels	.0542, .0539	.0400, .0442	.0489, .0528	.0657	.0763
	3 levels	.0589, .0547	.0359, .0430	.0399, .0522	.0694	.0776
	4 levels	.0696, .0570	.0413, .0436	.0442, .0535	.0802	.0880

Equations

$$f_{rec}=\arg\min_{f}\Big{\{}\|Af-b\|_{2}^{2}+\lambda\|Tf\|_{1}\Big{\}},$$

$$({\varphi }_k{)}_m=\{\begin{array}{ll}(-1{)}^k& \text{mbox}ifm=0\\ 0& \text{mbox}if1\le m<N-k\\ (-1{)}^{k+m+N}(\genfrac{}{}{0ex}{}{k}{N-m})& \text{mbox}ifN-k\le m<N\end{array}.$$

$$f\ast {\varphi }_k,$$

$$\Phi_{2}=\left(\begin{array}[]{c c c c c c}1&-2&1&0&\dots&0\\ 0&1&-2&1&\dots&0\\ 0&0&1&-2&\dots&0\\ \vdots&&\ddots&&\dots&\vdots\\ 1&0&\dots&&1&-2\\ -2&1&\dots&&0&1\end{array}\right).$$

$$f_{rec}=\arg\min_{f}\Big{\{}\|Af-b\|_{2}^{2}+\lambda\|\Phi_{k}f\|_{1}\Big{\}}.$$

$$f(x)\approx \sum _{m=0}^M{\alpha }_m(t)\frac{(x-t{)}^m}{m!},$$

$$({\varphi }_{k,j}{)}_m=$$

$$\{\begin{array}{ll}(-1{)}^k& \text{mbox}ifm=0\\ 0& \text{mbox}if1\le m\le N-j(k+1)\\ (-1{)}^{k+\lfloor \frac{N-m}{j}\rfloor }(\genfrac{}{}{0ex}{}{k}{\lfloor \frac{N-m}{j}\rfloor })& \text{mbox}ifN-j(k+1)<m<N\end{array}.$$

$$f\ast {\varphi }_{k,j},$$

$$\Phi_{2,2}=\left(\begin{array}[]{c c c c c c c c c}1&1&-2&-2&1&1&0&\dots&0\\ 0&1&1&-2&-2&1&1&\dots&0\\ 0&0&1&1&-2&-2&1&\dots&0\\ \vdots&&&\ddots&&&&\dots&\vdots\\ -2&-2&1&1&0&0&\dots&1&1\\ 1&-2&-2&1&1&0&\dots&0&1\end{array}\right).$$

$$f_{rec}=\arg\min_{f}\Big{\{}\|Af-b\|_{2}^{2}+\frac{\lambda}{\ell+1}\sum_{j=0}^{\ell}2^{-(j+k-1)}\|\Phi_{k,2^{j}}f\|_{1}\Big{\}}.$$

$$f=\sum _{t=0}^{2^{\ell }-1}⟨f,{\phi }_{k,\ell ,t}⟩{\phi }_{k,\ell ,t}+\sum _{j=\ell }^{n-1}\sum _{t=0}^{2^j-1}⟨f,{\psi }_{k,j,t}⟩{\psi }_{k,j,t}.$$

$${\tilde{\psi }}_{k,j,t}(x)=2^{j/2}{\psi }_k\left(2^j(x-t{2}^{-n})\right),t=0,1,\dots ,2^{n-1}.$$

$$f$$

$$=\sum _{t=0}^{2^{\ell }-1}⟨f,{\phi }_{k,\ell ,t}⟩{\phi }_{k,\ell ,t}+\sum _{j=\ell }^{n-1}{2}^{j-n}{\Psi }_{k,j}^T(f\ast {\psi }_{k,j,0}(-x)),$$

$$\sum _{j=\ell }^{n-1}\Vert {\Psi }_{k,j}f{\Vert }_1,$$

$$f\ast {\varphi }_{k,j}=F^{-1}\left(F(f)\cdot F({\varphi }_{k,j})\right),$$

$$F({\varphi }_{k,j}{)}_{\xi }=\frac{{\left(e^{\frac{i2\pi \xi j}{N}}-1\right)}^{k+1}}{e^{\frac{i2\pi \xi }{N}}-1},$$

$$F({\varphi }_{k,j}{)}_{\xi }=\sum _{m=0}^{N-1}({\varphi }_{k,j}{)}_m{e}^{\frac{-i2\pi \xi }{N}m}.$$

$$F({\varphi }_{k,j}{)}_{\xi }=\sum _{n=0}^{j(k+1)-1}(-1{)}^{k+\lfloor \frac{n}{j}\rfloor }\left(\genfrac{}{}{0ex}{}{k}{\lfloor \frac{n}{j}\rfloor }\right)e^{\frac{-i2\pi \xi }{N}(N-n)},$$

$$F({\varphi }_{k,j}{)}_{\xi }=\sum _{m=0}^k\left((-1{)}^{m+k}(\genfrac{}{}{0ex}{}{k}{m})\left[\sum _{\ell =0}^{j-1}{e}^{\frac{i2\pi \xi }{N}(jm+\ell )}\right]\right).$$

$$F({\varphi }_{k,j}{)}_{\xi }$$

$$=\sum _{\ell =0}^{j-1}{\left(e^{\frac{i2\pi \xi }{N}j}-1\right)}^k{e}^{\frac{i2\pi \xi }{N}\ell }.$$

$$p_{m,n}=\{\begin{array}{ll}1& \text{mbox}ifm=n\\ 1& \text{mbox}ifn=(m+j-1)\mathrm{mod}(N)+1\\ 0& \text{mbox}ifotherwise\end{array}.$$

$$p_{m,n}(j,k)=\{\begin{array}{ll}\left(\genfrac{}{}{0ex}{}{k+1}{\ell }\right)& \text{mbox}ifn=(m+j\ell -1)\mathrm{mod}(N)+1\\ 0& \text{mbox}ifotherwise\end{array},$$

$${\Phi }_{k,2^j}=P_j^{k+1}{P}_{j-1}^{k+1}\cdots P_1^{k+1}{\Phi }_k$$

$${\Phi }_{k,2^j}=P_j^{k+1}{\Phi }_{k,2^{j-1}}.$$

$$P_2$$

$$P_2^2$$

$$Rf(t,\theta )={\int }_{\Omega }f(x,y)\delta (t-(x,y)\cdot (\cos \theta ,\sin \theta ))dxdy,$$

$$f_{rec}=\arg \underset{f}{\min }\Vert Tf{\Vert }_1\text{s.t.}Af=b,$$

$$(-1{)}^p(\genfrac{}{}{0ex}{}{k}{p})=\sum _{j=0}^{\ell }(-1{)}^j\left(\genfrac{}{}{0ex}{}{k}{j}\right)\left(\genfrac{}{}{0ex}{}{k+1}{\ell -j}\right),$$

Full text

1 Introduction

Over the past couple of decades, $\ell_{1}$ regularization techniques such as total variation have become increasingly popular methods for image and signal denoising and reconstruction problems. Along with TV [28], a large variety of approaches for similar $\ell_{1}$ regularization approaches have been proposed for an array of problems. Signal and image recovery methods continue to attract a great deal of interest due to the wide variety of potential applications and ever increasing means of various sensing mechanisms to acquire data. To name a few, synthetic aperture radar (SAR) [46, 2], magnetic resonance imaging (MRI) [23, 24, 25], electron tomography[21, 33], and inpainting [32, 19] are all image recovery applications that have advanced in part due to $\ell_{1}$ regularization methods, and in each case the approach can be tailored to the challenges that the particular application poses. With many problems such as MRI and electron tomography, the challenge is often to acquire as little data as necessary due to possible damage of the subject being imaged or because of time constraints, driving the need for inverse methods that can achieve the absolute best results from very limited and noisy data.

The mathematical description of the general problem we are interested in is to recover a signal or image $f\in\mathbb{R}^{N}$, from noisy measurements $b$ of the form $b=Af+\epsilon$, where $A\in\mathbb{R}^{m\times N}$ is some sensing matrix that approximates the physical model of the particular problem. Then the $\ell_{1}$ regularized solution is given by

[TABLE]

where $T$ is some sparsifying linear transform and $\lambda$ is a parameter that balances the effects of the data and regularization terms. The appropriateness of this approach is that some prior knowledge of the signal suggests that $Tf$ is sparse, and that the formulation with the $\ell_{1}$ norm encourages such sparsity [12, 5, 6]. In many applications, some knowledge of the appropriate transform is available, particularly with images and for other signals, this knowledge is in the form of some “smoothness.”

In the case of TV, the sparsifying transform is given by $T:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N-1}$, where $(Tf)_{i}=f_{i+1}-f_{i}$. The general idea for this approach is that the signal $f$ is assumed to be piecewise constant with a few discontinuities, in which case $Tf$ is sparse. If this is not precisely true, this approach still effectively reduces unwanted oscillations at the cost of the well documented stair-casing effect [8, 3]. However, for more general piecewise smooth functions higher order TV (HOTV) regularization methods are effective [8, 4, 17], and they do not suffer from the stair-casing effects. In this case the transform maps $f$ to approximations of discrete derivatives of $f$, e.g. higher order finite differences of $f$.

Another popular choice for $T$ are wavelet transforms [36, 26, 23]. For instance, such a transform can be written as $T:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$, where $(Tf)_{j}=\langle f,\psi_{j}\rangle$ and $\psi_{j}$ are orthonormal so that $f=\sum_{j}\langle f,\psi_{j}\rangle\psi_{j}$. The idea here is that for appropriately smooth signals, most of the signal’s energy is captured in the few low frequency, larger scaled elements of the basis. Thus most of the coefficients can be neglected, and thus a sparse approximation of $f$ exists with respect to the basis.

1.1 Discussion and Contribution

The crux of general $\ell_{1}$ regularization methods is that recovering a signal with the most sparse representation, that is recovering the solution with the smallest so called $\ell_{0}$ norm, is often equivalent to its convexly relaxed variant of recovering the signal with the smallest $\ell_{1}$ norm, which is a field of study called compressed sensing (CS) [12, 5, 6]. Although convex $\ell_{1}$ optimization algorithms are useful in promoting sparsity, some small nonzero coefficients may still persist, an obvious sign that the assumptions needed for the exactness guarantees given by CS theory sometimes do not hold in practice. This observation is largely the original motivation our present work in developing a multiscale HOTV approach related to multiscale wavelet regularization.

Much work has been devoted to understanding and developing sparsity promoting regularization methods, which are related to our current work. Numerous variants of higher order TV methods have been proposed [8, 1, 17]. For example, in [1] the authors propose an edge detection operator that annihilates polynomials, which leads them to operators close to finite difference matrices. In [8] a combination of a TV regularizer with a quadratic second order regularizer is developed in the continuous domain to eliminate staircasing effects. Likewise, several authors have shown that using some combination of first and second order methods to be beneficial [37, 4, 35, 7]. Unfortunately, since there are multiple regularization terms these methods typically introduce additional parameters that need to be tuned. In terms of theory, it has been well documented that under certain conditions TV and HOTV are equivalent to reconstruction with splines [44, 38], i.e. the solution of such methods recovers a piecewise polynomial with a sparse set of jumps.

TV denoising in particular has several very interesting equivalences. It is well known that TV denoising and other more general first order denoising methods are equivalent to smoothing with a certain nonlinear diffusion models[34], a typical result of writing the equivalent Euler-Lagrange equations. Perhaps discussed less frequently and most related to the observations in our current work, TV denoising is equivalent to soft threshold denoising with the highest frequency basis elements of the Haar wavelets [39, 40], in particular with the so called cycle spinning [18]. In general however, the main difference between these methods is that with TV the smoothness analysis is limited to the finest scales, whereas wavelet regularizations promote function smoothness at multiple scales. A main contribution of this article is to expand further on the relationship between wavelets in $\ell_{1}$ regularization and those $\ell_{1}$ methods related to HOTV. In regards to extension of wavelets, a number of multidimensional generalizations have been invented including curvelets and shearlets [15, 20, 36], which are primarily used for sparse function approximation and improve the approximation rates in two and three-dimensions compared with their one-dimensional counterparts.

The method we develop here an alternative for HOTV regularization which we refer to as multiscale HOTV (MHOTV). In contrast to previous work, our approach considers combining both a multiscale approach and higher order TV methods for the class of image reconstruction problems. The motivation for such an approach is in observable sub par results due to the relaxation of the sparsity promotion through the $\ell_{1}$ norm, contrary to the aforementioned results with splines [44, 38]. In light of this, we determined this calls for analysis of the function behavior at multiple scales. As can be deduced, this multiscale strategy is similar to the treatment of wavelets, and we argue that our approach is indeed related to the use of Daubechies wavelets, with the main divergence coming in the orthogonality and/or frame conditions prescribed by the wavelets. Orthogonality may be unnecessary for general $\ell_{1}$ regularization techniques, although fundamental to thresholding denoising techniques, and the relaxation of this condition in our approach allows for better localization of the transform. In the development of MHOTV, we carefully address the computational concerns associated with our approach through the use of both the FFT and operator decompositions. We are able to show through several numerical examples that MHOTV provides a notable improvement to the current alternatives.

The organization of the remainder of the article is as follows. In section 2 we define the HOTV operators and the corresponding multiscale generalizations. We also motivate the approach via a numerical example, and make the connection with Daubechies wavelets. In section 2.1 we precisely define the MHOTV $\ell_{1}$ regularization model and give precise normalizations to deal with proper parameter selection. In section 3 we address the computational concerns associated with calculating MHOTV coefficients, devising two distinct ways that they can be calculated in an efficient manner. In section 4 we provide numerical results for 1-D and 2-D reconstruction problems, showing that MHOTV is an improvement to the original HOTV and the related Daubechies wavelets. Some proofs and definitions are provided in the appendix.

2 HOTV and Multiscale Generalizations

As an alternative to TV regularization, general order TV methods have been shown to be effective for $\ell_{1}$ regularization [8, 4, 33, 1]. The TV transform can be thought of as a finite difference approximation of the first derivative, thus annihilating a function in locations where the function is a constant, i.e. a polynomial of degree zero. Likewise, higher order TV transforms can be considered higher order finite difference approximations to higher order derivatives, thus annihilating the corresponding degree polynomials. With this in mind, we have the following definition:

Definition 1** (Finite Differences).**

Let $\phi_{k}\in\mathbb{R}^{N}$ be defined by

[TABLE]

Then for $f\in\mathbb{R}^{N}$, the periodic $k^{th}$ order finite difference of $f$ is given by

[TABLE]

where $*$ denotes the discrete convolution.

Remark 1**.**

The convolution in this definition (and in general) can be represented by multiplication with a circulant matrix $\Phi_{k}$, where each row of $\Phi_{k}$ holds a shifted version of $\phi_{k}$. An example of the matrix in the case $k=2$ is given in (3). Note that this definition uses a periodic extension of $f$ and can be ignored by dropping the last $k$ rows of $\Phi_{k}$.

[TABLE]

With this definition, the HOTV model can be said to recover

[TABLE]

Unfortunately for many real world imaging problems the equivalence between $\ell_{1}$ and $\ell_{0}$ may not hold in practice, yet the $\ell_{1}$ regularization still tends to encourage favorable solutions. In terms of the sparsity promoting transform, this means that the transform of the recovered function may not be truly sparse, but most of the values are instead relatively close to zero. For HOTV, this means that a local Taylor expansion of the recovered function will still contain some small nonzero higher order coefficients, yet essentially unobservable at the very local scale. In other words, at some point $t$, there exists a polynomial expansion of minimal degree of $f$ given by

[TABLE]

which holds for all $x$ within some small interval $I$ around the point $t$. Ideally a solution given by the order $k$ HOTV model recovers a solution so that the coefficients $\alpha_{m}(t)$ vanish for $m\geq k$. The $\ell_{1}$ model allows for these coefficients to remain, although very small, and the function still appears to essentially be a polynomial of degree less than $k$. However, when this behavior persists over many points at a larger scale, the result can be a function which looks more like a trigonometric polynomial rather than an algebraic one.

This phenomenon is demonstrated in Figure 1, where a piecewise polynomial of degree two was reconstructed from random noisy samples with 50% sampling111The number of samples is half the number of grid points. using TV and HOTV regularizations. The sampling matrix $A\in\mathbb{R}^{N/2\times N}$ is constructed so that a random 10% of its entries are set to be nonzero, where these nonzero values are uniformly distributed over $[0,1]$. The samples were corrupted with normaly distributed mean zero noise. Two different grid sizes are demonstrated, 256 and 1024, and it can be observed that these small oscillations become increasingly abundant with more grid points. However, in the bottom of the figure, the third order finite difference of the HOTV3 solution plotted in logarithmic scale shows that locally this oscillatory behavior results in almost exact low order polynomials, although very small amplitudes persist in the transformed domain and thus not truly sparse in the $\ell_{0}$ sense. Nevertheless, all regularization approaches should still be deemed useful, as evidenced by the least squares solution shown.

Due to this phenomena we propose a multiscale HOTV approach, which considers the regularization transform at multiple scales. The idea is that a larger stencil would penalize these oscillations even with the $\ell_{1}$ norm. As TV generalizes to the Haar wavelet by stretching and scaling of the elements, we propose the same with HOTV. To this end we give the following definition.

Definition 2** (Multiscale Finite Differences).**

Let $\phi_{k,j}\in\mathbb{R}^{N}$ be defined by

[TABLE]

Then for $f\in\mathbb{R}^{N}$, the periodic $k^{th}$ order finite difference of scale $j$ of $f$ is given by

[TABLE]

where $*$ denotes the discrete convolution.

Remark 2**.**

Again, this convolution can be represented as multiplication with a circulant matrix $\Phi_{k,j}$. An example of $\Phi_{k,j}$ in the case $k=2$ and $j=2$ is given in (7).

[TABLE]

2.1 MHOTV Reconstruction Model

We now present the general model for MHOTV reconstruction. Generally speaking, we still use the model presented in (1), where $A$ maps the unknown function $f$ to some perhaps noisy measurements given by $b$, from which we use to reconstruct $f$. Our sparsity promoting transforms are now given by the matrices $\Phi_{k,2^{j}}$, for $j=0,1,\dots,\ell$, where $\ell$ is the maximum scaling of the operator used and $k$ is the chosen order. Setting our maximum scaling to $\ell=0$ corresponds to the traditional HOTV approach. Although not completely necessary, we choose a dyadic scaling of the operators, similar to the treatment of wavelets. As with wavelets, we will show that this is convenient for computational purposes. Finally then our reconstruction model is given by

[TABLE]

The factor of $2^{-j}$ is a normalization that accounts for the increasing norms of each operator, which would otherwise weigh too heavily to the largest scaling operator 222This is akin to the dyadic scaling of the wavelet basis elements after the dyadic stretching.. The scaling of the parameter $\lambda$ by $\ell+1$ simplifies the selection of the parameter, which would otherwise need to be manually scaled by such a factor to account for the number of scales being used. By similar reasoning, the additional scaling by $2^{-k+1}$ is used to account for the order $k$ of the method [31].

2.2 Relationship to Daubechies Wavelets

Wavelets can be characterized as an orthonormal basis that is generated through a multiresolutional analysis [10, 26]. The multiresolutional analysis leads to the shifting and dyadic stretching and scaling of a single generating mother wavelet, analogous to our treatment of MHOTV by shifting and stretching of a single row or element of the matrices $\Phi_{k}$. From this very general characterization, there are a number of parameters in the design of the wavelets. For Daubechies wavelets the smoothness is characterized by the number of vanishing moments, i.e. the number of polynomial orders to which the wavelet is orthogonal. A wavelet with $k$ vanishing moments acts as a multiscale differential operator of order $k$. As a trade off, an increasing number of vanishing moments chosen for the wavelet basis results in an increase in the support of the wavelet functions, and Daubechies wavelets are designed to yield the orthonormal wavelet basis of minimal support given a selected number of vanishing moments [26].

To develop a basic mathematical description of a wavelet expansion, suppose we want to represent a uniform pixelated function with $2^{n}$ pixels on $[0,1]$ in terms of the wavelet basis. Then denoting our scaling function and mother wavelet with $k$ vanishing moments by $\varphi_{k}$ and $\psi_{k}$ respectively, we have the following orthonormal wavelet representation

[TABLE]

Here, $\psi_{k,j,t}(x)=2^{j/2}\psi_{k}\left(2^{j}x-t\right)$ and similarly for $\varphi_{k,j,t}$, i.e. shrinking and scaling of the of the generating wavelet functions. The parameter $\ell$ is a positive integer with $0\leq\ell\leq n$, and the value $n-\ell$ is said to be the number of scales in the wavelet expansion 333For $\ell=n$ it is understood that the second sum is removed.. With the representation in (9), the coefficients for the scaling functions in the first sum capture most of the energy of the signal, and the wavelet coefficients $c_{k,j,t}=\langle f,\psi_{k,j,t}\rangle$ vanish for values of $t$ where $f$ is a polynomial of degree $k-1$ over the support of $\psi_{k,j,t}$. For $\ell_{1}$ regularization, we only need to be concerned with regularization of the wavelet coefficients in (9), and thus the coefficients for the scaling functions in the first sum are not included in the regularization.

Connecting these ideas to HOTV, we see that these transforms are playing similar roles. Both are prescribed by the number of vanishing moments, or in the language of HOTV, the highest order polynomial that is annihilated by the approach. Furthermore, both are designed to yield minimal support given the number of vanishing moments. The crucial difference lies in the orthogonality condition prescribed by wavelets, which further increases the support of the wavelet elements. We emphasize again, that this condition is fundamental to compression and threshold denoising methods, but not necessarily useful with general image reconstruction problems.

Finally, one additional technique utilized for $\ell_{1}$ regularization and denoising as well is the use a wavelet frame by taking all possible shifts for each scaling of the wavelets, which is sometimes referred to as translational invariant cycle spinning [9, 42, 18]. This eliminates the lack of translation invariance of a wavelet basis that can otherwise result in unwanted artifacts near discontinuities. With this in mind, we may define the wavelet frame elements by

[TABLE]

Then the averaged wavelet frame representation of $f$ may be written as

[TABLE]

where $\Psi_{k,j}^{T}=(\tilde{\psi}_{k,j,0},\tilde{\psi}_{k,j,1},\dots,\tilde{\psi}_{k,j,2^{n}-1})$. Hence a wavelet approach promotes sparsity with respect to the vectors $f*\psi_{k,j,0}$, or equivalently with respect to $\Psi_{k,j}f$. Then a regularization norm in this setting takes the form

[TABLE]

which is analogous to our regularization norm in (8). For wavelets, the scalings are inherent to function definitions, and the dyadic stretching of the elements is indicated by $j$ as opposed to $2^{j}$.

The case when $\ell=n-1$ would be most closely related to the original HOTV, and for smaller values of $\ell$ the wavelets are more comparable to the MHOTV development in this article.

Since computing both MHOTV operators and wavelets coefficients are convolutional operations, we may visualize their corresponding filters in Fourier space, providing us another basis for comparison, which we have done in Figure 2. Each of these can be interpreted as high pass filters, where the higher levels pass increasingly lower frequencies. A very close similarity of the wavelet filters and MHOTV filters can be observed in Figure 2, providing a strong visual confirmation to our preceding discussion on the close relationship between the two.

3 Fast Calculation of MHOTV Operators

Calculation of traditional HOTV coefficients is a computationally inexpensive task, due to the sparsity of the matrix operator. However, with increasing dyadic scales the direct calculation increases exponentially. Due to this, in the proceeding section we develop two distinct approaches that show that these calculations can be carried out with linear increase in the flop count with respect to the number of scales used.

Fast computation of standard HOTV can be done in several ways. One can construct the sparse matrix $\Phi_{k}$ and perform matrix computations directly, a calculation with runtime of $kN$ flops. One could make use of other procedures, such as MATLAB’s “diff” command which requires the same flop count without storing the matrix. With MHOTV, these approaches become less appealing. With matrix construction, if one is using several scales, then several matrices need to be computed and stored, and the matrices become significantly less sparse for larger scales. The“diff” command cannot be implemented directly for larger scale HOTV operators.

Another alternative is to use the Fourier convolution theorem to perform the convolution operation via a product in Fourier space. For the traditional HOTV operators, this can be fairly slow compared with the matrix and “diff” approach, since the necessary two discrete Fourier transforms would require $\sim 2N\log_{2}N$ flops compared with the $kN$ flops for the alternative implementations. However, this method is relatively comparable for MHOTV, since the Fourier transforms only need to be computed once to determine the coefficients at all scales.

We outline two procedures for efficient calculation of MHOTV. First, we describe the Fourier approach, where we derive precise formulas for the MHOTV Fourier filters. Second, we describe an alternative efficient approach by decomposition of the MHOTV matrix operators.

3.1 Computation via Fourier Transforms

By the Fourier convolution theorem, the MHOTV operators can be computed as multiplications in Fourier space, i.e.

[TABLE]

where $F$ denotes the discrete Fourier transform. Although this can be numerically computed, it is a convenient to have an exact formula for the discrete Fourier transform of $\phi_{k,j}$. Moreover, analytic determination of $F(\phi_{k,j})$ allows us to generalize the MHOTV to fractional orders.

Proposition 1**.**

The DFT of the vector $\phi_{k,j}$ defined in (6) has an explicit expression given by

[TABLE]

for $\xi=0,1,\dots,N-1$.

Proof.

The expression for the $\xi^{th}$ Fourier coefficient in the DFT of $\phi_{k,j}$ is given by

[TABLE]

Notice that the terms $1\leq m\leq N-j(k+1)$ vanish by definition of $\phi_{k,j}$. For the latter terms, we make the substitution $n=N-m$ and flip the sum to give the expression

[TABLE]

where the term $n=0$ corresponds to $j=0$ and the following indices $n=1,2,\dots,m(k+1)-1$, correpsond to $j=N-1,N-2,\dots,N-m(k+1)+1$, respectively. Notice that we may drop the $N$ in the numerator of the exponential and that the values of $\phi_{k,j}$ repeat over strings of length $j$. Therefore each of these corresponding strings of exponential terms in (13) get the same weights, leading to the following sum:

[TABLE]

Here the inner sum represents the $j$ consecutive terms in (13) that receive the same weights from $\phi_{k,j}$, namely $(-1)^{m+k}{k\choose m}$. Switching the order of summation, we recognize the sum over $m$ as a binomial expansion leading to

[TABLE]

The remainder of the proof follows by elementary calculations. ∎

3.2 Fast Computation via Operator Decomposition

In this section, we give a decomposition for the matrix operator $\Phi_{k,2^{j}}$ and describe how this decomposition can be used for rapid calculation of MHOTV operators. The decomposition of $\Phi_{k,2^{j}}$ is given in the following theorem.

Theorem 1**.**

Let the matrix $P_{j}$ with entries $\{p_{m,n}\}_{m,n=1}^{N}$ be defined by

[TABLE]

Then the following holds:

The entries of $P_{j}^{k+1}$, which we denote by $\{p_{m,n}(j,k)\}_{m,n=1}^{N}$, are given by

[TABLE]

where it is implied $\ell$ is an integer satisfying $0\leq\ell\leq k+1$. 2. 2.

$\Phi_{k,2^{j}}$* has the decomposition*

[TABLE]

and therefore

[TABLE] 3. 3.

The equality in (17) holds for any rearrangement of the product of matrices.

The proof of this theorem is given in the appendix. The matrices $P_{2}$ and $P_{2}^{2}$ are shown below to illustrate the sparse structure of these operators:

[TABLE]

Proposition 2**.**

Direct calculation of $\Phi_{k,2^{j}}$ requires $2^{j}Nk$ flops. The same calculation using the decomposition in (17) requires $jN(k+1)+Nk$ flops. The same calculation using the Fourier method requires $2Nlog_{2}N+N$.

Proposition 2 is a direct result of Theorem 1, the Fourier convolution theorem combined with the FFT, and the flops required for the direct calculation. We assume that the FFT and inverse FFT can be computed in $N\log_{2}N$ flops, although the exact count is somewhat vague, depending on the precise algorithm and if $N$ is a power of 2. To compute the full set of operators, we can get away with less flops then adding the flops for each level. If we use the decomposition approach to calculate the operators as determined by (17), the associated computations are limited to that at the highest scale, since the intermediate scales are determined in this calculation as pointed out in (18). If we use the Fourier approach for calculating the coefficients in (8), only one foward FFT is required for the function $f$. Then the product of $F(f)$ and $F(\phi_{k,2^{j}})$ must be computed for each $j$, as well as the inverse FFT for each of these products. The observations lead to the following corollary.

Corollary 1**.**

Let $T$ be the matrix containing the complete set of $\ell+1$ operators involved in the MHOTV $\ell_{1}$ regularization norm, so that $T^{T}=[\Phi_{k,1}^{T},\Phi_{k,2}^{T},\dots,\Phi_{k,2^{\ell}}^{T}].$ Then calculating $T$ using the operator decomposition given in Theorem 1 requires $\ell N(k+1)+Nk$ flops. Calculating $T$ using the Fourier approach requires a total flop count of $(\ell+2)N\log_{2}N+(\ell+1)N$.

A few concluding remarks are in order.

Remark 3**.**

All of the results presented are for 1-D signals. For higher dimensions say 2-D, the operators can be applied along each row and column, and the flop count is only doubled, disregarding the likely increased number of indices.

Remark 4**.**

To solve (8), we use the well establised alternating direction method of multipliers (ADMM) [22, 14, 45]. This approach introduces splitting variables that allows one to split the objective functional into equivalent subproblems that can be solved relatively fast. Our algorithm can be downloaded at [29], and some of the simulations in the proceeding section can also be found there.

4 Numerical Experiments

4.1 Repeat of 1-D Simulations

To compare MHOTV and wavelet regularized reconstructions we repeat the numerical examples presented in Figure 1 with the same noisy data used for the HOTV reconstruction. The corresponding reconstruction with MHOTV and wavelets are presented in Figure 3. Recall that the measurements were collected at a 50% sampling rate and corrupted with normally distributed mean zero noise. For the multiscale HOTV and wavelets, three scaling levels were used. The selection of the regularization parameter $\lambda$ was set to the same value for each order for HOTV and the wavelets, where we used a similar normalization approach for the wavelets coefficients as presented in (8).

The results in Figure 3 were generated with orders 1, 2, and 3. The order is indicated with the numbers next to the approach in the legends, e.g. we denote the order $k=3$ MHOTV approach with MHOTV3. For a baseline comparison, the least squares solution is shown as well. Compared with the corresponding $1024$ reconstructions from HOTV in Figure 1, these solutions show clear improvements, particularly with the higher orders. As we expect, although the MHOTV1 and Haar wavelet coefficients are computed in a different manner, the resulting reconstruction are identical since the models are theoretically equivalent. They both exhibit the staircasing and noise effects in precisely the same locations. The higher order approaches also show many similar effects of the noisy features, exhibiting certain oscillatory features with the same general behavior in precisely the same locations. However, with the higher orders, these approaches are not equivalent and MHOTV provides regulatory information at finer scales due to the minimal support of the transform elements. The result appears to be a modest improvement in the resulting reconstructions.

Finally, in the bottom of the figure the third order finite difference of the MHOTV3 solution is plotted in logarithmic scale. Comparing this with the original HOTV3 finite difference in Figure 1, we observe that the solution exhibits much better sparsity with respect to this transform domain, as desired.

4.2 2-D Tomographic Simulations

In this section we investigate the regularization methods on the common 2-D tomographic image reconstruction problem [27]. The phantom test image is shown in Figure 4 (a). The data generated for tomography are 1-D line integrals of the image, well-known as Radon transform data. Formally, the Radon transform of a 2-D function or image $f$ is defined as

[TABLE]

where $\Omega$ is the image domain and $\delta$ is the Dirac delta function. As in many applications, the data collected for reconstruction are of the form known as parallel beam geometry. In this setting, the full knowledge of noisy $Rf(t,\theta)$ is known for some finite set of angles, $\theta.$444There is also a discretization over $t$, but it is small enough to ignore. In this numerical experiment, we use a total of 29 angles that are equally distributed across the full 180$\degree$ angular range, which are visualized as a sinogram in Figure 4 (b). Such a limited set of data is sometimes referred to as limited data tomography. Mean zero normally distributed noise was again added to the data values. Classically tomographic reconstruction from parallel beam geometry can be done by first transforming the data into Fourier space by the Fourier slice theorem, and then applying a chosen ramp filter to this data and the inverse Fourier transform. This direct approach, called filtered backprojection, is sensitive to noise and is shown in Figure 4 (c).

The problem can however be discretized and approximated by a set of linear equations $Af=b$ (see for instance [30] on pages 8-9 within section 1.5.), where $A$ is sparse matrix that is a discretized approximation of the Radon transform, $f$ is the vectorized image, and $b$ is a vector holding the data values. With this set up we can apply regularization models such as (1) and (8). We use a $512\times 512$ pixelated mesh for the image domain in this experiment. The results for applying these models with HOTV, MHOTV, and Daubechies wavelets all at orders one and three are shown in Figure 5. Each of the models are also supplemented with a nonnegativity constraint, $f\geq 0$, which is carried out with a projected gradient method. A baseline comparison obtained by a conjugate gradient least squares solver is also shown in the figure. To ensure accurate comparison between the methods by appropriate parameter selection and algorithmic convergence, the relative data errors defined by $\|Af-b\|_{2}/\|b\|_{2}$ are shown in the figure, and it confirms that each approach approximately fit the data equally well, with all of the errors contained within an interval of size 0.0129.

As has been observed previously [33], due to a number of reasons including undersampling, noise, fine details between the image features, and nature of the regularization, the order 1 solutions (TV) can leave the fine features under resolved, even though the underlying image is truly a piecewise constant that classical TV was originally designed to recover. Each of these order 1 images appear relatively similar, with the MHOTV and Daubechies approaches showing modest improvements in resolving some features. As in the 1-D case, the HOTV3 solution exhibits some small local oscillations that appear as noise in the image. However, this image, as well as the other order 3 approaches resolve the features notably more clear than the order 1 approaches. Both of the order 3 multiscale approaches appear less noisy than the HOTV order 3, while still maintaining a good approximation of the image features.

4.3 Quantitative Results

We performed two sets of simulations to compare the methods in a more quantitative manner. The first set of results presented here involved setting up 100 different test problems and then running all of our methods over each time for multiple noise levels, and the mean reconstruction error over all simulations is presented in Table 1, with the MHOTV resulting in the left of each column and Daubechies wavelets in the right of each column. It is important to note here, that the parameter $\lambda$ in 8 was optimized in every reconstruction to yield the solution that minimized the true error between the test signal and the reconstruction, making for objective comparisons. In order to set up each test problem, a 1D piecewise quadratic polynomial (presumably ideal for order 3) was randomly generated over a 1024 stencil, and the entries in sampling matrix $A\in\mathbb{R}^{1024\times 1024}$ and added noise to $b$ were randomly generated from a mean zero Gaussian distribution. Overall, these results show that MHOTV moderately outperforms Daubechies wavelets in each case, and remaining comparisons between the order and number of levels are generally consistent between MHOTV and the wavelets.

For the single level case (original TV and HOTV), the error generally increases for higher orders, contrary to the results in previous work [1]. Multiple scales show notable improvement for the higher orders, whereas they show a mild reduction in accuracy for order 1. The most benefit for both orders 2 and 3 is seen when using 3 levels, and order 2 actually outperforms order 3. Finally, using the fact that (12) gives us a way to compute fractional orders of the method, we present also the results from orders 1.5 and 2.5. These are notably worse than the integer orders, a testament to the fact that these fractional order derivates result in highly nonlocal differences 555To observe these nonlocal stencils, one can compute the inverse Fourier transform of (12) for fractional orders $k$..

In the second set of results presented here we ran a series of numerical simulations and measured the rate of successful recovery for each method as a function of the sampling rate. For each simulation we randomly generated a piecewise polynomial of specified maximal degree over a 1024 stencil. This function was randomly sampled at the specified sampling rate precisely as in the previous 1-D simulations in section 4.1, where the sampling rate is defined by the number of samples divided by the number of grid points. Each regularization procedure was then used for reconstruction, and the $\ell_{2}$ error between the true function and reconstructed functions is determined. If the error was less than $10^{-2}$, then the reconstruction was said to yield a successful recovery. This simulation was carried out for each sampling rate in 20 trials, and the fraction of those 20 trials that yielded successful recovery is set as our probability of success. In each case, the generated test functions had five discontinuities, and the location of the jumps were drawn randomly from a uniform distribution on the approximation interval.

No noise was added for these simulations, as this can make the likelihood of an exact recovery unlikely. Therefore, for this case our general $\ell_{1}$ model as a modification of (1) is given by

[TABLE]

and similarly for our specific MHOTV model in (8). This constrained data fitting problem is solved by reformulating as an unconstrained problem with an augmented Lagrangian function [16, 22].

The results for these simulations are shown in Figure 6. The results are separated in two ways, by the degree of the piecewise polynomial function that is sampled (varying along the rows) and the order of the regularization method (varying along the columns). In the first row are results for piecewise constant functions, in the second row are piecewise linear functions, and in the third row are piecewise quadratic functions. In all cases, the MHOTV yields the highest probability of success, regardless of the degree of the polynomial or order of the regularization, and the Daubechies wavelets success appears to generally lie somewhere between MHOTV and HOTV. The order 1 regularizations perform well only in the case of piecewise constant functions. On the other hand, the order 2 and 3 regularizations perform well for all function types, with order 2 again outperforming order 3 both with piecewise linear and quadratic signals.

5 Summary

HOTV circumvents the staircasing often observed in TV solutions and has been shown to be more effective for problems with fine features, where resolution can be improved by increasing the order of derivatives in the regularization term [33]. In some applications, however, high order derivatives promote solutions with spurious local oscillations, as shown in Figure 1. The MHOTV regularization we introduce in this work is shown to mitigate unwanted oscillations while maintaining the resolution power of high order regularization.

Although the theory for MHOTV reconstructions remains underdeveloped when compared to wavelets regularization [11, 43, 15, 20, 36, 13, 41], our experiments indicate that MHOTV can outperform wavelets regularization in practical applications. Figure 3, for instance, shows fewer spurious oscillations in the MHOTV reconstruction than for Daubechies wavelets penalization. A feature that can also be observed for the 2-D tomographic data. Moreover, our results show that MHOTV regularization requires fewer samples for successful reconstructions than for HOTV and wavelets. Computational efficiency is achieved by performing the transformation in Fourier space or by matrix decomposition, as derived in section 3. The associated matlab algorithms can be downloaded at [29], and some of the simulations in the proceeding section can also be found there.

Appendix A Proof of Theorem 1

Lemma 1**.**

Let $k,\ell\in\mathbb{Z}$ with $0\leq\ell\leq k$. Then we have the following Vandermonde-like identity:

[TABLE]

where $p=\ell/2$ for $\ell$ even and $p=(\ell-1)/2$ for $\ell$ odd.

Proof of lemma 1.

Consider the polynomial $p(x)=(1-x^{2})^{k}(1+x)$, which can be factored as $p(x)=(1-x)^{k}(1+x)^{k+1}$. Both representations can be expanded using the binomial sum giving

[TABLE]

by the first representation and

[TABLE]

by the second representation. Since (22) and (23) are equivalent for all $x$, the coefficients of any particular power of $x$ are equivalent, which is the equality we set out to prove. ∎

Proof of theorem 1.

Statement 3 is an immediate consequence of statement 2, since each matrix involved in the product is a convolution operator, and convolution operations are commutative and associative.

To prove statement 1, first observe that with increasing $m$, the nonzero entries in the rows of $P_{m}$ become increasingly spaced, and it easy to see that the general resulting product $P_{m}^{k+1}$ is essentially the same for each $m$ with different spacings between the nonzero entries. Thus it is enough to show statement 1 for $m=1$. In the case $k=1$, this calculation can be checked directly. So suppose 1 holds for some arbitrary $k$. Then we need to show that $(P_{1}P_{1}^{k+1})$ yields the desired result as defined by (16). It is fairly easy to see that the resulting entries of this product is simply the addition of two neighboring entries (modulo $N$) in $P_{1}^{k+1}$. Any such entries added together trivially yields the desired values, and the proper location of these values is also easy to confirm.

Similar arguments used for statement 1 also apply to statement 2. First, we can consider an inductive approach, over $m$, where we will need to show $\Phi_{k,2^{m+1}}=P_{m+1}^{k+1}\Phi_{k,2^{m}}$. Note that again due to the spacing of the entries of $P_{m}^{k+1}$, the argument for any arbitrary $m$ is parallel to that for $m=1$, with just different handling of the indices. Therefore the case for $m=1$ suffices for the inductive step, and the case for $m=1$ is an immediate consequence of the previous lemma. ∎

Appendix B Definitions

If $f,g\in\mathbb{R}^{N}$, then the convolution of $f$ and $g$ is given by

[TABLE]

where for indices of $g$ running outside of domain of $g$, a periodic extension of $g$ is assumed. The discrete Fourier transform (DFT) of $f$ is defined by

[TABLE]

and the inverse discrete Fourier transform (IDFT) of $f$ is given by

[TABLE]

Acknowledgements

This work is supported in part by the grants NSF-DMS 1502640, NSF-DMS 1522639 and AFOSR FA9550-15-1-0152.