On the Autoconvolution Equation and

Total Variation Constraints^{∗}

Gunter Fleischer^{†}  Rudolf Gorenflo^{‡}  Bernd Hofmann^{§}

1991 Mathematics Subject Classifications: 65 J20, 45 G10, 65 R30 Keywords: autoconvolution, ill-posed problem, discretization, constrained least-squares approach, bounded total variation

Abstract

This paper is concerned with the numerical analysis of the autoconvolution equation x ∗ x = y restricted to the interval [0, 1]. We present a discrete constrained least squares approach and prove its convergence in L^{p}(0, 1), 1 ≤ p < ∞, where the regularization is based on a prescribed bound for the total variation of admissible solutions. This approach includes the case of non-smooth solutions possessing jumps. Moreover, an adaption to the Sobolev space H^{1}(0, 1) and some remarks on monotone functions are added. The paper is completed by a numerical case study concerning the determination of non-monotone smooth and non-smooth functions x from the autoconvolution equation with noisy data y.

∗The research was supported by the Deutsche Forschungsgemeinschaft (DFG) Grant Ho 1454/3-2 and by the Alexander von Humboldt Foundation, Bonn.

†Faculty of Mathematics, Technische Universität Chemnitz–Zwickau, D–09107 Chemnitz, Germany, email:

fleische@mathematik.tu-chemnitz.de

‡Department of Mathematics and Computer Science, Free University of Berlin, Arnimallee 2-6, D-14195 Berlin, Germany, email: gorenflo@math.fu-berlin.de

§Faculty of Mathematics, Technische Universität Chemnitz–Zwickau, D–09107 Chemnitz, Germany, email:

hofmannb@mathematik.tu-chemnitz.de

1 Introduction

In the paper [8] by Gorenflo and Hofmann the nonlinear ill-posed autoconvolution equation

∫_{0}^{s} x(s − t) x(t) dt = y(s), 0 ≤ s ≤ 1,  (1.1)

on the finite interval [0, 1] has been analyzed. This autoconvolution problem can be written as an operator equation

F (x) = y  (1.2)

with the continuous nonlinear operator F : D(F) ⊂ X → Y defined by

[F (x)](s) := [x ∗ x](s) := ∫_{0}^{s} x(s − t) x(t)dt, 0 ≤ s ≤ 1,  (1.3)

and mapping between Banach spaces X and Y with norms ∥.∥_{X} and ∥.∥_{Y}, respectively, containing real functions on the interval [0, 1]. In [8] there have been discussed intrinsic properties of the autoconvolution operator F from (1.3) and conditions for its compactness, injectivity and weak closedness, in particular for the Hilbert space X = Y = L^{2}(0, 1). As a consequence the general theory of Tikhonov regularization became applicable to equation (1.1). The character of ill-posedness in this equation strongly depends on the solution point x and its local degree of ill-posedness. Applications of the autoconvolution equation arising in physics and in stochastics are also mentioned in [8].

On the other hand, we discussed in a recent paper (cf. [5]) including numerical results the case that x is considered as a function of the space L^{2}(−∞, ∞) possessing a support supp x ⊂ [0, 1], where the complete data function [x ∗ x](s) (0 ≤ s ≤ 2) is observable. In such a case Fourier transform techniques are applicable and yield some more insight into the behaviour of the autoconvolution equation. However, the knowledge of [x ∗ x](s) for s > 1 is not always realistic. Therefore, in the present paper we are going to investigate stable approximate discretized solutions to (1.1), where both the function x to be determined and the data function y that can be measured are restricted to arguments from the interval [0, 1].

The approximate solution of the autoconvolution equation (1.1) will be based for Y := L^{2}(0, 1) on the restriction of admissible solutions x to compact subsets of the domain D(F) with prescribed properties. Provided that F is injective the inverse operator F^{−1} becomes continuous. We will show in Section 2 that a compactification of the autoconvolution equation in X := L^{p}(0, 1) can be based on a prescribed upper bound c for the total variation T (x) of solutions x, which are in addition uniformly bounded below and above by positive constants a and b, respectively. This allows us to construct convergent discretized solutions also in the case of non-smooth solutions possessing jumps. In this context, we generalize the well-known descriptive regularization approach using the set of monotone functions uniformly bounded below and above as a compact subset in L^{p}(0, 1), 1 ≤ p < ∞ (cf. Section 4, [13] and [4]). The total variation bound c plays in our consideration the role of a regularization parameter. In Section 3, the ideas of Section 2 are extended to the Sobolev space case X := H^{1}(0, 1). A brief reference to the case of monotone functions is given in Section 4. The paper is completed by a case study presented in Section 5 that illustrates the theoretical assertions of Section 2. In this case study the behaviour of discretized least-squares solutions to the autoconvolution equation subject to uniform bounds of the total variation is investigated, where both the case of a smooth and of a non-smooth solution are reflected.

2 Discretizing the Autoconvolution Equation under Total Vari-

ation Constraints

Let us consider the autoconvolution operator (1.3) between the Banach spaces X := L^{p}(0, 1) for fixed 2 ≤ p < ∞ with norm ∥x∥_{L^{p}(0,1)} = (∫_{0}^{1} |x(t)|^{p} dt)^{1/p} and Y := L^{2}(0, 1). In this context, we define the sets

D_{ε}^{+} := {x ∈ L^{p}(0, 1) : x(t) ≥ 0 a.e. in [0, 1], ε = sup{τ : x(t) = 0 a.e. in [0, τ]}} (2.1)

and

R_{ε}^{+} := {y ∈ L^{2}(0, 1) : y(s) ≥ 0 a.e in [0, 1], ε = sup{χ : y(s) = 0 a.e. in [0, χ]}}. (2.2)

Then we have the following proposition which, because of L^{p}(0, 1) being densely embedded in L^{2}(0, 1), follows from [8, Theorem 1 and Lemma 6] and [5, Proposition 2.5]:

Proposition 2.1 The autoconvolution operator F : L^{p}(0, 1) → L^{2}(0, 1) from (1.3) is a continuous nonlinear operator for all 2 ≤ p < ∞. In the restricted case F : D_{0}^{+} ⊂ L^{p}(0, 1) → R_{0}^{+} ⊂ L^{2}(0, 1) the operator is injective, but the autoconvolution equation (1.2) is locally illposed in the sense of Definition 2.2 in all points x ∈ D_{0}^{+}.

