Niv Dayan

Martin Kjaer Svendsen

Matias Bjørling

Philippe Bonnet

Luc Bouganim
email: luc.bouganim@inria.fr

EagleTree: Exploring the Design Space of SSD-Based Algorithms

Solid State Drives (SSDs) are a moving target for system designers: they are black boxes, their internals are undocumented, and their performance characteristics vary across models. There is no appropriate analytical model and experimenting with commercial SSDs is cumbersome, as it requires a careful experimental methodology to ensure repeatability. Worse, performance results obtained on a given SSD cannot be generalized. Overall, it is impossible to explore how a given algorithm, say a hash join or LSM-tree insertions, leverages the intrinsic parallelism of a modern SSD, or how a slight change in the internals of an SSD would impact its overall performance. In this paper, we propose a new SSD simulation framework, named EagleTree, which addresses these problems, and enables a principled study of SSD-Based algorithms. The demonstration scenario illustrates the design space for algorithms based on an SSD-based IO stack, and shows how researchers and practitioners can use EagleTree to perform tractable explorations of this complex design space.

INTRODUCTION

Flash-based Solid State Drives (SSDs) offer the same block device interface as hard disk drives (HDDs). It is thus seemingly transparent for a database administrator to replace HDDs by SSDs. The problem is that SSDs do not respect the performance contract that has always been valid for HDDs: e.g., sequential IOs are no longer orders of magnitude faster than random IOs. Worse, there is no consistent performance contract that all SSDs adhere to. We have argued in [START_REF] Bjørling | The Necessary Death of the Block Device Interface[END_REF] that this absence of a performance contract is due to the high level of software and hardware complexity, throughout the portion of the IO stack hidden behind the block device interface. So, even if SSDs look like HDDs, they behave in a very different way. In fact, the behavior of a given SSD is hard to characterize as it depends on undocumented internal features at the hardware level, e.g., SSD geometry, underlying flash chips, and mainly at the software level, e.g., the Flash Translation Layer (FTL) embedded on the SSD controller. The fact that complex SSDs internals are hidden and undocumented raises problems, not just in terms of experimental methodology -as care must be taken to bring an SSD to a well defined state before running experiments in order to obtain reproducible results [START_REF] Bonnet | System Co-Design and Data Management for Flash Devices[END_REF] -but essentially in terms of system design, as it is impossible for practitioners and researchers to explore how SSD internals impact overall system performance. More specifically, we seek to explore the following questions that have received little attention so far: • How does SSD parallelism impact performance (or the dual question: how can an algorithm efficiently leverage SSD parallelism)? A flash-based SSD contains tens to hundreds of flash chips wired in parallel to the SSD controller through multiple channels. The SSD scheduler deals with parallelism across and within flash chips and must decide which IO should be scheduled, and for writes, where (i.e., on which LUN1 ) it should be done, and precisely when (with respect to other IOs). • How do garbage collection and wear leveling interfere with application IOs? Each update leaves an obsolete flash page (with a before image). Over time, obsolete flash pages accumulate, and are reclaimed through garbage collection (GC).

In addition, the FTL relies on wear leveling (WL) to distribute the erase count across flash blocks and mask bad blocks. Note that both GC and WL read live pages from a victim block and write those pages at other locations, before that block is erased. As a result, GC and WL interfere with the application's IOs, possibly compromising throughput and contributing to latency variability. Note that application IOs also interfere with each other, which raises issues of fairness that have not been properly addressed so far. • What is the impact of replacing the block layer by a communication interface, hopefully better suited to deal with the complexity of SSDs? We have argued for such a cross layer approach [START_REF] Bjørling | The Necessary Death of the Block Device Interface[END_REF], but we cannot rely on existing SSDs to experiment with its design.

In order to study these questions, we must resort to simulations where SSD internals can be manipulated with controlled experiments. Existing SSD simulators fall short of this requirement.

Early efforts such as FlashSim [START_REF] Kim | FlashSim: A Simulator for NAND Flash-based Solid-State Drives[END_REF] (open source) and SSDSim [START_REF] Agrawal | Design tradeoffs for SSD performance[END_REF] (Microsoft license) are inadequate at the hardware level: they fail to represent SSD parallelism and chip characteristics. While the more recent open source simulator NANDFlashSim [START_REF] Jung | NANDFlashSim: Intrinsic latency variation aware NAND flash memory system modeling and simulation at microarchitecture level[END_REF] improves the hardware and parallelism representation, it does not expose a design space at the controller level that allows for rapid experimentation with mapping, GC, WL or scheduling policies.

In this paper, we present EagleTree, an open source software simulation framework for SSD-based applications. EagleTree encompasses the complete IO stack: application, operating system, SSD controller and flash chips array. EagleTree exposes a large, complex design space (see Section 2), which can be explored in a tractable way. The demonstration scenario (see Section 3) illustrates the complexity of this design space and shows how researchers and practitioners can use EagleTree to explore interplays between parameters across the IO stack.

THE EAGLETREE SIMULATOR 2.1 Overview

