Aurélien Alfonsi
email: alfonsi@cermics.enpc.fr

Alexander Schied
email: schied@uni-mannheim.de

Capacitary measures for completely monotone kernels via singular control

Keywords: Singular control, verification argument, capacity theory, optimal order execution, transient price impact, infinite-dimensional Riccati differential equation

We give a singular control approach to the problem of minimizing an energy functional for measures with given total mass on a compact real interval, when energy is defined in terms of a completely monotone kernel. This problem occurs both in potential theory and when looking for optimal financial order execution strategies under transient price impact. In our setup, measures or order execution strategies are interpreted as singular controls, and the capacitary measure is the unique optimal control. The minimal energy, or equivalently the capacity of the underlying interval, is characterized by means of a nonstandard infinite-dimensional Riccati differential equation, which is analyzed in some detail. We then show that the capacitary measure has two Dirac components at the endpoints of the interval and a continuous Lebesgue density in between. This density can be obtained as the solution of a certain Volterra integral equation of the second kind.

1 Introduction and statement of results

Background

Let G : R + → R + be a function. The problem of minimizing the energy functional

E(µ) := 1 2 G(|t -s|) µ(ds) µ(dt)
over probability measures µ supported by a given compact set K ⊂ R plays an important role in potential theory. A minimizing measure µ * , when it exists, is called a capacitary measure, and the value Cap (K) := 1/E(µ * ) is called the capacity of the set K; see, e.g., [START_REF] Choquet | Theory of capacities[END_REF], [START_REF] Fuglede | On the theory of potentials in locally compact spaces[END_REF][START_REF] Landkof | Foundations of Modern Potential Theory[END_REF].

In this paper, we develop a control approach to determining the capacitary distribution µ * when K is a compact interval and G is a completely monotone function. In this approach, measures µ on K will be regarded as singular controls and E(µ) is the objective function. Our goal is to obtain qualitative structure theorems for the optimal control µ * and characterize µ * by means of certain differential and integral equations.

The intuition for this control approach, and in fact our original motivation, come from the problem of optimal order execution in mathematical finance. In this problem, one considers an economic agent who wishes to liquidate a certain asset position of x shares within the time interval [0, T ]. This asset position can either be a long position (x > 0) or a short position (x < 0). The order execution strategy chosen by the investor is described by the asset position X t held at time t ∈ [0, T ]. In particular, one must have X 0 = x. Requiring the condition X T + = 0 assures that the initial position has been unwound by time T . The leftcontinuous path X = (X t ) t∈[0,T ] will be nonincreasing for a pure sell strategy and nondecreasing for a pure buy strategy. A general strategy can consist of both buy and sell trades and hence can be described as the sum of a nonincreasing and a nondecreasing strategy. That is, X is a path of finite variation.

The problem the economic agent is facing is that his or her trades impact the price of the underlying asset. To model price impact, one starts by informally defining q dX t as the immediate price impact generated by the (possibly infinitesimal) trade dX t executed at time t. Next, it is an empirically well-established fact that price impact is transient and decays over time; see, e.g., [START_REF] Moro | Market impact and trading profile of hidden orders in stock markets[END_REF]. This decay of price impact can be described informally by requiring that G(t -s) dX s is the remaining impact at time t of the impact generated by the trade dX s . Here, G : R + → R + is a nonincreasing function with G(0) = q, the decay kernel. Thus, s<t G(t -s) dX s is the price impact of the strategy X, cumulated until time t. This price impact creates liquidation costs for the economic agent, and one can derive that, under the common martingale assumption for unaffected asset prices, these costs are given by C(X) := 1 2 [0,T ] [0,T ] G(|t -s|) dX s dX t (1) plus a stochastic error term with expectation independent of the specific strategy X; see [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF]. Thus, minimizing the expected costs amounts to minimizing the functional C(X) over all leftcontinuous strategies X that are of bounded variation and satisfy X 0 = x and X t+ = 0. This problem was formulated and solved in the special case of exponential decay, G(t) = e -ρt , by [START_REF] Obizhaeva | Optimal trading strategy and supply/demand dynamics[END_REF]. The general case was analyzed by [START_REF] Alfonsi | Order book resilience, price manipulation, and the positive portfolio problem[END_REF] in discrete time and by [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF] in the continuous-time setup we have used above. We refer to [START_REF] Alfonsi | Optimal execution strategies in limit order books with general shape functions[END_REF], Alfonsi & Schied (2010), [START_REF] Gatheral | Exponential resilience and decay of market impact[END_REF], [START_REF] Predoiu | Optimal execution in a general one-sided limit-order book[END_REF], Schied & Slynko (2011), and Gatheral & Schied (2011) for further discussions and additional references in the context of mathematical finance.

Clearly, the cost functional C(X) coincides with the energy functional E(ν X ) of the measure ν X (dt) := dX t . So finding an optimal order execution strategy is basically equivalent to determining a capacitary measure for [0, T ]. There is one important difference, however: capacitary measures are determined as minimizers of E(µ) with respect to all nonnegative measures µ on [0, T ] with total mass 1, while ν X may be a signed measure with given total mass ν X ([0, T ]) = -x. This difference can become significant if G(| • |) is only required to be positive definite in the sense of Bochner (which is equivalent to C(X) ≥ 0 for all X), because then minimizers of the unconstraint problem need not exist. It was first shown by [START_REF] Alfonsi | Order book resilience, price manipulation, and the positive portfolio problem[END_REF], and later extended to continuous time by [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF], that a unique optimal order execution strategy X * exists and that X * is a monotone function of t when G is convex and nonincreasing. This result has the important consequence that the constrained problem of finding a capacitary measure is equivalent to the unconstrained problem of determining an optimal order execution strategy.

In this paper, we aim at describing the structure of capacitary measures/optimal order execution strategies. To this end, it is instructive to first look at two specific examples in which the optimizer is known in explicit form. [START_REF] Obizhaeva | Optimal trading strategy and supply/demand dynamics[END_REF] find that for exponential decay, G(t) = e -ρt , the capacitary measure µ * has two singular components at t = 0 and t = T and a constant Lebesgue density on (0, T ):

µ * (dt) = 1 2 + ρT δ 0 (dt) + ρ 2 + ρT dt + 1 2 + ρT δ T (dt). ( 2
)
Numerical experiments show that it is a common pattern that capacitary measures for nonincreasing convex kernels have two singular components at t = 0 and T = 0 and a Lebesgue density on (0, T ). However, the capacitary measure for G(t) = max{0, 1 -ρt} is the purely discrete measure

µ * = 1 2 + N N i=0 1 - i N + 1 δ i ρ + δ T -i ρ ,
where N := ⌊ρT ⌋ (Gatheral et al. 2012, Proposition 2.14). So it is an interesting question for which nonincreasing, convex kernels G the capacitary measure µ * has singular components only at t = 0 and t = T and is (absolutely) continuous on (0, T ). It turns out that a sufficient condition is the complete monotonicity of G, i.e., G belongs to C ∞ ((0, ∞)) and (-1) n G (n) is nonnegative in (0, ∞) for n ∈ N. More precisely, we have the following result, which is in fact an immediate corollary to the main results in this paper.

Corollary 1. Suppose that G is completely monotone with G ′′ (0+) := lim t↓0 G ′′ (t) < ∞. Then the capacitary measure µ * has two Dirac components at t = 0 and t = T and is has a continuous Lebesgue density on (0, T ).

Statement of main results

Our main results do not only give the preceding qualitative statement on the form of µ * but they also provide quantitative descriptions of the Dirac components of µ * and of its Lebesgue density on (0, T ). To prepare for the statement of these results, let us first assume that G(0) = 1, which we can do without loss of generality. Then we recall that by the celebrated Hausdorff-Bernstein-Widder theorem (Widder 1941, Theorem IV.12a), G is completely monotone if and only it is the Laplace transform of a Borel probability measure λ on R + :

G(t) = e -ρt λ(dρ), t ≥ 0.
In particular, every exponential polynomial,

G(t) = d i=0 λ i e -ρ i t , (3)
with λ i , ρ i ≥ 0 and i λ i = 1 is completely monotone. Another example is power-law decay,

G(t) = 1 (1 + t) γ
for some γ > 0, which is a popular choice for the decay of price impact in the econophysics literature. We assume henceforth that G ′′ (0+) < ∞, which is equivalent to

ρ := ρ λ(dρ) < ∞ and ρ 2 := ρ 2 λ(dρ) < ∞. (4)
A crucial role will be played by the following infinite-dimensional Riccati equation for func-

tions ϕ : [0, ∞) × R 2 + → R, ϕ ′ (t, ρ 1 , ρ 2 ) + (ρ 1 + ρ 2 )ϕ(t, ρ 1 , ρ 2 ) = 1 2ρ ρ 1 + xϕ(t, ρ 1 , x) λ(dx) ρ 2 + xϕ(t, x, ρ 2 ) λ(dx)
(5) where ϕ ′ denotes the time derivative of ϕ, and the function ϕ satisfies the initial condition

ϕ(0, ρ 1 , ρ 2 ) = 1 for all ρ 1 , ρ 2 ≥ 0. ( 6
)
Remark 1. When writing (5) in the form ϕ ′ = F (ϕ) one sees that the functional F is not a continuous map from some reasonable function space into itself, unless λ is concentrated on a compact interval. For instance, it involves the typically unbounded linear operator ϕ → (ρ 1 + ρ 2 )ϕ. Therefore, existence and uniqueness of solutions to (5), ( 6) does not follow by an immediate application of standard results such as the Cauchy-Lipschitz/Picard-Lindelöf theorem in Banach spaces (Hille & Phillips 1957, Theorem 3.4.1) or more recent ones such as those in [START_REF] Teixeira | Strong solutions for differential equations in abstract spaces[END_REF] and the references therein. In fact, even in the simplest case in which λ reduces to a Dirac measure, the existence of global solution hinges on the initial condition; it is easy to see that solutions blow up when ϕ(0) is not chosen in a suitable manner.

We now state a result on the global existence and uniqueness of ( 5), ( 6). It states that the solution takes values in the locally convex space C(R 2 + ) endowed with topology of locally uniform convergence. For integers k ≥ 0, the space C k ([0, ∞); C(R 2 + )) will consist of all continuous functions ϕ : [0, ∞) → C(R 2 + ) which, when considered as functions ϕ : [0, ∞) → C(K) for some compact subset K of R 2 + , belong to C k ([0, ∞); C(K)).

Theorem 1. When G ′′ (0+) < ∞ the initial value problem (5), (6) admits a unique solution ϕ in the class of functions ϕ in C 1 ([0, ∞); C(R 2 + )) that satisfy an inequality of the form

0 ≤ ϕ(t, ρ 1 , ρ 2 ) ≤ c(1 + ρ 1 )(1 + ρ 2 ), ( 7
)
where c is a constant that may depend on ϕ and locally uniformly on t. Moreover, ϕ has the following properties.