Definition 2.2 We call the equation (1.2) locally ill-posed in x ∈ D(F) if, for arbitrarily small r > 0 and balls B_{r} := {x̃ ∈ X : ∥x̃ − x∥_{X} ≤ r}, there is an infinite sequence {x_{k}} ⊂ D(F) ∩ B_{r}(x) with

∥F (x_{k}) − F (x)∥_{Y} → 0, but ∥x_{k} − x∥_{X} ↛ 0 as k → ∞.  (2.3)

Otherwise the equation is called locally well-posed in x ∈ D(F).

To overcome the difficulties of ill-posedness of a problem under consideration one can restrict the domain D(F) to a subset, which is compact in the Banach space X.

For a real function x(t) (0 ≤ t ≤ 1) we denote by

T (x) := sup_{0≤t_{0}<t_{1}<...<t_{k−1}<t_{k}≤1} ∑_{i=1}^{k} |x(t_{i}) − x(t_{i−1})|  (2.4)

the total variation of the function x on [0, 1] and by T_{S}(x) the analogously defined total variation of x on a closed subinterval S ⊂ [0, 1]. Note that the supremum in formula (2.4) is to be taken over all possible finite grids of the form 0 ≤ t_{0} < t_{1} < ... < t_{k−1} < t_{k} ≤ 1 with an arbitrarily chosen integer k. We consider, for given positive constants a, b and c, where

0 < a < b,  (2.5)

the domain

D := {x : [0, 1] → [a, b], T (x) ≤ c, x left-continuous for t ∈ (0, 1], x right-continuous for t=0} .  (2.6)

For technical reasons we assume that the lower bound a is strictly positive (see the remark after formula (2.21)). Obviously we have D ⊂ L^{p}(0, 1) for all 1 ≤ p < ∞. The requirement of the left- and right-continuity for the functions x ∈ D is reasonable, since a function of bounded variation has due to [12, Corollary 2, Chap. VIII, § 3] only a countable set of discontinuity points, namely jumps. Therefore, the left limit lim_{t→t_{0}−0} x(t) exists in all points of the interval (0, 1]. In the continuity points t_{0} this limit coincides with the value x(t_{0}). In all other points let be the values of x defined by x(t_{0}) := lim_{t→t_{0}−0} x(t). That means, with respect to L^{p}(0, 1)-elements we consider the representative, which is left-continuous in every point t ∈ (0, 1]. Moreover let x(0) := lim_{t→0+0} x(t), i.e. we consider no jumps at t = 0.

Lemma 2.3 The domain D from (2.5) – (2.6) is a compact subset of L^{p}(0, 1), 1 ≤ p < ∞, and we have D ⊂ D_{0}^{+}.

The proof of compactness of D is based on Helly’s theorem (cf. e.g. [12, Chap. VIII, §4]). For the proof ideas we refer to [4, Lemma 4.2]. On the other hand, note that Lemma 2.3 is a corollary of Theorem 2.5 in the paper [1] of Acar and Vogel. Namely, the set D from (2.5) – (2.6) is bounded with respect to the bv-norm

∥x∥_{BV [0,1]} := ∥x∥_{L^{1}[0,1]} + T (x).  (2.7)

Based on Lemma 2.3 providing compactness the following well-known Lemma of Tikhonov will allow us to prove stability results.

Lemma 2.4 Let F : D(F) ⊂ X → Y be a continuous and injective operator between the Banach spaces X and Y with a compact domain D(F). We denote by x^{∗}, for given righthand side y^{∗} ∈ F (D(F)), the unique solution of the operator equation (1.2). Then for a family of approximate solutions x_{η} ∈ D(F) the convergence of residual norms

∥F (x_{η}) − F (x^{∗})∥_{Y} → 0 as η → 0  (2.8)

implies the convergence of the approximate solutions

∥x_{η} − x^{∗}∥_{X} → 0 as η → 0.  (2.9)

A slightly modified version of this theorem and its proof can be found in Baumeister’s book [2, p. 18].

In order to obtain numerical approximate solutions, in the sequel we are going to discretize the autoconvolution equation (1.1) – (1.3), where the restriction of F to the compact subset D from (2.5) – (2.6),

F : D ⊂ L^{p}(0, 1) → L^{2}(0, 1),  (2.10)

is used. Similar to the discretization methods in [7] and [11], where also a total variation constraint is essential, we subdivide the interval [0, 1] into n subintervals I_{i} of the uniform length h := 1/n, where

I_{i} := ((i − 1)h, ih] (i = 1, ..., n).

For simplicity we set T_{i}(x) := T_{[(i−1)h,ih]}(x) for x ∈ D. Moreover, let

t_{j} := \frac{h}{2} + (j − 1)h (j = 1, ..., n)

denote the midpoints and

s_{i} := ih (i = 1, ..., n)

the right endpoints of such intervals.

To discretize the nonlinear integral equation (1.1), for all i, j = 1, 2, ..., n the values x(t_{j}) and y(s_{i}) will be approximated by some x_{j} and y_{i}, respectively. A discrete autoconvolution operator

\underline{F} : IR^{n} → IR^{n}  (2.11)

can be defined by

\underline{F} (\underline{x}) := (∑_{j=1}^{i} hx_{i−j+1}x_{j})_{i=1}^{n} , \underline{x} = (x_{1}, ..., x_{n})^{T} .  (2.12)

In its discrete form the autoconvolution equation then reads as

\underline{F} (\underline{x}) = \underline{y}, \underline{y} = (y_{1}, ..., y_{n})^{T} ,  (2.13)

or as

∑_{j=1}^{i} hx_{i−j+1}x_{j} = y_{i}, (i = 1, ..., n).  (2.14)

The realistic situation that the given data are noisy can be included. Instead of the exact data y_{i} for the right-hand side we will use perturbed data ŷ_{i}, where

∥ˆ\underline{y} − \underline{y}∥_{2} ≤ δ  (2.15)

and δ is a fixed upper bound for the noise of the data vector ˆ \underline{y} = (ŷ_{1}, ..., ŷ_{n})^{T}. Here we have used the scaled Euclidean norm