EagleTree is an open source simulator (available at https://github.com/ClydeProjects/EagleTree). It does not only simulate an SSD but also (a subset of) the OS and applications utilizing it. The first advantage of this approach is that EagleTree is an entire system operating in virtual time. It is thus possible to conduct large and complex design-space explorations, involving hundreds of experiments, in a tractable way. A second advantage is that we can experiment with cross-layer designs.

Figure 1: EagleTree architecture and design space exploration

EagleTree is composed of four layers, from the bottom-up: the hardware, the SSD controller, the operating system, and the application layers (see Figure 1). With EagleTree we can study questions related to an individual layer, or we can conduct crosslayer studies. For example:

At the hardware layer: Where are the hardware bottlenecks? What is the impact of the underlying flash chips characteristics? How should we use advanced commands (e.g. copybacks, pipelining), and what trade-offs is their usage subject to?

At the SSD controller layer: What is the impact of the mapping strategy? When should we trigger garbage collection or wear leveling? How should we schedule application reads and writes (1) relative to each other and (2) relative to internal IOs? What is the impact of different interleaving strategies? What is the best usage for RAM or for battery-backed RAM? What is the impact of applications' IOs queue size?

At the OS scheduler level: What is the best scheduling strategy (e.g., FIFO, CFQ, priorities)? How many outstanding IOs should be submitted to the SSD? What useful meta-information may be transmitted to the SSD in case of cross layer optimizations?

At the application layer: How can an algorithm leverage SSD internal parallelism? How should we submit synchronous and asynchronous IOs?

Finally, considering a cross-layer approach: How can we leverage an open interface to cooperate with the OS and the SSD? What is the impact on application algorithms (e.g., indexing, extent-based space allocation, join algorithms, external sorting algorithms)? How should work be divided across layers (e.g., scheduling at OS level or SSD level or both)?

The obvious drawback of our approach is that the entire system is simulated, and thus, far away from a deployable system. Our goal is that EagleTree should be useful to answer some basic questions, but we recognize that the answers that we obtain will have to be confronted to actual system components. We are exploring an approach where the fully simulated environment of EagleTree is complemented by (1) a software-based SSD simulator (based on the hardware and controller layers of EagleTree) connected to a real OS and real applications --the problem is then to store the manipulated data in RAM (otherwise it may slowdown the execution), and to map the virtual IO time to real time, and (2) a hardware-based SSD simulator using the OpenSSD hardware [START_REF]EagleTree game window[END_REF]. We think that these approaches, more realistic but less agile [START_REF] Saxena | Getting Real: Lessons in Transitioning Research Simulations into Hardware Systems[END_REF], are best suited in a second phase, once EagleTree has been used to get a basic understanding of the key trade-offs.

Design and Implementation

Each EagleTree layer consists of configurable parameters and customizable policies. From the bottom-up:

Hardware: EagleTree allows users to set up every hardware parameter of the simulated SSD: basic flash chip timings (i.e., to send a command, transfer data on a channel, read, write or erase), flash chip and SSD geometry. Moreover, EagleTree allows specifying the flash chip type (i.e., SLC or MLC) and its support for advanced commands. Finally, EagleTree includes a memory manager used to track the amount of RAM and battery-backed RAM used for the controller's metadata and IO buffers. All these parameters are variables that can be set, viewed and updated with ease. Predefined configurations are provided based on existing SSDs and flash chip datasheets.

SSD Controller:

The SSD controller is responsible for orchestrating mapping, garbage-collection, wear leveling modules and scheduling.

The mapping scheme supports the virtualization of the address space, mapping logical addresses onto physical ones. For reads, the mapping scheme must look up the physical address corresponding to the logical address of the incoming IO. For writes, the mapping scheme imposes constraints on which physical address a given IO might be bound to. As a result, the mapping scheme potentially restricts the scheduling policy. For now, we have considered the most flexible schemes i.e., pagebased mappings: the well-known DFTL [START_REF] Gupta | DFTL: a flash translation layer employing demand-based selective caching of page-level address mappings[END_REF] and a page-based mapping scheme where the entire mapping is kept in RAM.

For page-mapping FTLs, garbage-collection should fulfill the following requirements. First, it is desirable to wait as long as possible before performing garbage-collection. Doing so maximizes the number of invalid pages across the SSD, thereby ensuring that victim blocks have few live pages. On the other hand, GC must not occur so late that the FTL actually runs out of available space for incoming writes. Second, it is desirable to maintain free space on every LUN to maximize the flexibility for where writes can be made. The default GC module strives to fulfill these goals by triggering GC so that a given number of blocks (GC Greediness parameter) are always free on each LUN. The default wear leveling module keeps track of (1) the ages of all blocks, (2) a timestamp for each block marking the time in which it was last erased, (3) the average length of time it takes a block to be erased, and (4) the current time. Using this information, the WL module can identify particularly young blocks that have not been erased for a very long time, and can target them for static wear leveling. Moreover, EagleTree contains dynamic wear-leveling strategies that allow maintaining several free blocks of different ages in each LUN. The overall goal is to associate hot data with young blocks and cold data with old blocks. Temperature detection for pages can be done by (1) assuming the pages migrated in static wear-leveling are cold, and everything else is hot, or (2) using a temperature detection mechanism for each page such as the one described in [START_REF] Park | Hot Data Identification for Flash-based Storage Systems Using Multiple Bloom Filters[END_REF], which we have implemented, or (3) using information about the temperature of data coming through an open interface from the application.

EagleTree supports parallelism among channels and operation interleaving within a channel thus making the IO scheduler a central part in EagleTree. With regards to IO scheduling, the problem is roughly the following. Given the state of the flash chip array and a queue of pending IOs from various sources (e.g. application, garbage-collection, mapping, etc.), of various types (e.g. read, write, erase, copy-back), and that have been waiting in the queue for different lengths of time, which IO should be executed next and where? EagleTree provides a modular framework for exploring scheduling policies. For example, it is possible to:

• Control when and where to trigger internal operations, such as garbage-collection, wear-leveling; • Implement priority schemes that differentiate among IOs based on their sources, types, and waiting times; • Set deadlines for different IO types, and control the way in which overdue IOs are handled as relative to other IOs;

• Record and exploit information about logical address patterns (e.g. sequential vs random, hot vs cold); • Control the aggressiveness of interleaving and copy-back operations;

Obviously, these defaults strategies can be overwritten or complemented. In addition, other modules can be added to the SSD controller, e.g., a write-buffering module that uses battery-backed RAM to temporarily store data before it is written on flash pages.

Open Interface: EagleTree takes a departure from the traditional block device interface by basing communication between the OS and the SSD on an extensible messaging framework that allows the operating system and SSD to communicate as peers. Users are able to create new types of messages between the SSD and the OS conveying any amount of information or instructions. Using this framework, it is possible to build arbitrarily complex communication protocols between the SSD and OS. Investigating what we stand to gain from more explicit communication between the OS and SSD is the subject of our ongoing work. For now, we hope to inspire and excite the reader by sketching a few examples and their potential impact on performance.

• Priorities: the OS can communicate to the SSD the priority of an IO. The SSD can take this into account by offering the IO special treatment in terms of scheduling. • Update-locality: the OS can inform the SSD which pages share update-locality. The SSD can then write these pages so as to minimize subsequent garbage-collection. • Temperatures: the OS can inform the SSD whether the page being written is likely to be updated soon. The SSD can use this to benefit wear-leveling and garbage-collection efficiency.

OS Scheduler:

The Operating System manages IO requests incoming from multiple simulated concurrent threads. It maintains a pool of pending IOs from each thread and decides, based on a customizable scheduling policy, which IOs to issue next to the SSD. This policy can take into account the IO type (e.g. read/write/trim), its priority, the dispatching thread, etc. The default scheduling strategy is FIFO. Once the SSD has completed executing an IO, it interrupts and notifies the OS. The OS then activates the thread that dispatched the IO. The thread can respond by issuing additional IOs.

Threads:

The Thread layer is a programming framework that gives users absolute control over the workload. Users are able to extend an abstract thread class by providing a definition for two methods: init() and call_back(). The init() method is called by the OS when the thread is initialized, and the call_back() method is triggered every time an IO originating from that thread completes. Within each of these methods, it is possible to issue any number of messages to the OS. A user can create any number of threads with arbitrarily complex behaviors. For example, we have implemented threads simulating the behavior of a file system as well as a thread that follows the IO pattern of Grace hash Join.

Running Experiments

EagleTree contains an experimental suite API, which consists of experiment templates. An experiment template takes (1) an SSD parameter or policy (2) a strategy for how to vary it in an experiment, and (3) a workload definition. It runs an experiment and produces a comprehensive amount of visual statistical output. This includes graphs showing how performance metrics (e.g., throughput, latency, latency variability) evolved with respect to the given parameter or policy, as well as graphs showing how various metrics evolved across time in the experiments, and massive visual traces showing exactly how every IO was handled throughout the simulator components. EagleTree contains other useful features to allow conducting controlled, repeatable experiments. For example, it is possible to attach statistics gathering objects to an individual thread to measure its performance. It is also possible to create dependencies among threads. This latter feature is particularly useful for bringing the SSD to a well-defined state. This can typically be done by starting thread(s) that write over the entire logical address space sequentially and/or randomly [START_REF] Bouganim | uFLIP: Understanding flash I/O patterns[END_REF] and then triggering the experiment workload once the preparation threads finished, and measuring performance only for the experiment workload.

DEMONSTRATION SCENARIO

The purpose of the demonstration is to show through experiments, that: (1) the design space of SSD-based algorithms is vast; that (2) opening the interface is interesting but still increases that design space; that (3) EagleTree can help to explore that design space; and finally that (4) interesting solutions are sometime counter-intuitive. The demonstration will be centered on scheduling issues and will roughly follow the following outline.

Layered Tour: This part uses the demonstration GUI shown on Figure 2. We will first introduce the approach and EagleTree architecture. Then, attendees will choose configuration parameters (e.g. hardware setup, controller and OS policies, application workload) then run the simulator, and observe live results in terms of numerical performance metrics, traces, and graphical outputs (lower part of the GUI and right part of Figure 2). We will also show pre-computed graphs using longer experiments to show the full power of EagleTree's experimental suite. We will pay particular attention to the impact of scheduling policies on performance, and explain why prioritizing between application reads and writes is not always easy. We will introduce the challenge of scheduling internal operations as non-obtrusively as possible.

Open Interface Appetizers: In a second phase, we will open the block device interface. By unlocking the red lock shown on the figure, we will enable the attendee to choose from a range of possible extensions to the block device interface (e.g., explicit communication of locality, data temperature), guiding the attendee such that large variation in performance can be observed. Game: We will close the demonstration proposing a game to the attendee. Using the second GUI (see Figure 3) and knowledge acquired during the presentation, the user will have to guess the optimal combination of scheduling policies given a subset of the SSD scheduling design space. The attendee's objective will be to maximize throughput for a given workload while balancing mean latency and latency variability between different types of IOs. The user who draws nearest to the optimal configuration at each demo session will win an EagleTree Tshirt.

Figure 2 :

 2 Figure 2: EagleTree demonstration main window (on the left) and example of code and results

Flash chips are decomposed into logical units, denoted LUNs that constitute the minimum granularity of parallelism. This notion of LUN, introduced in the Open Nand Flash Interface (ONFI) standard, abstracts away the notions of packages, chips and dies.Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Articles from this volume were invited to present their results at The 39th International Conference on Very Large Data Bases, August

26th -30th 2013, Riva del Garda, Trento, Italy. Proceedings of the VLDB Endowment, Vol. 6, No. 12 Copyright 2013 VLDB Endowment 2150-8097/13/10... $ 10.00.


Michel Kieffer
email: kieffer@lss.supelec.fr

Eric Walter
email: walter@lss.supelec.fr

Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS

Keywords: Confidence regions, Interval analysis, Nonlinear system identification

In parameter estimation, it is often desirable to supplement the estimates with an assessment of their quality. A new family of methods proposed by Campi et al. for this purpose is particularly attractive, as it makes it possible to obtain exact, nonasymptotic confidence regions under mild assumptions on the noise distribution. A bottleneck of this approach, however, is the numerical characterization of these confidence regions. So far, it has been carried out by gridding, which provides no guarantee as to its results and is only applicable to low dimensional spaces. This paper shows how interval analysis can contribute to removing this bottleneck.

Introduction

When a vector p of parameters of some approximate mathematical model is estimated from a noisy data vector y, this is usually via the minimization of some cost function J (p), for instance J (p) = yy m (p) 2 [START_REF] Campi | Non-asymtotic confidence sets for the parameters of linear transfer functions[END_REF] , where y m (p) is the vector of model outputs, assumed here to be a deterministic function of p and • 2 is a (possibly weighted) ℓ 2 norm. Then p = arg min p J (p). Even if y and y m ( p) are reassuringly similar, it would be naive to consider p as the final answer to the estimation problem. One should instead attempt to attach some quality tag to p by assessing the reliability of the numerical values thus obtained. A key issue is drawing conclusions that are as little prejudiced as possible, and the approaches recently proposed by Campi et al. for this purpose [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF][START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF] are particularly attractive. A difficulty with these approaches, however, is the numerical characterization of these confidence regions. So far, it has been carried out by gridding, which provides no guarantee as to its results and is only applicable to low dimensional spaces.

The aim of this paper is to show how interval analysis can contribute removing this bottleneck, by provid-ing guaranteed results, as well as results in high dimensional spaces. The approaches Leave-out Sign-dominant Correlated Regions (LSCR) [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF] and Sign-Perturbed Sums (SPS) [START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF] recently proposed by Campi et al. are briefly recalled in Section 2. Section 3 shows how interval analysis can be used to characterize the exact confidence regions defined by LSCR and SPS in a guaranteed way. Examples are in Section 4 and conclusions in Section 5.

LSCR and SPS

The most striking feature of LSCR and SPS is that they avoid a large number of the usual assumptions about the noise corrupting the data yet provide an exact characterization of parameter uncertainty in non-asymptotic conditions. Both require only the noise samples to be independently distributed with distributions symmetric with respect to zero. For LSCR, arbitrary noise can even be dealt with if the regressors are random, independently identically and symmetrically distributed, and independent of the noise [START_REF] Campi | Non-asymtotic confidence sets for the parameters of linear transfer functions[END_REF][START_REF] Granichin | The nonasymptotic confidence set for parameters of a linear control object under an arbitrary external disturbance[END_REF].

LSCR [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF] looks for a region Θ to which the parameter vector p * of the model assumed to have generated the data belongs with a specified probability. Let ε t (p) be a prediction error, such that ε t (p * ) is a realization of the noise corrupting the data at time t. It may, for instance, be the difference between some output measurement y t , and the corresponding model output y m t (p). The procedure for computing one such confidence region is as follows:

(1) Select two integers r 0 and q 0. (2) For t = 1 + r, . . . , k + r = n, compute

c ε t-r,r (p) = ε t-r (p) ε t (p) , (1)
(3) Compute s ε i,r (p) = k∈Ii c ε k,r (p) , i = 1, . . . , m, where I i is a subset of a set I of indexes and the collection G of these subsets I i , i = 1, . . . , m, forms a group under symmetric difference, i.e., (I i ∪ I j ) -

(I i ∩ I j ) ∈ G.
(4) Find the set Θ ε r,q such that at least q of the s ε i,r (p) are larger than 0 and at least q are smaller than 0.

The probability that p * belongs to Θ ε r,q is 1-2q/m. The shape and size of Θ ε r,q depend not only on the values given to q and r but also on the group G and its number of elements m. A procedure for generating a group of appropriate size is suggested in [START_REF] Gordon | Completely separating groups in subsampling[END_REF].

The set Θ ε r,q may be defined more formally as Θ ε r,q = Θ ε,1 r,q ∩ Θ ε,2 r,q , with, for j = 1, 2,

Θ ε,j r,q = p ∈ P such that m i=1 τ ε,j i (p) q , (2)
where P is the prior domain for p and where τ ε,j i (p) = 1 if (-1) j s ε i,r (p) 0 and τ ε,j i (p) = 0 else. The set Θ ε,1 r,q contains all values of p ∈ P such that at least q of the functions s ε i,r (p) are smaller than 0, whereas Θ ε,2 r,q contains all values of p ∈ P such that at least q of the functions s ε i,r (p) are larger than 0.

When the model studied is driven by an input u t , one may obtain a similar confidence region by substituting c u t-s,s (p) = u t-s (p) ε t (p) for (1) in the procedure above, thus replacing autocorrelations by intercorrelations. One then computes a set Θ u s,q , again such that Pr p * ∈ Θ u s,q = 1 -2q/m.

The fact that the set Θ ε r,q (or Θ u r,q ) obtained by this approach is exact does not mean that its volume is minimal, and the resulting confidence region may turn out to be much too large to be useful. One may then intersect several such regions. For a given value of q and m, assume that n ε confidence regions Θ ε r,q and n u confidence regions Θ u s,q have been obtained for n ε values of r and n u values of s. The probability that p * belongs to the intersection Θ of these (n ε + n u ) regions then satisfies Pr (p * ∈ Θ) 1 -(n ε + n u )2q/m. The price to be paid for taking the intersection of several confidence regions is that the probability that p * belongs to the resulting confidence region is no longer known exactly, as only a lower bound for this probability is available. SPS [START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF] also provides a confidence region to which p * belongs with a specified probability, by exploiting the symmetry of the noise distribution and the independence between noise samples. It is designed for linear regression, where y t = ϕ T t p * + w t , t = 1, . . . , n, with ϕ t a known regression vector that does not depend on the unknown parameters. It computes an exact confidence region for p * around the least-squares estimate p, which is the solution to the normal equations n t=1 ϕ t y tϕ T t p = 0. For a generic p, define

s 0 (p) = n t=1 ϕ t y t -ϕ T t p , (3)
and the sign-perturbed sums

s i (p) = n t=1 α i,t ϕ t y t -ϕ T t p , (4)
where i = 1, . . . , m -1 and α i,t are independent and identically distributed (i.i.d.) random signs, so α i,t = ±1 with equal probability, and

z i (p) = s i (p) 2 2 , i = 0, . . . , m -1. (5)
A confidence region Σ q is obtained as the set of all values of p such that z 0 (p) is not among the q largest values of (z i (p))

m-1 i=0 . In [START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF], it has been shown that p * belongs to Σ q with exact probability 1q/m. Σ q may be defined more formally as

Σ q = p ∈ P such that m-1 i=1 τ i (p) q (6)
where

τ i (p) = 1 if z i (p) -z 0 (p) > 0 and τ i (p) = 0
else. This is justified by the fact that if m-1 i=1 τ i (p) q, then one has τ i (p) = 1 for at least q out of the m -1 functions τ i (p). As a consequence, there are at least q functions z i (p) such that z i (p) > z 0 (p) and z 0 (p) is not among the q largest values of (z i (p)) m-1 i=0 .

Guaranteed characterization

In LSCR and SPS, one has to characterize a set or an intersection of sets defined as

Ψ q = p ∈ P such that m i=1 τ i (p) q , (7)
where τ i (p) is some indicator function τ i (p) = 1 if f i (p) 0 and τ i (p) = 0 else, with f i (p) depending on the model structure, the measurements, and the parameter vector p.

Set inversion

Characterizing Ψ q may be alternatively formulated as a set-inversion [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] problem

Ψ q = P ∩ τ -1 ([q, m]), with τ (p) = m i=1 τ i (p) , (8)
which may be efficiently solved via interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF][START_REF] Moore | Interval Analysis[END_REF] using the SIVIA algorithm [START_REF] Jaulin | Applied Interval Analysis[END_REF]. For that purpose, inclusion functions for the τ i 's and consequently for the f i 's are required. SIVIA recursively partitions P into boxes (vectors of intervals) proved to belong to Ψ q , boxes proved to have no intersection with Ψ q , and undetermined boxes for which no conclusion can be obtained. SIVIA bisects undetermined boxes until their width is less than some precision parameter ε.