(a) ϕ is strictly positive.

(b) ϕ is symmetric: ϕ(t, ρ 1 , ρ 2 ) = ϕ(t, ρ 2 , ρ 1 ) for all (ρ 1 , ρ 2 ) ∈ R 2 + . (c) 1 = ϕ(t, ρ, x) λ(dx) = ϕ(t, x, ρ) λ(dx) for all ρ ≥ 0. (d) ϕ ∈ C 2 ([0, ∞); C(R 2 + )). (e) For every t, the kernel ϕ(t, •, •) is nonnegative definite on L 2 (λ), i.e., f (x)f (y)ϕ(t, x, y) λ(dx) λ(dy) ≥ 0 for f ∈ L 2 (λ). (8)
(f) The functions ϕ(t, ρ 1 , ρ 2 ) and ϕ ′ (t, ρ 1 , ρ 2 ) satisfy local Lipschitz conditions in (ρ 1 , ρ 2 ), locally uniformly in t.

We can now explain how to use singular control in approaching the minimization of E(µ) or C(X). To this end, using order execution strategies X = (X t ) will be more convenient than using the formalism of the associated measures µ(dt) = dX t because of the natural dynamic interpretation of X. Henceforth, a [0, T ]-admissible strategy will be a leftcontinuous function (X t ) of bounded variation such that X T + = 0. Our goal is to minimize the cost functional C(X) defined in (1) over all [0, T ]-admissible strategies with fixed initial value X 0 = x. Clearly, this problem is not yet suitable for the application of control techniques since C(X) depends on the entire path of X. We therefore introduce the auxiliary functions

E X t (ρ) := [0,t) e -ρ(t-s) dX s , for ρ ≥ 0. ( 9
)
These functions will play the role of state variables that are controlled by the strategy X.

Lemma 1. For any [0, T ]-admissible strategy X, the function E X t (ρ) is uniformly bounded in ρ and t. Moreover,

C(X) = [0,T ) E X t (ρ) λ(dρ) dX t + 1 2 t≤T (∆X t ) 2 , ( 10
)
where ∆X t := X t+ -X t denotes the jump of X at t.

Proof. Clearly, |E X t (ρ)| ≤ X var , where X var denotes the total variation of X over [0, T ]. To obtain (10), we integrate by parts to get

C(X) = [0,T ) [0,t) G(t -s) dX s dX t + G(0) 2 t≤T (∆X t ) 2 .
Now we write G(t -s) as e -ρ(t-s) λ(dρ) and apply Fubini's theorem. .

The form (10) of our cost functional is now suitable for the application of control techniques. To state our main result, we let ϕ be the solution of our infinite-dimensional Riccati equation as provided by Theorem 1 and we define

ϕ 0 (t) := ϕ(t, 0, 0) and ψ(t, ρ) := xϕ(t, x, ρ) λ(dx) (11)
Theorem 2. Let X * be the unique optimal strategy in the class of [0, T ]-admissible strategies with initial value

X 0 = x. Then C(X * ) = x 2 2ϕ 0 (T ) . ( 12
)
Moreover, X * has jumps at t = 0 and t = T of size

∆X * 0 = ∆X * T = - ψ(T, 0) 2ρϕ 0 (T ) x
and is continuously differentiable on (0, T ). The derivative x(t) = d dt X * t is the unique continuous solution of the Volterra integral equation

x(t) = f (t) + t 0 K(t, s)x(s) ds, (13)
where, for

Θ(t, ρ) := ρ + ψ(t, ρ) ψ(t, 0) x 2 ϕ(t, x, 0) λ(dx) -x 2 ϕ(t, x, ρ) λ(dx) + ρ 2 , ( 14
)
the function f and the kernel K(•, •) are given by

f (t) = ∆X * 0 2ρ e -ρt Θ(T -t, ρ) λ(dρ), K(t, s) = 1 2ρ e -ρ(t-s) Θ(T -t, ρ) λ(dρ). (15)
Let us recall that we know in addition from Theorem 2.20 in [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF]

) that t ∈ [0, T ] → X *
t is monotone. The identity (12) immediately yields the following formula for the capacity of a compact interval.

Corollary 2. If G ′′ (0+) < ∞, the capacity of a compact interval [a, b] is given by Cap ([a, b]) = 2ϕ 0 (b -a).

Computational aspects

In general, the Riccati equation ( 5), ( 6) cannot be solved explicitly. A closed-form solution exists, however, when G is an exponential polynomial as in (3), i.e., when λ has a discrete support. Let us assume that λ(dx)

= d i=0 λ i δ ρ i (dx), with ρ 0 = 0 < ρ 1 < • • • < ρ d , λ i ≥ 0, and d i=0 λ i = 1. All the input that is needed in Theorem 2 are the values ϕ ij (t) := ϕ(t, ρ i , ρ j ), for 0 ≤ i, j ≤ d.
By Theorem 1, ϕ(t) is a symmetric matrix that solves the following matrix Riccati equation:

ϕ ′ = -ϕM (3) ϕ -ϕM (4) + M (1) ϕ + M (2) , with M (3) ij = -1 2ρ λ i ρ i λ j ρ j , M (4) ij = - λ i ρ i ρ j 2ρ + δ ij ρ i , M (1) = -(M (4) ) ⊤ and M (2) ij = ρ i ρ j
2ρ . According to [START_REF] Levin | On the matrix Riccati equation[END_REF], the solution of this equation is given by

ϕ(t) = (R (1) (t)1 + R (2) (t))(R (3) (t)1 + R (4) (t)) -1 ,
where 1 ij = 1 and

R(t) = R (1) (t) R (2) (t) R (3) (t) R (4) (t) = exp t M (1) M (2) M (3) M (4) .
In the special cases d = 1 and d = 2, the solution of the Riccati equation ( 5), ( 6) becomes even easier and, to some extend, becomes explicit. We demonstrate this here for d = 1:

Example 1. In the case d = 1, G is of the form G(t) = λ + (1 -λ)e -ρt for some λ ∈ [0, 1)
and some ρ > 0. Clearly, we can set λ := 0 without changing the optimization problem. Then ρ = ρ 1 = ρ, and (5) becomes

ϕ ′ 00 = ρ 2 ϕ 2 01 , ϕ ′ 01 + ρϕ 01 = ρ 2 (1 + ϕ 11 )ϕ 01 , ϕ ′ 11 + 2ρϕ 11 = ρ 2 (1 + ϕ 11 ) 2 .
For the initial condition ϕ kl (0) = 1, it has the unique solution ϕ 11 ≡ ϕ 01 ≡ 1 and ϕ 00 (t) = 1 + ρt/2. The condition (53) thus reduces to 0 = X t + E 1 (t) 1 + ρ(T -t) , which easily yields (2) as unique solution.

Given the solution ϕ of the Riccati equation, we can approximate the continuous time strategy by a discrete one as follows (x i will denote the trading size at time iT /N ).

• We first set x 0 = ψ(T,0) 2ρϕ 0 (T ) and E 0 (ρ ℓ ) = x 0 , 0 ≤ ℓ ≤ d. • Suppose that 1 ≤ i < N and that x i-1 and E i-1 (ρ ℓ ) have been computed. Then, we set thanks to (53):

x i = 1 - i-1 j=0 x j -E i-1 (ρ)e -ρT /N θ(T -iT /N, ρ)λ(dρ), E i (ρ ℓ ) = E i-1 (ρ ℓ )e -ρ ℓ T /N + x i .
• Set x N = 1 -i-1 j=0 x j . Alternatively, we could have approximated the minimization of the cost (1) by the following discrete problem. Let M i,j = G |i -j| T N , 0 ≤ i, j ≤ N , and consider

minimize 1 2 x T M x over x ∈ R N +1 s.t N i=0 x i = 1. ( 16
)
The solution of this problem is obviously given by

1 1 T M -1 1 M -1 1, where 1 i = 1 for 0 ≤ i ≤ N .
The study of the convergence when N → +∞ of this discrete solution to the continuous one given by Theorem 2 is beyond the scope of the paper, but has been already proved in some cases (e.g. Alfonsi et al. (2010) and [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF]) and is clear heuristically.

We have computed and drawn the solutions given by both methods in Figure 1 for T = 1, N = 50, and λ(dρ) = 0.1δ 0 (dρ) + 0.2δ 1 (dρ) + 0.2δ 3 (dρ) + 0.2δ 5 (dρ) + 0.2δ 7 (dρ) + 0.1δ 10 (dρ). They are already rather close together for N = 50, and they merge when N → +∞. Let us discuss briefly the time complexity of the two methods. The one given by ( 16) gets very slow when N gets large since it involves the inversion of a N × N matrix. Instead, when λ has a discrete support, the matrix Riccati equation can be solved quickly and the algorithm above has a O(N ) time complexity, which is much faster. However, this is no longer true when λ does not have discrete support. In that case, we have to approximate λ by a discrete measure, which means that we have to increase d. Doing so, will slow down the algorithm based on the Riccati equation. A rigorous treatment of the convergence rate and time complexity of both algorithms is beyond the scope of this paper and is left for future research.

Proofs

Proof of Theorem 1

Let us write (5) in the form ϕ ′ (t) = F λ (ϕ(t)), where

F λ (f )(ρ 1 , ρ 2 ) = -(ρ 1 +ρ 2 )f (ρ 1 , ρ 2 )+ 1 2ρ ρ 1 + xf (ρ 1 , x) λ(dx) ρ 2 + xf (x, ρ 2 ) λ(dx) . (17) Lemma 2. Suppose that λ is supported by the compact interval [0, ρ max ]. Then (5), (6) admits a unique solution ϕ ∈ C 1 ([0, ∞); C(R 2 + )).
Moreover, ϕ has the properties (a), (b), and (c) in the statement of Theorem 1.

Proof. Let J ⊂ R + be any compact interval containing [0, ρ max ]. Then F λ defined in (17) maps C(J × J) into itself. Moreover, F λ is Lipschitz continuous with respect to the sup-norm on every bounded subset of C(J × J). Hence, the Cauchy-Lipschitz/Picard-Lindelöf theorem in Banach spaces implies the existence of a unique local solution ϕ J ∈ C 1 ([0, t J ); C(J × J)) for some maximal time t J > 0 (Hille & Phillips 1957, Theorem 3.4.1). We will show below that Figure 2: Comparison of the approximated optimal strategies given by the Riccati method and N = 50, with the optimal continuous one X * t (computed with N = 1000) .

t J = ∞. Then, if J ⊃ J is another compact interval, the restriction of ϕ J (t) to J must coincide with ϕ J (t) due to the uniqueness of solutions. This consistency then implies the existence and uniqueness of solutions ϕ ∈ C 1 ([0, ∞); C(R 2 + ))). Moreover, the uniqueness of solutions and the fact that both (5) and ( 6) are symmetric in ρ 1 and ρ 2 implies that ϕ(t, ρ 1 , ρ 2 ) = ϕ(t, ρ 2 , ρ 1 ) for all (ρ 1 , ρ 2 ), which is property (b) in Theorem 1.