∥\underline{z}∥_{2} := (∑_{i=1}^{n} hz_{i}^{2})^{\frac{1}{2}}

for \underline{z} ∈ IR^{n}. For our further investigations we introduce the restriction operators

R : D ⊂ L^{p}(0, 1) → IR^{n} and Q : F (D) ⊂ L^{2}(0, 1) → IR^{n}

by

(R(x))_{j} := x(t_{j}) (j = 1, ..., n)  (2.16)

and

(Q(y))_{i} := y(s_{i}) (i = 1, ..., n),  (2.17)

as well as the extension operators E_{1} : IR^{n} → L^{p}(0, 1) and E_{2} : IR^{n} → L^{2}(0, 1) by

(E_{1}(\underline{x}))(t) := x_{j} (t ∈ I_{j}, j = 1, ..., n), (E_{1}(\underline{x}))(0) := x_{1}  (2.18)

and

(E_{2}(\underline{y}))(s) := y_{i} (s ∈ I_{i}, i = 1, ..., n), (E_{2}(\underline{y}))(0) := y_{1}.  (2.19)

We are searching now for an optimal solution vector

\underline{x}^{opt} = (x_{1}^{opt}, ..., x_{n}^{opt})^{T}

solving the discrete least-squares problem

∥\underline{F} (\underline{x}) − ˆ \underline{y}∥_{2} → min, subject to \underline{x} ∈ M,  (2.20)

where M is defined as

M := \{\underline{x} ∈ IR^{n} : 0 < a ≤ x_{i} ≤ b (i = 1, ..., n), ∑_{i=1}^{n−1} |x_{i+1} − x_{i}| ≤ c\} .  (2.21)

There exist solutions of (2.20), since M is compact in IR^{n} and ∥\underline{F} (\underline{x}) − ˆ \underline{y}∥_{2} : IR^{n} → IR^{1} is a continuous functional possessing a minimum over M. The condition 0 < a ≤ x_{i} ≤ b is more restrictive than the discretized version of x ∈ D_{0}^{+}. We require this stronger condition, because we want M to be a compact subset of IR^{n}.

For the vectors η := (δ, h)^{T}, \underline{x}^{opt} ∈ M and ˆ \underline{y} we define the piecewise constant function x_{η} ∈ D by

x_{η}(t) := E_{1}(\underline{x}^{opt})(t) (0 ≤ t ≤ 1).  (2.22)

and the piecewise constant function y_{δ} by

y_{δ}(s) := E_{2}(ˆ\underline{y})(s) (0 ≤ s ≤ 1).

Lemma 2.5 If we define the operator F_{η} : L^{p}(0, 1) → L^{2}(0, 1) by the formula

[F_{η}(x)](s) := ∑_{j=1}^{i} ∫_{I_{j}} x(s_{i} − t)x(t)dt (s ∈ I_{i}),  (2.23)

then we have the equation

∥F_{η}(ξ) − ζ∥_{L^{2}(0,1)}^{2} = ∥\underline{F} (\underline{ξ}) − \underline{ζ}∥_{2}^{2}  (2.24)

for all ξ := E_{1}(\underline{ξ}), where \underline{ξ} := (ξ_{1}, .., ξ_{n})^{T} ∈ IR^{n} and all ζ := E_{2}(\underline{ζ}) ∈ L^{2}(0, 1), where \underline{ζ} := (ζ_{1}, .., ζ_{n})^{T} ∈ IR^{n}.

Proof:

∥F_{η}(ξ) − ζ∥_{L^{2}(0,1)}^{2} = ∫_{0}^{1}([F_{η}(ξ)](s) − ζ(s))^{2}ds
= ∑_{i=1}^{n} ∫_{I_{i}} (∑_{j=1}^{i} ∫_{I_{j}} ξ(s_{i} − t)ξ(t)dt − ζ(s))^{2} ds
= ∑_{i=1}^{n} h (∑_{j=1}^{i} hξ_{i−j+1}ξ_{j} − ζ_{i})^{2}
= ∥\underline{F} (\underline{ξ}) − \underline{ζ}∥_{2}^{2}.

This proves the lemma Lemma 2.6 Let x ∈ D from (2.5) – (2.6). Then we have the estimation

∥F (x) − F_{η}(x)∥_{L^{2}(0,1)} ≤ 2hb^{2} + 2hbc.

Proof: We write

∥F (x) − F_{η}(x)∥_{L^{2}(0,1)} = (∑_{i=1}^{n} ∫_{I_{i}} (∫_{0}^{s} x(s − t)x(t)dt − ∫_{0}^{s_{i}} x(s_{i} − t)x(t)dt)^{2} ds)^{\frac{1}{2}} . (2.25)

Then we can estimate the expression in the inner parentheses by

|∫_{0}^{s} x(s − t)x(t)dt − ∫_{0}^{s_{i}} x(s_{i} − t)x(t)dt|

≤ |∫_{s_{i−1}}^{s} x(s − t)x(t)dt| + |∫_{0}^{s_{i−1}} x(s − t)x(t)dt − x(s_{i} − t)x(t)dt| + |∫_{s_{i−1}}^{s_{i}} x(s_{i} − t)x(t)dt|
≤ hb^{2} + ∑_{j=1}^{i−1} ∫_{I_{j}} |x(s − t) − x(s_{i} − t)||x(t)|dt + hb^{2}
≤ b ∑_{j=1}^{i−1} ∫_{I_{j}} |x(s − t) − x(s_{i} − t)|dt + 2hb^{2}.  (2.26)

Now we substitute u := s_{i} − t, du := −dt. For a fixed point t ∈ (s_{j−1}, s_{j}] = I_{j} we obtain u ∈ (s_{i−j}, s_{i−j+1}] = I_{i−j+1} and in view of −h ≤ s − s_{i} ≤ 0

s − s_{i} + u ∈ (s_{i−j−1}, s_{i−j+1}] = I_{i−j} ∪ I_{i−j+1}.

Moreover, we can estimate (2.26) by

b ∑_{j=1}^{i−1} ∫_{I_{j}} |x(s − t) − x(s_{i} − t)|dt + 2hb^{2}