Contractors for guaranteed characterization

Indetermination often results from range overestimation by inclusion functions. As a consequence, boxes have to be bisected by SIVIA many times to allow one to conclude on the position of the resulting boxes with respect to Ψ q . This may entail an intractable computational complexity, even for a moderate dimension of p.

Contractors [START_REF] Jaulin | Applied Interval Analysis[END_REF] partly address this issue. Consider a setinversion problem where one has to characterize the set

X = [x] ∩ f -1 (Y) , (9)
with f :

D ⊂ R n → R m , Y ⊂ R m ,
C f ,Y ([x]) ⊂ [x] such that [x] ∩ X = C f ,Y ([x]) ∩ X, so no part of X in [x]
is lost. It allows parts of the candidate box [x] that do not belong to X to be eliminated, without bisection.

In the problems considered here, the role of x is taken by p. The fact that the function τ introduced in ( 8) is not differentiable forbids the use of most classic contractors [START_REF] Jaulin | Applied Interval Analysis[END_REF]. The new contractor proposed here assumes that the functions f i are differentiable and is implemented in two steps. First, a set of m possibly overlapping subboxes of [p] is built, trying to remove all values of p from [p] such that f i (p) < 0, i = 1, . . . , m, see Sections 3.2.1 and 3.2.2. Second, the union of all non-empty intersections of at least q of these boxes is computed to get a possibly contracted box, see Section 3.2.3.

Contractor for LSCR and SPS

The first step uses the centered form of f i , which, for some m ∈ [p], may be written as

