Continuous-Domain Formulation of Inverse Problems for Composite Sparse-Plus-Smooth Signals

We present a novel framework for the reconstruction of 1D composite signals assumed to be a mixture of two additive components, one sparse and the other smooth, given a finite number of linear measurements. We formulate the reconstruction problem as a continuous-domain regularized inverse problem with multiple penalties. We prove that these penalties induce reconstructed signals that indeed take the desired form of the sum of a sparse and a smooth component. We then discretize this problem using Riesz bases, which yields a discrete problem that can be solved by standard algorithms. Our discretization is exact in the sense that we are solving the continuous-domain problem over the search space specified by our bases without any discretization error. We propose a complete algorithmic pipeline and demonstrate its feasibility on simulated data.


Introduction
In the traditional discrete formalism of linear inverse problems, the goal is to recover a signal c 0 ∈ R N based on some measurement vector y ∈ R M .These measurements are typically acquired via a linear operator H ∈ R M×N that models the physics of our acquisition system (forward model), so that Hc 0 ≈ y.The recovery is often achieved by solving an optimization problem that aims at minimizing the discrepancy between the measurements Hc of the reconstructed signal c and the acquired data y.This data fidelity is measured with a suitable convex loss function E : R M × R M → R, the prototypical example being the quadratic error E(x, y) = 1 2 x−y 2 2 .A regularization term is often added to the cost functional, which yields the optimization problem arg min where R is the regularization functional, L specifies a suitable transform domain, and λ > 0 is a tuning parameter that determines the strength of the regularization.The use of regularization can have multiple motivations: 1. to handle the ill-posedness of the inverse problem, which occurs when different signals yield identical measurements; 2. to favor certain types of reconstructed signal (e.g., sparse or smooth) based on our prior knowledge; 3. to improve the conditioning of the inverse problem and thus increase its numerical stability and robustness to noise.

Continuous-Domain Formulation
Until now, we have focused on the discrete setting, as it constitutes the vast majority of the inverse-problem literature for the obvious reason of computational feasibility.However, most real-world signals are inherently continuous.Therefore, when feasible, to formulate the inverse problem in the continuous domain is a natural and desirable objective.
In this work, we adapt the discrete approach of (2) to 1D continuous-domain composite signals by solving an optimization problem of the form min s1,s2 where s 1 , s 2 are the two components of the signal s = s 1 + s 2 : R → R, ν = (ν 1 , . . ., ν M ) : s → ν(s) ∈ R M is the continuous-domain linear forward model, and L 1 , L 2 are suitable continuously defined regularization operators.Typical choices are L i = D N0,i for i ∈ {1, 2}, where D is the derivative operator and N 0,i the order of the derivative.The regularization norm • M is the total-variation (TV) norm for measures, which is the continuous counterpart of the discrete ℓ 1 norm [20][21][22].We refer to this term as the generalized TV (gTV) regularizer due to the presence of the operator L 1 .Finally, • L2 is the usual norm over the space L 2 (R) of signals with finite energy; we refer to the corresponding term as the generalized Tikhonov (gTikhonov) regularizer, which promotes smoothness in combination with the operator L 2 .

Representer Theorems and Discretization
A classical way of discretizing a continuous-domain problem is to reformulate it as a finite-dimensional one by relying on a representer theorem that gives a parametric form of the solution.Prominent examples include representer theorems for problems formulated over reproducing-kernel Hilbert spaces (RKHS), which are foundational to the field of machine learning [23,24].As demonstrated in [25,Theorem 3], the minimization problem (3) over the component s 2 (with a fixed s 1 ) -i.e., gTikhonov regularization -falls into this category: the representer theorem states that there is a unique solution of the form where the additional component p 2 lies in the null space of L 2 (i.e., L 2 {p 2 } = 0), h m is a (typically quite smooth) kernel function that is fully determined by the choice of ν m and L 2 , and a m,2 ∈ R are expansion coefficients.Therefore, to solve the continuous-domain problem, one need only optimize over the a m,2 coefficients and the null-space component p 2 which lives in a finite-dimensional space.This leads to a standard finite-dimensional problem with Tikhonov regularization.Concerning the minimization over the component s 1 (gTV regularization), several representer theorems give a parametric form of a sparse solution in different settings [21,[26][27][28][29].The case of our exact setting is tackled by [25,Theorem 4], which states that there is an L 1 -spline solution of the form where where δ is the Dirac impulse), K is the number of atoms of s 1 which is bounded by K ≤ (M − N 0,1 ), N 0,1 being the dimension of the null space of L 1 , and p 1 lies in the null space of L 1 .For example, when L 1 = D N0,1 , s 1 is a piecewise polynomial of degree (N 0,1 − 1) with smooth junctions at the knots x k .These representer theorems have paved the way for various exact discretization methods.In the gTikhonov case, one can optimize over the a m,2 coefficients in (4) directly [25].For the gTV case (5), grid-based techniques using a well-conditioned B-spline basis [30] as well as grid-free techniques [31] have been proposed.