= b ∑_{j=1}^{i−1} ∫_{I_{i−j+1}} |x(s − s_{i} + u) − x(u)|du + 2hb^{2}
≤ hb ∑_{j=1}^{i−1} (T_{i−j}(x) + T_{i−j+1}(x)) + 2hb^{2}
≤ hbT (x) + hbT (x) + 2hb^{2}
≤ 2hbc + 2hb^{2}.

Finally we substitute this estimation into equation (2.25). This yields the assertion of the lemma

Lemma 2.7 Under the assumptions stated above we have

∥F (x_{η}) − F (x^{∗})∥_{L^{2}(0,1)} ≤ 4hb^{2} + 6hbc + 2δ → 0 as η → 0.  (2.27)

Proof (for similar ideas see also [6]): From the triangle inequality we obtain

∥F (x_{η})−F (x^{∗})∥_{L^{2}(0,1)} ≤ ∥F (x_{η})−F_{η}(x_{η})∥_{L^{2}(0,1)}+∥F_{η}(x_{η})−y_{δ}∥_{L^{2}(0,1)}+∥y_{δ}−y∥_{L^{2}(0,1)}. (2.28)

The right-hand side of (2.28) consists of three terms which we want to estimate one by one: Due to Lemma 2.6 for the first term it holds

∥F (x_{η}) − F_{η}(x_{η})∥_{L^{2}(0,1)} ≤ 2hb^{2} + 2hbc (x_{η} ∈ D).

To estimate the second term of (2.28) we define \underline{x}^{∗} := R(x^{∗}) as the vector of the function values of the exact solution x^{∗} of the autoconvolution equation (1.1) in the midpoints of the intervals I_{i}. Since we have \underline{x}^{opt} as the least-squares solution of (2.20), the residual norm of \underline{x}^{∗} cannot be smaller than the residual norm of \underline{x}^{opt}. Furthermore, we can apply Lemma 2.5 with ξ := x_{η} and ζ := y_{δ}. This yields

∥F_{η}(x_{η}) − y_{δ}∥_{L^{2}(0,1)} = ∥\underline{F} (\underline{x}^{opt}) − ˆ \underline{y}∥_{2} ≤ ∥\underline{F} (\underline{x}^{∗}) − ˆ \underline{y}∥_{2}.

Using the identity

F_{η}(x) = E_{2}(Q(F (x))) (x ∈ D),

this allows us to estimate further as follows:

∥\underline{F} (\underline{x}^{∗}) − ˆ \underline{y}∥_{2} ≤ ∥\underline{F} (\underline{x}^{∗}) − Q(F (x^{∗}))∥_{2} + ∥Q(F (x^{∗})) − ˆ \underline{y}∥_{2}
= ∥F_{η}(E_{1}(R(x^{∗}))) − E_{2}(Q(F (x^{∗})))∥_{L^{2}(0,1)} + ∥\underline{y} − ˆ \underline{y}∥_{2}

= (∑_{i=1}^{n} ∫_{I_{i}} (∑_{j=1}^{i} ∫_{I_{j}}(x̃(s_{i} − t)x̃(t) − x^{∗}(s_{i} − t)x^{∗}(t))dt)^{2} ds)^{\frac{1}{2}} + δ
≤ (∑_{i=1}^{n} ∫_{I_{i}} (∑_{j=1}^{i} ∫_{I_{j}} |x̃(s_{i} − t)||x̃(t) − x^{∗}(t)| + |x^{∗}(s_{i} − t) − x̃(s_{i} − t)||x^{∗}(t)|dt)^{2} ds)^{\frac{1}{2}} + δ
≤ (∑_{i=1}^{n} ∫_{I_{i}} (∑_{j=1}^{i} ∫_{I_{j}} 2bT_{j}(x^{∗})dt)^{2} ds)^{\frac{1}{2}} + δ ≤ 2hbc + δ,

where x̃ := E_{1}(R(x^{∗})). The last inequalities essentially used Lemma 2.5 with ξ = E_{1}(R(x^{∗})) = x̃ and ζ = E_{2}(Q(F (x^{∗}))), respectively. Note that we have x̃(t) = x^{∗}(t_{j}) for t ∈ I_{j} and thereby |x̃(t) − x^{∗}(t)| ≤ T_{j}(x^{∗}). Taking into account |y_{i} − ŷ_{i}| ≤ δ and the identity

∥E_{2}(\underline{y})∥_{L^{2}(0,1)} = ∥\underline{y}∥_{2},

which can easily be proved, we hence can estimate the third term of (2.28) as follows (cf. Lemma 2.6):

∥y_{δ} − y∥_{L^{2}(0,1)} ≤ ∥y − E_{2}(Q(y))∥_{L^{2}(0,1)} + ∥E_{2}(Q(y)) − y_{δ}∥_{L^{2}(0,1)}
= ∥F (x^{∗}) − E_{2}(Q(F (x^{∗})))∥_{L^{2}(0,1)} + ∥E_{2}(Q(y)) − E_{2}(Q(y_{δ}))∥_{L^{2}(0,1)}
= ∥F (x^{∗}) − F_{η}(x^{∗})∥_{L^{2}(0,1)} + ∥Q(y) − Q(y_{δ})∥_{2} ≤ 2hb^{2} + 2hbc + δ.

Finally we can add the three terms and obtain by (2.28) the inequality (2.27). Evidently, the right-hand side of (2.27) tends to zero as h and δ both tend to zero. This proves the lemma

By the result of Lemma 2.7 we can apply Lemma 2.4 to prove in L^{p}-spaces the convergence of approximate solutions to the autoconvolution equation under total variation constraints.

Theorem 2.8 Consider the autoconvolution problem (1.1) – (1.3) with D(F) := D from (2.5)-(2.6) and denote by x^{∗} ∈ D, for given right-hand side y^{∗} ∈ F (D(F)), the unique solution of the autoconvolution equation. Then the family of approximate solutions x_{η} according to (2.22) converges to the solution x^{∗} of (1.2):

∥x_{η} − x^{∗}∥_{L^{p}(0,1)} → 0 as η → 0 for all 1 ≤ p < ∞.  (2.29)