[f i,c ] ([p]) = f i (m) + ([p] -m) T [g i ] ([p]) (10) = f i (m) + np j=1 ([p j ] -m j ) [g i,j ] ([p]) , (11)
where g i is the gradient of f i and [g i ] ([p]) is the natural inclusion function for g i , see [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Using [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF], we build a contractor C fi,[0,∞[ for the set of all values of p ∈ [p] such that f i (p) 0, as follows. ). Some of them may be empty, in which case, [p ′ i ] = ∅ indicates that there is no p ∈ [p] such that f i (p) 0.

With the k-th component [p k ] of [p], when 0 / ∈ [g i,k ] ([p]), C fi,[0,∞[ associates the contracted interval p ′ i,k = [p k ] ∩ ([f i,c ] ([p]) ∩ [0, ∞[) -f i (m) - np j=1,j =k ([p j ] -m j ) [g i,j ] ([p]) / [g i,k ] ([p]) + m k . ( 12
) When 0 ∈ [g i,k ] ([p]), C fi,[0,∞[ leaves [p k ] unchanged, i.e., [p ′ i,k ] = [p k ].

Contractor for SPS

We take advantage of the fact that the functions s i (p), i = 0, . . . , m are affine in p to reduce the number of occurrences of p in their formal expression, and thus the pessimism of the corresponding inclusion functions. Equation ( 3) is rewritten as

s 0 (p) = b 0 -A 0 p (13) with b 0 = n t=1 ϕ t y t and A 0 = n t=1 ϕ t ϕ T t .
Similarly, (4) is rewritten as

s i (p) = b i -A i p (14) with b i = n t=1 α i,t ϕ t y t and A i = n t=1 α i,t ϕ t ϕ T t .
Equations ( 5), ( 13), ( 14), and the fact that the A i 's are symmetric imply that

z i (p) -z 0 (p) = p T A 2 i -A 2 0 p -2 b T i A i -b T 0 A 0 p + b T i b i -b T 0 b 0 . (15)
The matrices A 2 i -A 2 0 are symmetric and may thus be diagonalized as

A 2 i -A 2 0 = U T i D i U i
, where U i is an orthonormal matrix (i.e., such that U T i = U -1 i ), and D i = diag d i,1 , . . . , d i,np is a diagonal matrix. With the change of variables π = U i p, (15) becomes

z i (p) -z 0 (p) = π T D i π -2β T i π + γ i , (16)
where

β T i = b T i A i -b T 0 A 0 U T i and γ i = b T i b i - b T 0 b 0 .
Then, provided that d i,j = 0 for j = 1, . . . , n p , (16) can be rewritten as

z i (p) -z 0 (p) = np j=1 d i,j π j - β i,j d i,j 2
+ γ i - np j=1 β 2 i,j d i,j . (17) Let [π] = U i [p]. A contractor for [π j ] is obtained from (17) as follows π ′ j = [π j ] ∩ ± 1 d i,j ([z i ] ([p]) -[z 0 ] ([p])) ∩ [0, ∞[ - np k=1 k =j d i,k [π k ] - β i,k d i,k 2 -γ j + np k=1 β 2 i,k d i,k 1 2 + β i,j d i,j . (18)
The contractor

C zi-z0,[0,∞[ for [p] is then such that [p ′ i ] = C zi-z0,[0,∞[ ([p]) = [p] ∩ U T i [π ′ ] . ( 19
)
When n is large enough, and provided that the ϕ t 's have been well designed, it is very unlikely that A 2 i -A 2 0 is rank deficient. If this occurs, (17) and (18) have to be rewritten distinguishing the zero and nonzero d i,j 's.

Proposition 1 Provided that d i,j = 0 for j = 1, . . . , n p , for all [p ′ i ], i = 1, . . . , m -1, built using (18) and (19), [p ′ i ] ⊂ [p] and [p ′ i ] ∩ (z i -z 0 ) -1 ([0, ∞[) = [p] ∩ (z i -z 0 ) -1 ([0, ∞[) . (20) Proof [p ′ i ] ⊂ [p] is true by construction. It remains to be proved that [p] ∩ (z i -z 0 ) -1 ([0, ∞[) ⊂ [p ′ i ] ∩ (z i -z 0 ) -1 ([0, ∞[). Let p 0 ∈ [p] ∩ (z i -z 0 ) -1 ([0, ∞[) and π 0 = U i p 0 ∈ [π]. To prove that p 0 ∈ [p ′ i ] ∩ (z i -z 0 ) -1 ([0, ∞[), it suffices to prove that π 0 ∈ [π ′ ]. By definition of p 0 , one has np j=1 d i,j π 0 j - β i,j d i,j 2 + γ i - np j=1 β 2 i,j d i,j = z i p 0 -z 0 p 0 0. ( 21
)
Since π 0 ∈ [π], after some manipulations of (21), one gets for j = 1, . . . , n p

π 0 j ∈ [π j ] ∩ ± 1 d i,j ([z i ] ([p]) -[z 0 ] ([p])) ∩ [0, ∞[ - np k=1 k =j d i,k [π k ] - β i,k d i,k 2 -γ j + np k=1 β 2 i,k d i,k 1 2 + β i,j d i,j . Thus π 0 j ∈ π ′ j ,

Building a q-relaxed intersection

During the second step, the contractor builds a box [p ′ ] enclosing the q-relaxed intersection [6, 7]

P = q j∈{1,...,m-1} p ′ j = J⊂[1,...,m-1] card(J) q j∈J p ′ j , (22)
of the boxes in L = {[p ′ 1 ] , . . . , [p ′ m ]}, i.e., the union of all intersections of at least q boxes in L, such that [p ′ ] satisfies

P ⊂ [p ′ ] ⊂ [p] . ( 23
)
Proposition 2 For any box [p ′ ] satisfying (23), one has

Ψ q ∩ [p ′ ] = Ψ q ∩ [p],
with Ψ q as defined in [START_REF] Jaulin | Set-membership localization with probabilistic errors[END_REF].

Proof Assume that there exists p 0 ∈ [p] such that

p 0 ∈ Ψ q ∩ [p] but p 0 / ∈ Ψ q ∩ [p ′ ]. Since p 0 ∈ Ψ q ∩ [p], we have p 0 ∈ Ψ q . According to (7), m i=1 τ i (p 0 )
q. There are thus at least q functions τ i such that τ i (p 0 ) 1. Assume, without loss of generality, that τ

1 (p 0 ) 1, . . . , τ q (p 0 ) 1. Since τ i (p 0 ) 1, i = 1, . . . , q, by definition of C fi,[0,∞[ , one has p 0 ∈ [p ′ i ], i = 1, . . . , q and p 0 ∈ i=1,...,q [p ′ i ]. By definition of P and [p ′ ], p 0 ∈ i=1,...,q [p ′ i ] ⊂ P ⊂ [p ′ ],
which contradicts the initial assumption. ⋄

Evaluating the q-relaxed intersection

Consider m scalar intervals [p 1 ] , . . . , [p m ]. Algorithm 1, which formalizes a computation carried out on an example in [START_REF] Jaulin | Set-membership localization with probabilistic errors[END_REF], builds the smallest interval containing the

1 [p] = ∅;
2 Reindex the boxes [p i ] in such a way that

p 1 p 2 • • • p m ; 3 For i = q to m 4 if m j=1 p i ∈ [p j ]
q then p = p i ; break;

5 Reindex the boxes [p i ] in such a way that

p 1 p 2 • • • p m 6 For i = q to m 7 if m j=1 (p i ∈ [p j ]) q then p = p i ; break; Algorithm 1. [p] = q relaxed intersection ([p1] , . . . , [pm])
union of all intersections of q intervals with a complexity O (m log m). This is the smallest interval containing P as defined by ( 22

Examples

All computations were carried out with Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF], the interval-analysis toolbox in Matlab, on an Intel Core i7 at 3.7 GHz with 8 GB RAM. The computations required by SIVIA and the q-relaxed intersection, which form the major part of the computational burden, could be speeded up considerably in C++. The source code is available at http://www.l2s.supelec.fr/perso/kieffer-0.

y t = α (p * ) (exp (λ 1 (p * ) t) -exp (λ 2 (p * ) t)) + w t , where α (p) = k 21 / (k 01 -k 12 + k 21 ) 2 + 4k 12 k 21 , λ 1,2 (p) = -1 2 ((k 01 + k 12 + k 21 ) ± α (p) /k 21
) and the w t 's are realizations of i.i.d. N 0, σ 2 variables, for t = 0, T, . . . , (n -1) T . The variance of the measurement noise is σ 2 = 10 Here, the set Θ ε r,q has been characterized using LSCR with r = 1 and q = 3, which corresponds to a 90 % confidence region, see Figure 2 Table 1 shows the evolution of the computing times for the gridding approach and SIVIA. The increase is quadratic with 1/ε for gridding. It is slower with SIVIA, since only undetermined boxes are further bisected when ε decreases. An additional advantage of SIVIA is that the results provided are guaranteed. T .

To evaluate the performance of the proposed technique for a large number of parameters, FIR models with n a = 20 random parameters in [-2, 2] na are generated. Then, n = 512, 1024, 2048, 4096, and 8192 noise-free data points are generated applying a random i.i.d. sequence u t of ±1, which is the D-optimal input under the constraint that the input has to remain in [-1, 1]. White Laplacian noise is then added to these data. The standard deviation of the noise is set up to get SNRs from 10 dB to 40 dB.

Our aim is to characterize a 95 % confidence region with SPS. A possible choice is m = 255 and q = 13. The initial search box in parameter space is taken as P = -10 4 , 10 4 20 . Getting accurate inner and outer approximations using unions of non-overlapping boxes is hopeless with n a = 20 parameters, because of the curse of dimensionality. Our aim is instead to provide a guaranteed outer-approximation of the confidence region 1 . For that purpose, the contractor of Section 3.2.2 is applied once to P (iterations are useless). A box estimate is obtained in 5 s in average, whatever the value of n. This is because the computational complexity of the contractor of Section 3.2.2 is mainly determined by m and n a . On a log-log scale, maximum width seems linear in the SNR and in the number of samples. In the latter case, the slope is about -1/2, consistent with what is observed when maximum-likelihood estimation is carried out assuming an additive Gaussian noise, although this hypothesis on the noise is neither true nor assumed here. 1 The computed least-square estimate belongs to the confidence region, as showed in [START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF]. It thus forms a (point) inner approximation.

Conclusions

Interval analysis provides tools to evaluate guaranteed inner and outer-approximations of non-asymptotic confidence regions defined by LSCR and SPS. Symbolic manipulations are particularly useful to design more efficient contractors and struggle against the curse of dimensionality.

  Due to the pessimism of centered forms on large boxes, the contractor C fi,[0,∞[ becomes efficient only when [p] is small enough. Considering the m functions f i and applying all the contractors C fi,[0,∞[ , i = 1, . . . , m, to [p], one obtains m possibly contracted boxes [p ′ 1 ] = C f1,[0,∞[ ([p]) , . . . , [p ′ m ] = C fm,[0,∞[ ([p]

  which completes the proof. ⋄ Applying the contractors C zi-z0,[0,∞[ , i = 1, . . . , m, to [p], as in Section 3.2.1, one obtains m possibly contracted boxes C z1-z0,[0,∞[ ([p]) , . . . , C zm-z0,[0,∞[ ([p]).

  ) in the scalar case. At Steps 4 and 7 of Algorithm 1, (p ∈ [p j ]) = 1 if p ∈ [p j ] and (p ∈ [p j ]) = 0 otherwise. The extension to boxes is obtained by applying Algorithm 1 componentwise.

4. 1

 1 Nonlinear estimation Consider the two-compartment model described by Figure 1. Only the content of the second compartment is observed. The model parameters to be estimated are p = (k 01 , k 12 , k 21 ) T . The data are generated by this model for p * = (1, 0.25, 0.5) T . They satisfy

- 4 .

 4 The sampling period is T = 0.2 s, and n = 64. To facilitate illustration, only k 01 et k 12 are estimated. The value k * 21 of k 21 is assumed known. The prediction errors are ε t (p) = y ty m t (p) , with y m t (p) = α (p) (exp (λ 1 (p) t)exp (λ 2 (p) t)), for t = 0, T, . . . , (n -1) T .

  Fig. 1. Two-compartment model 90 % confidence region LSCR -SIVIA

Fig. 2 .

 2 Fig. 2. Gridding (left) and paving (right) of search space obtained using LSCR in the two-compartment model case for the characterization of the confidence region Θ ε r,q ; the two bottom subfigures are obtained by zooming in on one of the connected components of Θ ε r,q

  . The initial search set in parameter space is P = [0, 5]×[0, 5]×[0.5, 0.5]. The top left part of Figure 2 represents the result obtained in 284 s by gridding as in [1] with a grid step-size ε = 0.0025. The top right part of Figure 2 has been obtained by SIVIA with ε = 0.0025 in 175 s. The top right part of Figure 2 proves that the confidence region consists of two disconnected subsets, a consequence of the lack of global identifiability of the model (exchanging the values of k 01 and k 12 does not change the model output).

  Figure 2 (bottom part) zooms on one of the two confidence subsets, which turns out to contain the actual value of the unknown parameters, although this is not guaranteed, of course.

  Figure 3 represents the width of the largest component of the resulting outer box as a function of the SNR and of the number of data points.

Fig. 3 .

 3 Fig. 3. Width of the largest component of the outer box resulting from a single application of the contractor of Section 3.2.2 as a function of the number of data points and of the SNR.

  , . . . , a na-1 ) T and u t = 0 for t 0. For t = 1, . . . , n, the w t 's are i.i.d. noise samples. In linear regression form, (24) becomes y t = ϕ T t p * + w t with ϕ T t = (u t , . . . , u t-na+1 ) and p * = a * 0 , . . . , a *

	ε	0.1	0.025 0.01 0.0025 0.001
	Gridding (s) 0.26	2.9	18	284	1750
	SIVIA (s)	15	57	93	175	400
	Table 1
	Computing times for the example of Section 4.1
	4.2 FIR model
	Consider now the system
		y t = y m t (p * ) + w t ,		(24)
	with the FIR model y m t (p) =	na-1 i=0 a i u t-i , where
	p = (a 0 na-1

This work has been partly supported by the ANR CPP.


Elsa Dupraz

Aline Roumy

Michel Kieffer

Source Coding with Side Information at the Decoder and Uncertain Knowledge of the Correlation

This paper considers the problem of lossless source coding with side information at the decoder, when the correlation model between the source and the side information is uncertain. Four parametrized models representing the correlation between the source and the side information are introduced. The uncertainty on the correlation appears through the lack of knowledge on the value of the parameters.

For each model, we propose a practical coding scheme based on non-binary Low Density Parity Check Codes and able to deal with the parameter uncertainty. At the encoder, the choice of the coding rate results from an information theoretical analysis. Then we propose decoding algorithms that jointly estimate the source vector and the parameters. As the proposed decoder is based on the Expectation-Maximization algorithm, that is very sensitive to initialization, we also propose a method to produce first a coarse estimate of the parameters.

I. INTRODUCTION

The problem of lossless source coding with side information at the decoder has been well investigated when the correlation model between the source X and the side information (SI) Y is perfectly known. Slepian and Wolf showed that this case induces no loss in performance compared to the conditional setup, i.e., the setup where the side information is also known at the encoder [START_REF] Slepian | Noiseless coding of correlated information sources[END_REF]. Following this principle, several works, see, e.g., [START_REF] Puri | PRISM: A new robust video coding architecture based on distributed compression principles[END_REF], [START_REF] Stankovic | On code design for the Slepian-Wolf problem and lossless multiterminal networks[END_REF], [START_REF] Xiong | Distributed source coding for sensor networks[END_REF], propose practical coding schemes for the Slepian-Wolf (SW) problem. Most of them are based on channel codes [START_REF] Stankovic | Design of Slepian-Wolf codes by channel code partitioning[END_REF], and particularly Low Density Parity Check (LDPC) codes [START_REF] Liveris | Compression of binary sources with side information at the decoder using LDPC codes[END_REF], [START_REF] Matsuta | Universal Slepian-Wolf source codes using Low-Density Parity-Check matrices[END_REF]. This approach allows to leverage on many results on LDPC codes for the code construction and optimization [START_REF] Li | Density evolution for nonbinary LDPC codes under Gaussian approximation[END_REF], [START_REF] Richardson | The capacity of Low-Density Parity-Check codes under message-passing decoding[END_REF] even if there is a need to adapt the developed algorithms to the case of SW coding [START_REF] Chen | The equivalence between Slepian-Wolf coding and channel coding under density evolution[END_REF].

Nonetheless, most of these works assume perfect knowledge of the joint distribution P (X, Y ).

In [START_REF] Jalali | A universal scheme for Wyner-Ziv coding of discrete sources[END_REF], it is shown that the performance of the SW coding scheme remains the same if P (X) is unknown. Here we consider the case where the correlation channel P (Y |X) is uncertain because it is in general more difficult to obtain in practical situations. In this way, [START_REF] Sgarro | Source coding with side information at several decoders[END_REF] considers the case where P (Y |X) is given to the decoder but not perfectly known at the encoder. Here we assume that P (Y |X) is uncertain at both the encoder and the decoder. A usual solution to address this problem is to use a feedback channel [START_REF] Aaron | Wyner-Ziv coding of motion video[END_REF], [START_REF] Eckford | Rateless Slepian-Wolf Codes[END_REF], [START_REF] Varodayan | Rate-adaptive codes for distributed source coding[END_REF], or to allow interactions between the encoder and the decoder [START_REF] Yang | Interactive encoding and decoding for one way learning: Near lossless recovery with side information at the decoder[END_REF]. The advantage of the feedback channel is that the rate is adapted to the true characteristics of the source. However, a feedback channel can be difficult to implement in many practical situations such as sensor networks. Moreover, the feedback channel is in general used by the decoder to ask for additional packets to the encoder or to stop the transmission.

Each time a new packet is received, the decoder processes again to try to reconstruct the source from all the received packets. This can result in huge decoding delays.

When no feedback is allowed, several practical solutions based on LDPC codes and proposed February 1, 2013 DRAFT for channel coding may be adapted to the SW problem. When hard decoding is performed, as proposed in [START_REF] Liu | Weighted symbol-flipping decoding algorithm for nonbinary LDPC codes with flipping patterns[END_REF], [START_REF] Richardson | Design of capacity-approaching irregular Low-Density Parity-Check codes[END_REF] for channel coding, only symbol values are used, at the price of an important loss in performance. An alternative solution is the min-sum decoding algorithm proposed in [START_REF] Chen | Near optimum universal belief propagation based decoding of Low-Density Parity Check codes[END_REF], [START_REF] Savin | Min-Max decoding for non binary LDPC codes[END_REF] for channel coding, respectively for binary and non-binary sources. The min-sum algorithm uses soft information for decoding, but does not require the knowledge of the correlation channel. However, for the min-sum to be as efficient as the soft decoding algorithm, a coefficient is required for the initialization of the algorithm and needs to be chosen carefully.

Unfortunately, its value depends on the correlation channel.

In many applications, it is possible to restrict the correlation channel model to a given class (e.g., binary symmetric, Gaussian, etc.) due to the nature of the problem. Consequently, in this paper, we introduce four signal models. Each model assumes that the correlation channel belongs to a given class and is parametrized by some unknown parameter vector. For two of the models, the correlation channel between source symbols (X n , Y n ) is parametrized by an unknown parameter π n , varying from symbol to symbol. One of these two models assumes the knowledge of a prior distribution P Π (π n ) for π n . The case where no prior on π n is known corresponds to arbitrarily varying sources [START_REF] Ahlswede | Coloring hypergraphs: A new approach to multi-user source coding-1[END_REF], [START_REF] Berger | The source coding game[END_REF]. For the two other models, the correlation channel is parametrized by an unknown parameter θ, fixed for the sequence {(X n , Y n )} +∞ n=1 but allowed to vary from sequence to sequence. This corresponds to universal source coding [START_REF] Gallager | Source coding with side information and universal coding[END_REF].

The distinction between the two models is also in the knowledge of a prior for θ. The distinction between varying parameters π n and a fixed parameter θ has been proposed earlier in [START_REF] Piantanida | Outage behavior of discrete memoryless channels under channel estimation errors[END_REF] in the case of channel coding.

The coding scheme we propose is based on non-binary LDPC codes. Hard and min-sum LDPC decoding are not able to exploit the knowledge of the structure of the class. Therefore, the sumproduct LDPC decoding algorithm is considered. From an analysis of the performance bounds, we explain for each model how to choose the coding rate and the LDPC coding matrix. Then, we show that the classical sum-product LDPC decoding algorithm can be used for only one February 1, 2013 DRAFT model. For the three other models, we propose a decoding algorithm that jointly estimates the source vector and the parameters. As the method is based on the EM (Expectation Maximization) algorithm [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF], which is very sensitive to initialization, we also propose a method to obtain first coarse estimates of the values of the parameters.

The paper is organized as follows. In Section II, the four signal models we consider are described formally. Section III explains how to choose the coding rates and to design the LDPC coding matrices. Section IV proposes a decoding method adapted to each model. Finally, Section V presents simulation results.

II. SIGNAL MODEL

The source X to be compressed and the SI Y available at the decoder produce sequences of symbols {X n } +∞ n=1 and {Y n } +∞ n=1 , respectively. X and Y denote the source and SI discrete alphabets. Bold upper case letters, e.g., X N 1 = {X n } N n=1 , denote random vectors, whereas bold lower case letters, x N 1 = {x n } N n=1 , represent their realizations. When it is clear from the context that the distribution of a random variable X n does not depend on n, the index n is omitted.

The goal of this section is to model the uncertainty on the correlation channel P (Y |X). Each of the four proposed models consists of a family of parametric distributions. In every case, the source distribution P (X) is assumed perfectly known and does not depend on the uncertain parameters. The first two models allow parameter variations from symbol to symbol.

Definition 1. (DP-Source).

A Dynamic with Prior source (X, Y ), or DP-Source, produces a

sequence of independent symbols {(X n , Y n )} +∞ n=1 drawn ∀n from P (X n , Y n ) that belongs to a family of distributions {P (X, Y |Π = π) = P (X)P (Y |X, Π = π)} π∈Pπ parametrized

by a

The four models defined in this section were also introduced with different names in two papers [START_REF] Dupraz | Source coding with side information at the decoder: Models with uncertainty, performance Bounds, and practical coding schemes[END_REF], [START_REF] Dupraz | Practical coding scheme for universal source coding with side information at the decoder[END_REF], of the same authors.

M-Source was for DP-Source, WPM-Source for DwP-Source, P-Source for SP-Source, WP-Source for SwP-Source. The names were changed for the sake of clarity.

{(X n , Y n )} +∞ n=1 .
The SwP-source, completely determined by P θ and {P (X, Y |θ)} θ∈P θ , is stationary but nonergodic [START_REF] Gallager | Information theory and reliable communication[END_REF]Section 3.5]. The only difference between the SP-and SwP-Sources lies in the definition of θ (no distribution for θ is specified in the SwP-Model). Note that both the encoder and the decoder are aware of the model characteristics given in Definitions 1 to 4.

In the SW setup, the infimum of achievable rates for our models are given by 1) for the DP-Source [START_REF] Slepian | Noiseless coding of correlated information sources[END_REF],

R = H(X|Y ) (1)
where H(X|Y ) is calculated from P (X = x|Y = y) = π∈Pπ P (π), P (X = x|, Y = y, π).

2) for the DwP-Source [START_REF] Ahlswede | Coloring hypergraphs: A new approach to multi-user source coding-1[END_REF], R = sup

P (X,Y )∈Conv({Pπ(X,Y )} π∈Pπ ) H(X|Y ) (2)
where Conv({P π (X, Y )} π∈Pπ ) is the convex hull of the elements of {P π (X, Y )} π∈Pπ ,

3) for the SP-Source [START_REF] Han | Information-spectrum methods in information theory[END_REF]Theorem 7.3.4],

R = P Θ -ess. sup H(X|Y, Θ = θ), (3)
where P Θ -ess. sup is the essential sup (the sup on the support of the distribution) with respect to the prior distribution P Θ , 4) for the SwP-Source [START_REF] Csiszar | Linear codes for sources and source networks: Error exponents, universal coding[END_REF],

R = sup θ∈P θ H(X|Y, θ) . ( 4
)
February 1, 2013 DRAFT

We see that for the DwP-Model, the SP-Model and the SwP-Model, the infimum of achievable rates are given by worst cases defined on the set of values the parameters may take (SP-and SwP-Models), or on the convex hull of this set of values (DwP-Model).

III. ENCODING

The coding schemes we propose are based on LDPC codes for SW coding. As suggested by [START_REF] Liveris | Compression of binary sources with side information at the decoder using LDPC codes[END_REF], [START_REF] Matsuta | Universal Slepian-Wolf source codes using Low-Density Parity-Check matrices[END_REF], LDPC codes initially introduced for channel coding can also be used for SW coding, after adaptation of the coding process and the decoding algorithm. In channel coding, LDPC codes were proposed for binary-input channels [START_REF] Gallager | Low-Density Parity Check Codes[END_REF] and generalized to non-binary input channels in [START_REF] Davey | Low Density Parity Check codes over GF(q)[END_REF]. The adaptation to the SW setup is described in [START_REF] Liveris | Compression of binary sources with side information at the decoder using LDPC codes[END_REF] for the binary case. In this paper, we propose a generalization of this adaptation to the non-binary case. This section describes the encoding part and introduces the involved notations. Note that the encoding part is as in the binary case, except that, now, the encoding operations are performed in GF(q). There are more differences in the decoding part.

We assume that the source symbols X are discrete and belong to GF(q). The SW coding of a source vector x of length N is performed by producing a vector s = H T x of length M < N .

The matrix H is sparse, with non-zero coefficients uniformly distributed in GF(q)\{0}. In the following, ⊕, ⊖, ⊗, ⊘ are the usual operators in GF(q). In the bipartite graph representing the dependencies between the random variables of X and S, the entries of X are represented by Variable Nodes (VN) and the entries of S are represented by Check Nodes (CN). The set of CN connected to a VN n is denoted N (n) and the set of VN connected to a CN m is denoted

N (m).
The sparsity of H is determined by the VN degree distribution λ(x) = i≥2 λ i x i-1 and the CN degree distribution ρ(x) = i≥2 ρ i x i-1 with i≥2 λ i = 1 and i≥2 ρ i = 1. In SW coding, the rate r(λ, ρ) of a code is given by r(λ, ρ) = M N = i≥2 ρ i /i i≥2 λ i /i . In order to perform the encoding of a source vector X, one needs to choose properly the coding rate and to design the LDPC coding matrix, i.e., to impose good degree distributions February 1, 2013 DRAFT (λ(x), ρ(x)) [START_REF] Li | Density evolution for nonbinary LDPC codes under Gaussian approximation[END_REF], [START_REF] Richardson | Design of capacity-approaching irregular Low-Density Parity-Check codes[END_REF]. The performance analysis of Section II suggests the following approach.

For the DP-Source, the LDPC coding matrix is designed for the known distribution P (X|Y ).

For the three other models, the LDPC coding matrix is designed for the worst cases defined by ( 2)-( 4).

IV. DECODING ALGORITHM

This section introduces LDPC-based decoding algorithms capable of dealing with the uncertainty on the value of the parameters of the models. For the DP-Source, the decoding algorithm is the sum-product LDPC decoder adapted to SW coding. For the other sources, the LDPC decoding algorithm cannot be used directly because of the lack of knowledge on the parameters. We thus propose to jointly estimate the encoded source sequence X N 1 and the unknown parameters. This joint estimation is performed with an EM algorithm [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF]. A method producing a first coarse estimate of the parameters is also presented to properly initialize the EM algorithm.

A. DP-Source: Standard LDPC decoding

In [START_REF] Liveris | Compression of binary sources with side information at the decoder using LDPC codes[END_REF] the standard sum-product LDPC decoding algorithm has been adapted to SW coding of binary sources with perfect correlation channel knowledge. This section generalizes the adaptation of the decoding algorithm to non-binary SW coding. Indeed, in the SW case, one needs to take into account both the probability distribution of X and of the received codeword S. For the DP-Source, the conditional distribution is perfectly determined as

P (X n = k|Y n = y n ) = π∈Pπ P (π)P (X n = k|Y n = y n , π) . ( 5
)
The sum-product decoder performs an approximate Maximum A Posteriori (MAP) estimation of x from the received codeword s and the observed side information y. The messages exchanged in the dependency graph are vectors of length q. The initial messages for each VN n are denoted

m (0) (n, y n ), with components m (0) k (n, y n ) = log P (X n = 0|Y n = y n ) P (X n = k|Y n = y n ) , k = 0 . . . q -1 . ( 6
)
February 1, 2013 DRAFT

The messages from CN to VN are computed with the help of a particular Fourier Transform (FT), denoted F(m). Denoting r the unit root associated to GF(q), the i-th component of the FT is given by [START_REF] Li | Density evolution for nonbinary LDPC codes under Gaussian approximation[END_REF] as F i (m) = q-1 j=0 r i⊗j e -m j / q-1 j=0 e -m j .

At iteration ℓ, the message

m (ℓ) (m, n, s m ) from CN m to VN n is m (ℓ) (m, n, s m ) = A[s m ]F -1   n ′ ∈N (m)\n F W H n ′ m m (ℓ-1) (n ′ , m, y n ′ )   (7)
where

sm = ⊖s m ⊘ H n,m , H n ′ m = ⊖H n ′ ,m ⊘ H n,m and W [a] is the q × q matrix such that W [a] k,n = δ(a ⊗ n ⊖ k), 0 k, n q -1. A[k] is a q × q matrix that maps a vector message m into a vector message l = A[k]m with l j = m j⊕k -m k . Note that A[k]
does not appear in the channel coding version of the algorithm and is specific to SW coding. The derivation of ( 7) is

shown in the appendix. At a VN n, a message m (ℓ) (n, m, y i ) is sent to the CN m and an a posteriori message m(ℓ) (n, y n ) is computed. They both satisfy m (ℓ) (n, m, y n ) = m ′ ∈N (n)\m m (ℓ) (m ′ , n, s m ′ ) + m (0) (n, y n ) , (8)
m(ℓ) (n, y n ) = m ′ ∈N (n) m (ℓ) (m ′ , n, s m ′ ) + m (0) (n, y n ) . (9)
From ( 9), each VN n produces an estimate x ℓ) or if ℓ = L max , the maximum number of iterations.

(ℓ) n = arg max k m(ℓ) k (n, y n ) of x n . The algorithm ends if s = H T x (
When the conditional distribution P (Y |X) is uncertain, the previously described decoding algorithm cannot be applied directly, because the initial messages (6) cannot be evaluated accurately.

B. SwP-Source: EM algorithm

We first consider the SwP-Source and then extend the proposed algorithm to the cases of the DwP-and SP-Sources. For the SwP-Source, one needs the actual value of the parameter vector θ because the sum-product LDPC decoder requires the knowledge of the conditional distribution February 1, 2013 DRAFT P (X|Y ). The EM algorithm is thus used to estimate jointly the source sequence X and the parameter θ. A method to produce coarse estimates of the parameters is also described.

1) Joint estimation of θ and x:

The joint estimation of the source vector x and of the parameter θ from the observed vectors y and s is performed via the EM algorithm [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF]. Knowing some estimate θ (ℓ) obtained at iteration ℓ, the EM algorithm maximizes, with respect to θ,

Q(θ, θ (ℓ) ) = E X|y,s,θ (ℓ) [log P (y|X, s, θ)] (10)
= x∈GF(q) n P (x|y, s, θ (ℓ) ) log P (y|x, s, θ) (11)
= N n=1 q-1 k=0 P (X n = k|y n , s, θ (ℓ) ) log P (y n |X n = k, θ) .
Solving this maximization problem gives the update equations detailed in Lemma 1. For simplicity, the correlation model between X and Y is assumed to be additive, i.e., there exists a random The update equation for the EM algorithm is [START_REF] Toto-Zarasoa | Maximum Likelihood BSC Parameter Estimation for the Slepian-Wolf Problem[END_REF]

θ (ℓ+1) = 1 N N n=1 |y n -p (ℓ) n | ( 12
)
where p

(ℓ) n = P (X n = 1|y n , s, θ (ℓ) ).
Let (X, Y ) be a SwP-Source that generates symbols in GF(q). Let the correlation channel be such that Y = X ⊕ Z, where Z is a random variable in GF(q), and P (Z = k) = θ k . The update equations for the EM algorithm are ∀k ∈ GF(q), θ k,n = P (X n = k|y n , s, θ (ℓ) ).

(ℓ+1) k = N n=1 P (ℓ) yn⊖k,n N n=1 q-1 k ′ =0 P (ℓ) yn⊖k ′ ,n (13)
Proof: The binary case is provided by [START_REF] Toto-Zarasoa | Maximum Likelihood BSC Parameter Estimation for the Slepian-Wolf Problem[END_REF]. In the non-binary case, the updated estimate is obtained by maximizing [START_REF] Dupraz | Practical coding scheme for universal source coding with side information at the decoder[END_REF] taking into account the constraints 0 ≤ θ k ≤ 1 and q-1 k=0 θ k = 1.

Note that P (ℓ)

k,n = P (X n = k|y n , s, θ (ℓ) ) in ( 13) can be estimated with a sum-product algorithm that assumes that the true parameter is θ (ℓ) .

2) Initialization of the EM algorithm:

We now propose an efficient initialization of the EM algorithm valid for irregular codes and for sources X and Y taking values in GF(q). This generalizes the method proposed in [START_REF] Toto-Zarasoa | Maximum Likelihood BSC Parameter Estimation for the Slepian-Wolf Problem[END_REF] for regular and binary codes. The rationale is to derive a Maximum Likelihood (ML) estimate of θ from a function u = H T x ⊕ H T y of the observed data H T x and y.

a) The BSC with irregular codes: In this case, each binary random variable U m is the sum of random variables of Z. Although each Z n appears in several sums, the following assumption is made in this section.

Assumption 1. Each U m is obtained from i.i.d. random variables Z (m) j .
The validity of this assumption depends on the choice of the matrix H and is not true in general.

Although it produces an approximate solution, this choice may lead to a reasonable initialization for the EM algorithm. Furthermore, the number of terms in the sum for U m depends on the degree of the CN m. One can use the CN degree distribution ρ(x) as a probability distribution for these degrees, or decide to take into account the knowledge of the CN degrees. Both cases lead to a probability model for the U m and enable to obtain an ML estimate for θ, as described in the two following lemmas.

(θ) = 1 2 -1 2 dc j=2 ρ j (1 -2θ) j , ∀θ ∈ [0, 1 2 ].
Proof: The random variables U m are independent (sums of independent variables). They are identically distributed because the J m and the Z (m) j are identically distributed. α = P (U = 1) = dc j=2 ρ j P (U = 1|J = j). Then, from [START_REF] Toto-Zarasoa | Maximum Likelihood BSC Parameter Estimation for the Slepian-Wolf Problem[END_REF],

P (U = 1|J = j) = j i=1,i odd j i θ i (1 -θ) j-i
and from [12, Section 3.8],

P (U = 1|J = j) = 1 2 -1 2 (1 -2θ) j . Thus α = f (θ). The ML estimate α of α given u is α = 1 M M m=1 u m . Finally, as f is invertible for θ ∈ 0, 1 2 , then from [19,
Theorem 7.2], the ML estimate of θ is given by θ = f -1 ( α).

L(θ) = dc j=2 N 1,j (u) log 1 2 - 1 2 (1 -2θ) j + dc j=2 N 0,j (u) log 1 2 + 1 2 (1 -2θ) j ( 14
)
where N 1,j (u) and N 0,j (u) are the number of symbols in u obtained from the sum of j elements and respectively equal to 1 and 0.

Proof: The random variables U m are independent (sums of independent variables). Therefore, the likelihood function satisfy L(θ) = log P (u|θ) = M m=1 log P (u m |j m , θ). Then, as in the proof of Lemma 2, we obtain [START_REF] Gallager | Source coding with side information and universal coding[END_REF].

The method of Lemma 2 is simpler to implement than the one of Lemma 3 but does not take into account the actual matrix H, at the price of a small loss in performance. b) The non-binary discrete case: Only the case of a regular code is presented here, but the method can be generalized to irregular codes (see the previous section). Assumption 1 also holds in this case. Now, the probability mass function of Z is given by θ = [θ 0 . . . θ q-1 ] with θ k = P (Z = k) ∀k ∈ GF(q). Now, each U m is the sum of symbols of Z, weighted by the coefficients contained in H. A first solution does not exploit the knowledge of these coefficients, but uses the fact that the non-zero coefficients of H are distributed uniformly in GF(q)\{0} (Lemma 4). A second solution takes into account the knowledge of the coefficients (Lemma 5). take their values in GF(q). Denote

θ k = P (Z = k), α k = P (U = k) and assume that θ = [θ 0 . . . θ q-1 ] and α = [α 0 . . . α q-1 ] are unknown.
Then the random variables of U are independent and the parameters satisfy α = f (θ), with

f (θ) = n 1 ,...,n q-1 d c n 1 , . . . , n q-1 1 q -1 dc F -1 q-1 j=0 (F (W [j]θ))) n j ( 15
)
where the sum is over all the possible combinations of integers n 1 , . . . , n q-1 such that 0 ≤ n k ≤ d c and q-1 k=1 n k = d c and dc n 1 ,...,n q-1 is a multinomial coefficient.

Denote θ and α the ML estimates of θ and α, obtained from an observed vector u, with

α k = N k (u) M where N k (u) is the number of occurrences of k in the vector u. Then, if f is invertible, θ = f -1 ( α).
Proof: The random variables U m are independent (sums of independent variables). Then,

α k = P (U = k) = {h j } dc j=1 P ({h j } dc j=1 )P (U = k|{h j } dc j=1
) in which the sum is on all the possible combinations of coefficients {h j } dc j=1 . This can be simplified as α k = n 1 ,...,n q-1 P (N 1 = n 1 , . . . , N q-1 = n q-1 )P (U = k|n 1 , . . . , n q-1 ) where n k is the number of occurrences of k in February 1, 2013 DRAFT {h j } dc j=1 . One has P (N 1 = n 1 , . . . , N q-1 = n q-1 ) = dc n 1 ,...,n q-1 1 q-1 dc . Then, the vector denoted P U|n 1 ,...,n q-1 = [P (U = 0|n 1 , . . . , n q-1 ) . . . P (U = q -1|n 1 , . . . , n q-1 )]

can be expressed as P U|n 1 ,...,n q-1 = F -1

q-1 j=1 (F (W [j]θ))) n j . Therefore, α = [α 0 , . . . , α q-1 ] = n 1 ,...,n q-1 d c n 1 , . . . , n q-1 1 q -1 dc F -1 q-1 j=1 (F (W [j]θ))) n j . (17) The ML estimates α k of α k are α k = N k (u)
M . Finally, if f is invertible, then from [19, Theorem 7.2], the ML estimate of θ is given by θ = f -1 ( α).

Lemma 5. Let U be a length M random vector with entries in GF(q) such that each U m is the

sum of d c i.i.d. random variables, i.e., U m = dc j=1 h (m) j Z (m) j . The Z (m) j are independent ∀j, m,

and identically distributed random variables taking their values in GF(q). The values of the coefficients h

(m) j are known and belong to GF(q)\{0}. Denote

θ k = P (Z = k), α k = P (U = k)
and assume that θ = [θ 0 , . . . , θ q-1 ] and α = [α 0 , . . . , α q-1 ] are unknown. Then the random variables of U are independent and the ML estimate θ from an observed vector u maximizes

L(θ) = M m=1 log F -1 um dc j=1 F(W [h (m) j ]θ) (18)
under the constraints 0 ≤ θ k ≤ 1 and q-1 k=0 θ k = 1.

Proof:

The random variables U m are independent (sums of independent variables). The ML estimate θ is the value that maximizes the likelihood function given by

L(θ) = log P (u|θ, {h (m) j } dc,M j=1,m=1 ) (19)
= M m=1 log P (u m |θ, {h (m) j } dc j=1 ) (20)
under the constraint that 0 ≤ θ k ≤ 1 and q-1 k=0 θ k = 1. The second equality (20) comes from the independence assumption. Following the steps of Lemma 4, we show that (20) becomes

L(θ) = M m=1 log F -1 um dc j=1 F(W [h (m) j ]θ) .
February 1, 2013 DRAFT

C. The DwP-Source

The DwP-Source is non-stationary. Consequently, if one assumes a stationary model such that

P (X n = k|Y n = m) = α k,m (21)
and tries to produce an estimate α(n) k,m from observed sequences (x, y) of length n, then the sequence of estimates α(n) k,m does not necessarily converge as n goes to infinity. However, such an estimate is well defined for a fixed length n. Thus, we apply the procedure defined for the SwP-Source to get α(n) k,m from y and u.

D. SP-Source: MAP with EM

For the SP-Source, the distribution P Θ (θ) is available and one can perform the MAP estimation of Θ. Then, the EM equation ( 10) for the MAP estimation becomes [3]

Q(θ, θ (ℓ) ) = E X|y,s,θ (ℓ) [log P (X|y, s, θ)] + log P Θ (θ) . (22)
Knowing some estimate θ (ℓ) of θ at iteration ℓ, one has to maximize [START_REF] Liveris | Compression of binary sources with side information at the decoder using LDPC codes[END_REF] with respect to θ to obtain θ (ℓ+1) . As for the SwP-Source, the LDPC decoding algorithm initialized with θ (l) provides an approximate version of P (X|y, s, θ (ℓ) ), required to perform the MAP estimation of θ (l+1) .

A coarse estimation of θ can be obtained from u = H T x + H T y as

θ (0) = arg max θ∈P θ log P Θ (θ) + log P (u|H, θ) (23)
in order to initialize the EM algorithm. In the binary case and from the assumptions of Lemma 3 this corresponds to maximizing

L MAP (θ) = log P Θ (θ) + dc j=2 N 1,j (u) log 1 2 - 1 2 (1 -2θ) j + dc j=2 N 0,j (u) log 1 2 + 1 2 (1 -2θ) j (24)
with respect to θ. In the non-binary case and from the assumptions of Lemma 5 this corresponds to maximizing

L MAP (θ) = log P Θ (θ) + M m=1 log F -1 um dc j=1 F(W [h szm,j ]θ) (25)
under the constraints 0 ≤ θ k ≤ 1 and q-1 k=0 θ k = 1.

However, this approach does not fully exploit the density over θ but only its mode, because a hard value of θ is estimated at each iteration and used for the following iterations. To deal with this problem, one could think of using Variational Bayesian Expectation Maximization (VBEM) [START_REF] Beal | Variational Algorithms for Approximate Bayesian Inference[END_REF]. Unfortunately, the VBEM equations are intractable for most of the distributions, particularly in the discrete case. The discrete additive model considered here is not a conjugate exponential model, for which a tractable implementation exists.

V. SIMULATIONS

The performance of the initialization techniques obtained in Lemmas 2 to 5 are first compared.

Then, we evaluate the joint estimation methods proposed for the various models introduced in Section II. In the following, the random variables are either binary or non-binary with values in GF( 4). The correlation model is such that there exists a random variable Z with Y = X ⊕ Z, and X is distributed uniformly.

A. Performance of the initialization techniques (SwP-Model)

The binary case is considered first. Z is such that P (Z = 1) = θ, θ ∈ P θ = [0, 0.18]. We choose a code λ(x) = 0.4295x + 0.2750x 2 + 0.0745x 8 + 0.1150x 9 + 0.0035x 11 + 0.0930x 14 + 0.0095x 15 and ρ(x) = 0.2187x 5 + 0.7760x 6 + 0.0053x 7 , designed for the worst possible parameter θ = 0.18 and obtained from a code optimization realized with a differential evolution algorithm [START_REF] Storn | Differential evolution -a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. The rate of this code is R = 0.75 bit/symbol. The initialization methods of Lemmas 2 and 3 are evaluated and compared through two experiments. Indeed, the models defined in the formulations of the lemmas are supposed to represent the behavior of the LDPC encoding using Assumption 1. In this section, we want to determine whether this assumption is meaningful.

First, we wish to evaluate the performance of the estimation methods on simulated codewords, i.e., generated at random from the models as they are defined in the formulations of then applied to each realization to evaluate the MSE. This gives the two superposed upper curves of Figure 1. As before the two methods give the same performance. However, we observe a loss of a factor 10 in MSE compared to the ideal case, due to the fact that the entries of U are not independent. Nevertheless, the performance seems sufficient for the initialization of the EM algorithm.

For the non-binary case, the probability distribution of Z is given by θ = [θ 0 , . . . , θ 3 ] where P r(Z = k) = θ k . The set P θ is such that ∀θ ∈ P θ , θ 0 ≥ 0.76. We choose a code λ(x) = 0.413x + 0.375x 2 + 0.012x 4 and ρ(x) = x, giving R = 1.6 bit/symbol. In this case, the code was tuned for the worst case θ = [0.76, 0.08, 0.08, 0.08]. We do not perform code optimization here because it is not defined, neither for the non-binary SW setup (see [START_REF] Chen | The equivalence between Slepian-Wolf coding and channel coding under density evolution[END_REF] for the binary SW setup), neither for the models we consider (see [START_REF] Li | Density evolution for nonbinary LDPC codes under Gaussian approximation[END_REF] for density evolution for non-binary codes with Gaussian approximation). However, a code with variable node degree 2 exhibit good performance in the non-binary case [START_REF] Poulliat | Design of regular (2, dc)-LDPC codes over GF(q) using their binary images[END_REF] and the check node degree distribution is chosen to February 1, 2013 DRAFT obtain the rate R = 1.6 bit/symbol. The experiments described in the binary case are repeated for the methods proposed in Lemmas 4 and 5. The parameter estimates are now obtained from a projected gradient descent. Figure 2 shows the MSE of the two cases obtained by averaging over 1000 vectors of length 10000, generated from θ = [0.82 0.06 0.06 0.06] The conclusions of the binary case hold also in this setup and in the following, the method of Lemma 4 is used since it is less complex.

B. Complete coding scheme for the SwP-and SP-Sources

The performance of the complete scheme is now evaluated, in the non-binary case. The parameter θ of the distribution of Z can take the following forms:

θ (0) = [1 -3α, α, α, α] or θ (1) = [α, 1 -3α, α, α] or θ (2) = [α, α, 1 -3α, α] or θ (3) = [α, α, α, 1 -3α], where α ∈ [0, 0.08].
The worst case in terms of rate is given by α = 0.08 and corresponds to four distinct elements of P θ .

The case of the SwP-Source is treated first, and four setups are compared. In each setup, 100 source vectors of length 1000 are generated. For each vector, a θ (i) is selected uniformly at random and α is sampled uniformly at random in [0, 0.08]. The considered LDPC code is the one defined in Section V-A for the non-binary case. We set 20 iterations for the LDPC decoder and 3 iterations for the EM algorithm (when required). Three performance criteria are considered:

the coding rate, given by the choice of the code, the mean error calculated for each vector as

1000 n=1 1 Xn = Xn 1000
where 1 is the indicator function, and the mean decoding time. All the results are presented in Table I.

For the genie-aided setup, θ is given to the decoder.

For the second setup, a learning sequence of 200 source symbols is first transmitted by the encoder, to enable the decoder to produce a ML estimate of θ. This results in a rate increase because the learning sequence is transmitted without the help of the SI, at H(X) = 2 bits/symbol. Note that the mean decoding time is smaller than in the genie-aided setup, because the LDPC The third setup corresponds to the method presented in the paper. Coarse estimate of θ obtained from Lemma 4 initializes the EM algorithm. Here P θ consists of disjoint sets. The first idea is to try to estimate θ successively with the four possible θ (j) , and to keep the θ(j)

that maximizes L( θ(j) ), introduced in [START_REF] Jalali | A universal scheme for Wyner-Ziv coding of discrete sources[END_REF]. Unfortunately, we see that this solution only works when j = 0. This suggests that the method performs well only if z has a subsequent number of null entries. For this reason, the following approach is applied. For j = 0 . . . 3, we compute ũ(j) = H T x⊖H T y ⊖H T j where j is a column vector full of j. Then an estimate θ(j) is produced from each ũ(j) . The corresponding log-likelihood value L( θ(j) ) is calculated and we keep the θ(j) that maximizes L( θ(j) ). In this case, the method performs well.

In the fourth setup, the EM algorithm is initialized at random. We first assume that the θ (0)

is the true distribution, and initialize α at random in [0, 0.08]. If the EM algorithm does not converge with this assumption, the procedure is repeated, assuming that the distribution is θ (1) , and so on until convergence or until θ (3) has been tested. We see that this method increases the decoding time and produces poor performance.

For the SP-Source, the same model is considered and the prior distribution on θ is such that P r(θ (j) ) = p j , p θ = [p 0 . . . p 3 ] = [0.5, 0. initialized at random, are tested again over 100 source vectors of length 1000. In the genieaided setup, θ is given to the decoder. In the learning-sequence case, a MAP estimate of θ is produced from a learning sequence of length 200. For the method described in the paper, the initialization method is performed as before by trying to estimate θ successively with the four possible θ (j) . But now, the selection is performed by keeping the θ(j) that maximizes L MAP ( θ(j) ) = L( θ(j) ) + log p j . For the EM initialized at random, the successive assumptions θ (j) are tested as before, beginning with the θ (j) of highest probability. The results are presented in Table II. We see that there is no change compared to the SwP-setup, except that the EM random is slightly faster.

C. Comparison to a solution with feedback

In this section we compare our no-feedback coding approach with a 1-bit feedback transmission for a source generated from the SwP-Model of Section V-B.. The 1-bit feedback is sent by the receiver to the encoder to ask for additional packets or stop the transmission. The goal is to save rate by avoiding sending data at the worst rate as in the no-feedback method. However, it results in multiple decoding and thus potentially large delays. It is therefore of interest to study the rate/decoding delay tradeoff. Only an evaluation of the achievable rate and estimated mean-time decoding are provided, because it is sufficient to compare the advantages and the February 1, 2013 DRAFT drawbacks of the solution with feedback. In the solution with feedback we consider, when the decoder cannot decode with the received codeword, it requests more check equations via the feedback channel. Each time it receives new equations, the decoder tries to reconstruct the source vector with the use of a sum-product LDPC decoder.

Denote n the length of the source vector and α = 0.08, consider K rate levels R 1 , . . . , R k , . . . , R K associated to K intervals

I 1 = [0, α/K] . . . I K = [α(K -1)/K, α].
The coding system processes as follows. The encoder first sends nR 1 symbols to the decoder. The decoder tries to reconstruct the source, assuming the true parameter is α/K and considering successively the four assumptions θ (0) to θ (3) . If it fails, it sends a request via the feedback channel and the encoder sends n(R 2 -R 1 ) new symbols. The decoder then tries to reconstruct the source from the nR 2 received symbols, assuming the true parameter is 2 × α/K. The process continues until the source vector has been decoded. Note that here, it is assumed that the I k are small enough to allow the decoder to perform good with a parameter that is not exactly the true one.

Five setups are compared, in terms or achievable rate (R) and of estimated mean decoding time (T ). The results are shown in Figure 3. Denote t the decoding time of one LDPC decoder iteration and N it the required number of iterations. In the following, we set K = 8, t = 0.2s, N it = 20. Remark that the entropy of a source of parameter α is the same whatever θ (i) and denote h(α) = H(X|Y, α). In each case, we assume that a code or a sequence of codes reaching the entropy can be constructed.

For the solution with feedback, assume that we can construct a sequence of codes such that R 1 = h(α/K), . . . R K , = h(α) and achieving small probability of error respectively for α ∈ I 1 , . . . , α ∈ I K . Thus, α ∈ I k , R k = h(kα/K). We also assume that the delay induced by the feedback is negligible compared to the decoding time. Then, for α ∈ I k the mean decoding time is estimated as T k = 2 × t × N it × k where the coefficient 2 comes from the assumption that θ (i) is selected uniformly at random. In the curve of Figure 3, the circles represent the (R k , T k )

and the line can be seen as an interpolation when K increases. For the genie-aided setup, the rate is dimensioned for the worst case, i.e., R = h(α) and an approximation of the mean decoding time is calculated as T = N it . For the setup with learning sequence, assuming a sequence length representing a fraction 1/5 of the total length, R = 4/5h(α) + 1/5H(X) and T = N it . For the coding scheme described in the paper, R = h(α) and we approximate T = 2 × N it , assuming that 2 iterations of the EM algorithm are required, and that the θ (i) is correctly retrieved at the initialization step. For the coding scheme with EM initialized at random, R = h(α) and we approximate T = 3 × 2 × N it , assuming that 3 iterations of the EM algorithm are required, and that θ (i) is selected uniformly at random. We see that a solution with feedback can induce a considerable rate gain but at the price of huge decoding delays. Remark that the choice of the parameter K is important. If K decreases, the size of the I k increases which reduces the mean decoding time. On the other hand, as for α ∈ I k , the effective coding rate is R k , the rate needed to decode for α can increase.

D. DwP-Source

The solution proposed for the DwP-Source is now evaluated in the non-binary case. The distribution of Z is given by π = [1 -3π, π, π, π], where π can take four different values:

π 0 = 0, π 1 = 0.02, π 2 = 0.04 and π 3 = 0.06, defining four states. Indeed, for the DwP-Model, the worst case is taken on the convex hull of the possible distributions.

We now consider source vectors of length 1000 and fix a block length m. For each block of length m in a vector, a probability distribution for the states is generated uniformly at random.

The values m = 1, 100 and 500 are tested. A fourth setup is tested for m generated uniformly at random in {1, ..., 1000}. More precisely, a first value m 0 is generated, a second one m 1 is generated after m 0 symbols and so on. The method proposed for the SwP-Source is then applied with the same code over 100 realizations for each m. The complete decoding technique described for the SwP-Source is used: coarse estimate of the parameter from Lemma 4 followed by EM algorithm. The decoder is provided the information that the set of possible distributions for Z is on the form In this case, from (2), the infimum of achievable rates is given by R = 2 bits/symbol = H(X) bits/symbol, because the infimum of achievable rates is given by the worst cases on the convex hull of the set of possible distributions. Therefore, the SI cannot be exploited.

VI. CONCLUSION

This paper introduced four signal models modeling the uncertainty on correlation channel between the source and the SI. Practical coding schemes based on non-binary LDPC codes were proposed for the SW setup and for the four models. Simulation results exhibit good performance in terms of probability of error, rate, or decoding delay, compared to the solution with a learning sequence or the solution with an EM algorithm initialized at random. Future works will be on the design of good degree distributions for our models with nonbinary symbols, and on the extension to the lossy case. We will also investigate correlation model selection, i.e., the choice of one of the four source correlation models and of the structure of the family distribution for the model. denoted m (ℓ-1) (n ′ , m, y n ′ ). Thus, the k-th component of the CN message m to VN n (7) is log P X n = 0|s m , {m (ℓ-1) (n ′ , m, y n ′ )} n ′ ∈N (m)\n P X n = k|s m , {m (ℓ-1) (n ′ , m, y n ′ )} n ′ ∈N (m)\n [START_REF] Richardson | Design of capacity-approaching irregular Low-Density Parity-Check codes[END_REF] We first detail the impact of the operator ⊗ on a message, to study the term ⊖(H m,n ′ ⊘ H m,n ) ⊗ x n ′ in [START_REF] Puri | PRISM: A new robust video coding architecture based on distributed compression principles[END_REF]. Given a random variable Z taking its values in GF(q), and with a probability vector p = [P (Z = 0), . . . , P (Z = q -1)] T , the probability vector of a ⊗ Z satisfies q = [P (a ⊗ Z = 0), . . . , P (a ⊗ Z = q -1)] T = W [a]p, where the matrix W [a]

has been defined just after equation [START_REF] Csiszar | Linear codes for sources and source networks: Error exponents, universal coding[END_REF]. Similarly, l = log P (a⊗Z=0) P (a⊗Z=0) , . . . , log P (a⊗Z=0)

P (a⊗Z=q-1)
is obtained from m = log P (Z=0) P (Z=0) , . . . , log P (Z=0) P (Z=q-1) as l = W [a]m. Therefore, in [START_REF] Richardson | Design of capacity-approaching irregular Low-Density Parity-Check codes[END_REF], we need W H n ′ m m (ℓ-1) (n ′ , m, y n ′ ), ∀n ′ ∈ N (m)\n, where H n ′ m = ⊖H n ′ ,m ⊘ H n,m .

We now detail the impact of the operator ⊖ on a message to deal with ⊖ n ′ ∈N (m)\n (H m,n ′ ⊘ H m,n ) ⊗ x n ′ . The probabilities of a sum of random variables in GF(q) can be evaluated with the help of a particular Fourier transform [START_REF] Goupil | FFT-based BP decoding of general LDPC codes over Abelian groups[END_REF]. From [START_REF] Li | Density evolution for nonbinary LDPC codes under Gaussian approximation[END_REF], the ith component of the Fourier transform applied on a message vector m is F i (m) = q-1 j=0 r i⊗j e -m j / q-1 j=0 e -m j and the k-th component of its inverse is F -1 k (f ) = log q-1 i=0 f i / q-1 i=0 r -i⊗k f i .

Finally, the term s m ⊘ H m,n , specific to SW coding, is taken into account. Denote Γ a random variable taking its values in GF(q) and m = log P (Γ=0) P (Γ=0) , . . . , log P (Γ=0) P (Γ=q-1) . The message vector l = log P (a⊕Γ=0) P (a⊕Γ=0) , . . . , log P (a⊕Γ=0) P (a⊕Γ=q-1) corresponding to a ⊕ Γ is obtained as l = A[a]m. Setting a = s n ⊘ H m,n gives the final message vector m (ℓ) (n, m, y n ).

  variable Z such that Y = X ⊕Z and θ parametrizes the distribution of Z. The Binary Symmetric correlation Channel (BSC) of unknown transition probability θ = P (Y = 1|X = 0) = P (Y = 0|X = 1) is a special case, where Z is a binary random variable such that P (Z = 1) = θ. Lemma 1. Let (X, Y ) be a binary SwP-Source. Let the correlation channel be a Binary Symmetric channel (BSC) with parameter θ = P (Y = 0|X = 1) = P (Y = 1|X = 0), θ ∈ [0, 1].

February

Lemma 2 .

 2 Let U be a binary random vector of length M . Each U m is the sum of J m identically distributed binary random variables Z (m) j , i.e., U m = Jm j=1 Z (m) j , where the Z (m) j are independent February 1, 2013 DRAFT ∀j, m. {J m } M m=1 are i.i.d. random variables taking their values in {2, . . . , d c } with known probability P (J = j) = ρ j . Denote θ = P (Z = 1), α = P (U = 1) and assume that θ and α are unknown. Then their ML estimates θ and α from an observed vector u satisfy α = 1 M M m=1 u m and θ = f -1 ( α), where f is the invertible function f

Lemma 3 .

 3 Let U be a binary random vector of length M . Each U m is the sum of j m identically distributed binary random variables Z (m) j , i.e., U m = jm j=1 Z The values of j m are known and belong to {2, . . . , d c }. Denote θ = P (Z = 1) and assume that θ is unknown. Then the entries of U are independent and the ML estimate θ from an observed vector u is the argument of the maximum of

February

Lemma 4 .

 4 Let U be a length M random vector with entries in GF(q) such that each U m is the sum of d c i.i.d. products of random variables, i.e., U m = dc j=1 H j are identically distributed random variables, mutually, and individually independent ∀j, m.The H (m) j are uniformly distributed in GF(q)\{0}. The Z (m) j

FebruaryFig. 1 .FebruaryFig. 2 .

 12 Fig. 1. MSE of the estimators for the binary case.

[ 1 -

 1 3θ, θ, θ, θ] and θ ∈ [0, 0.06]. The results are presented in Table III. Mean error and mean decoding time are measured. We see that smaller m give smaller decoding time, because with a small m, the source behavior is better approximated by the SwP-Model. Now consider a DwP-Model with four states, giving four possible distributions for Z: π (0) = [0.91, 0.03, 0.03, 0.03], π (1) = [0.03, 0.91, 0.03, 0.03], π (2) = [0.03, 0.03, 0.91, 0.03], π (3) = February 1, 2013 DRAFT [0.03, 0.03, 0.03, 0.91].

  |X, θ)} θ∈P θ parametrized by a vector θ. The vector θ takes its value in a set P θ that is either discrete or continuous. The source symbols X and Y take their values in the discrete sets X and Y, respectively. Moreover, the parameter θ is fixed for the sequence

	Definition 4. (SwP-Source). A Static without Prior source (X, Y ) (SwP-Source) produces a
	sequence of independent symbols {(X n , Y n )} +∞ n=1 drawn from a distribution belonging to a family
	{P (X, Y |θ) = P (X)P (Y
	The following models consider a time-invariant parameter vector.
	Definition 3. (SP-Source) A Static with Prior source (X, Y ) (SP-Source) produces a sequence of
	independent symbols {(X n , Y n )} +∞ n=1 drawn from a distribution belonging to a family {P (X, Y |Θ =
	θ) = P (X)P (Y |X, Θ = θ)} θ∈P θ parametrized by a random vector Θ. The random vector Θ,
	with distribution P Θ (θ), takes its value in a set P θ that is either discrete or continuous. The
	source symbols X and Y take their values in the discrete sets X and Y, respectively. Moreover,
	the realization of the parameter θ is fixed for the sequence {(X n , Y n )} +∞ n=1 .
	The SP-source, determined by P θ , P Θ (θ), and {P (X, Y |Θ = θ)} θ∈P θ , is stationary but non-
	ergodic [13, Section 3.5].
	February 1, 2013	DRAFT

February 1, 2013 DRAFT random vector Π n . The {Π n } +∞ n=1 are i.i.d. with distribution P (Π) and take their values in a discrete set P π . The source symbols X n and Y n take their values in the discrete sets X and Y, respectively. The DP-Source, completely determined by P π , P (Π), and {P (X, Y |Π = π)} π∈Pπ , is stationary and ergodic, see [13, Section 3.5]. Definition 2. (DwP-Source). A Dynamic without Prior source (X, Y ), or DwP-Source, produces a sequence of independent symbols {(X n , Y n )} +∞ n=1 drawn ∀n from P (X n , Y n ) that belongs to a family of distributions {P (X, Y |π) = P (X)P (Y |X, π)} π∈Pπ parametrized by a vector π n . Each π n takes its values in a discrete set P π . The source symbols X n and Y n take their values in the discrete sets X and Y, respectively.

The DwP-Source, determined by P π and {P (X, Y |π)} π∈Pπ , is non-stationary and non-ergodic

[START_REF] Gallager | Information theory and reliable communication[END_REF] Section 3.5]

. The only difference between the DP-and DwP-Sources lies in the definition of the parameters π n . In the DwP-Source, no distribution for π n is specified, either because its distribution is not known or because π n is not modeled as a random variable.

  [START_REF] Poulliat | Design of regular (2, dc)-LDPC codes over GF(q) using their binary images[END_REF], 0.125, 0.125] and α is distributed uniformly on [0, 0.08]. The four setups: genie-aided, learning sequence, method described in the paper, EM

	Setup	Err Time (s) Rate (bit/symb.)
	Genie-aided	< 10 -5	5.4	1.6
	Learn. Seq.	< 10 -5	4.2	1.7
	EM	< 10 -5	9.1	1.6
	EM random 8.4 × 10 -3	41.0	1.6
	February 1, 2013			DRAFT

February 1, 2013 DRAFT

ACKNOWLEDGMENTS

The authors would like to thank Valentin Savin for the helpful discussions and advices.

Part of this paper will be presented at the Data Compression Conference (DCCPart of this work was supported by the ANR-09-VERS-019-02 grant (ARSSO project).

APPENDIX

In this Appendix, we detail the derivation of the update rule (7) at a CN for the SW problem, when the LDPC code is non binary and the decoder is the sum-product algorithm. This update rule derives from the parity check equation at CN m, given by n ′ ∈N (m) H m,n ′ ⊗ x n ′ = s m , that can be restated as

The update rule at a CN, and for the sum-product algorithm, consists in computing the reliability information on the variable x n as a function of the reliability information on the variables x n ′ , February 1, 2013 DRAFT