Our Contribution
In this work, we show that the representer theorems presented in the previous section can be combined into a composite one when dealing with Problem (3).More specifically, we prove that there exists a solution to (3) of the form s * 1 = s * 1 + s * 2 such that s * 1 is of the form (5) and s * 2 is of the form (4): a "sparse plus smooth" solution.Building on this representation, we propose an exact discretization scheme.Both components s i for i ∈ {1, 2} are expressed in a suitable Riesz basis as , where c i [k] are the coefficients to be optimized.This leads to an infinite-dimensional optimization problem reminiscent of the infinite-dimensional compressed sensing framework of Adcock and Hansen [32].
To solve this infinite-dimensional problem numerically, we cast it as a finite-dimensional problem under some mild assumptions.This requires a careful handling of the boundaries of our interval of interest.In our implementation, we choose basis functions , where β L is the B-spline for the operator L. B-splines are popular choices of basis functions [33][34][35], in large part due to their minimalsupport property.Indeed, β Li has finite support when L i = D N0 , and it is the shortest-support generating function of the space of uniform L i splines [36].We show that optimizing over the spline coefficients leads to a discrete problem similar to (2) of the form min where . This discretization is exact in the sense that it is equivalent to the continuous problem (3) when each component s i lies in the space generated by the basis functions {ϕ i,k } k∈Z .This is a consequence of our informed choice of these basis functions ϕ i,k .Moreover, the short support of the B-splines leads to well-conditioned H i matrices and, thus, to a computationally feasible problem.

Related Works
The use of multiple regularization penalties is quite common in the litterature.However, in most cases, each penalty is applied to the full signal instead of a component-wise [37][38][39][40][41][42].A prominent example of such an approach is the elastic net [43], which is widely used in statistics.The spirit of these approaches is however quite different from ours: the reconstructed signal is encouraged to satisfy different priors simultaneously.
Conversely, in (2), each component satisfies different priors independently from the others, which will give very different results.The model of Meyer [44] and its generalization by Vese and Osher [45,46] follow the same idea as Problem (3), with the important difference that they use calculus-of-variation techniques to solve it.There is a connection as well with the Mumford-Shah functional [47], which is commonly used to segment an image in piecewise-smooth regions.The main difference lies in the fact that the optimization is not performed over the different components of the signal, but over the region boundaries.Another difference is that these models assume that one has full access to the noisy signal over a continuum, whereas (3) assumes that we only have access to some discrete measurements specified by the forward model ν.

Outline
In Section 2, we give the necessary mathematical preliminaries to formulate our optimization problem.In Section 3, we formulate the continuous-domain problem and present our representer theorem, which is our main theoretical result.In Section 4, we explain our discretization strategy, which relies on the selection of suitable Riesz bases.Finally, in Section 6, we present experiments on simulated data.