Proof: In the case p ≥ 2 based on Lemma 2.7 the Lemma 2.4 immediately yields the convergence property (2.29), since the autoconvolution operator F : D ⊂ L^{p}(0, 1) → L^{2}(0, 1) is continuous and injective. Furthermore, D is a compact subset in L^{p}(0, 1) because of Lemma 2.3. For 1 ≤ p < 2 the norm ∥ · ∥_{L^{p}(0,1)} is ’weaker’ than the norm ∥ · ∥_{L^{2}(0,1)}. This ensures the convergence condition (2.29) also in this case

By using the method of Tikhonov regularization in Hilbert spaces X and Y the minimizers x_{α} of the auxiliary extremal problems

∥F (x) − y∥_{Y}^{2} + α∥x∥_{X}^{2} → min, subject to D(F)  (2.30)

with the regularization parameter α > 0 are exploited to find stable approximate solutions of an ill-posed operator equation (1.2). The smaller the regularization parameter α is chosen, the ’closer’ the original and the auxiliary problem are related, but the more instable and highly oscillating the solution of the auxiliary problem will become. In general, α has to be selected such that an appropriate trade-off between stability and approximation is realized. In our compactification approach using upper bounds c of the total variation the inverse value \frac{1}{c} plays a comparable role. In fact, if we consider small values \frac{1}{c}, then highly oscillating functions with large total variation values are admissible. On the other hand, for small values c the solutions obtained cannot oscillate very much, and the approximate solutions will be computed in a more stable way. However, if c is selected too small, then it may occur that the (unknown) exact solution is not an element of the set D. In such a case we would ’overregularize’ the autoconvolution equation. By controlling the upper bound c of total variation we are able to suppress oscillations. Compared to the frequently used compactification in L^{p} by using monotonicity constraints and lower and upper bounds for the function values (see Section 4) the approach of this section allows us to handle a more comprehensive class of (also non-monotone) functions. A numerical case study presented in Section 5 will illustrate the theoretical results of this section and some specific effects of the discretized solution of the autoconvolution equation under total variation constraints.

In the case p = ∞ we cannot assert convergence under our assumption of bounded total variation. If the solution x^{∗} has a jump point, then ∥x_{η} − x^{∗}∥_{L^{∞}(0,1)} → 0 as η → 0 is not true in general.

3 The Sobolev Space Case

In [8] it was already mentioned that the operator F of autoconvolution according to (1.3) mapping from X := L^{2}(0, 1) into the space Y := L^{2}(0, 1) is non-compact, but it becomes a compact operator if we change the problem to the Sobolev space X := H^{1}(0, 1) ≅ W_{2}^{1}(0, 1) of functions x with a quadratically integrable generalized derivative x′ and norm

∥x∥_{H^{1}(0,1)} = (∫_{0}^{1} |x(t)|^{2} dt + ∫_{0}^{1} |x′(t)|^{2} dt)^{1/2} .  (3.1)

In both cases the autoconvolution equation is locally ill-posed everywhere. But for compact operators F, we have in general a stronger form of ill-posedness. If our pairs of spaces X and Y are Hilbert spaces, following the concept of [8] (see also [9, Sect. 2.2.2]) we can express the local degree of ill-posedness μ (0 ≤ μ ≤ ∞) of the autoconvolution equation in a solution point x^{∗} by the decay rate of the singular value sequence σ_{1} ≥ σ_{2} ≥ ... ≥ σ_{i} ≥ ... > 0 tending to zero as i → ∞ of the Fréchet derivative F′(x^{∗}) in the form

μ := sup{ν : σ_{i} = O(i^{−ν}) as i → ∞},  (3.2)

where this linear operator given by F′(x^{∗})h = 2h ∗ x^{∗} is compact. Since the compact embedding operator from H^{1}(0, 1) into L^{2}(0, 1) has a sequence of singular values κ_{1} ≥ κ_{2} ≥ ... ≥ κ_{i} ≥ ... > 0 tending to zero with a rate κ_{i} ∼ 1/i as i → ∞, for the Sobolev space X := H^{1}(0, 1) under consideration in this section the ill-posedness degree grows at least by one (cf. Hofmann and Tautenhahn [10]) compared to the L^{2}(0, 1) case of Section 2. Thus, for a compactification in H^{1}(0, 1) ’stronger’ restrictions on the admissible solutions x are necessary. However, our aim in this section is also stronger, namely to obtain convergence of approximate solutions x_{η} to x^{∗} in the H^{1}(0, 1)-norm (3.1).

Here we consider, for given constants a_{1}, a_{2}, b_{1}, b_{2} and c with

0 < a_{1} < b_{1}, a_{2} < b_{2},  (3.3)

the domain

D := \{x : [0, 1] → [a_{1}, b_{1}], {∃x′ : [0, 1] → [a_{2}, b_{2}], \atop T (x′) ≤ c,} x′ left-continuous for t ∈ (0, 1], x′ right-continuous for t=0\} , (3.4)

where the function x′(t) (0 ≤ t ≤ 1) a.e. in [0, 1] coincides with a derivative of x(t) in the classical sense. Obviously we get D ⊂ H^{1}(0, 1) and hence every function x ∈ D with D from (3.3) – (3.4) is continuous. In analogy to Lemma 2.3 we have in the Sobolev space case:

Lemma 3.1 The domain D from (3.3) – (3.4) is a compact subset of H^{1}(0, 1) with D ⊂ D_{0}^{+}.

In contrast to the L^{p}-case the restriction of the total variation, here T (x′) ≤ c, is only needed to show the compactness of the domain D. It has no relevance for the convergence of the images F (x_{η}) of approximate solutions x_{η} to F (x^{∗}) in L^{2}(0, 1) as η tends to zero.

The discretization of the autoconvolution problem (1.1) – (1.3), where the operator F from (1.3) maps in the form

F : D ⊂ H^{1}(0, 1) → L^{2}(0, 1)  (3.5)

and where the domain D is defined by (3.3) – (3.4) will be performed similar to the L^{p}(0, 1) case. However, piecewise constant functions are not in H^{1}(0, 1). Therefore, we use continuous piecewise linear approximate functions. Here, let (in contrast to Section 2)

t_{j} := jh (j = 0, ..., n)

denote the n + 1 nodes subdividing the interval [0, 1], and again I_{j} = ((j − 1)h, jh]. Furthermore, the x_{j} again denote approximate values of x(t_{j}). As the discrete autoconvolution operator we introduce here:

\underline{F} : IR^{n+1} → IR^{n},  (3.6)

where \underline{F} (\underline{x}) = (z_{1}, ..., z_{n})^{T} and for i = 1, 2, ..., n :

z_{i} = ∫_{0}^{ih} (E_{1}(\underline{x})) (ih−t) (E_{1}(\underline{x})) (t) dt = ∑_{j=1}^{i} \frac{h}{6} (2x_{i−j}x_{j} + x_{i−j+1}x_{j} + x_{i−j}x_{j−1} + 2x_{i−j+1}x_{j−1}) .

(3.7)

By E_{1} : IR^{n+1} → H^{1}(0, 1) we denote in contrast to Section 2 the operator of piecewise linear interpolation according to

(E_{1}(\underline{x}))(t) := \frac{t − jh}{h}(x_{j} − x_{j−1}) + x_{j} (t ∈ I_{j}, j = 1, ..., n).  (3.8)

For noisy data (see (2.15)) we search for a minimizer

\underline{x}^{opt} = (x_{0}^{opt}, x_{1}^{opt}, ..., x_{n}^{opt})^{T}

of the least-squares problem (2.20) with M from

M := \{\underline{x} ∈ IR^{n+1} : {0 < a_{1} ≤ x_{i} ≤ b_{1} (i = 0, ..., n), \atop ha_{2} ≤ x_{i} − x_{i−1} ≤ hb_{2} (i = 1, ..., n),} ∑_{i=1}^{n−1} |x_{i+1} − 2x_{i} + x_{i−1}| ≤ hc\} .

(3.9)

With the same arguments as before it follows that (2.20) is solvable. The choice of \underline{F} is due to the fact that we have to guarantee the validity of formula (2.24) with F_{η} from (2.23).

By setting for the approximate solution

x_{η} := E_{1}(\underline{x}^{opt}),  (3.10)

where η = (δ, h)^{T}, we also have x_{η} ∈ D with D according to (3.3)-(3.4). Moreover, it can be shown that as in Lemma 2.7 we have ∥F (x_{η}) − F (x^{∗})∥_{L^{2}(0,1)} → 0 for η → 0. The proof dealing with the H^{1}(0, 1) approximation of functions by linear splines is omitted here. Using again Lemma 2.4 with X := H^{1}(0, 1) and Y := L^{2}(0, 1) we obtain:

Theorem 3.2 Consider the autoconvolution problem (1.1) – (1.3) with D(F) := D from (3.3) – (3.4) and denote by x^{∗} ∈ D, for given right-hand side y^{∗} ∈ F (D(F)), the unique solution of the autoconvolution equation. Then the family of approximate solutions x_{η} converges to the solution x^{∗} of (1.2):

∥x_{η} − x^{∗}∥_{H^{1}(0,1)} → 0 as η → 0.  (3.11)

4 Monotonicity Constraints

In this section we deal with solutions of the autoconvolution equation subject to the set of monotone and uniformly bounded functions considered as a particular subset of the functions possessing a bounded total variation.

First we consider the domain

D := {x : 0 ≤ x(t) ≤ b, t ∈ [0, 1], x non-increasing}  (4.1)

forming a compact subset in L^{p}(0, 1), 1 ≤ p < ∞. Then the operator F from (2.10) is also injective, since D ⊂ D_{0}^{+} ∪ {0} and x(t) = 0 (0 ≤ t ≤ 1) is the only function of D according to (4.1) with x(0) = 0. The discretization of this monotonicity case is completely the same as given in Section 2 for the total variation case with the exception of the fact that we have to introduce

M := \{\underline{x} ∈ IR^{n} : 0 ≤ x_{n} ≤ ... ≤ x_{1} ≤ b\}.  (4.2)

replacing (2.21). Since each monotone function is of bounded variation, we obtain the convergence results of Section 2 with c = b and a = 0.

Now we change to the case of non-decreasing solutions, where

D := {x : 0 ≤ x(t) ≤ b, t ∈ [0, 1], x non-decreasing}  (4.3)

and

M := \{\underline{x} ∈ IR^{n} : 0 ≤ x_{1} ≤ ... ≤ x_{n} ≤ b\}.  (4.4)

The set D from (4.3) is also compact in L^{p}(0, 1), but the injectivity of F fails (cf. [8]). Because of that we have to distinguish two cases:

On the one hand let y ∈ R_{0}^{+}, i.e. y(s) > 0 if s > 0. Then the corresponding solution x^{∗}(t) is uniquely determined from y a.e. in [0, 1] and ∥F (x_{η}) − F (x^{∗})∥_{L^{2}(0,1)} → 0 for η → 0 also implies ∥x_{η} − x^{∗}∥_{L^{p}(0,1)} → 0, since Tikhonov’s lemma (see Lemma 2.4) in fact only needs the local injectivity condition F (x) = F (x^{∗}) (x ∈ D) ⟹ x = x^{∗}.

On the other hand, let y ∈ R_{ε}^{+} for ε > 0, i.e. y(s) = 0 if s ∈ [0, ε]. As shown in [8], in such a case the autoconvolution operator F is non-injective and it holds:

x^{∗}(t) = \{ \begin{matrix}0 & a.e. in [0, \frac{ε}{2}] \\ uniquely determined & a.e. in [\frac{ε}{2}, 1 − \frac{ε}{2}] \\ arbitrarily non-negative & in [1 − \frac{ε}{2}, 1] \end{matrix} .

Consequently, we have x^{∗} ∈ D+_{\frac{ε}{2}}. Since the values x^{∗}(t) do not depend on y for t ∈ [1 − \frac{ε}{2}, 1], we cannot expect any information about the solution in this subinterval from the data. Therefore, it makes sense to solve the equation (1.1) only on the interval [\frac{ε}{2}, 1 − \frac{ε}{2}]. We will show that this case is reducible to the already treated case y ∈ R_{0}^{+}. Because of this we define the operator F_{ε} : L^{p}(\frac{ε}{2}, 1 − \frac{ε}{2}) → L^{2}(ε, 1) as

[F_{ε}(x)](s) := ∫ s− \frac{ε}{2}_{\frac{ε}{2}} x(s − t)x(t)dt.

Then we have [F (x)](s) = [F_{ε}(x)](s) for \frac{ε}{2} ≤ s ≤ 1 − \frac{ε}{2}. By using the transformations

t̃:= \frac{t − \frac{ε}{2}}{1 − ε}, s̃:= \frac{s − ε}{1 − ε},

and

x̃(t̃) := x((1 − ε)t̃+ \frac{ε}{2}) = x(t), ỹ(s̃) := y((1 − ε)s̃ + ε) = y(s), F̃_{ε}(x̃) := F_{ε}(x),

we obtain an operator F̃_{ε} : L^{p}(0, 1) → L^{2}(0, 1) defined by

[F̃_{ε}(x̃)](s̃) := (1 − ε) ∫ s̃ ; 0 x̃(s̃ − t̃)x̃(t̃)dt̃.

Then we get x̃ ∈ L^{p}(0, 1) if x ∈ L^{p}(0, 1), and instead of (1.2) we have to solve the equation F̃_{ε}(x̃) = ỹ now. From y ∈ R_{ε}^{+} and x ∈ R+_{\frac{ε}{2}} it follows that ỹ ∈ R_{0}^{+} and x̃ ∈ D_{0}^{+}, respectively. Hence we have F̃_{ε}(x̃) = (1 − ε)F (x̃) for all x̃ ∈ D_{0}^{+}. Therefore, we can proceed as in the injective case and compute converging approximate solutions x̃_{η}. Then we transform back to the interval [\frac{ε}{2}, 1 − \frac{ε}{2}] and obtain approximate solutions with satisfactory properties on this interval, where the performed linear transformation retains the monotonicity. Finally we extend the solution by zero on the interval [0, \frac{ε}{2}). On the other remaining subinterval (1− \frac{ε}{2}, 1] the solution can be extended arbitrarily provided that the monotonicity requirement is satisfied. Unfortunately, the value of ε is unknown if only discrete noisy data are given. In some situations, however, this value can be estimated and the transformation procedure becomes applicable.

5 Numerical Examples

In the concluding section we present some case studies on the behaviour of approximate discrete least-squares solutions to the autoconvolution equation (1.1) from noisy data, where we follow the approach of Section 2.

The first study is devoted to the case of a continuous, but non-monotone exact solution. We use the example

x^{∗}(t) = −3t^{2} + 3t + \frac{1}{4} = 1 − 3 (t − \frac{1}{2})^{2}  (0 ≤ t ≤ 1)  (5.1)

with the right-hand side

y^{∗}(s) = \frac{3}{10}s^{5} − \frac{3}{2}s^{4} + s^{3} + \frac{3}{4}s^{2} + \frac{1}{16}s (0 ≤ s ≤ 1)  (5.2)

and obtain a := 0.25 ≤ x^{∗}(t) ≤ b := 1, T (x^{∗}) = 1.5 and x^{∗} ∈ D with D from (2.5) - (2.6). The noisy data ˆ \underline{y} were generated by adding normally distributed pseudorandom numbers with zero mean and standard deviation σy_{i} (σ fixed) to the discrete values y_{i} of (5.2). We used varying values c as upper bounds for the total variation of the discretized solutions.

The nonlinear optimization problem (2.20) was numerically solved by a Gauss-Newton code. In the case of unacceptable Gauss-Newton steps this code uses the Marquardt method. The theory of this procedure is due to [3, pp. 348-368] (for the algorithm see [3, pp. 369-383]). We used penalty terms to handle the constraints of D. In all figures presented below the solid lines give the exact solution x^{∗} according to (5.1), whereas the lines with small circles express the approximate solutions x_{η} such that every circle corresponds to a grid point of discretization.

In the Figures 1 and 2 we compare approximate solutions x_{η} in the case of unperturbed data (σ = 0) using n = 50 grid points and different bounds c for the total variation. For an appropriate choice c = 1.5 associated with the really arising total variation level, the approximate solution is very good in the noiseless case (see Figure 1), whereas an underestimated value c = 0.8 < T (x^{∗}) corresponds to an overregularized solution (see Figure 2), which is much too ’flat’ compared to the function x^{∗} to be determined.

Now we turn to the case of noisy data. For all computations in the context of the Figures 3 – 6 a per mille noise level σ = 10^{−3} was used. We begin with a situation (see Figure 3), where the total variation bound was omitted (c = ∞). Then the set M of admissible discrete solutions contains strongly oscillating vectors. Especially for t from the right halfinterval of [0, 1] the quality of the approximate solutions may be very bad in that situation.

The Figure 4 illustrates in a rather convincing manner the utility of the total variation approach presented above in Section 2 for handling noisy data. In particular, the approximation quality of x_{η} in Figure 4 with c = 1.5 at the right end of the interval is much better than in Figure 3. We can motivate this right-end effect as follows: By the autoconvolution of a function x(t) (0 ≤ t ≤ 1) the values x(t) for small t influence the function values y(s) in some sense more than the values x(t) with t close to 1. Namely, x(t) only influences y(s) for s > t. As a consequence, the reconstruction of x(t) from y is more stable for smaller t, since then the function y(s) = [x∗x](s) (0 ≤ s ≤ 1) has collected more information about the value x(t) to be determined. In the case of overregularization (c is selected too small compared to T (x^{∗})), this phenomenon may cause large reconstruction errors specifically at the right end of the interval [0, 1] (see Figure 5 with c = 1.0).

We should mention that the analysis of the problem based on Lemma 2.4 does not provide any rate of convergence for the solution error ∥x_{η} − x^{∗}∥_{L^{p}(0,1)} depending on ∥F (x_{η}) − F (x^{∗})∥_{L^{2}(0,1)}. On the other hand, Lemma 2.7 shows that the order of magnitude for the discrepancy norm ∥F (x_{η})−F (x^{∗})∥_{L^{2}(0,1)} corresponds with the maximum max(hb̄^{2}, hb̄c̄, δ), where h := 1/n, b̄ := sup_{{x∈D \atop t∈[0,1]}} x(t) and c̄:= sup_{{x∈D \atop t∈[0,1]}} T (x). For sufficiently large n this discrepancy norm is dominated by the noise level δ, or in our case study by the value σ. So we can see comparing the Figures 4 and 6 that different discretization levels n = 50 and n = 25 yield approximate solutions with nearly the same accuracy provided that the noise level (σ = 10^{−3}) does not change.

To a second study we have been motivated by numerical experiments carried out by Kutsche in her thesis [11]. There it was shown that the constraint of bounded variation is very useful in the L^{1}-approximation of piecewise continuous solutions to Abel integral equations. We now will demonstrate the effects of using the least-squares method under total variation constraints in the case of non-smooth functions possessing jumps. Therefore we consider as the exact solution the step function

x^{∗}(t) = \{ \begin{matrix}0.5 & if 0 ≤ t ≤ 0.5 \\ 0.25 & if 0.5 < t ≤ 0.8 \\ 0.75 & if 0.8 < t ≤ 1 \end{matrix}  (5.3)

with the right-hand side

y^{∗}(s) = \{ \begin{matrix}0.25t & if 0 ≤ t ≤ 0.5 \\ 0.125 & if 0.5 < t ≤ 0.8 \\ 0.5t − 0.275 & if 0.8 < t ≤ 1 \end{matrix} .  (5.4)

The exact solution x^{∗} is discontinuous, non-monotone but a function of bounded variation. Its total variation can easily computed as T (x^{∗}) = 0.75. The function x^{∗} is bounded, positive and left-continuous on the whole interval [0, 1]. Therefore the requirements of the set D from (2.5) - (2.6) are fulfilled.

We will now compare the approximate solutions of this example for different choices of the parameter c. Let the number of discretization points n = 50 and the value σ = 10^{−2} of noise be constant throughout this study. Then we are able to control the solution by changing the parameter c.

In the Figures 7 – 10 the graphs of both the numerical solution x_{η} and the exact solution x^{∗} (bold line) are drawn as piecewise constant functions. In our first example (Figure 7) we computed the solution without any total variation restriction (c = ∞). The solution is – as in the first example – rather bad and highly oscillating. However, it is to mention that the jumps of x^{∗} are reconstructed relatively good in this case. In Figure 8 the situation c = T (x^{∗}) = 0.75 is illustrated. Here the solution is much smoother than in the unconstrained case, but the points with maximal approximation errors are now the jumps at t = 0.5 and t = 0.8. In these points the approximate solution is ’oversmoothed’. This depends on the fact that the smoothing effect of regularization acts uniformly on the whole interval [0, 1], but the character of a jump function does not correspond with this property. Therefore the jumps are blurred by that choice of c. Moreover, the ’right-end effect’ discussed above is superposed and leads to growing errors near t = 1.

Finally, in Figure 9 with c = 0.5 and Figure 10 with c = 2 we demostrate two more situations. If c underestimates the value T (x^{∗}), then the effect of blurring the jumps is still more pronounced. On the other hand, the admissible oscillation level grows if c overestimates T (x^{∗}). In that case, however, the location of jumps can be determined rather precise. That means, if one supposes that the exact solution is a step function, then it is recommended to choose c not too small. This allows some oscillations around the exact solution whose amplitudes are small if c is not chosen much too large.

References

[1] Acar, R.; Vogel, C.R.: Analysis of total variation penalty methods. Inverse Problems 10 (1994), 1217-1229.

[2] Baumeister, J.: Stable Solution of Inverse Problems. Vieweg, Braunschweig 1987.

[3] Dennis, J.E.; Gay, D.M.; Welsch, R.E.: An adaptive nonlinear least–squares algorithm. ACM Trans. Math. Softw. 7 (1981), 348-368.

[4] Engl, H.W.; Hofmann, B.; Zeisel, H.: A decreasing rearrangement approach for a class of ill-posed nonlinear integral equations. Journal of Integral Equations and Applications 5 (1993), 443-463.

[5] Fleischer, G; Hofmann, B.: On inversion rates for the autoconvolution equation. Inverse Problems 12 (1996), 419-435.

[6] Gellrich, C.; Hofmann, B.: A study of regularization by monotonicity. Computing 50 (1993), 105-125.

[7] Gorenflo, R.: Computation of rough solutions of Abel integral equations. In: Inverse and Ill-Posed Problems (Eds.: H.W. Engl, C.W. Groetsch), Academic Press, Boston, 1987, 195-210.

[8] Gorenflo, R.; Hofmann, B.: On autoconvolution and regularization. Inverse Problems 10 (1994), 353-373.

[9] Hofmann, B.: Regularization for Applied Inverse and Ill-Posed Problems. B.G. Teubner, Leipzig 1986.

[10] Hofmann, B.; Tautenhahn, U.: On ill-posedness measures and space change in Sobolev scales. Paper submitted.

[11] Kutsche, C.: Produktquadraturverfahren für nicht-glatte Lösungen Abelscher Integralgleichungen erster und zweiter Art. Dissertation, Freie Universität Berlin, FB Mathematik und Informatik, Berlin 1994.

[12] Natanson, I.P.: Theorie der Funktionen einer reellen Veränderlichen. Akademie-Verlag, Berlin, 1969.

[13] Tikhonov, A.N; Goncharsky, A.V.; Stepanov, V.V.; Yagola, A.G.: Regularizing algorithms and a priori information (in Russian). Nauka, Moscow, 1983.

[figure: 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1; 1.1; 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1]

Figure 1: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 0, n = 50, c = 1.5

[figure: 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1; 1.1; 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1]

Figure 2: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 0, n = 50, c = 0.8 Figure 3: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−3}, n = 50, c = ∞

[figure: 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1; 1.1; 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1]

Figure 4: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−3}, n = 50, c = 1.5 Figure 5: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−3}, n = 50, c = 1.0,

(inappropriate initial values used in Gauss-Newton-method)

[figure: 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1; 1.1; 0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 0  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9  1

Figure 6: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−3}, n = 25, c = 1.5 Figure 7: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−2}, n = 50, c = ∞

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 0  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9  1

Figure 8: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−2}, n = 50, c = 0.75 Figure 9: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−2}, n = 50, c = 0.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 0  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9  1

Figure 10: Least-squares solution of \underline{F} (\underline{x}) = ˆ \underline{y}, σ = 10^{−2}, n = 50, c = 2.0