We now fix an interval J ⊃ [0, ρ max ]. Before proving that t J = ∞, we will show that

ϕ J (t, ρ, x) λ(dx) = 1 for ρ ∈ J and t < t J . ( 18
)
This will then establishes property (c) in the statement of Theorem 1 for t ∈ [0, t I ). Then we will use ( 18) to derive some estimates on ϕ J that will yield ϕ J > 0 and t J = ∞.

To prove (18), we let I(t, ρ) := ϕ J (t, ρ, x) λ(dx) and ψ J (t, ρ) := xϕ J (t, ρ, x) λ(dx). We have

I ′ (t, ρ) + ρI(t, ρ) + ψ J (t, ρ) = 1 2ρ ρ + ψ J (t, ρ) ρ + xI(t, x) λ(dx) . ( 19
)
This is a (non-homogeneous) affine ODE of the form

I ′ (t) = b(t) + A(t)I(t)
, where the operator

(A(t)f )(ρ) = -ρf (ρ) + 1 2ρ ρ + ψ J (t, ρ) xf (x) λ(dx)
is a continuous map from [0, δ] into the space of bounded linear operators on C(J) for each δ < t J . Hence this ODE admits a unique solution in C 1 ([0, δ]; C(J)) with initial condition I(0, ρ) = 1. But ( 19) is solved by I(t, ρ) = 1, which which establishes (18).

For the next step, we let

t 0 := inf t ∈ [0, t J ) min ρ 1 ,ρ 2 ∈J ϕ J (t, ρ 1 , ρ 2 ) < 0 .
Since ϕ J is a continuous map from [0, t J ) into C(J × J) and ϕ J (0) = 1, we must have t 0 > 0.

Due to (18) we have on [0, t 0 ) that

ρ 1 ρ 2 2ρ ≤ ϕ ′ J (t, ρ 1 , ρ 2 ) + (ρ 1 + ρ 2 )ϕ J (t, ρ 1 , ρ 2 ) ≤ (ρ 1 + ρ max )(ρ 2 + ρ max ) 2ρ . ( 20
)
When defining φJ (t, ρ 1 , ρ 2 ) := e t(ρ

1 +ρ 2 ) ϕ J (t, ρ 1 , ρ 2 ), (21)
the preceding inequality can be rewritten as

ρ 1 ρ 2 2ρ • e t(ρ 1 +ρ 2 ) ≤ φ′ J (t, ρ 1 , ρ 2 ) ≤ (ρ 1 + ρ max )(ρ 2 + ρ max ) 2ρ • e t(ρ 1 +ρ 2 ) .
Integrating these inequalities yields that for 0

≤ t < t 0 ϕ J (t, ρ 1 , ρ 2 ) ≥ e -t(ρ 1 +ρ 2 ) + ρ 1 ρ 2 (1 -e -t(ρ 1 +ρ 2 ) ) 2ρ(ρ 1 + ρ 2 ) > 0 ( 22
)
with the convention 1-e -(ρ 1 +ρ 2 )t

ρ 1 +ρ 2 = t for ρ 1 = ρ 2 = 0. Hence ϕ J (t, ρ 1 , ρ 2 ) ≤ e -t(ρ 1 +ρ 2 ) + (ρ 1 + ρ max )(ρ 2 + ρ max ) 2ρ(ρ 1 + ρ 2 ) (1 -e -t(ρ 1 +ρ 2 ) ). ( 23
)
Inequality ( 22) ensures that t 0 ≥ t J . Both inequalities ( 22) and ( 23) ensure the solution ϕ J (t) does not explode in finite time, which by standard arguments yields that t J = +∞. This proves the global existence of solutions as well as property (a) in the statement of Theorem 1.

The preceding lemma works only for measures λ that are concentrated on a finite interval. To obtain solutions for more general measures λ, we need to find upper bounds that are independent of ρ max . To this end, we first derive such bounds for the function ψ(t, ρ) defined in (11). By Lemma 2, this function is well-defined whenever λ has compact support, and it follows from dominated convergence together with ( 22) and ( 23

) that ψ ∈ C 1 ([0, ∞); C(R + )) and that ψ ′ (t, ρ) = xϕ ′ (t, ρ, x) λ(dx).
Lemma 3. Under the assumptions of Lemma 2, we have

0 < ψ(t, ρ) ≤ ρ 2 ρ for all ρ ≥ 0. ( 24
)
Proof. The lower bound in ( 24) is clear from ϕ > 0. To prove the upper bound, we suppose by way of contradiction that there exist t, ρ, and ε > 0 such that ψ(t, ρ) ≥ ε + ρ 2 /ρ. Then there must be a compact interval J ⊃ [0, ρ max ] such that

τ ε := inf t ≥ 0 max ρ∈J ψ(t, ρ) ≥ ρ 2 ρ + ε
is finite. Since ψ(0, ρ) = ρ and ρ 2 ≤ ρ 2 , the time τ ε must also be strictly positive. Moreover, there exists

ρ ε ∈ J such that max ρ∈J ψ(τ ε , ρ) = ψ(τ ε , ρ ε ) = ρ 2 ρ + ε.
Then τ ε is the first time at which the function t → ψ(t, ρ ε ) reaches a new maximum, and so ψ ′ (τ ε , ρ ε ) ≥ 0. Integrating (5) with respect to ρ 1 λ(dρ 1 ) and evaluating at ρ 2 = ρ ε gives