Operators and Splines
The crucial elements of our formulation are the regularization operators L 1 and L 2 .In this section, we specify which type of operators are suitable in our framework.
Let S ′ (R) denote the space of tempered distributions, defined as the dual of the Schwartz space S(R) of infinitely smooth functions on R whose successive derivatives are rapidly decaying.Let F be the generalized Fourier transform with f F {f }.Then, it is a standard result in distribution theory that the frequency response L = F {L{δ}} of an ordinary differential operator L : S ′ (R) → S ′ (R) is a slowly increasing smooth function L : R → C [48, Chapter 7, §5].Moreover, for any f ∈ S ′ (R), we have that F {L{f }} = L f ∈ S ′ (R).Next, we require L to be spline-admissible in the sense of Definition 1.
Definition 1 (Spline-admissible operator).A continuous LSI operator L : S ′ (R) → S ′ (R) is spline-admissible if it verifies the following properties: • there exists a function of slow growth ρ L : R → R (the Green's function of L) that satisfies L{ρ L } = δ; The prototypical example of a spline-admissible operator is the multiple-order derivative L = D N0 for N 0 ≥ 1.Its causal Green's function is the one-sided power function ρ L = x N 0 −1 + (N0−1)!, where x + = max(0, x).The null space of L is the space of polynomials of degree less than N 0 .
A spline-admissible operator L specifies the family of L-splines provided in Definition 2.
Definition 2 (Nonuniform L-spline).Let L be a spline-admissible operator in the sense of Definition 1.A nonuniform L-spline with K knots where a k ∈ R is the amplitude of the kth singularity.The weighted sum of Dirac impulses in (7) is known as the innovation of s.The spline s can equivalently be written as where p ∈ N L .
For example, the operator L = D N0 leads to the well-known polynomial splines, which are piecewise polyomials of degree (N 0 − 1) and of differentiability class C N0−2 .

Sparse Component
The other crucial elements of our framework are the native spaces for each component.Let M(R) be the space of bounded Radon measures, which is known by the Riesz-Markov theorem [49,Chapter 6] to be the continuous dual of C 0 (R).The latter is the space of continuous functions vanishing at infinity, which is a Banach space when equipped with the supremum norm • ∞ .The sparsity-promoting regularization norm • M is defined for a tempered distribution w ∈ S ′ (R) as Practically, the two critical features of the • M norm are the following: 1. it generalizes the L 1 norm in the sense that w M = w L1 for any w ∈ L 1 (R); Accordingly, the native space for s 1 in (3) is defined as which is a Banach space when equipped with the direct-sum topology.It is also the largest space for which the regularization is well-defined.We refer to [50] for technical details on the construction of M L1 (R).

Smooth Component
The regularization norm • L2 for the smooth component s 2 in ( 3) is defined over the Hilbert space L 2 (R).
The corresponding native space of the smooth component s 2 is the Hilbert space

Boundary Conditions
Finally, to ensure the well-posedness of our optimization problem, boundary conditions need to be introduced for one of the two native spaces.Let N 0 = N L1 ∩ N L2 be the intersection of the null spaces.We introduce a biorthogonal system (φ 0 , p 0 ) for N 0 in the sense of [21,Definition 3].An example of a valid choice is given in Appendix .2.The search space with boundary conditions φ 0 is then given by

Continuous-Domain Inverse Problem
Now that the relevant spaces have been introduced, we present in Theorem 1 the optimization task that we use to reconstruct sparse-plus-smooth composite signals.This representer Theorem gives a parametric form of a solution of our optimization problem; the proof is given in Appendix .1.
Theorem 1 (Continuous-domain representer theorem).Let E : R M × R M → R be a nonnegative, coercive, proper, convex, and lower-semicontinuous functional.Let L 1 , L 2 be spline-admissible operators in the sense of Definition 1 and let ν = (ν 1 , . . ., ν M ) be a linear measurement operator composed of the M linear functionals , where N ν is the null space of ν (well-posedness assumption).Then, for any λ 1 , λ 2 > 0, the optimization problem has a solution (s * 1 , s * 2 ) ∈ S with the following components: for some K 1 ≤ (M − N 0,1 ), where p 1 ∈ N L1 , and a 1,k , x k ∈ R; where Moreover, for any pair of solutions A pleasing outcome of Theorem 1 is that it combines Theorems 3 and 4 of [25] into one.There is, however, an added technicality due to the boundary conditions φ 0 .The latter are necessary to ensure the well-posedness of Problem (13).Otherwise, for any (s * 1 , s * 2 ) ∈ S and p ∈ N 0 , we would have that(s * 1 + p, s * 2 − p) ∈ S which would imply that S is unbounded.Note, however, that these conditions do not restrict the search space, since

Exact Discretization
In order to discretize Problem (13), we restrict the search spaces M L1,φ0 (R) and H L2 (R).The standard approach to achieve this is to choose a sequence of appropriate basis functions {ϕ i,k } k∈Z that span the reconstruction spaces for i ∈ {1, 2} that are subject to the constraints V 1 (R) ⊂ M L1,φ0 (R) and V 2 (R) ⊂ H L2 (R).These continuous spaces are linked to discrete spaces V i (Z), the choices of which will be made explicit in (20) and (27).More precisely, there is a one-to-one mapping between them using the basis functions ϕ i,k .

Riesz Bases and B-Splines
For numerical purposes, a desirable property is that our basis functions satisfy the Riesz property.Riesz bases are highly important concepts in that they generalize orthonormal bases, while leaving more flexibility for other desirable properties such as short support [51].
Definition 3 (Riesz basis).A sequence of functions {ϕ k } k∈Z with ϕ k ∈ L 2 (R) is said to be a Riesz basis if there exist constants 0 < A ≤ B such that, for any c ∈ ℓ 2 (Z), we have that Popular examples of Riesz bases are B-spline bases, which are introduced in Definition 4.
Definition 4 (B-spline).The B-spline for a spline-admissible operator L is characterized by a finite-differencelike filter (d L [k]) k∈Z , and is defined as The criteria for choosing a valid filter d L for a general class of operators L are given in [52,Theorem 2.7].
The best-known example of a B-spline is the polynomial B-spline for the operator L = D N0 , whose filter A key feature of B-splines is that they are the L-splines with the shortest support or, when finite support is impossible, with the fastest decay.Moreover, by [52,Theorem 2.7], for a valid B-spline β L as specified by Definition 4, the sequence of functions {β L (• − k)} k∈Z forms a Riesz basis in the sense of Definition 3.
It is clear from (18) that the innovation of the B-spline is a sum of Dirac impulses given by

Choice of Basis Functions
We now present and discuss our choice for the basis functions ϕ 1,k and ϕ 2,k .

Sparse Component
For the sparse component, we choose basis functions with the digital-filter space is the largest possible native reconstruction space [53,Equation (22)].The choice of the basis functions ϕ 1,k is guided by the following considerations: • they generate the space of uniform L 1 splines.This conforms with Theorem 1, which states that the component s * 1 is an L 1 -spline; • they enable exact computations in the continuous domain.In particular, we have that L • the Riesz-basis property of B-splines leads to a well-conditioned system matrix, which is paramount in numerical applications.
B-splines are the only functions that satisfy all these properties.Based on these criteria, B-splines are thus optimal.

Smooth Component
At first glance, the most natural choice for ϕ 2,k is to select the basis functions suggested by (15) in Theorem 1: h m for 1 ≤ m ≤ M and a basis of N L2 , which yield a finite number M + N 0,2 of basis functions.However, this approach runs into the following hitches: • the basis functions h m are typically increasing at infinity, which contradicts the Riesz-basis requirement and leads to severely ill-conditioned optimization tasks [25]; • depending on the measurements operator ν, h m may lack a closed-form expression.
We therefore focus on these criteria, in a spirit similar to [54].The ϕ 2,k are chosen to be regular shifts of a generating function ϕ 2 , with ϕ 2,k = ϕ 2 (• − k) such that {L 2 {ϕ 2,k }} k∈Z forms a Riesz basis in the sense of Definition 3. Contrary to ϕ 1,k , these requirements allow for many choices of ϕ 2,k .In order to perform exact discretization, one then only needs to compute the following autocorrelation filter.
Proposition 1 (Autocorrelation filter for the smooth component).Let ϕ 2 be a generating function such that ϕ 2,k = ϕ 2 (• − k) form a Riesz basis.Then, the following two items hold: • the inner product L 2 {ϕ 2,k }, L 2 {ϕ 2,k ′ } L2 only depends on the difference (k − k ′ ).We can thus introduce the autocorrelation filter for any k, k ′ ∈ Z; Proof.The first item is proved with a simple change of variable in the integral that defines the inner product.
The second item is derived by observing that, for any c 2 , we have

Formulation of the Discrete Problem
The autocorrelation filter introduced in Proposition 1 enables us to discretize Problem (13) in an exact way in the V i (R) spaces.
Proposition 2 (Riesz-basis Discretization).Let ϕ i,k be chosen as specified in Section 4.2 for i ∈ {1, 2}, k ∈ Z, and The cost function is given by where d L1 is the finite-difference-like filter from Definition 4, ρ is defined in Proposition 1, and •, • ℓ2 is the inner product over ℓ 2 (Z).Then, Problem (23) is equivalent to a restriction of the search spaces M L1,φ0 (R) and H L2 (R) to the spaces V 1 (R) and V 2 (R) defined in (16), respectively, so that in the sense that there exists a bijective linear mapping Proof.By plugging the expansions s i = k∈Z c i [k]ϕ i,k into the cost function J , using the linearity of ν, we get the data-fidelity term of (24).Using (19), we readily deduce that L [30,Equation (25)].As for the second regularization term, we observe that using (22) for the last step.This proves the equivalence between Problems ( 25) and ( 23), up to the specified mapping which is indeed a bijective linear mapping due to the Riesz-basis properties of {ϕ 1,k } k∈Z and {ϕ 2,k } k∈Z .

Practical Implementation
We now discuss how to solve our discretized problem (23) in practice, which involves recasting it as a finitedimensional problem.

Finite Domain Assumptions
To solve problem (23) numerically in an exact way, we must make assumptions that will enable us to restrict the problem to a finite interval of interest.The first item is fulfilled for common one-dimensional operators L 1 such as ordinary differential operators [55] or rational operators [56].
The second assumption is natural and is often fulfilled in practice, for instance in imaging with a finite field of view.The support length T then roughly corresponds to the number of grid points in the interval of interest.Note that, for simplicity, we only consider integer grids.However, the finesse of the grid can be tuned at will by adjusting T and rescaling the problem over the interval of interest.

Choice of Basis Functions ϕ 2,k
For our implementation, we make a specific choice of basis functions ϕ 2,k for the second component, since many different choices satisfy the requirements of Section 4.2.We choose the L * 2 L 2 B-spline basis and ϕ 2,k = ϕ 2 (• − k), where L * 2 denotes the adjoint operator of L 2 .In addition to the items discussed in Section 4.2, this choice has the following advantages: • the generator ϕ 2 has a simple explicit expression that does not depend on the measurement operator ν; • the autocorrelation filter ρ also has a simple expression, as will be shown in Proposition 4; • in the special case of the sampling operator ν m = δ(• − x m ), where the x m are the sampling locations, this choice conforms with (15) in Theorem 1 since s * 2 is then an L * 2 L 2 -spline.Note, however, that we do not exploit the knowledge that s * 2 has knots at the sampling locations x m .

Formulation of the Finite-Dimensional Problem
Our choice of basis functions together with the assumptions in Section 5.1 enable us to restrict Problem (23) to the interval of interest I T .More precisely, we introduce the indices m i , M i ∈ Z for i ∈ {1, 2}; the range [m i . . .M i ] corresponds to the set of indices k for which Supp(ϕ i,k ) ∩ I T = ∅, so that the basis function ϕ i,k affects the measurements.Hence, the number of active basis functions (i.e., the number of spline coefficients Finally, we introduce the native digital-filter space which is a valid choice because V 2 (R) ⊂ H L2 (R).Indeed, we can verify that, for any c 2 ∈ V 2 (Z), the function . This is due to the finite support of both (d L2 * c 2 ) and b 1/2 , where the filter b 1/2 and the decomposition g = b 1/2 * d L2 are introduced in Proposition 4 in Appendix .3.(20), our choice of V 2 (Z) in (27) is not the largest valid space: there exist larger vector spaces such that V 2 (R) ⊂ H L2 (R).However, the support restriction implies that for any has a finite support.This is a desirable property both for simplicity of implementation and because it conforms with Theorem 1, since s * 2 in (15) also satisfies this property.Our specific choice of support for (d L2 * c 2 ) is guided by boundary considerations and will be justified in the proof of Proposition 3.
The restriction to a finite number of active spline coefficients leads to finite-dimensional system and regularization matrices.The system matrices are of the form The regularization matrix for the sparse component, denoted by , is of the form The second component requires a careful handling of the boundaries in order to achieve exact discretization.This leads to a more complicated expression for the associated regularization matrix, which is given in (61) in Appendix .3.Finally, we introduce the matrix A ∈ R N0×N1 associated to the boundary condition functionals φ 0 .Our choice of boundary condition functionals φ 0 is presented in Appendix .2.With this choice, the constraint φ 0 ( k∈Z c 1 [k]ϕ 1,k ) = 0 leads to N 0 linear constraints on the coefficients c 1 = (c 1 [m 1 ], . . ., c 1 [M 1 ]), which can be written in matrix form as Ac 1 = 0.In common cases, these constraints simply lead to the N 0 first coefficients of c 1 to be set to zero, which thus reduces the dimension of the optimization problem.
These matrices enable an exact discretization of Problem (23), as shown in Proposition 3, the proof of which being given in Appendix .4.

Proposition 3 (Recasting as a finite problem). Let ϕ
, and let the assumptions in Section 5.1 be satisfied.Then, Problem (23) is equivalent to the optimization problem where the cost function is given by The matrices H i and L i for i ∈ {1, 2} are defined in (28), (29), and (61).This equivalence holds in the sense that there exists a bijective linear mapping from S d to S between their solution sets.
The combination of Propositions 2 and 3 allows us to solve the continuous-domain infinite-dimensional problem ( 25) by finding a solution (c * 1 , c * 2 ) ∈ S of the finite-dimensional problem (30).We obtain the corresponding solution of by extending these vectors to digital filters (c * 1 , c * 2 ) ∈ S d (this extension is unique as specified by Proposition 3), which yields the continuous-domain reconstruction

Sparsification Step
Although problem (30) can be solved using standard solvers such as the alternating-direction method of multipliers (ADMM), there is no guarantee that such solvers will yield a solution of the desired form specified by Theorem 1, i.e., s * 1 is an L 1 -spline with fewer than (M − N 0,1 ) knots, and s * 2 is a sum of M kernel functions and a null space element.This is a particularly relevant observation for the first component since, at fixed second component s * 2 , only extreme-point solutions s * 1 of Problem ( 13) take the prescribed form [21].This problem can be alleviated by computing a solution (c * 1 , c * 2 ) to Problem (30), and then finding an extreme point of the solution set c extr 1 2 ), which leads to a solution (c extr 1 , c * 2 ) of the prescribed form.This is achieved by recasting the problem as a linear program and using the simplex algorithm [57] to reach an extreme-point solution [25,Theorem 7].

Experimental Validation
We now validate our reconstruction algorithm in a simulated setting.

Grid Size
We rescale the problem by a factor T so that the interval of interest I T is mapped into [0, 1].We tune the finesse of the grid (and the dimension of the optimization task) by varying T , which amounts to varying the grid size h = 1/T in the rescaled problem.

Ground Truth
We generate a ground-truth signal s GT = s GT 1 + s GT 2 .The sparse component s GT 1 is chosen to be an L 1 -spline of the form (8) with few jumps, for which gTV is an adequate choice of regularization, as demonstrated by (14) in our representer theorem.For the smooth component s GT 2 , we generate a realization of a solution s 2 of the stochastic differential equation L 2 s 2 = w, where w is a Gaussian white noise with standard deviation σ 2 by following the method of [58].The operator L 2 then acts as a whitening operator for the stochastic process s 2 .The reason for this choice is the connection between the minimum mean-square estimation of such stochastic processes and the solutions to variational problems with gTikhonov regularization L 2 s 2 2 L2 [23,59,60].

Forward Operator
Our model is the Fourier-domain cosine sampling operator of the form ν 1 (s) = 1 0 s(t)dt (DC term) and for 2 ≤ m ≤ M , where the sampling pulsations ω m are chosen at random within the interval (0, ω max ], and the phases θ m are chosen at random within the interval [0, 2π).Notice that ν m is a Fourier-domain measurement of the restriction of s to the interval of interest [0, 1], in conformity with the finite-domain assumption in Section 5.1.
For the data-fidelity term, we use the standard quadratic error E(x, y) = 1 2 x − y 2 2 .

Comparison with Non-Composite Models
We now validate our new sparse-plus-smooth model against more standard non-composite models.More precisely, for i ∈ {1, 2} we solve the regularized problems arg min with regularizers R 1 (f ) = L 1 {f } M (sparse model with native space X 1 = M L1 (R)) and R 2 (f ) = L 2 {f } L2 (smooth model with native space X 2 = H L2 (R)).We discretize these problems using the reconstruction spaces V i (R) described in this paper (without restricting V 1 (R) with the boundary conditions φ 0 ).The sparse model thus amounts to an ℓ 1 -regularized discrete problem which we solve using ADMM, while the smooth model has a closed-form solution that can be obtained by inverting a matrix.For this comparison, we choose regularization operators L 1 = D and L 2 = D 2 with M = 50 Fourierdomain measurements (cosine sampling with ω max = 100).We generate the ground-truth signal according to Section 6.1.2,with K 1 = 5 jumps whose i.i.d.Gaussian amplitudes have the variance σ 2 1 = 1 for s GT 1 .For the smooth component s GT 2 , we generate a realization of a Gaussian white noise w with the variance σ 2 2 = 100, such that L 2 {s GT 2 } = w.The measurements are corrupted by some i.i.d.Gaussian white noise n ∈ R M so that y = ν(s GT ) + n.We set the signal-to-noise ratio (SNR) between ν(s GT ) and n to be 50 dB.For all models, we use the grid size h = 1/2 7 .The regularization parameters are selected through a grid search to maximize the SNR of the reconstructed signal with respect to the ground truth.
The results of this comparison are shown in Figure 2. As expected, due to the fact that our sparse-plussmooth signal model matches the ground truth, our reconstructed signal yields a higher SNR (21.46 dB) than the sparse-only (21.07 dB) and smooth-only (18.17 dB) models.Moreover, our reconstruction is qualitatively much more satisfactory.As can be observed in the zoomed-in section, the sparse-only model is subject to a staircasing phenomenon in the smooth regions of the ground-truth signal, a well-known shortcoming of totalvariation regularization.Our reconstruction does not suffer from this phenomenon and is remarkably accurate in the smooth regions.In fact, in our reconstruction, most of the error with respect to the ground truth comes from a lack of precision in the localization of the jumps due to gridding, which is costly in terms of SNR but does not affect much the visual impression.Finally, the smooth-only model fails both visually and in terms of SNR, due to its inability to represent sharp jumps.0 0.2 0.4 0.6 0.8 Ground truth Sparse + smooth Sparse Smooth Smooth : SNR = 18.17 dB with λ = 10 −11 .

Conclusion
We have introduced a continuous-domain framework for the reconstruction of multicomponent signals.It assumes two additive components, the first one being sparse and the other being smooth.The reconstruction is performed by solving a regularized inverse problem, using a finite number of measurements of the signal.
The form of a solution to this problem is given by our representer This form justifies the choice of the search space in which we discretize the problem.Our discretization is exact, in the sense that it amounts to solving a continuous-domain optimization problem restricted to our search space.The discretized problem is then solved using our ADMM-based algorithm, which we validate on simulated data. .

Preliminaries
We extend the biorthogonal system (φ 0 , p 0 ) for N 0 to the biorthogonal systems ( φ1 , p1 ) and ( φ2 , p2 ) for N L1 and N L2 , respectively, where φi = φ 0 φ i and pi = p 0 p i for i ∈ {1, 2}.It is known from [50,Theorem 4] that any function s 1 ∈ M L1 (R) has the unique decomposition where is the pseudo-inverse operator of L 1 for the biorthogonal system ( φ1 , p1 ) [50,Section 3.2].Using this decomposition, we can equip the space M L1 (R) with the norm Finally, an element s 1 ∈ M L1 (R) is in the restricted search space M L1,φ0 (R) if and only if c0 = 0.
Similarly, for any s 2 ∈ H L2 (R), there is a unique decomposition where c 0 = φ 0 (s 2 ) ∈ R N0 , c 2 = φ 2 (s 2 ) ∈ R N0,2−N0 , and h ∈ L 2 (R).Consequently, the associated norm for the space H L2 (R) is defined as Existence of a Solution The first step is to prove that (13) has a minimizer.We do so by reformulating the problem as the minimization of a weak*-lower semicontinuous functional over a weak*-compact domain.We then prove the existence by relying on the generalized Weierstrass theorem.
We denote the cost at the trivial point (0, 0) as J 0 = J (0, 0) = E(0, y).Adding the constraint J (s 1 , s 2 ) ≤ J 0 does not change the solution set of the original problem, as it must hold for any minimizer of (13).So, from now on, we assume that the cost functional is upper-bounded by J 0 .This readily implies that The coercivity of E(•, y) implies the existence of a constant C 1 > 0 such that E(z, y) ≤ J 0 ⇒ z 2 ≤ C 1 .Together with (40), this yields Moreover, since ν is weak*-continuous over M L1 (R), it is also continuous.This is due to the fact that a Banach space (in this case, the predual of M L1 (R)) is isometrically embedded in its double dual [49].Moreover, by assumption, ν is continuous over H L2 (R).Hence, there exists a second constant C 2 > 0 such that Now, by taking and, together with ( 41) and ( 42), we deduce that By using the triangle inequality and the two bounds ( 46) and ( 44), we have Finally, the well-posedness assumption in Theorem 1 ensures the existence of a constant B > 0 such that Hence, by taking q = φ 1 (s 1 ) T p 1 + φ 0 (s 2 ) T p 0 + φ 2 (s 2 ) T p 2 (49) and by applying the Inequality (48), we have that Therefore, the original problem ( 13) is equivalent to the constrained minimization problem min where . The cost functional in (51), which is the same as in (13), is weak* lower-semicontinuous.Moreover, the constraint cube is weak*-compact in the product topology due to the Banach-Anaoglu theorem.Hence, (51) reaches its infimum, and so does (13).
Form of the Solution Let (s 1 , s2 ) be a solution of (13) and consider the minimization problem min Unser et al. have shown in [21] that (52) has a minimizer s * 1 of the form (14).One can also readily verify that (s * 1 , s2 ) is a minimizer of the original problem.Similarly, one can consider the minimization problem min It is known from [25,Theorem 3] that (53) has a minimizer s * 2 of the form (15). Again, (s * 1 , s * 2 ) is a solution of the original problem, which matches the form specified by Theorem 1.
Uniqueness of the Second Component To prove the final statement of Theorem 1, let us consider two arbitrary pairs of solutions ( f1 , f2 ) and ( f1 , f2 ) of Problem ( 13) and let us denote by J min their minimal cost value.The convexity of the cost functional yields that, for any α ∈ (0, 1) and The optimality of ( f1 , f2 ) and ( f1 , f2 ) implies that (54) must be an equality.In particular, we must have that Now, due to the strict convexity of L 2 {•} 2 L2 , we deduce that L 2 { f2 − f2 } = 0, and hence that ( f2 − f2 ) ∈ N L2 .This implies that all solutions have the same second component up to a term in the null space of L 2 .
.2 Choice of Boundary Condition Functionals φ 0 We discuss here our choice of the boundary-condition functionals φ 0 for certain common choices of operators L i .We focus on multiple-order derivative operators L i = D N0,i , although the discussion remains valid for the more general class of rational operators [56], which, to the best of our knowledge, is the largest class of spline-admissible operators that satisfy the first assumption in Section 5.1.The null spaces N Li are thus the spaces of polynomials of degree smaller than N 0,i .We assume for now that we have N 0,1 ≤ N 0,2 , in which case we have N L1 ⊂ N L2 and thus N 0 = N L1 and N 0 = (D 1 − 1).Then, for any ǫ > 0, the functionals φ 0 = 1 ǫ rect • ǫ for N 0 = 1 and for N 0 > 1, where rect(t) = 1 for −1/2 ≤ t < 1/2 and 0 elsewhere, are valid choices of a biorthogonal system matched to the basis p 0 = 1, (•), . . ., (•) N 0 −1 (N0−1)!* δ(• − ǫ 2 ) of N 0 .Indeed, one can easily verify that this choice satisfies the biorthonormality relation φ 0,i , p 0,j = δ i−j (Kronecker delta).Moreover, we have φ 0,i ∈ X L1 (the predual of M L1 (R)), which implies that (p 0 , φ 0 ) is indeed a valid biorthogonal system of N 0 [50, Proposition 5].The fact that φ 0,i ∈ X L1 is proved in [61] for the case N 0 = 2; this proof can readily be extended to higher orders.
The boundary conditions (56), along with a choice of ǫ such that ǫ h is arbitrarily small, is numerically equivalent to φ 0 (f ) = (f (0), . . ., f (N0−1) (0 + )), where f (N0−1) (0 + ) is the right limit of f (N0−1) at 0. It can easily be shown in this case that, for s which leads to the constraints Ac 1 = (c 1,1 , . . ., c 1,N0 ) = 0 in Problem (30).This choice simplifies the optimization task by reducing the dimension of the problem, whereas other boundary conditions could lead to more complicated linear constraints and would make the optimization task more difficult.
So far, we have assumed that N 0,1 ≤ N 0,2 since this condition is always satisfied for the most common case L 1 = D.However, if N 0,1 > N 0,2 , then it is more convenient to apply the boundary conditions φ 0 to the second component, which leads to the same simple boundary conditions Ac 2 = (c 2,1 , . . ., c 2,N0 ) = 0.By default, we implicitly consider the more common case N 0,1 ≤ N 0,2 throughout the paper and thus impose the boundary conditions on the first component. .

Factorization of the Autocorrelation Filter
To specify the regularization matrix for the second component, we must first express in a convenient form the autocorrelation filter ρ defined in Proposition 1.This is done in Proposition 4, which gives the expression of ρ and its "square root" g for the choice of basis function ϕ 2 = β L * 2 L2 made in Section 5.2.Proposition 4 (Factorization of the autocorrelation filter).Let the assumptions in Section 5.1 be satisfied, and let Then, the basis {ϕ 2,k } k∈Z forms a Riesz basis as required in Section 4.2, and the autocorrelation filter ρ defined in Proposition 1 is of the form where b Proof.We have that ×HL 2 denotes the duality product between H L2 (R) and its dual H ′ L2 (R), and the last line results from the symmetry of ρ and b.
Next, we prove that b is positive-semidefinite.Indeed, for any finitely supported filter c, we have that where we have used the property  Table 1: Relevant filters and their supports (C = 2− √ 3 6 ).
We summarize in Table 1 the different filters and their mutual relations.Without loss of generality, we take the filters d Li for i ∈ {1, 2} to be causal, which leads to causal B-splines.These filters will be useful for the definition of the regularization matrix L 2 .

1 . 2 .
The operators L i for i ∈ {1, 2} admit a B-spline with finite support, which implies that the filters d Li (introduced in Definition (18)) and ρ (Proposition 4) have finite support.Without loss of generality, the support of d Li is chosen to be [0 . . .D i − 1] (of length D i > 0), which leads to causal B-splines β Li .The measurement functionals ν m are supported in an interval I T = [0, T ], where T ∈ N.