ψ ′ (τ ε , ρ ε ) + ρ ε ψ(τ ε , ρ ε ) + ρ 2 ϕ(τ ε , ρ, ρ ε ) λ(dρ) = 1 2ρ ρ ε + ψ(τ ε , ρ ε ) ρ 2 + ρψ(τ ε , ρ) λ(dρ) .( 25
)
Since ϕ(τ ε , ρ, ρ ε ) λ(dρ) = 1, the Cauchy-Schwarz inequality (or, alternatively, Jensen's inequality) implies that ρ 2 ϕ(τ ε , ρ, ρ ε ) λ(dρ) ≥ ψ(τ ε , ρ ε ) 2 . Moreover, the definition of ρ ε and the fact that λ is supported on J yield that ρψ(τ ε , ρ) λ(dρ) ≤ ρψ(τ ε , ρ ε ). Plugging these two inequalities into (25) leads to

ψ ′ (τ ε , ρ ε ) ≤ 1 2ρ ρ ε + ψ(τ ε , ρ ε ) ρ 2 + ρψ(τ ε , ρ ε ) -ρ ε ψ(τ ε , ρ ε ) -ψ(τ ε , ρ ε ) 2 = ρ ε ρ 2 2ρ + ρ 2 2ρ - ρ ε 2 ψ(τ ε , ρ ε ) - 1 2 ψ(τ ε , ρ ε ) 2 =: p(ψ(τ ε , ρ ε )),
where p(•) is a polynomial function of degree two. It has the two roots -ρ ε ≤ 0 and ρ 2 /ρ > 0. Therefore p(x) < 0 for x > ρ 2 /ρ and in turn 0

> p(ψ(τ ε , ρ ε )) = ψ ′ (τ ε , ρ ε ), which contradicts the fact that ψ ′ (τ ε , ρ ε ) ≥ 0.
Lemma 4. Under the assumptions of Lemma 2, we have

e -t(ρ 1 +ρ 2 ) + ρ 1 ρ 2 (1 -e -t(ρ 1 +ρ 2 ) ) 2ρ(ρ 1 + ρ 2 ) ≤ ϕ(t, ρ 1 , ρ 2 ) ≤ exp(-(ρ 1 + ρ 2 )t) + (ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ(ρ 1 + ρ 2 ) 1 -exp(-(ρ 1 + ρ 2 )t) , (26)
-(ρ 1 + ρ 2 ) - (ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ ≤ ϕ ′ (t, ρ 1 , ρ 2 ) ≤ (ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ . ( 27
)
Proof. The ODE (5) can be rewritten as

ϕ ′ (t, ρ 1 , ρ 2 ) + (ρ 1 + ρ 2 )ϕ(t, ρ 1 , ρ 2 ) = 1 2ρ ρ 1 + ψ(t, ρ 1 ) ρ 2 + ψ(t, ρ 2 ) . ( 28
)
Defining φ as in ( 21) and using the upper bound in (24) thus yields that

ρ 1 ρ 2 2ρ • e t(ρ 1 +ρ 2 ) ≤ φ′ (t, ρ 1 , ρ 2 ) ≤ (ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ • e t(ρ 1 +ρ 2 ) . ( 29
)
Arguing as in the final step of the proof Lemma 2 now yields (26). By plugging (26) back into (28) and using once again (24), we obtain (27).

Lemma 5. For all R, T > 0 there exist constants L 2 , L 2 ≥ 0 depending only on R, T , ρ, and ρ 2 such that for all t ∈ [0, T ] and

ρ 1 , ρ 1 ρ 2 ∈ [0, R], |ϕ(t, ρ 1 , ρ 2 ) -ϕ(t, ρ 1 , ρ 2 )| ≤ L 1 |ρ 1 -ρ 1 | (30) and |ϕ ′ (t, ρ 1 , ρ 2 ) -ϕ ′ (t, ρ 1 , ρ 2 )| ≤ L 2 |ρ 1 -ρ 1 |. (31)
Proof. We consider ρ 1 , ρ 1 , ρ 2 ≥ 0 and define

∆ρ 1 := ρ 1 -ρ 1 , ∆ϕ(t) = ϕ(t, ρ 1 , ρ 2 ) -ϕ(t, ρ 1 , ρ 2
) and ∆ψ(t) = ψ(t, ρ 1 ) -ψ(t, ρ 1 ).

By subtracting the equation ( 28) satisfied by ϕ(t, ρ 1 , ρ 2 ) from the corresponding one satisfied by ϕ(t, ρ 1 , ρ 2 ), we get

∆ϕ ′ (t) + ϕ(t, ρ 1 , ρ 2 )∆ρ 1 + (ρ 2 + ρ 1 )∆ϕ(t) = 1 2ρ ρ 2 + ψ(t, ρ 2 ) ∆ρ 1 + ∆ψ(t) . ( 32
)
This equation is a linear non-homogeneous ODE for ∆ϕ(t) and, since ∆ϕ(0) = 0, solved by

∆ϕ(t) = t 0 1 2ρ (ρ 2 + ψ(s, ρ 2 )) -ϕ(s, ρ 1 , ρ 2 ) ∆ρ 1 + 1 2ρ ρ 2 + ψ(s, ρ 2 ) ∆ψ(s) e -(ρ 1 +ρ 2 )(t-s) ds, Since |ψ(s, ρ 2 )| ≤ ρ 2
/ρ, we get with (26) and sup α≥0

1-e -αt α = t that |∆ϕ(t)| ≤ 1 2ρ ρ 2 + ρ 2 ρ t 0 |∆ρ 1 |+ |∆ψ(s)| ds + 1+ ( ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ T t 0 |∆ρ 1 | ds. (33)
Now, we have that

|∆ψ(t)| = (ϕ(t, ρ 1 , x) -ϕ(t, ρ 1 , x))x λ(dx) ≤ |ϕ(t, ρ 1 , x) -ϕ(t, ρ 1 , x)|x λ(dx) ≤ ρ 2 ρ t 0 |∆ρ 1 | + |∆ψ(s)| ds + ρ + ρ 1 + ρ 2 ρ ρ 2 ρ T t 0 |∆ρ 1 | ds.
For the last inequality, we have used Fubini's theorem and (33). Now, Gronwall's Lemma gives:

|∆ψ(t)| ≤ |∆ρ 1 |t ρ 2 ρ + ρ + ρ 1 + ρ 2 ρ ρ 2 ρ T exp ρ 2 ρ t . (34)
Plugging this back into (33), we get the existence of a constant L 1 , which depends only on R, T , ρ, and ρ 2 , such that

|∆ϕ(t)| ≤ L 1 |∆ρ 1 |. (35)
Finally, using ( 34) and ( 35) in ( 32) and recalling the locally uniform bounds ( 26) and ( 24) on ϕ and ψ gives (31).

Now drop the assumption that λ is supported on a compact interval and aim at proving existence and uniqueness of solutions in this general case. To this end, we take a sequence R n ↑ ∞ for which λ([0, R 1 ]) > 0 and define

f n := 1 λ([0, R n ]) ½ [0,Rn] and dλ n = f n dλ, (36)
so that each λ n satisfies the assumptions of Lemma 2. By ϕ n we denote the corresponding solution of ( 5), ( 6) provided by that lemma. For each n ≥ 1, we have

ρ n := ρ λ n (dρ) ≥ R 1 0 ρ λ(dρ) =: ρ 0 and ρ 2 n := ρ 2 λ n (dρ) ≤ ρ 2 λ([0, R 1 ]) =: ρ 2 0 . (37)
Hence, Lemma 4 yields that for each n,

e -t(ρ 1 +ρ 2 ) + ρ 1 ρ 2 (1 -e -t(ρ 1 +ρ 2 ) ) 2ρ 0 (ρ 1 + ρ 2 ) ≤ ϕ n (t, ρ 1 , ρ 2 ) ≤ 1 + (ρ 1 + ρ 2 0 ρ 0 )(ρ 2 + ρ 2 0 ρ 0 ) 2ρ 0 1 -e -(ρ 1 +ρ 2 )t ρ 1 + ρ 2 (38)
and

-(ρ 1 + ρ 2 ) - (ρ 1 + ρ 2 0 ρ 0 )(ρ 2 + ρ 2 0 ρ 0 ) 2ρ 0 ≤ ϕ ′ n (t, ρ 1 , ρ 2 ) ≤ (ρ 1 + ρ 2 0 ρ 0 )(ρ 2 + ρ 2 0 ρ 0 ) 2ρ 0 . (39)
Similarly, Lemma 5 yields that for all R, T > 0 there is a constant L ≥ 0 such that for all n

|ϕ n (t, ρ 1 , ρ 2 ) -ϕ n (t, ρ 1 , ρ 2 )| ≤ L|ρ 1 -ρ 1 | for all t ∈ [0, T ], ρ 1 , ρ 1 , ρ 2 ∈ [0, R], (40)
The inequalities (39), ( 40) and the Arzela-Ascoli theorem imply that the sequence (ϕ n ) n∈N is relatively compact in the class of continuous functions on [0, T ] × [0, R] 2 whenever T , R > 0, and hence admits a convergent subsequence in that class. By passing to a subsequence arising from a diagonalization argument if necessary, we may assume that there exists a continuous function ϕ : [0, ∞) × R 2 + → R + such that ϕ n → ϕ locally uniformly. The uniform bounds (38) and dominated convergence imply that

ψ n (t, ρ) := xϕ n (t, ρ, x) λ n (dx) = xϕ n (t, ρ, x)f n (x) λ(dx) (41) -→ xϕ(t, ρ, x) λ(dx) = ψ(t, ρ) (42)
locally uniformly in (t, ρ). Hence,

F λn (ϕ n (t))(ρ 1 , ρ 2 ) → F λ (ϕ(t))(ρ 1 , ρ 2 ), locally uniformly in (t, ρ 1 , ρ 2 ), where F λn is defined through (17). Since ϕ ′ n = F λn (ϕ n ), we conclude that ϕ ′ n → F λ (ϕ) locally uniformly in [0, ∞) × R 2
+ . Moreover, we have for each n that

ϕ n (t, ρ 1 , ρ 2 ) -1 = t 0 ϕ ′ n (s, ρ 1 , ρ 2 ) ds.
The lefthand side of this equation converges to ϕ(t, ρ 1 , ρ 2 ) -1, whereas the righthand side converges to t 0 F λ (ϕ(s))(ρ 1 , ρ 2 ) ds. This proves that ϕ solves (5) and that ϕ ∈ C 1 ([0, ∞); C(R 2 + )).

Remark 2. By sending R 1 to infinity in (38) we get that the solution ϕ constructed above satisfies the bounds

e -t(ρ 1 +ρ 2 ) + ρ 1 ρ 2 (1 -e -t(ρ 1 +ρ 2 ) ) 2ρ(ρ 1 + ρ 2 ) ≤ ϕ(t, ρ 1 , ρ 2 ) ≤ 1 + (ρ 1 + ρ 2 ρ )(ρ 2 + ρ 2 ρ ) 2ρ 1 -e -(ρ 1 +ρ 2 )t ρ 1 + ρ 2 . ( 43
)
From ( 24), ( 41), and the lower bound in ( 43) we get moreover that

0 < ψ(t, ρ) ≤ ρ 2 ρ . ( 44
)
Now we turn to prove the uniqueness of solutions in the class of functions ϕ ∈ C 1 ([0, ∞); C(R 2 + )) satisfying a bound of the form (7) To this end, let ϕ 1 and ϕ 2 be two solutions in that class and set

δ(t, ρ 1 , ρ 2 ) = ϕ 2 (t, ρ 1 , ρ 2 ) -ϕ 1 (t, ρ 1 , ρ 2 ).
We will show that δ(t) 2

L 2 ( λ⊗ λ)

= 0 for all t, whenever λ be a positive finite Borel measure of the form λ = λ + µ, where µ is a positive finite Borel measure with compact support. Taking, for instance, µ as the Lebesgue measure on [0, R] will then imply that δ(t, ρ 1 , ρ 2 ) = 0 for ρ 1 , ρ 2 ∈ [0, R]. So this will give the uniqueness of solutions.

Let us define F λ as

F λ (f )(ρ 1 , ρ 2 ) = 1 2ρ ρ 1 + xf (ρ 1 , x) λ(dx) ρ 2 + xf (x, ρ 2 ) λ(dx) . ( 45
)
Lemma 6. We have

F λ (ϕ i (t)) ∈ L 2 ( λ ⊗ λ) and F λ (ϕ 1 (t)) -F λ (ϕ 2 (t)) 2 L 2 ( λ⊗ λ) ≤ C( δ(t) 2 L 2 ( λ⊗ λ) + δ(t) 4 L 2 ( λ⊗ λ) ), ( 46
)
where C a positive constant that depends only on ρ and ρ 2 .

Proof. For simplicity, we will drop the argument t throughout the proof. We may write

F λ (ϕ 2 ) -F λ (ϕ 1 ) = 1 2ρ xδ(ρ 1 , x) λ(dx) ρ 2 + xϕ 1 (x, ρ 2 ) λ(dx) + xδ(x, ρ 2 ) λ(dx) ρ 1 + xϕ 1 (ρ 1 , x) λ(dx) + xδ(ρ 1 , x) λ(dx) xδ(x, ρ 2 ) λ(dx) .
Thus,

( F λ (ϕ 2 ) -F λ (ϕ 1 )) 2 ≤ 3 4ρ 2 xδ(ρ 1 , x) λ(dx) ρ 2 + xϕ 1 (x, ρ 2 ) λ(dx) 2 + xδ(x, ρ 2 ) λ(dx) ρ 1 + xϕ 1 (ρ 1 , x) λ(dx) 2 + xδ(ρ 1 , x) λ(dx) xδ(x, ρ 2 ) λ(dx) 2 .
Now we integrate this inequality with respect to λ(dρ 1 ) λ(dρ 2 ). The two first terms can be analyzed in the same way. First, we observe that (ρ 2 + ρ 2 ρ ) 2 λ(dρ 2 ) is finite. Then we note that

xδ(ρ 1 , x) λ(dx) 2 ≤ ρ 2 δ(ρ 1 , x) 2 λ(dx) ≤ ρ 2 δ(ρ 1 , x) 2 λ(dx). ( 47
)
Hence,

xδ(ρ 1 , x) λ(dx) 2 λ(dρ 1 ) ≤ ρ 2 δ 2 L 2 ( λ⊗ λ) .
Thus, the two first terms can be bounded by

C 0 δ 2 L 2 ( λ⊗ λ)
, where C 0 is a constant that only depends on ρ and ρ 2 . Using once again (47), we get that the third term can be bounded from above by

C 1 δ 4 L 2 ( λ⊗ λ)
, where the constant C 1 depends only on ρ and ρ 2 . Now we differentiate δ 2 and integrate over [0, t]:

δ(t, ρ 1 , ρ 2 ) 2 = -2 t 0 (ρ 1 + ρ 2 )δ(s, ρ 1 , ρ 2 ) 2 ds +2 t 0 δ(s, ρ 1 , ρ 2 ) F λ (ϕ 2 (s))(ρ 1 , ρ 2 ) -F λ (ϕ 1 (s))(ρ 1 , ρ 2 ) ds ≤ 2 t 0 δ(s, ρ 1 , ρ 2 ) F λ (ϕ 2 (s))(ρ 1 , ρ 2 ) -F λ (ϕ 1 (s))(ρ 1 , ρ 2 ) ds.
We now integrate w.r.t. λ(dρ 1 ) λ(dρ 2 ) and get by using the Cauchy-Schwarz inequality,

δ(t, •) 2 L 2 ( λ⊗ λ) ≤ 2 t 0 δ(s, •) L 2 ( λ⊗ λ) F λ (ϕ 2 (s)) -F λ (ϕ 1 (s)) L 2 ( λ⊗ λ) ds.
By continuity of t → δ(t, •) L 2 ( λ⊗ λ) , we know that for each T > 0 there is a constant

K such that δ(t, •) L 2 ( λ⊗ λ) ≤ K when t ∈ [0, T ]. Thus, we get from Lemma 6 that δ(t, •) 2 L 2 ( λ⊗ λ) ≤ C(1 + K 2 ) t 0 δ(s, •) 2 L 2 ( λ⊗ λ) ds, which in turn gives that δ(t, •) 2 L 2 ( λ⊗ λ)
= 0 on [0, T ] by Gronwall's Lemma. This concludes the proof of uniqueness. Now we turn to proving the properties (a) through (f) in Theorem 1. Property (a) (strict positivity) can be proved just as in the case of a compactly supported measure λ in Lemma 2. Property (b) (symmetry) is already clear. Property (c) ( ϕ(t, ρ, x) λ(dx) = 1) follows from the corresponding property of the approximating functions ϕ n , the uniform bounds (38), and dominated convergence.

Property (d) states that ϕ ∈ C 2 ([0, ∞); C(R 2 + )). By dominated convergence and the bound (39), which also holds for ϕ ′ in place of ϕ ′ n , we get that ψ(t, ρ) belongs to C 1 ([0, ∞); C(R + )). Thus, our ODE gives ϕ ′ ∈ C 1 ([0, ∞); C(R 2 + )), which proves property (d). We now prove property (e). It is clearly enough to prove it when f : R + → R is a bounded measurable function with compact support. To this end, let φ(t, ρ 1 , ρ 2 ) := e t(ρ 1 +ρ 2 ) ϕ(t, ρ 1 , ρ 2 ). Then φ belongs to

C 1 ([0, ∞); C(R 2 + )) and φ′ (t, ρ 1 , ρ 2 ) = 1 2ρ • e t(ρ 1 +ρ 2 ) (ρ 1 + ψ(t, ρ 1 ))(ρ 2 + ψ(t, ρ 2 )).
That is, φ′ (t, ρ 1 , ρ 2 ) = g(t, ρ 1 )g(t, ρ 2 ) for a function g. Thus,

f (x 1 )f (x 2 ) φ′ (t, x 1 , x 2 ) λ(dx 1 ) λ(dx 2 ) = f (x)g(t, x) λ(dx) 2 ≥ 0. Since φ(t) = 1 + t 0 φ′ (s) ds, we find that φ(t) is nonnegative definite. Finally, with f (x) = e -tx f (x), f (x 1 )f (x 2 )ϕ(t, x 1 , x 2 ) λ(dx 1 ) λ(dx 2 ) = f (x 1 ) f (x 2 ) φ(t, x 1 , x 2 ) λ(dx 1 ) λ(dx 2 ) ≥ 0.
This establishes property (e) in Theorem 1. Finally, property (f) (the local Lipschitz property for ϕ and ϕ ′ ) follows just as in Lemma 5. This concludes the proof of Theorem 1.

Proof of Theorem 2

The strategy in the proof of Theorem 2 is to use a verification argument and based on guessing the optimal costs V (T, E(•), x) for liquidating x shares over [0, T ] with additional and arbitrary initial data E(•). The result of our guess is formula (48) below. We explain its heuristic derivation in Appendix A.

Let ϕ be a solution of the infinite-dimensional Riccati equation ( 5), ( 6). This solution gives rise to a family of linear operators Φ

t : L 2 (λ) → L 2 (λ) ∩ C(R + ) defined by Φ t f (ρ) = f (x)ϕ(t, x, ρ) λ(dx), f ∈ L 2 (λ). By (7), t → Φ t f is a continuous map into both L 2 (λ) and C(R + ) for each f ∈ L 2 (λ)
. By the inequality (39), which also hold for ϕ in place of

ϕ n , t → Φ t f is a continuously differentiable map into both L 2 (λ) and C(R + ) for each f ∈ L ∞ (λ). For t ≥ 0, E(•) ∈ L 2 (λ)
, and x ∈ R, we define

V (t, E(•), x) := 1 2 1 ϕ 0 (t) x -(Φ t E)(0) 2 -E, Φ t E , (48)
where •, • denotes the usual inner product in L 2 (λ). For t ∈ [0, T ] and a [0, T ]-admissible strategy X we define

C X t := [0,t) E X s dλ dX s + 1 2 s<t (∆X s ) 2 + V (T -t, E X t , X t ).
By Lemma 1, the first two terms on the right correspond to the cost accumulated by the strategy up to time t. Moreover,

V (0, E X T , X T ) = 1 2 (X T ) 2 -E X T dλ • X T = 1 2 (∆X T ) 2 + E X T dλ • ∆X T ,
due to the requirement X T + = 0. This gives C X T = C(X). Our goal is thus to show the following verification lemma: dC X t ≥ 0 with equality if and only if X = X * for a certain strategy X * . This will identify X * as the optimal strategy and V (T, E(•), x) as the optimal cost for liquidating x shares over [0, T ] with additional initial data E(•) at time t = 0. In the formalism of potential theory, V (T, 0, -1) will then be the minimal energy of a probability measure on [0, T ].

Lemma 7. For every [0, T ]-admissible strategy X, C X t is absolutely continuous in t and

dC X t dt = 1 2 ψ(T -t, 0) ϕ 0 (T -t) X t -(Φ T -t E X t )(0) + E X t (ρ) ρ + ψ(T -t, ρ) λ(dρ) 2 (49) for a.e. t ∈ [0, T ].
Proof. Recall the following integration by parts formula for leftcontinuous functions α t , β t of locally bounded variation:

α t β t -α s β s = [s,t) α r dβ r + [s,t) β r dα r + r∈[s,t)
∆α r ∆β r .

It follows that t → E X t (ρ) is of bounded variation and

E X t (ρ) -E X s (ρ) = X t -X s -ρ t s E X r (ρ) dr (50)
as well as

E X t (ρ 1 )E X t (ρ 2 ) -E X s (ρ 1 )E X s (ρ 2 ) = [s,t) E X r (ρ 1 ) + E X r (ρ 2 ) dX r - t s (ρ 1 + ρ 2 )E X r (ρ 1 )E X r (ρ 2 ) dr + r∈[s,t) (∆X r ) 2 . ( 51
)
Therefore,

ϕ(T -t, ρ 1 , ρ 2 )E X t (ρ 1 )E X t (ρ 2 ) -ϕ(T -s, ρ 1 , ρ 2 )E X s (ρ 1 )E X s (ρ 2 ) = - t s ϕ ′ (T -r, ρ 1 , ρ 2 )E X r (ρ 1 )E X r (ρ 2 ) dr - t s ϕ(T -r, ρ 1 , ρ 2 )E X r (ρ 1 )E X r (ρ 2 )(ρ 1 + ρ 2 ) dr + [s,t) ϕ(T -r, ρ 1 , ρ 2 ) E X r (ρ 1 ) + E X r (ρ 2 ) dX r + r∈[s,t) ϕ(T -r, ρ 1 , ρ 2 )(∆X r ) 2 .
We have already observed in the proof of Lemma 1 that |E X r (ρ)| is uniformly bounded in r ∈ [0, T ] and ρ ≥ 0 by the total variation of X. Hence we may integrate both sides of the preceding identity with respect to λ(dρ 1 ) λ(dρ 2 ) to obtain, with the symmetry of ϕ and the notation E X t := ρE X t (ρ), that

E X t , Φ T -t E X t -E X s , Φ T -s E X s = - t s E X r , Φ ′ T -r E X r dr -2 t s E X r , Φ T -r E X r dr + 2 [s,t) 1, E X r dX r + r∈[s,t) (∆X r ) 2 .
By a similar reasoning we obtain

Φ T -t E X t (0) -Φ T -s E X s (0) = - t s Φ ′ T -r E X r (0) dr - t s Φ T -r E X r (0) dr + X t -X s and X t -Φ T -t E X t (0) 2 -X s -Φ T -s E X s (0) 2 = 2 t s X r -Φ T -r E X r (0) Φ ′ T -r E X r (0)+Φ T -r E X r (0) dr.
Using these formulas, we can now compute

C X t -C X s = 1 2 t s ϕ ′ 0 (T -r) ϕ 2 0 (T -r) X r -Φ T -r E X r (0) 2 dr + t s X r -Φ T -r E X r (0) ϕ 0 (T -r) Φ ′ T -r E X r (0) + Φ T -r E X r (0) dr + 1 2 t s E X r , Φ ′ T -r E X r dr + t s E X r , Φ T -r E X r dr
Therefore, C X t is absolutely continuous on [0, T ] and has the derivative

dC X t dt = 1 2 ϕ ′ 0 (T -t) ϕ 2 0 (T -t) X t -Φ T -t E X t (0) 2 + X t -Φ T -t E X t (0) ϕ 0 (T -t) Φ ′ T -t E X t (0) + Φ T -t E X t (0) + 1 2 E X t , Φ ′ T -t E X t + E X t , Φ T -t E X t for a.e. t.
To further analyze the preceding formula, we take an extra point ∆. We let λ := λ + δ ∆ and extend E X t and ϕ to functions on {∆} ∪ [0, ∞) by putting

E X t (∆) := 1 ϕ 0 (T -t) X t -(Φ T -t E X t )(0) ,
t ≥ 0, ϕ(t, ∆, ρ) = ϕ(t, ρ, ∆) := ϕ(t, 0, ρ), ϕ(t, ∆, ∆) := ϕ(t, 0, 0) = ϕ 0 (t).

(52)

We furthermore define the function

f (x) = x if x ∈ [0, ∞), 0 if x = ∆,
and we extend the definition of E X via E X t (x) = f (x)E X t (x) for x ∈ {∆} ∪ [0, ∞). Finally, we set for g ∈ L 2 (λ), Φ t g(x) = ϕ(t, x, y)g(y)λ(dy) and one easily checks that Φ :

L 2 (λ) → L 2 (λ).
With this notation, we get

dC X t dt = 1 2 E X t , Φ ′ T -t E X t L 2 (λ) + E X t , Φ T -t E X t L 2 (λ) = 1 2 E X t , Φ ′ T -t E X t L 2 (λ) + E X t , Φ T -t E X t L 2 (λ) + Φ T -t E X t , E X t L 2 (λ) = 1 2 E X t (x)E X t (y) ϕ ′ (T -t, x, y) + (f (x) + f (y))ϕ(T -t, x, y) λ(dx) λ(dy) = 1 2 E X t (x)E X t (y) f (x) + ψ(T -t, x) f (y) + ψ(T -t, y) λ(dx) λ(dy) = 1 2 E X t (x) f (x) + ψ(T -t, x) λ(dx) 2 ,
where we have used the Riccati equation ( 5) and the notation (11) in the fourth step. This proves the assertion.

It follows from the above that a [0, T ]-admissible strategy X * with X

* 0 = x satisfies C(X * ) = V (T, 0, x) ≤ C(X)
for all other [0, T ]-admissible strategies X with X 0 = x if dC X * t /dt vanishes for a.e. t ∈ [0, T ]. Using (49), we write this latter condition as

0 = X * t + E X * t (ρ)θ(T -t, ρ) λ(dρ) for a.e. t, (53)
where

θ(τ, ρ) = ϕ 0 (τ )(ρ + ψ(τ, ρ)) ψ(τ, 0) -ϕ(τ, ρ, 0). Then θ(τ, ρ) λ(dρ) = ϕ 0 (τ )2ρ ψ(τ, 0) -1. ( 54
)
Plugging this and ( 50) into (53) yields that for a.e. t

X * t = X * 0 1 - ψ(T -t, 0) ϕ 0 (T -t)2ρ + ψ(T -t, 0) ϕ 0 (T -t)2ρ t 0 ρE X * s (ρ)θ(T -t, ρ) λ(dρ) ds. (55)
Thus, the leftcontinuous function X * t coincides with an absolutely continuous function for a.e. t ∈ [0, T ]. It follows that these two functions coincide for every t ∈ (0, T ]. Thus, E X * is continuous on (0, T ] by (50), which in turn implies via (55) that X * is continuously differentiable throughout (0, T ).

When taking the limit t ↓ 0 in (55), we get

X * 0+ = x 1 - ψ(T, 0) ϕ 0 (T )2ρ , which gives ∆X * 0 = - ψ(T, 0) 2ρϕ 0 (T ) x. (56)
That ∆X * 0 = ∆X * T follows from Remark 2.10 in [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF]. Since ψ(t, ρ) is continuously differentiable in t, θ(t, ρ) is also continuously differentiable in t. Differentiating (53) with respect to t ∈ (0, T ) yields

0 = d dt X * t + dE X * t (ρ) dt θ(T -t, ρ) λ(dρ) -E X * t (ρ)θ ′ (T -t, ρ) λ(dρ) = 2ρϕ 0 (T -t) ψ(T -t, 0) • d dt X * t -E X * t (ρ) θ ′ (T -t, ρ) + ρθ(T -t, ρ) λ(dρ),
where we have used ( 50) and (54) in the second step. This gives

d dt X * t = ψ(T -t, 0) 2ρϕ 0 (T -t) E X * t (x) θ ′ (T -t, x) + xθ(T -t, x) λ(dx). (57)
Then, we have

ϕ(T -t, ρ 1 , ρ 2 ) = [t,T ] [t,T ] G(|r -s|) dB s (ρ 1 ) dB r (ρ 2 ) = ϕ(T -t, ρ 2 , ρ 1 ), ϕ(T -t, ρ 1 , ρ 2 ) λ(dρ 1 ) = [t,T ] G(r -t)dB r (ρ 2 ) = 1.
Moreover, we observe that ϕ(T -t, 0, ρ) = -B t (ρ) since B T (ρ) = 0. This gives in particular that Z t = -E X t (ρ)ϕ(T -t, 0, ρ)λ(dρ). Besides, we have

[t,T ] [t,T ] G(|s -r|) dZ r dZ s = [t,T ] A(s) dZ s = E X t (ρ 1 ) [t,T ] e -ρ 1 (s-t) dZ s λ(dρ 1 ) = E X t (ρ 1 )E X t (ρ 2 ) [t,T ] e -ρ 1 (s-t) dB s (ρ 2 ) λ(dρ 1 ) λ(dρ 2 ) = E X t (ρ 1 )E X t (ρ 2 )ϕ(T -t, ρ 1 , ρ 2 ) λ(dρ 1 ) λ(dρ 2 ).
Thus, we obtain altogether that

C t (X) = 1 2 X t -E X t (ρ)ϕ(T -t, 0, ρ)λ(dρ) 2 ϕ(T -t, 0, 0) - E X t (ρ 1 )E X t (ρ 2 )ϕ(T -t, ρ 1 , ρ 2 ) λ(dρ 1 ) λ(dρ 2 ) .
This suggests the definition

V (t, E(.), x) := 1 2 X -E(ρ)ϕ(T -t, 0, ρ)λ(dρ) 2 ϕ(t, 0, 0) - E(ρ 1 )E(ρ 2 )ϕ(t, ρ 1 , ρ 2 ) λ(dρ 1 ) λ(dρ 2 ) , so that we have C(X) = [0,t) E X s (ρ) λ(dρ) dX s + 1 2 s<t (∆X s ) 2 + V (T -t, E t (.), X t ). (62)
Thus, we have obtained the formula given by ( 48), and it remains to explain why ϕ should solve a Riccati equation. To do so, we consider an arbitrary strategy (X s , s ∈ [0, T ]) and consider the cost (62), which is the cost of the strategy that is equal to X on [0, t) and optimal on [t, T ]. To make the dependence on t explicit, we denote this cost by C t (X). To simplify things, we will focus on the particular case (3) of a discrete measure λ andd i=0 λ i = 1. With this choice, V only depends on E(ρ i ), 0 ≤ i ≤ d. We introduce the following notations:

(dx) = d i=0 λ i δ ρ i (dx), with ρ 0 = 0 < ρ 1 < • • • < ρ d , λ i ≥ 0,
ϕ ij (t) = ϕ(t, ρ i , ρ j ), 0 ≤ i, j ≤ d, V (t, E 0 , . . . , E d , X) = 1 2 1 ϕ 00 (t) X - d i=0 λ i E i ϕ 0i (t) 2 - d i=0 d j=0 λ i λ j E i E j ϕ ij (t) , E X t = d i=0 λ i E X t (ρ i ), C t (X) = [0,t) E X s dX s + 1 2 s<t (∆X s ) 2 + V (T -t, E X t (ρ 0 ), . . . , E X t (ρ d ), X t ).
Lemma 9. We have ∆C t (X) = 0 for all t ∈ [0, T ].

Proof. Note that ∆X t = ∆E X t (ρ). Hence it is clear that ∆C t (X) = 0 if ∆X t = 0. Now suppose that ∆X t = 0. Then

∆C t (X) = d i=0 λ i E X t (ρ i )∆X t + 1 2 (∆X t ) 2 + ∆ C t (X). ( 63
)
On the other hand, we have

V (t, E 0 + δ, . . . , E d + δ, X + δ) = 1 2 1 ϕ 00 (t) X + δ - d i=0 λ i (E i + δ)ϕ 0i (t) 2 - d i=0 d j=0 λ i λ j (E i + δ)(E j + δ)ϕ ij (t) = V (t, E 0 , . . . , E d , X) -δ d j=0 d i=0 λ i λ j E i ϕ ij (t) - δ 2 2 d j=0 d i=0 λ i λ j ϕ ij (t) = V (t, E 0 , . . . , E d , X) -δ d i=0 λ i E i - δ 2 2 .
Here we have used the facts that i λ i = 1 and i λ i ϕ ij = 1. Putting everything together yields the assertion.

We can now focus on infinitesimal variations, and we denote V t := ∂V /∂t, V i := ∂V /∂E i and V X := ∂V /∂X. We have, when ∆X t = 0,

dC t (X) = E X t dX t -V t dt + d i=0 V i dE X t (ρ i ) + V X dX t = E X t + d i=0 V i + V X dX t -V t + d i=0 ρ i E X t (ρ i )V i dt
By simple calculations, we get E X t + d i=0 V i + V X = 0, and our expression simplifies to

dC t (X) = -V t + d i=0 ρ i E X t (ρ i )V i dt.
Let us now calculate V t :

V t = - ϕ ′ 00 2ϕ 2 00 X - d j=0 λ j E j ϕ 0j 2 - d i=0 λ i E i X - d j=0 λ j E j ϕ 0j ϕ ′ 0i ϕ 00 - 1 2 d i,j=0 λ i λ j E i E j ϕ ′ ij .
To simplify computations, we define λ -1 := 1 and E -1 := X -d j=0 λ j E j ϕ 0j ϕ 00 (64)

as well as ϕ -1i := ϕ 0i and ϕ -1-1 = ϕ 00 . Then

V t = - 1 2 d i,j=-1 λ i E i λ j E j ϕ ′ ij .
With ρ -1 := ρ 0 = 0, we get V i = -λ i λ -1 E -1 ϕ -1i + d j=0 λ j E j ϕ ij = -λ i d j=-1 λ j E j ϕ ij and therefore

d i=0 ρ i E i V i = - d i,j=-1 ρ i + ρ j 2 λ i E i λ j E j ϕ ij .
Altogether, we obtain

dC t (X) = 1 2 d i,j=-1 λ i E X t (ρ i )λ j E X t (ρ j ) ϕ ′ ij + (ρ i + ρ j )ϕ ij dt,
where E X t (ρ -1 ) is defined according to (64).

Thus, we arrive at the following quadratic form 1 2 d k,l=0 E k E l λ k λ l (ϕ ′ kl + (ρ k + ρ l )ϕ kl ). Since we should have dC t (X) ≥ 0 with dC t (X) = 0 for the optimal strategy, this quadratic form should be nonnegative with rank one. Indeed, since the control X is of dimension one, it would not be possible to make dC t (X) = 0 if the rank of the quadratic form were higher than two. We are now going to write the conditions that ensures that this quadratic form is nonnegative with rank one. To do so, We introduce the new coordinates (∆ 0 , . . . , ∆ d ) such that

E 0 = ∆ 0 , E 1 = ∆ 0 + ∆ 1 , • • • E d = ∆ 0 + ∆ d .
In these coordinates, our quadratic form becomes 1 2 ∆ 2 0 λ 2 0 (ϕ ′ 00 + 2ρ 0 ϕ 00 ) +

d l=1 ∆ 0 (∆ 0 + ∆ l )λ 0 λ l (ϕ ′ 0l + (ρ 0 + ρ l )ϕ 0l ) + 1 2 d k=1 d l=1 (∆ 0 + ∆ k )(∆ 0 + ∆ l )λ k λ l (ϕ ′ kl + (ρ k + ρ l )ϕ kl )
After some calculations, we get that the coefficient for ∆ 2 0 , ∆ 0 ∆ l , ∆ k ∆ l and ∆ 2 l (for 1 ≤ k, l ≤ d) are respectively ρ, λ l ρ l + d k=0 λ k ρ k ϕ kl , λ k λ l (ϕ ′ kl + (ρ k + ρ l )ϕ kl ) and 1 2 λ 2 l (ϕ ′ ll + 2ρ l ϕ ll ). Thus, the matrix Q for the quadratic form has coefficients

Q 00 = ρ, Q 0l = Q l0 = λ l 2 ρ l + d k=0 λ k ρ k ϕ kl , Q kl = Q lk = 1 2 λ k λ l (ϕ ′ kl + (ρ k + ρ l )ϕ kl ) if k, l ≥ 2, k = l, Q ll = 1 2 λ 2 l (ϕ ′ ll + 2ρ l ϕ ll ).
Since Q is of rank one, the determinant of the matrices

Q 00 Q 0l Q l0 Q ll , Q 00 Q 0l Q k0 Q kl

Figure 1 :

 1 Figure1: Comparison of the approximated optimal strategies (x i , 0 ≤ i ≤ N) obtained with (16) and with the method based on the Riccati equation (slightly shifted to the right) .

  must vanish for l = 1, . . . , d and k < l. That gives, respectively, kl + (ρ k + ρ l )ϕ kl = 1 2ρ ρ l + d i=0 λ i ρ i ϕ il ρ k + d j=0 λ j ρ j ϕ kj ,

The first author acknowledges the support of the "Chaire Risques Financiers" of Fondation du Risque. The second author gratefully acknowledges support by Deutsche Forschungsgemeinschaft DFG

We now want to simplify (57). To this end, we use the notation ψ(t) := xψ(t, x) λ(dx). and the formulas ϕ ′ (t, ρ, 0) = -ρϕ(t, ρ, 0) + 1 2ρ (ρ + ψ(t, ρ))ψ(t, 0), ϕ ′ 0 (t) = 1 2ρ ψ(t, 0) 2 ψ ′ (t, ρ) = -x 2 ϕ(t, x, ρ) λ(dx) -ρψ(t, ρ) + 1 2ρ (ρ + ψ(t, ρ))(ρ 2 + ψ(t)) (58)

Then a tedious computation shows that

where Θ(t, ρ) is as in ( 14). Therefore, (57) becomes

Now we have

Plugging this formula back into (59) and using Fubini's theorem yields that x(t) = d dt X * t solves the Volterra integral equation ( 13). This is a Volterra integral equation of the second kind with continuous kernel K(t, s) and continuous function f (t). It hence admits a unique continuous solution x(•) (Linz 1985, Theorem 3.1). Conversely, given such a solution x(•), we can define a [0, T ]-admissible strategy X * via (56) and x(t) = d dt X * t . Then X * satisfies (55) for t = 0+ as well as (57) for t > 0. Integrating (57) and reversing the steps made above in deriving ( 57) from ( 55) shows that X * satisfies (55) for t ∈ (0, T ], and so X * is optimal. This concludes the proof of Theorem 2.

A Heuristic derivation of the value function

We want to explain here how it is possible to guess the value function V (t, E(•), x) introduced in (48). We start our discussion by deriving a formula for the costs C(X) of a strategy (X s , s ∈ [0, T ]) that is arbitrary on [0, t) and optimal on [t, T ]. We set

G(|s -r|) dX s dX r .

The rightmost integral can be written as

The first-order condition of optimality thus reads

where ν is a suitable Lagrange multiplier (compare Theorem 2.11 in [START_REF] Gatheral | Transient linear price impact and Fredholm integral equations[END_REF]).

Lemma 8. Suppose that R and R are functions with finite variation such that R T = R T and

which implies the assertion in view of the fact that G is strictly positive definite.

Now suppose that we have auxiliary functions with finite variation B t (ρ) such that B T (ρ) = 0 and

We also define

Therefore,

G(|s -r|) dB r (0) for t ≤ s ≤ T .

Lemma 8 hence implies that X s + Z s = νB s (0) for t ≤ s ≤ T . Hence, we get

From these identities we get

G(|s -r|) dZ r dZ s , since the first double integral is equal to -B t (0). Now we define

e -ρ 1 (r-t) dB r (ρ 2 ), ρ 1 , ρ 2 ≥ 0.

which gives precisely the Riccati equation. Thus, equation ( 65) holds for 1 ≤ k, l ≤ d. In fact, the choice of E 0 = ∆ 0 is arbitrary. Had we chosen E i = ∆ i for some i > 0, we had obtained (65) for k, l = i. Therefore (65) holds in fact for all k, l = 0, . . . , d.


Mehrdad Yaghoobi
email: yaghoobi@ieee.org

Nam ‡ Sangnam
email: sangnam.nam@inria.fr

Rémi Gribonval
email: remi.gribonval@inria.fr

Mike E Davies
email: mike.davies@ed.ac.uk

NOISE AWARE ANALYSIS OPERATOR LEARNING FOR APPROXIMATELY COSPARSE SIGNALS

Keywords: Sparse Approximation, Analysis Framework, Cosparse Signal Model, Analysis Operator Learning, Douglas-Rachford Splitting

This paper investigates analysis operator learning for the recently introduced cosparse signal model that is a natural analysis complement to the more traditional sparse signal model. Previous work on such analysis operator learning has relied on access to a set of clean training samples. Here we introduce a new learning framework which can use training data which is corrupted by noise and/or is only approximately cosparse. The new model assumes that a p-cosparse signal exists in an epsilon neighborhood of each data point. The operator is assumed to be uniformly normalized tight frame (UNTF) to exclude some trivial operators. In this setting, an alternating optimization algorithm is introduced to learn a suitable analysis operator.

INTRODUCTION

The linear sparse model is a suitable low-dimensional model for many natural signals. A discrete signal y ∈ R m is called sparse when we can represent y using a vector x ∈ R n , with few non-zero elements, and a linear overcomplete transform Φ ∈ R m×q m < q, called the dictionary, as follows y = Φx.

(

) 1
This has also been called the synthesis parsimony model [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. As Φ is overcomplete, the model is underdetermined and for many reasons a convex sparsity penalty, i.e. ℓ1(x) := x 1 = P |xi|, has been used to recover the sparse coefficients x.

The formulation ( 1) is an ideal model which often does not precisely fit to the real world signals either because the signal is noisy or is only approximately sparse. A standard approach to compensate the effect of noise and/or model mismatch, is to assume that the noise is Gaussian and the signal is within an ǫ neighborhood of the observation. The linear constraint can then be relaxed to give the following convex program, called Basis Pursuit Denoising (BPDN) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF],

min x x 1 s. t. y -Φx 2 ≤ ǫ. (2)
A better approximation is sometimes possible by subsequently projecting y onto the span of the columns of Φ in the support of x: a technique known as debiasing.

This work is supported by EU FP7, FET-Open grant number 225913 and EPSRC grant EP/F039697/1. MED acknowledges support of his position from the Scottish Funding Council and their support of the Joint Research Institute with the Heriot-Watt University as a component part of the Edinburgh Research Partnership.

Dictionary Learning

In many instances the appropriate dictionary Φ is not known and must be learnt from training data Y = [yi]i∈I ∈ R m×L [START_REF] Tosic | Dictionary learning[END_REF]. One approach to synthesis dictionary learning is to solve (2) in a matrix form and apply a constraint on the dictionary to exclude trivial solutions, as follows:

min X,Φ∈D X 1 s. t. Y -ΦX 2 ≤ ǫ. (3)
Such an optimization has been dealt with in [START_REF] Yaghoobi | Dictionary learning for sparse representations: A pareto curve root finding approach[END_REF]. As Φ does not appear in the objective, many alternating optimization methods have difficulties with this formulation and they use a Lagrangian multiplier version of (3), see for example [START_REF] Yaghoobi | Dictionary learning for sparse approximations with the majorization method[END_REF][START_REF] Gribonval | Dictionary identificationsparse matrix-factorisation via ℓ1 minimisation[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF][START_REF] Geng | On the local correctness of l1 minimization for dictionary learning[END_REF] and references therein.

Cosparse and Approximately Cosparse Models

An alternative low-dimensional signal model is the analysis parsimony model [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. In this setting, there exists an overcomplete analysis linear operator Ω ∈ R N×m , which maps the signals y ∈ R m to low-dimensional subspaces of a higher-dimensional analysis space R N [START_REF] Nam | The cosparse analysis model and algorithms[END_REF]. In other words, z = Ωy, where z has many zero elements, i.e. z 0 = Np, where p is the number of zero components. The integer p, which is called the cosparsity of y [START_REF] Nam | The cosparse analysis model and algorithms[END_REF], has an important role in the analysis of signal recovery in linear inverse problems. Similar to the synthesis framework, most of the real world signals are not exactly cosparse: some noise corrupts the observations and there are also some model mismatches. It suggests an approximate cosparse model like z ≈ Ωy where z is cosparse and the accuracy of the approximation depends on the noise/model mismatch level. In other words, the model assumes that some cosparse signal e y, such that z := Ωe y has many zero elements, is within an ǫ neighborhood of y, i.e. e yy 2 ≤ ǫ [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF].

Analysis Operator Learning (AOL)

Unlike the synthesis framework, there is very little work on adapting the analysis model to a set of sample data Y. Recent works that have begun to investigate this problem are [START_REF] Peyre | Learning analysis sparsity priors[END_REF][START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF][START_REF] Ophir | Sequential minimal eigenvalues an approach to analysis dictionary learning[END_REF]. The main difficulty in the analysis operator learning is avoiding trivial solutions emerging from the learning process. [START_REF] Ophir | Sequential minimal eigenvalues an approach to analysis dictionary learning[END_REF] randomly cycles through the training samples and [START_REF] Peyre | Learning analysis sparsity priors[END_REF] implicitly constraints the operator update, to indirectly avoid such trivial solutions. Alternatively a constrained optimization approach has been presented in [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] that avoids most trivial solutions. One restriction of the AOL method of [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] is that it requires noiseless co-sparse data that exactly follow the model z = Ωy. Removing this restriction is the goal of the present paper.

CONSTRAINED ANALYSIS OPERATOR LEARNING

Cosparse signals, by definition, have a large number of zero elements in the analysis space. One proposal for adapting an operator to a set of signals Y is to maximize the number of zeros of ΩY or equivalently minimized the number of non-zero elements. The algorithm presented in [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] uses this approach along with the ℓ1 convex surrogate function for measuring sparsity. If we do not apply a constraint on the operator then trivial solutions such as Ω = 0 will invariably exist. It is explained in [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] why some simple constraints like row norm, row norm plus full-rank and the tight frame constraints do not resolve the issue. Instead, a constraint C called Uniform Normalized Tight Frame (UNTF) was proposed to simply exclude such trivial solutions. The constrained AOL can then be reformulated as,

min Ω ΩY 1 s. t. Ω ∈ C (4)
where,

C = {Ω ∈ R N×m : Ω T Ω = I, ∀i ωi 2 = c}, (5)
and ωi is the ith row of Ω. As • 1 is continuous and C is not convex, (4) may have many separate local optima. A projected sub-gradient based iterative method is presented in [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] to find a local minimum.

For convenience the procedure is summarized in Algorithm 1, where sgn is an entrywise version (i.e. {sgn(A)}i j = sgn(Ai j )) of the sign function which is defined by sgn(a) = sgn(a) when a = 0 and sgn(a) = [-1, 1] otherwise. The notations PUN (resp. PT F ) refer to projections onto the set of uniformly normalized (resp. tight frame) matrices, cf [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF]. The simulations in [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF] show that when the signals are cosparse enough, i.e. many zero elements in each Ωyi, generic analysis operators can be recovered using this framework.

Noise Aware AOL (NAAOL) Formulation

When the signal y is noise corrupted or only approximately cosparse an approximate signal b y, also called the denoised signal, can be found using an analysis denoising formulation as follows,

min b y Ωb y 1 + λ 2 b y -y 2 2 , (6)
where parameter λ regulates the level of noise/model missmatch. Such a noisy formulation motivates us to further extend the AOL framework to consider noise and approximate cosparsity. This leads to an optimization problem with two parameters b Y and Ω as follows,

min

Ω, b Y Ω b Y 1 + λ 2 b Y -Y 2 F s. t. Ω ∈ C. (7)
We use the constraint C defined in [START_REF] Yaghoobi | Dictionary learning for sparse approximations with the majorization method[END_REF], since the new formulation still has the same set of trivial solutions. Problem ( 7) is a non-convex optimization problem on the intersection of two manifolds, which is generally difficult to solve. A practical algorithm to find a local minimum of ( 7) is presented next.

Noise Aware AOL Algorithm

It is clear that the NAAOL problem can be simplified to the AOL problem, if b Y is known. Alternatively, if we know the optimal operator Ω, b Y can easily be found using a convex program. However, the size of the problem makes such an optimization challenging. The main difficulties are that ℓ1 is not differentiable and the existence of

Algorithm 1 Projected Subgradient Based AOL (Operator Update) Input: Y, Kmax, Ωin, stepsize η, threshold ǫ ≪ 1 1: initialization: k = 1, Ω [0] = 0, Ω [1] = Ωin 2: while ǫ ≤ Ω [k] -Ω [k-1] F and k ≤ Kmax do 3: ΩG = ∂f (Ω [k] ) = sgn(Ω [k] Y)Y T 4
:

Ω [k+1] = PUN n PT F n Ω [k] -η ΩG oo 5: k = k + 1 6: end while 7: output: Ωout = Ω [k-1]
Ω inside the ℓ1 penalty does not allow us to use conventional methods for solving problems like [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. We here use the Douglas-Rachford Splitting (DRS) technique to efficiently solve [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF] when Ω is fixed, here it is also called the alternating direction method of multipliers (ADMM), see [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF] for a brief overview. This technique has been used for approximation with the Total Variation (TV) penalty and analysis sparse approximation in [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF]

min b Y, b Z b Z 1 + λ 2 b Y -Y 2 F s. t. b Z = Ω b Y. (8)
The Augmented Lagrangian (AL) [START_REF] Rockafellar | Augmented lagrange multiplier functions and duality in nonconvex programming[END_REF] method is applied to solve [START_REF] Geng | On the local correctness of l1 minimization for dictionary learning[END_REF]. In the Lagrangian multiplier method we use the dual parameter B ∈ R N×L and add a penalty term,

D B, Ω b Y -b Z E .
In the AL, we also add an extra quadratic penalty related to the constraint and derive the new objective ψ( b Y, b Z, B) as follows,

ψ( b Y, b Z, B) = b Z 1 + λ 2 b Y -Y 2 F + γ D B, Ω b Y -b Z E + γ 2 Ω b Y -b Z 2 F = b Z 1 + λ 2 b Y -Y 2 F + γ 2 B + Ω b Y -b Z 2 F - γ 2 B 2 F
where 0 < γ ∈ R + is a constant parameter. According to the duality property, the solution of maxB minb Z, c Y ψ( b Y, b Z, B) coincides with the solution of [START_REF] Geng | On the local correctness of l1 minimization for dictionary learning[END_REF]. Using the DRS method, we iteratively optimize a convex/concave surrogate objective ψs( b Y and concave with respect to B. We can thus iteratively update each of the parameters, while keeping the rest fixed, as follows, b

Y, b Z, B, B [n] ) = ψ( b Y, b Z, B) -B -B [n] 2 F ,
Y [n+1] = (λI + γ Ω T Ω) -1 (λY + γ Ω T ( b Z [n] -B [n] )) b Z [n+1] = S 1 γ n Ω b Y [n+1] + B [n] o B [n+1] = B [n] + (Ω b Y [n+1] -b Z [n+1] ), (9)
where Sα, with an α > 0 , is the entrywise soft-threshold operator defined by Sα(β) = βα sgn(β) if |β| ≥ α and 0 otherwise.

Algorithm 2 Noise Aware Analysis Operator Learning (NAAOL)

Input: Y, Kmax, Ωin, η, ǫ ≪ 1, λ, γ 1: initialization: k = 1, Ω [0] = 0, Ω [1] = Ωin 2: while ǫ ≤ Ω [k] -Ω [k-1] F and k ≤ Kmax do 3:
Operator update from Algorithm 1.

4:

while not converged do 5:

Parameters update with (9).

Note that the update formula for b

Y involves a matrix inversion, which is computationally expensive. Here Ω is a tight frame and the matrix inversion is thus significantly simplified using the fact that Ω T Ω is identity. In this case, the operator (λI + γ Ω T Ω) -1 is simply a scaling with 1 λ+γ . We iterate [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] for a number of iterations or until the parameters cease to change significantly. Although the convergence of the iterative updates (9) can individually be investigated, it can also be deduced using the fact that it is a special case of DRS, which converges under mild conditions. Such an iterative algorithm only finds the solution of ( 8) and needs to be combined with the Operator Update (OU) of Algorithm 1. A pseudocode for the full NAAOL is presented in Algorithm 2.

SIMULATION RESULTS

For our experiments we used a set of face images, which are centred and cropped [START_REF] Lee | Acquiring linear subspaces for face recognition under variable lighting[END_REF]. Such images can be modeled as approximately piecewise smooth signals. A pseudo-random admissible Ω [0] ∈ R 128×64 has been used as an initial analysis operator and a training set of size L = 131072 of 8 × 8 image patches (m = 64) from 13 different faces, have been used to learn the operators (N = 128).

We applied both of the AOL and the NAAOL algorithms to demonstrate how much the cosparsities of the training samples increase using the noise aware formulation. AOL was iterated Kmax = 4000 times and NAAOL iterated for Kmax = 10 iterations with λ = γ = 0.5, while the inner-loop, i.e. Algorithm 1, was iterated 100 times. A plot of the analysis coefficients for an example y along with its corresponding cosparsities, with three different Ω, are presented in Figure 1. The initial operator Ω0 = Ω [0] has been used with y in (a). Not surprisingly, the signal is not cosparse with this arbitrary operator (p = 0). In (b), the same plot is drawn using the learned operator with AOL. Although some coefficients are small, most are not zero, and p = 1. In the last plot, we have shown the analysis coefficients for b y using the learned operator with NAAOL. It is clear that the cosparsity has been increased significantly (from p = 1 to p = 27). We have further plotted the cosparsities of the first 256 training samples y's using the learned operator found by AOL and corresponding approximations b y's, which are found by NAAOL, in Figure 2. This figure also shows, the operator learning (7) using the noise aware formulation [START_REF] Gribonval | Dictionary identificationsparse matrix-factorisation via ℓ1 minimisation[END_REF] results in much greater cosparsity.

The aim of second experiment is to denoise the image using a learned operator and a TV-type operator. We keep using the settings of the first experiment. The learned operator and the finite difference operator can now be used to denoise a corrupted version of another face from the database, using [START_REF] Gribonval | Dictionary identificationsparse matrix-factorisation via ℓ1 minimisation[END_REF]. The original face is shown in Figure 3 (a) and the noisy version with additive i.i.d Gaussian noise, is shown in (b). Denoising was performed using two different regularization settings: (λ = γ) = 0.3, 0.1. The bottom two rows show the denoised images using the TV-type operator and the learned dictionary. We can visually conclude that the two operators successfully denoise the corrupted images with some slight differences. The results with the learned operators are smoother (this is mostly visible on a screen rather than a printed copy of the paper).

As the initial goal was to increase the cosparsities of the signals, we have also shown the cosparsities of different patches of the selected face image. The horizontal axis presents the index number of the patches. To compare the cosparsity using these operators, we have plotted their differences and its average in the bottom plot. Negative values here demonstrates the cases when the finite difference operator is a better operator than the learned operator. The average, which is indeed negative, is plotted as a horizontal line. As a result, although the learned operator performs reasonably well, the finite difference still provides cosparser signals.

CONCLUSION

This paper presents a novel framework to learn an analysis operator which can be used for low-dimensional modeling of sample data. The new framework has the feature that it can account for the effects of noise and model mismatch within the training data. This fact improves the robustness of the previous AOL [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF]. The proposed framework is based on alternatively updating the operator and the approximations of signals. The operator can be updated using the previously proposed algorithm for AOL while a new fast algorithm for solving the signal approximation has been proposed. Two preliminary experiments are presented. The result of the first experiment demonstrates the ability of the new algorithm to compensate for model mismatch in the proposed framework. In the second experiment, we demonstrated the operator learning in a denoising application. Although the finite difference operator, which is known to be very successful for denoising of piecewise constant images, still performs better than the learned operator, the proposed operator learning scheme shows a promising result. For the signals where we do not know of any suitable operator, the proposed technique may provide some insight into this difficult low-dimensional modeling problem.

  where B [n] is the current estimation of B. The fixed points of the iterative updates of ψs( b Y, b Z, B, B [n] ) are the same as ψ( b Y, b Z, B), as the extra term B-B [n] 2 F vanishes in any fixed points. ψs( b Y, b Z, B, B [n] ) is convex with respect to b Z and b

Fig. 1 .

 1 Fig. 1. Signals in the analysis space R N , N = 128. The coefficients with almost zero magnitude, i.e. less than 0.01, are indicated with stars. The cosparsity in each case is: (a) p = 0, (b) p = 1, and (c) p = 27.

Fig. 2 .

 2 Fig. 2. The cosparsities of y (bottom plot) and b y (top plot) respectively with the operators learned with AOL and NAAOL.

Fig. 3 .

 3 Fig. 3. Face image denoising. Top row: original face (left), noisy face (right). Denoising results using (6). Middle row: λ = 0.3 using the learned analysis operator (left) and the finite difference operator (right). Bottom row: same as middle row with λ = 0.1.

Fig. 4 .

 4 Fig. 4. Cosparsities of image paches (a)-(b) and a comparison (c).

  [START_REF] Goldstein | The split Bregman method for l1-regularized problems[END_REF] 1 . We here only present a simple version of the DRS technique, tailored for this problem. A new parameter b Z := Ω b Y is defined, such that the optimization with respect to b Z is easier. Now the problem is a constrained (but still convex) program with two parameters b Z and b Y as follows,

[START_REF] Goldstein | The split Bregman method for l1-regularized problems[END_REF] derives the formulation by incorporating the Bregman distance. However, it has been shown that the new method, called Alternating Bregman Splitting method, is the same as DRS applied to the dual problems[START_REF] Setzer | Split Bregman Algorithm, Douglas-Rachford Splitting and Frame Shrinkage[END_REF].