\documentclass[12pt]{./styles/osuthesis}
\usepackage[top=1.1in,bottom=1.25in,left=1.7in,right=1.4in]{geometry}
%\usepackage{graphicx}
%\usepackage[draft]{graphicx}
%\usepackage{subfigure}
%\usepackage{verbatim}
%\usepackage{enumerate}
%\usepackage{amsmath}
%\usepackage{amscd}
\usepackage{caption}
\usepackage{color}
%\usepackage{amsxtra}
%\usepackage{amstext}
\usepackage{amssymb}
%\usepackage{stmaryrd}
%\usepackage{hyperref}
%\usepackage{mathrsfs}
%\usepackage{latexsym}
%\usepackage{newlfont}
%\usepackage{indentfirst,color}
%\usepackage{doublespace}
\usepackage{amsmath}
\usepackage{float}
\usepackage{graphicx,amssymb,amstext,amsmath}
\usepackage{setspace,chngpage}
\usepackage{cite}
\usepackage{amsmath,eucal}
\usepackage{subfigure}
\usepackage{wrapfig}
\usepackage{algorithm}
\usepackage{algorithmic}
\setcounter{tocdepth}{1} % set how depth the table of content will list. 1 for section; 2 for subsection
\makeatletter
\newcommand{\rmnum}[1]{\romannumeral #1}
\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}
\makeatother
\newcommand{\sw}[1]{\color{black}#1}
\newtheorem{remark1}{Remark}
\newcommand{\fc}[1]{\textcolor{black}{#1}}

\ifdefined\showcomment
  \newcommand{\ths}[1]{{\textcolor{green}{#1}}}
  \newcommand{\scc}[2]{{\textcolor{green}{#2}} {\textcolor{red}{(#1)}}}
  \newcommand{\comm}[1]{\emph{\textcolor{red}{\textbf{SC:}#1}}}
  \newcommand{\com}[1]{{\textcolor{blue}{#1}}}
 \newcommand{\delete}[1]{\sout{#1}}
 \newcommand{\replace}[2]{{\textcolor{green}{#1}} \sout{#2}}
\else
  \newcommand{\ths}[1]{#1} % remove comment
  \newcommand{\scc}[2]{#2} % remove comment
  \newcommand{\comm}[1]{} % remove comment
  \newcommand{\com}[1]{#1} % remove comment
  \newcommand{\delete}[1]{} % remove comment
\newcommand{\replace}[2]{#1}
\fi

\begin{document}
\begin{center}
%\vspace[*]{3 cm}
%\headsep = 1.5cm
\thispagestyle{empty} °°
\\
UNIVERSITY OF OKLAHOMA\\
\vspace{0.2cm}
GRADUATE COLLEGE\\
\vspace{2.5cm}
ADAPTIVE DISTRIBUTED SOURCE CODING BASED ON\\
%\linebreak[1]
\vspace{0.2cm}
BAYESIAN INFERENCE\\
%\linebreak[1]
%\linebreak[1]
%\linebreak[1]
%USING BELIEF PROPAGATION\\
\vspace{3.5cm}
A DISSERTATION\\
\vspace{0.2cm}
SUBMITTED TO THE GRADUATE FACULTY\\
\vspace{0.2cm}
in partial fulfillment of the requirements for the\\
\vspace{0.2cm}
Degree of\\
\vspace{0.2cm}
DOCTOR OF PHILOSOPHY\\
\vspace{2.5cm}
By\\
\vspace{0.4cm}
\singlespace
FENG CHEN\\
Norman, Oklahoma\\
2015
\end{center}

%%Signature Page
\newpage
\thispagestyle{empty}
%\headsep = 5cm
\begin{center}
\msp
ADAPTIVE DISTRIBUTED SOURCE CODING BASED ON\\ 
BAYESIAN INFERENCE
\linebreak[1]
\linebreak[1]
\linebreak[1]
A DISSERTATION APPROVED FOR THE\\
SCHOOL OF ELECTRICAL AND COMPUTER ENGINEERING
\linebreak[1]
\linebreak[1]
\linebreak[1]
\linebreak[1]
\linebreak[1]
\linebreak[1]
\linebreak[1]
\linebreak[1]
%\linebreak[1]
BY
\linebreak[1]
\linebreak[1]
\end{center}
%\raggedright
\begin{flushright}
\rule{6.5cm}{0.02em}\\
\vspace{-0.15in}
    Dr. Samuel Cheng, Chair\\
%\linebreak[1]
%\linebreak[1]
\vspace{0.2in}
\rule{6.5cm}{0.02em}\\
\vspace{-0.15in}
       Dr. Pramode Verma\\
%\linebreak[1]
%\linebreak[1]
\vspace{0.2in}
\rule{6.5cm}{0.02em}\\
\vspace{-0.15in}
  Dr. Nikola Petrov\\
\vspace{0.2in}
%\linebreak[1]
%\linebreak[1]
\rule{6.5cm}{0.02em}\\
\vspace{-0.15in}
  Dr. Choon Yik Tang\\
\vspace{0.2in}
%\linebreak[1]
%\linebreak[1]
\rule{6.5cm}{0.02em}\\
\vspace{-0.15in}
  Dr. Gregory MacDonald

\end{flushright}
%%copyright page
\newpage
\singlespace
\thispagestyle{empty}
\headsep = 20.1cm
\begin{center}
\copyright\ Copyright by FENG CHEN 2015\\
All Rights Reserved.
\end{center}

\newpage
\doublespace
\headsep = 0in
\pagenumbering{roman}\setcounter{page}{4}
\begin{center}
{\large \bf Acknowledgments}
\end{center}

I would like to take this opportunity to express my deepest gratitude to my adviser and mentor Dr. Samuel Cheng for his hard work and his support during my academic career. During the last three years, Dr. Cheng offered me an amazing research environment, and in this environment I progressed from the status of an amateur researcher to the status of a real researcher. Without his guidance and support, I would not be able to complete this dissertation. I would also like to express my gratitude to Dr. Verma, Dr. Petrov, Dr. Tang, and Dr. MacDonald for serving as my committee members and for guiding me during my dissertation work.

I would like to thank all of the members of Dr. Cheng's group and my other classmates at OU-Tulsa TCOM: Shuang Wang, Lijuan Cui, Amin, Nafise, Lei Yang, Dong Han, etc.

I also want to thank all of my former mathematics and physics teachers. These teachers prompted and encouraged my love for math and science from the time that I was 7 years old. 

I would like to offer my sincere gratitude to my family members. My parents raised me and gave me my preliminary education. I would also like to thank my wife, Huiling Wu, who encouraged me to complete this dissertation, and who gave birth to my dear son Michael.

%\doublespacing
\tableofcontents
\listoftables
\listoffigures
\newpage
%\pagenumbering{roman}\setcounter{page}{\value{placeholder}}
\renewcommand{\thepage}{\roman{page}}
\setcounter{page}{\value{placeholder}}
\begin{center}
{\large \bf Abstract}
\end{center}

\noindent Distributed Source Coding (DSC) is an important topic for both in information theory and communication. DSC utilizes the correlations among the sources to compress data, and it has the advantages of being simple and easy to carry out. In DSC, Slepian-Wolf (S-W) and Wyner-Ziv (W-Z) are two important problems, which can be classified as lossless compression and loss compression, respectively. Although the lower bounds of the S-W and W-Z problems have been known to researchers for many decades, the code design to achieve the lower bounds is still an open problem. 

This dissertation focuses on three DSC problems: the adaptive Slepian-Wolf decoding for two binary sources (ASWDTBS) problem, the compression of correlated temperature data of sensor network (CCTDSN) problem and the streamlined genome sequence compression using distributed source coding (SGSCUDSC) problem. For the CCTDSN and SGSCUDSC problems, sources will be converted into the binary expression as the sources in ASWDTBS problem for encoding. The Bayesian inference will be applied to all of these three problems. To efficiently solve these Bayesian inferences, message passing algorithm will be applied. For a discrete variable that takes a small number of values, the belief propagation (BP) algorithm is able to implement the message passing algorithm efficiently. However, the complexity of the BP algorithm increases exponentially with the number of values of the variable. Therefore, the BP algorithm can only deal with discrete variable that takes a small number of values and limited continuous variables. For the more complex variables, deterministic approximation methods are used. These methods, such as the variational Bayes (VB) method and expectation propagation (EP) method, can efficiently incorporated into the message passing algorithm.

A virtual binary asymmetric channel (BAC) channel was introduced to model the correlation between the source data and the side information (SI) in ASWDTBS problem, in which two parameters are required to be learned. The two parameters correspond to the crossover probabilities that are $0\rightarrow 1$ and $1 \rightarrow 0$. Based on this model, a factor graph was established that includes LDPC code, source data, SI and both of the crossover probabilities. Since the crossover probabilities are continuous variables, the deterministic approximate inference methods will be incorporated into the message passing algorithm. The proposed algorithm was applied to the synthetic data, and the results showed that the VB-based algorithm achieved much better performance than the performances of the EP-based algorithm and the standard BP algorithm. The poor performance of the EP-based algorithm was also analyzed.   

For the CCTDSN problem, the temperature data were collected by crossbow sensors. Four sensors were established in different locations of the laboratory and their readings were sent to the common destination. The data from one sensor were used as the SI, and the data from the other $3$ sensors were compressed. The decoding algorithm considers both spatial and temporal correlations, which are in the form of Kalman filter in the factor graph. To deal with the mixtures of the discrete messages and the continuous messages (Gaussians) in the Kalman filter region of the factor graph, the EP algorithm was implemented so that all of the messages were approximated by the Gaussian distribution. The testing results on the wireless network have indicated that the proposed algorithm outperforms the prior algorithm.

The SGSCUDSC consists of developing a streamlined genome sequence compression algorithm to support alternative miniaturized sequencing devices, which have limited communication, storage, and computation power. Existing techniques that require a heavy-client (encoder side) cannot be applied. To tackle this challenge, the DSC theory was carefully examined, and a customized reference-based genome compression protocol was developed to meet the low-complexity need at the client side. Based on the variation between the source and the SI, this protocol will adaptively select either syndrome coding or hash coding to compress variable lengths of code subsequences. The experimental results of the proposed method showed promising performance when compared with the state of the art algorithm (GRS).
\newpage
\setcounter{placeholder}{\value{page}}
%\msp    %========  single space text.
        %========  For the final version, this command should not be used.
%====================== main  body of the dissertation ==========================
\doublespacing
\pagenumbering{arabic}\setcounter{page}{1}

 \chapter{INTRODUCTION}\label{chap1}

In a wireless sensor network, all of the sensors may send their own data to a common destination. However, due to the limited capacity of the communication channel and the limited memory, the data must first be compressed so that it can be sent to the common destination. \replace{One major difference between wireless sensor networks and other mobile networks is that the observations among sensors are highly correlated. This means that there is a significant amount of redundancy among the observations.}{Because all of the sensors are installed in the same network, there may exist some correlations among their data, which may help compress the sensors' data.} Distributed Source Coding (DSC) theory is concerned with the separate compressions of these correlated sources at each sensor and the joint decompression at the destination.  

In 1973, Slepian and Wolf published their classical DSC paper~\cite{SlepianW:73}. In this paper, Slepian and Wolf proved a very surprising result, that is, generally speaking, it is possible to have no performance loss in separate encoding compared to the case in which joint encoding is allowed. Wyner and Ziv later considered a lossy version of asymmetric DSC in \cite{WynerZ:1976,Wyner:78}. Because of the contributions of Slepian-Wolf and Wyner-Ziv, we usually refer to lossless DSC as an S-W problem and lossy source coding with side information as a W-Z problem. 

Since most of the problems in this dissertation can be classified as S-W coding or W-Z problems, these two codings will be briefly introduced first. Then, the contributions of this dissertation will be summarized. Finally, the organization of this dissertation will be outlined. 

\delete{DSC started as an information-theoretical problem in the renowed 1973 paper of Slepian and Wolf~\cite{SlepianW:73}, and their coding theory is called Slepian-Wolf (SW) coding.  SW coding is a lossless data compression technique, which encodes each source separately and decodes jointly. By this, the computation burdens are shifted from the encoded sides to the decoded side. A large number of SW coding schemes have been proposed\cite{PradhanR:99,AaronG:02,AaronRZ:03,RebolloZ:03,Ortega:05, VarodayanAG:05,xu2006layered,kubasov2007optimal,brites2008correlation, wang2011noise,krithivasan2011distributed,cui2011adaptive,CSVT2011, mauniversality,pan2012low,barcelo2012enhanced,cui2012adaptive, ma2013zero,deligiannis2014maximum,vaezi2014distributed}.

In 1976, Wyner and Ziv (WZ) considered a lossy version of the asymmetric SW problem, where one source is available at the decoder as side information~\cite{WynerZ:76}. WZ source coding is usually realized by quantization followed by SW coding of quantization induces based on channel coding\cite{AaronRZ:03,aaron2002wyner,St_louis_workshop:03,YangCXZ:03, RebolloZ:03,Aaron2004,aaron2004wyner,brites2006improving,ascenso2006content, xu2006layered,kubasov2007optimal,brites2008correlation,CSVT2011,wang2012onboard, huang2012cross,brites2012learning,chen2013compression,deligiannis2014maximum}. Low-Density Parity-Check (LDPC)-based method has been commonly used to realize SW coding\cite{WZcoding2003}. Noticeably, by using LDPC codes, it is possible to compress a joint binary source very close to the SW limit (i.e., the joint entropy). However, some data are correlated closely not only in spatial domain but also in temporal domain. Actually in real-world applications such as sensor networks, the value measured by each sensor is usually not independent with temporal neighbored values.  If we could make use of these correlations, the performance of distributed source coding (DSC) should be improved.}

Before continuing, some notations should be classified:
\begin{itemize}
\item the uppercase letters $X, Y, \dots$ denote the data source or random variables;
\item the lowercase letters $x, y, \dots$ denote one sample from $X, Y, \dots$;
\item $\mathcal{X}, \mathcal{Y}, \dots$ denote the alphabets;
\item $x^n$ denotes one n-length sequence $\{x_1,\dots,x_{n}\}$;
\item $\hat{x}$ denotes the decoded word (estimation) of $x$.
\end{itemize}

\section{S-W Coding}
\begin{figure}[!htb]
\includegraphics[width=0.8\textwidth]{./figures/S-W-Region}
\caption{Lossless Coding Region}
\label{fig:S-W-Region}
\end{figure}
\begin{definition}
Given two correlated joint sources $(X, Y)$, where $X \in \mathcal{X}$ and $Y \in \mathcal{Y}$, we will encode the $n$-length i.i.d sequences of $(X, Y)$  with rate $(R_X, R_Y)$ with the maps
\begin{equation}
\begin{split}
f_1: \mathcal{X}^n \rightarrow \{1, \cdots, 2^{nR_x}\}\\
f_2: \mathcal{Y}^n \rightarrow \{1, \cdots, 2^{nR_y}\}
\end{split}
\end{equation}
and decode the codeword with the map
\begin{equation}
\begin{split}
g: \{1, \cdots, 2^{nR_x}\} \times \{1, \cdots, 2^{nR_y}\}
 \rightarrow \mathcal{X}^n \times \mathcal{Y}^n.
\end{split}
\end{equation}
\delete{Suppose the following scenario: given two correlated i.i.d sequences $X \in \mathcal{X}$ and $Y \in \mathcal{Y}$ as the plaintext,  where $\mathcal{X}$ and $\mathcal{Y}$ are two alphabets, we will compress these two sequences to the codewords with rates $R_X$ and $R_Y$ respectively so that sending to the same decoder separately. The decoder will decode the codewords from $X$ and $Y$ back to the symbols in the alphabets $\hat{\mathcal{X}}$ and $\hat{\mathcal{Y}}$ respectively, i.e.
\begin{equation}
\begin{split}
\mbox{Encode} \ \ f: \mathcal{X} \times \mathcal{Y} \rightarrow \{1, \cdots, 2^{nR_x}\} \times \{1, \cdots, 2^{nR_y}\}, \\
\mbox{Decode} \ \ g: \{1, \cdots, 2^{nR_x}\} \times \{1, \cdots, 2^{nR_y}\}
 \rightarrow \hat{\mathcal{X}} \times \hat{\mathcal{Y}}.
\end{split}
\end{equation}}

If the decoding errors of these sequences are arbitrarily small, i.e.,
\begin{equation}
\begin{split}
\lim_{n \rightarrow \infty} \frac{1}{n}\sum_{i=1}^{n} D(x_i,\hat{x}_i)=0,\\ 
\lim_{n \rightarrow \infty} \frac{1}{n}\sum_{i=1}^{n} D(y_i,\hat{y}_i)=0,
\end{split}
\end{equation}
where $D(\cdot, \cdot)$ is the distortion function, then we designate these lossless codings as S-W codings.
\end{definition}

It turns out that the achievable rate region (see Figure \ref{fig:S-W-Region}) of the S-W coding is given by\cite{IT}
\begin{equation}
\begin{split} \label{SWConditions}
&R_X\ge H(X|Y) \\
&R_Y\ge H(Y|X) \\
&R_X + R_Y \ge H(X,Y),
\end{split}
\end{equation}
where $H(\cdot)$ is the entropy operation.

\begin{figure}[!htb]
  \centering
  \subfigure[ ]{\label{fig:Non-asymmetric_S-W}\includegraphics[width=0.4\textwidth]{./figures/Non-asymmetric_DSC.eps}}
  \subfigure[ ]{\label{fig:asymmetric_S-W}\includegraphics[width=0.4\textwidth]{./figures/asymmetric_DSC.eps}}
  \caption{Two types of S-W coding: (a) Non-asymmetric S-W coding; (b) asymmetric S-W coding.}
  \label{fig:S-W_Coding}
\end{figure}
Depending on whether the source $Y$ will be compressed or not, the S-W coding can be classified as two schemes: non-asymmetric S-W coding or asymmetric S-W coding(see Figure \ref{fig:S-W_Coding}). For asymmetric S-W coding, the source $Y$ is uncompressed, which indicates $R_Y = H(Y)$. Hence, to satisfy (\ref{SWConditions}), $R_X$ should be $R_X \ge H(X,Y) - R_Y = H(X,Y) - H(Y) = H(X|Y)$.
\begin{remark}
The proof of the S-W Theorem~\cite{SlepianW:73} is non-constructive, in that Slepian and Wolf did not indicate how to implement DSC efficiently. The research on practical algorithms for DSC was stagnant until the publication of the work by Pradhan and Ramchandran in 1999~\cite{PradhanR:99}. They rediscovered an early work by Wyner in which he suggested the use of channel codes for asymmetric S-W coding \cite{Wyner:74}. 
\end{remark}

\begin{remark}
With few exceptions, such as \cite{grangetto2007dac}, the majority of S-W code designs are based on this channel coding formulation. Practical syndrome-based schemes for S-W coding using channel codes have been studied in \cite{PradhanR:99,AaronG:02,AaronRZ:03,RebolloZ:03,Ortega:05, VarodayanAG:05,xu2006layered,kubasov2007optimal,brites2008correlation, wang2011noise,krithivasan2011distributed,cui2011adaptive,CSVT2011, mauniversality,pan2012low,barcelo2012enhanced,cui2012adaptive, ma2013zero,deligiannis2014maximum,vaezi2014distributed}.
\end{remark}

\section{W-Z Coding}
\delete{
Practical approaches to W-Z coding have been investigated in \cite{AaronRZ:03,aaron2002wyner,St_louis_workshop:03,YangCXZ:03, RebolloZ:03,Aaron2004,aaron2004wyner,brites2006improving,ascenso2006content, xu2006layered,kubasov2007optimal,brites2008correlation,CSVT2011,wang2012onboard, huang2012cross,brites2012learning,chen2013compression,deligiannis2014maximum}. The majority of these works are based on the source-channel coding framework, where sources at different encoders are quantized by source coding and the resulting quantization indices are further compressed using channel-coding based S-W coding. For example, in DISCUS \cite{PradhanR:99}, two sources codes (scalar quantization and trellis coded quantization \cite{MarcellinF:90} and two channel codes (scalar coset code and trellis-based coset code \cite{Ungerboeck:82}) are used in source-channel coding for the W-Z problem, resulting in four combinations.}

The lossy decoding version of Figure \ref{fig:asymmetric_S-W} is called the "W-Z coding" problem. \delete{Let sources $X$, $Y$ and $\hat{X}$ take values in the finite sets $\mathcal{X}$, $\mathcal{Y}$ and $\hat{\mathcal{X}}$ respectively, and $\{(X_n, Y_n)\}_{n =1}^N$ denotes the sequence of the sources $(X, Y)$.}For a given distortion $d$, the W-Z coding problem tries to find the coding method that achieves the distortion requirement with a minimum rate, i.e., for an $n$ block encoding scheme
\begin{equation}
\begin{split}
&\{f^{*}, g^{*}\} = \arg\min_{f, g}\{R(f,g)\},\\
&s.t.\\
&f: {\mathcal{X}}^n \rightarrow \{1,\dots, 2^{nR} \}\\
&g: {\mathcal{Y}}^n \times \{1, \dots, 2^{nR} - 1\} \rightarrow {\mathcal{X}}^{n}\\
&\lim\sup_{n \rightarrow \infty} \mathbb{E}[D(x^n, g(f(x^n)))] \le d
\end{split},
\end{equation}
where $f$ and $g$ are the encoding and decoding mappings, respectively, $\mathbb{E}[\cdotp]$ denotes the expectation operation, and $D(\cdotp,\cdotp)$ denotes the distortion function. The minimum rate $R$ is usually written as $R_{Y}^{*}(d)$

In \cite{WynerZ:1976,Wyner:78}, the authors obtained that 
\begin{equation}
R_{Y}^{*}(d) = \inf_{p \in \mathcal{M}(d)} [I(X;Z) - I(Y;Z)],
\end{equation}
$Z \in \mathcal{Z}$, where $\mathcal{Z}$ is arbitrary finite set, and $\mathcal{M}(d)$ is the set of $p(x, y, z)$ which satisfies
\begin{equation}
\begin{split}
&\sum_{z \in \mathcal{Z}}{p(x,y,z)}=p(x,y), \\
&Y,Z \mbox{ are conditionally independent given }X.
\end{split}
\end{equation}

The W-Z problem is usually implemented by quantizing the source data followed by the S-W coding.
\section{Contributions}
This dissertation will focus on how to apply the Bayesian inference framework to solve parts of the S-W and W-Z problems. The Bayesian model can be conveniently expressed by a factor graph, in which the message passing algorithm can be implemented to solve the model. Generally, Belief Propagation (BP) is an efficient message passing algorithm. However, for complex situations, such as those involving continuous random variables, it is generally infeasible to apply the BP algorithm. To solve this problem, the deterministic approximate methods will be incorporated into the message passing algorithm so that the factor graph can be efficiently decoded.

The contributions of this dissertation can be summarized as follows: 
\begin{itemize}
\item Research the properties of different deterministic approximate algorithms, such as VB and EP for DSC decoding.
\item Proposed an adaptive joint S-W decoding algorithm for binary sources. The correlation of the binary sources can be viewed as a virtual binary asymmetrical channel (BAC). Many of prior works  can be view as special cases of this model.

%\item Apply the proposed message passing algorithm to solve the generalized binary coding problem \textemdash binary asymmetrical coding (BAC). Many of the prior works, such as BSC estimation and Z-Channel estimation, can be viewed as special cases of this model.   
\item Build a temperature network and propose a new factor graph for this network. The proposed algorithm achieves a much better performance on the network data.
\item First apply the DSC to the genome data compression and decompression. The proposed algorithms have been applied to two genome sequence datasets - the Arabidopsis Information Resource (TAIR) \cite{websiteTAIR} and the Institute for Genomic Research (TIGR) \cite{websiteTIGR}. The proposed algorithm achieves a better performance than that of the GRS algorithm. 
\end{itemize}

Among all of above contributions, the DSC for temperature network and the DSC for genome data have already been published in \cite{chen2013compression,wang2014streamlined}.

\section{Organization}
This dissertation is organized as follows: The basic concepts of Bayesian inference, a factor graph, message passing algorithms and deterministic approximation inference  will be reviewed in Chapter \ref{chap:BayesianInference}. The BAC model will be studied in Chapter \ref{chap:BAC}. Chapter \ref{chap:WZ} is about DSC for wireless temperature networks. In Chapter \ref{chapter:genome}, the DSC for the genome data algorithm will be proposed and tested on two datasets. Chapter \ref{chap:conclusion} provides a conclusion for this dissertation.

\chapter{BAYESIAN INFERENCE}\label{chap:BayesianInference}
A Bayesian technique considers all possible modes when inferring an outcome. In this chapter, the general idea of Bayesian inference will be explained. Then, the prior selection criteria for the Bayesian inference will be introduced. Finally, the factor graph and the inference methods in the factor graph will be introduced. 
\section{Introduction}
First, let's classify some notations.
\begin{itemize}
\item $p(\cdot)$ is the probability distribution function.
\item The upper case letters $X$, $Y$, and $Z$ denote random variables or the variable nodes of the factor graph, and the lower case letters $x$, $y$, and $z$ denote their corresponding values.
\item $\theta$ is the parameter of the distribution model. 
\item $\mathbf{x}=\{x_1, x_2, \dots, x_N\}$ denotes the set of all observation data.
\end{itemize}

Suppose one wants to estimate the parameter of the model $P(x|\theta)$. From the point of view of classical probability, there exists a ground truth value of $\theta$. Given an observation set $\mathbf{x}$, one common approach to estimate $\theta$ is to maximize the likelihood (MLE) function $p(\mathbf{x}|\theta)$, i.e. 

\begin{equation}
\hat{\theta}_{MLE} = \arg\max_{\theta} p(\mathbf{x}|\theta). 
\end{equation}

Differing from classical probability, the parameter $\theta$, under the Bayesian framework, is treated as a random variable. According to Bayes' theorem, the posterior distribution of $\theta$ can be written as
\begin{equation}
p(\theta|\mathbf{x}) = \frac{p(\theta)p(\mathbf{x}|\theta)}{p(\mathbf{x})}, \label{PosteriorDistributionForm0}
\end{equation}
where $p(\theta)$ is the prior distribution, and $p(\mathbf{x})$ is the constant and called as evidence, which can be calculated as
\begin{equation}
p({\mathbf{x}}) = \int p(\theta) p(\mathbf{x}|\theta) d\theta.
\end{equation}
Therefore, one has the expression
\begin{equation}
\mbox{posterior} \propto \mbox{prior} \times \mbox{likelihood}. \label{PosteriorDistributionForm1}
\end{equation}

The prior reflects the existing knowledge of $\theta$ before the data $\mathbf{x}$ are sampled. As a matter of fact, the form of (\ref{PosteriorDistributionForm0}) is suitable for sequentially updating: When the new data are obtained, the previous posterior can be treated as the prior and then multiplied with the likelihood of the new data as~\eqref{PosteriorDistributionForm0} to update the posterior.

By maximizing a-posterior (MAP) as equation (\ref{PosteriorDistributionForm0}), one can obtain another estimate of $\theta$, i.e.,
\begin{equation}
\begin{split}
\hat{\theta}_{MAP} &= \arg\max_{\theta} p(\theta|{\mathbf{x}})=  \arg\max_{\theta} \frac{p({\mathbf{x}}|\theta)p(\theta)}{p(\mathbf{x})} \\
&=\arg\max_{\theta} p({\mathbf{x}}|\theta)p(\theta)
\end{split}
\end{equation}

Both $p(x|\theta_{MLE})$ and  $p(x|\theta_{MAP})$ can be used to make an inference on $x$. For a large number of the training set $\mathbf{x}$, $\theta_{MLE}$ and $\theta_{MAP}$ tend to converge to the same value; but for a limited number of the training set, $\theta_{MLE}$ may achieve better performance if the prior is appropriate. 

\delete{For example, if we toss a coin twice and get all heads, then can obtain $p(head|\theta_{MLE}) = 1$. This predication is obviously rather extreme. But we if we set the prior to $p(head) =p(tail)= 0.5$ by life experience or physical model, then $p(head|\theta_{MLE}) = 0.75$ for the same tossing result.}

As with the MAP method, the Bayesian inference techniques utilize a "posterior". However, the watershed between Bayesian and non-Bayesian inference is not drawn on whether a prior is used. A Bayesian inference does not first find the optimum parameter. Instead, it integrates all of the model parameters out, based on the belief (probability) of occurrence of each model. More precisely, the inference probability of ${\tilde{x}}$ is
\begin{equation}
p(\tilde{x}|\mathbf{x}) = \int p(\tilde{x}|\theta)p(\theta|\mathbf{x})d\theta.
\end{equation}

There are two important factors that determine the application scope of Bayesian inference: One is the prior selection, and another one is the problem of how to efficiently compute the integral (over the parameter) of Bayesian inference.

In the following sections, two prior selection methods will be introduced, and then some of the inference methods will be explained in detail.

\section{Prior Selection}
In this section, two priors will be introduced: the conjugate prior and reference prior. The conjugate prior has the advantage of reducing the computational complexity of the posterior. The reference prior is an "uninformative prior" that has a maximum distance from the posterior under Kullback-Leibler divergence.

\subsection{Conjugate Prior}
The conjugate prior has the exact form as the posterior. For example, for the $k$ dimension multinomial distribution
\begin{equation}
p(x|\pmb{\theta})=\frac{\Gamma({\sum_i}x_i+1)}{\prod_i^k\Gamma(x_i+1)} \prod_{i=1}^{k}{\theta_i^{x_i}},
\end{equation}
its conjugate prior Dirichlet distribution has the form
\begin{equation} \label{Equ:BetaFunction}
\mbox{Dir}(x|\pmb{\alpha}) = \frac{1}{\mbox{B}(\mathbf{\pmb{\alpha}})} \prod_{i=1}^{k}{x_i^{\alpha_i-1}}, 
\end{equation} 
where $\mbox{B}(\cdot)$ is the Beta function, which can be written in terms of the gamma function
\begin{equation} \label{Equ:BetaFunction1}
\mbox{B}(\pmb{\alpha})=\frac{\prod_{i=1}^k{\Gamma(\alpha_i)}}{\Gamma{(\sum_{i=1}^K{\alpha_i})}}.
\end{equation}

Since the prior and the posterior have the same form of distribution, to compute the posterior, one only needs to update the parameters of this distribution.

\subsection{Reference Prior}\label{subsec:reference_prior}
Reference priors were firstly introduced by Jose Bernardo \cite{Bernardo79RefPrior}. The idea of reference priors is that the prior should not strongly influence the posterior distribution, and such a prior is called the \textit{uninformative prior}~\cite{bayesian_choice}. Divergences, such as the Kullback-Leibler (KL) divergence, can be used to measure the influence of the priors on the posteriors. Suppose the random variable $X$ is parameterized by ${\Theta}$, and $T\left(X\right)$ is the sufficient statistic of $X$.A reference prior $p\left({\theta}\right)$ should be chosen that has a maximum distance to $p\left(\theta|t\right)$ under KL divergence
\begin{equation}
\int p\left(\theta|t\right) \log \frac{p\left(\theta | t\right)}{p\left(\theta\right)}\;d\theta,
\end{equation}
whose expectation over the distribution of $T$ can be written as
\begin{equation} \label{MutualInformation}
\begin{split}
I\left(\Theta;T\right)&=\int p\left(t\right) \int  p\left(\theta | t\right) \log \frac{p\left(\theta | t\right)}{p\left(\theta\right)}\;d\theta\;dt\\
&= \int \int  p\left(\theta , t\right) \log \frac{p\left(\theta | t\right)}{p\left(\theta\right)}\;d\theta\;dt
\end{split},
\end{equation}
which is the mutual information between $\theta$ and $t$. Under (\ref{MutualInformation}), one wants to choose $p^{*}\left(\theta\right)$ such that
\begin{equation}
p^{*}\left(\theta\right)=\arg\max\limits_{p\left(\theta\right)} I\left(\Theta,T\right).
\end{equation}

For one-dimensional parameters, the reference priors are equal to Jeffreys priors, which have the form
\begin{equation}
p\left({\theta}\right) \propto I\left({\theta}\right)^{\frac{1}{2}}, \label{Equ:Reference_Prior}
\end{equation}
where $I\left({{\theta}}\right)$ is the Fisher information, which can be expressed as
\begin{equation}
I\left({{\theta}}\right)=-\mathbb{E}_{{\theta}}\left[\frac{d^2\log p\left(\mathbf{x}|{\theta}\right)}{d{{\theta}}^2}\right]
\end{equation}.

\section{Inference Methods} \label{sec:Inference_methods}
A factor graph can model the correlation of a large set of variables in an intuitive manner. The variables in the factor graph can be either hidden variables or observable variables, and they will be connected with the factors, which take those variables as the arguments. Based on the factor graph, the message passing algorithm can be implemented, by which the beliefs of all the hidden variables can be derived. In this section, the factor graph will be introduced, and then the message passing algorithm in the factor graph will be derived.

\subsection{Factor Graph}\label{sec2:factor_graph}
\begin{wrapfigure}{r}{0.45\textwidth}
\includegraphics[width=2.5in]{./figures/ExamplesOfFactorGraph.eps}
\vspace{-0.1in}
\caption{An Example of a Factor Graph.}
\vspace{-0.2in}
\label{fig:ExampleOfFactorGraph}
\end{wrapfigure}
A factor graph represents the factorization of the joint distribution and is a bipartite graph composed of two different types of nodes: variable nodes and factor nodes. As one type of graphical model, the correlations among the variables can be easily observed. In a nutshell, variable nodes represent random variables in the problem, whereas factor nodes connecting immediately related variable nodes explain the correlations among them. Each factor node is assigned with a factor function that contains the variables connected to the factor node as arguments. The only assumption of a factor graph is that the joint probability of all the variables should be equal to the product of all the factor functions.

Figure \ref{fig:ExampleOfFactorGraph} shows an example of a factor graph. In this factor graph, there are $4$ variables $\{x_1, x_2, x_4, x_4\}$, and $3$ factors $\{f_{a}, f_{b}, f_{c}\}$. From the factor graph, all of the factors are multiplied to obtain the joint distribution $p(x_1, x_2, x_3, x_4)$ as
\begin{equation}
p(x_1, x_2, x_3, x_4) = \frac{f_a(x_1, x_2)f_b(x_2, x_3, x_4)f_c(x_1,x_3)}{Z},
\end{equation}
where $Z$ is the normalization constant, and can be computed as
\begin{equation}
Z= \sum_{x_1, x_2, x_3, x4} {f_a(x_1, x_2)f_b(x_2, x_3, x_4)f_c(x_1,x_3)}.
\end{equation}

\begin{table*}
\caption{Belief Propagation Algorithm.}
\vspace{0.1 in}
\centering
%\small
\begin{tabular}{p{5in}} \hline 
{\textbf{Inputs:}} $f_1(\mathbf{x}_1)$, $f_2(\mathbf{x}_2)$, \dots, $f_M(\mathbf{x}_M)$.\\
{\textbf{Initialization:}} initialize each message $m_{i\rightarrow a}(x_i)$ $\forall (a, i)$. \\
\textbf{Repeat until all $m_{a \rightarrow i}(x_i)$ converge or reach maximum iterations:}
\begin{itemize}
\item update all the messages $m_{a \rightarrow i}(x_i)$ via Equation (\ref{BP:f2v});
\item update all the messages $m_{i \rightarrow a}(x_i)$ via Equation (\ref{BP:v2f});
\end{itemize}
{\textbf{Output:}} the beliefs $b(x_i)$ $\forall i$ via Equation (\ref{BP:belief}).
\\ \hline 
\end{tabular}
\label{tab:belief_propagation}
\end{table*} 

\subsection{Belief Propagation Algorithm} \label{subsection:BP_algorithm}
In Bayesian inference, one usually needs to obtain the distributions of some specific random variables, which requires one to summarize/integrate the joint distribution over other random variables. However, due to the complex form of joint distribution, direct summations/integrations are usually infeasible. Belief propagation is an efficient algorithm that can obtain the beliefs (marginal distribution) of the random variables by implementing a sum-product algorithm in the factor graph.

\begin{figure}[!htb]
  \centering
  \subfigure[ ]{\label{fig:BP_F2V}\includegraphics[width=0.48\textwidth]{./figures/BP_F2V}}
  \subfigure[ ]{\label{fig:BP_V2F}\includegraphics[width=0.48\textwidth]{./figures/BP_V2F}}
  \subfigure[ ]{\label{fig:BP_Belief}\includegraphics[width=0.38\textwidth]{./figures/BP_Belief}}
  \caption{(a)message from factor to variable; (b) message from variable to factor; (c) variable belief.}
  \label{fig:BP}
\end{figure} 

In the factor graph, the belief propagation algorithm defines two types of message flows: One is the message from the factor to the variable, and another is the message from the variable to the factor. Suppose there are $N$ random variables $\{X_1, \dots, X_N\}$ and $M$ factors $\{f_1({\mathbf{x}}_1), \dots, f_M({\mathbf{x}}_M)\}$, where ${\mathbf{x}}_a\subset\{x_1, \dots, x_N\}$. For the BP algorithm, each message in the factor graph will be updated according to the following rules:
the message from the variable to the factor
\begin{equation}
m_{i \rightarrow a} \propto \prod_{b\in N(i) \backslash a}m_{b \rightarrow i}(x_i),  \label{BP:v2f}
\end{equation}
and the message from the factor to the variable
\begin{equation}
m_{a \rightarrow i} \propto \sum_{\mathbf{x}_a \backslash x_i}\left( f_a(\mathbf{x}_a) \prod_{j \in N(a) \backslash i} {m_{j\rightarrow a}(x_j)}\right),\label{BP:f2v}
\end{equation}
where the notation $N(i)\backslash a$ denotes the set of all neighbors of node $i$ excluding factor node $a$, and $\sum_{\mathbf{x}_a \backslash x_i}$ denotes a sum over all the variables in $\mathbf{x}_a$ except $x_i$. One uses $i$, $j$ as the index of set $\{1, \dots, N\}$ and $a$, $b$ as the index of set $\{1,\dots,M\}$. Moreover, the belief of $x_i$ is
\begin{equation}
b(x_i) \propto \prod_{a \in N(i)} m_{a \rightarrow i} (x_i). \label{BP:belief}
\end{equation}

Details of the BP algorithm are listed in Table \ref{tab:belief_propagation}.
 
\subsection{Variational Bayes (VB)} \label{subsection:VB_algorithm}
The BP algorithm can only compute the discrete variables, or continuous variables with limited forms such as the Gaussian distribution. Even for discrete random variables, the computation complexity will increase exponentially with the size of variable's the parameters. To cope with these complex situations, one can approximate these true distributions. 

Generally, there are two types of approximate methods: the random approximate method and the deterministic approximate method. The random approximate method usually utilizes sampling methods, such as the Gibbs sampling method~\cite{bishop2006pattern} and Metropolis-Hastings~\cite{metropolis1953equation} algorithms. Although the sampling method may theoretically approximate the true distribution with arbitrarily high precision as the number of samples increases, the complexity is usually much higher than that of the deterministic approximate method. In addition, the designs of the sampling methods are very tricky. For example, the Metropolis-Hastings method samples from a self-constructive Markov network, which determines the convergence speed and precision. However, there is no general design method for this Markov network. The deterministic approximate methods commonly have lower complexity. In addition, one may only concentrate on the main aspect of the inference variable, and the deterministic approximate method will meet our demands with high speed.

Based on above reasons, deterministic approximate methods will be adopted to deal with the complex distribution in the factor graph of this dissertation. Two of the deterministic approximate methods will be utilized in this dissertation:  variational Bayes (VB) and expectation propagation (EP).

Suppose one wants to make inferences on the hidden random variable $Z$ given $X$, and there are $N$ observations of $X$, which are $\mathbf{x}=\{x_1,\dots, x_N\}$.  Due to the complex form of the posterior $p(z|x)$, one will find a favorable distribution $q(z)$ to approximate $p(z|x)$ so that make inferences. Generally, one obtains the following form
\begin{equation}
 \mbox{KL}(q||p) = \ln p(x) + \mathcal{L}(q(z)), \label{KL-Divergence}
\end{equation}
where
\begin{equation}
\mbox{KL}(q||p) = \int q(z)\ln \left(\frac{q(z)}{p(z|x)}\right)dz,\label{Minimiazation:KL}
\end{equation}
\begin{equation}
\mathcal{L}(q) = \int q(z) \ln \left(\frac{q(z)}{p(x,z)}\right)dz.\label{Minimization:LQ}
\end{equation}

$\mbox{KL}(q||p)$ is the Kullback-Leibler divergence operation, which measures the similarity of the distribution between $q$ and $p$. $\mbox{KL}(q||p)$ is always non-negative, and is equal to $0$ when $q(x) = p(x)$. Comparing with the form of the $\mbox{KL}$ divergence, it is known that $\mathcal{L}(q)$ is the $\mbox{KL}$ divergence between $q(z)$ and $p(x,z)$, which indicates $\mathcal{L}(q)\ge 0$ and $\mbox{KL}(q||p) \ge \ln p(x) = \mbox{Constant}$. Hence, $\ln p(x)$ is the lower bound of $\mbox{KL}(q||p)$.

The optimization (minimization) function is~\eqref{KL-Divergence}. However, if one knows the exact posterior $p(z|x)$ in~\eqref{Minimiazation:KL}, then one can make inferences on $Z$ directly. Instead, one can minimize (\ref{Minimization:LQ}). In the following, let's suppose $z$ is a vector, i.e. $\mathbf{z}$, and $q(\mathbf{z})$ can be factorized as
\begin{equation}
q(\mathbf{z}) = \prod_{i=1}^M q_i(\mathbf{z}_i),\label{equ:formq}
\end{equation}
where $\mathbf{z}_i \subset \mathbf{z}$ and $i \in \{1,\dots,M\}$. 
Substituting~\eqref{equ:formq} into~\eqref{Minimization:LQ} and minimizing $\mathcal{L}(q)$ with respect to each $q(\mathbf{z}_j)$, one can obtain the optimal solution $q^{*}(\mathbf{z}_j)$ as
\begin{equation}
 q^{*}(\mathbf{z}_j) = \frac{\exp\left(\mathbb{E}_{i\ne j}[\ln p(\mathbf{x}, \mathbf{z})]\right)}{\int \exp\left(\mathbb{E}_{i\ne j}[\ln p(\mathbf{x}, \mathbf{z})]\right) d\mathbf{z}_j},
\end{equation}
where
\begin{equation}
\mathbb{E}_{i\ne j}[\ln p(\mathbf{x}, \mathbf{z})] = \int {\ln p(\mathbf{x}, \mathbf{z})\prod_{i \ne j}q_i d\mathbf{z}_i}. \label{Exp:ln}
\end{equation}

One can see from~\eqref{Exp:ln} that the solutions of all the $q^{*}(\mathbf{z}_j)$ are coupled with the others. To solve this issue, one can iteratively run through each of the $q^{*}(\mathbf{z}_j)$ until they all converge. In fact, this iteration algorithm is guaranteed to coverage~\cite{bishop2006pattern}.

\subsection{Expectation Propagation (EP)} \label{subsection:EP_algorithm}

\begin{table*}
\caption{Expectation Propagation Algorithm.}\label{Algorithm:EP}
\vspace{0.1 in}
\centering
%\small
\begin{tabular}{p{5in}} \hline 
{\textbf{Inputs:}} $t_1(\mathbf{z})$, $t_2(\mathbf{z})$, \dots, $t_N(\mathbf{z})$\\
{\textbf{Initialization:}} initialize $\tilde{t}_1(\mathbf{z})$, $\tilde{t}_2(\mathbf{z})$, \dots, $\tilde{t}_N(\mathbf{z})$; initialize $q(\mathbf{z})$ according to (\ref{EP:q_factorization}).\\
\textbf{Repeat until all $\tilde{t}_i$ converge or reach maximum iterations:}
\begin{itemize}
\item Choose a $\tilde{t}_i$ to refine;
\item Obtain the distribution $q^{\mathbf{\setminus} i}(\mathbf{z})$:
\begin{equation}
q^{\setminus i}(\mathbf{z}) \propto \frac{q(\mathbf{z})}{\tilde{t}_i(\mathbf{z})};
\end{equation}
\item Approximate $q^{\setminus i}(\mathbf{z})\tilde{t}_i(\mathbf{z})$ with new $q(\mathbf{z})$;
\item Compute normalizer $Z_i=\int t_i(\mathbf{z})q^{\setminus i}(\mathbf{z}) d{\mathbf{z}}$;
\item Update $\tilde{t}_i = Z_i q(\mathbf{z})/q^{\setminus i}(\mathbf{z})$ 
\end{itemize}
{\textbf{Output:}} $q(\mathbf{z})$.
\\ \hline 
\end{tabular}
\label{tab:expectation_propagation}
\end{table*}

The joint distribution of $\mathbf{z}$ and $N$ i.i.d. observations $\mathbf{x}=\{\mathbf{x}_1, \dots, \mathbf{x}_N\}$ is given by
\begin{equation}
p({\mathbf{x}},\mathbf{z}) = p(\mathbf{z})\prod_{i=1}^N p(\mathbf{x}_i|\mathbf{z}). \label{EP:JointDstribution}
\end{equation}

Let $t_0(\mathbf{z})=p(\mathbf{z})$ and $t_i(\mathbf{z})=p({\mathbf{x}}_i|\mathbf{z})$. Then (\ref{EP:JointDstribution}) can be written as 
\begin{equation}
p({\mathbf{x}},\mathbf{z}) =\prod_{i=0}^N t_i(\mathbf{z}).
\end{equation}

The EP algorithm approximates the posterior $p(\mathbf{z}|D)$ with $q(\mathbf{z})$ under 
\begin{equation}
\mbox{KL}(p({\mathbf{x}}|\mathbf{z})||q(\mathbf{z})) = \mbox{KL}\left(\frac{\prod_{i=0}^N t_i(\mathbf{z})}{p({\mathbf{x}})}||q(\mathbf{z})\right), \label{EP:KL}
\end{equation}
which has the reverse form of~\eqref{Minimiazation:KL}. However, the form of (\ref{EP:KL}) is intractable. To make~\eqref{Minimiazation:KL} tractable, the EP algorithm factorizes $q\left(\mathbf{z}\right)$ as
\begin{equation}
q\left(\mathbf{z}\right) = \frac{\prod_{i=1}^N \tilde{t}_i(\mathbf{z})}{\int \prod_{i=1}^N \tilde{t}_i(\mathbf{z})}, \label{EP:q_factorization}
\end{equation}
and let each $\tilde{t}_i(\mathbf{z})$ approximate each $t_i(\mathbf{z})$, and an approximate iteration algorithm will be listed in Table \ref{tab:expectation_propagation}~\cite{Minka:01}.

The commonly used distribution $\tilde{t}_i$ is the exponential family distribution
\begin{equation}
f(\mathbf{z}|\theta)=h(\mathbf{z})g(\mathbf{\theta})\exp\left({\mathbf{\theta}}^{T} \mathbf{u}(\mathbf{z})\right).
\end{equation} 

Once $q(\mathbf{z})$ is obtained, one can make inferences on $Z$ and compute $p(\mathbf{x})$ as
\begin{equation}
p({\mathbf{x}}) \approx \int {q(\mathbf{z})} d\mathbf{z}.
\end{equation}

\chapter{ADAPTIVE SLEPIAN-WOLF DECODING FOR TWO BINARY SOURCES}\label{chap:BAC} 
\section{Introduction}
The DSC with two sources problem can be modeled as a channel coding problem: The two sources of the DSC can be viewed as the input and output of a channel; The correlation of the two sources can be viewed as the property of the channel. \delete{The two binary sources DSC problem is usually modeled as the Binary symmetrical channel (BSC), which can be parameterized by the crossover probability. However, the practical correlation of two sources may be well modeled by the BSC.} In this chapter, the correlation of the two binary sources will be modeled as the Binary asymmetrical channel (BAC) problem, which has a wide application scope. 

\delete{ For binary channel communications, one can capture two transmission errors: $0$ is received as a $1$, and vice versa \cite{IT}. Suppose the probability of the two errors are $pr(1\rightarrow 0)$ and $pr(0 \rightarrow 1)$. If $pr(1\rightarrow 0)$ equals $pr(0\rightarrow 1)$, this channel is called binary symmetric channel (BSC). BSC is one of the most common channels in research community. However, the practical channel may not satisfy this symmetric condition. In this chapter, I will focus on the binary asymmetric channel (BAC), which does not assume any relation between $pr(1\rightarrow 0)$ and $pr(0\rightarrow 1)$. The BSC can be viewed as a special case of BAC.}    

\subsection{Related Work}
Since the DSC problem can be modeled as a channel coding problem, the algorithms of the channel problem can be directly applied to the DSC problem. In~\cite{SchonbergRP:04,pradhan2005gcc}, the DSC problem was modeled as a virtual BSC channel. To estimate the parameter (correlation in the DSC) in the BSC channel, various methods have been proposed. In  \cite{fang2009crossover,zia2007distributed,toto2011maximum}, the crossover probability of the BSC is assumed to be constant in the context of the source coding. \cite{cui2011adaptive,cuionline,wang2011noiseGlobeCOM} make further progress and assume that the crossover probability changes over time. However, if the correlation of the two sources does not satisfy the symmetrical condition, the BSC model will lead to a poor result. To be more applicable in practice, the symmetrical condition is not assumed in this section, and a model of the correlation of the two sources as the binary asymmetrical channel (BAC) is used. 

In the rest of this chapter, an adaptive DSC algorithm for two binary sources modeled as BAC is proposed. First, the BAC model is introduced. Then, the proposed algorithms are explained in detail; Finally, the experimental results and an analysis of the proposed algorithms are given.

\delete{
Compared the prior works~\cite{fang2009crossover,zia2007distributed,toto2011maximum,cui2011adaptive,cuionline,wang2011noiseGlobeCOM}, the proposed method is more adaptive and practical.}

\section{Problem Formulation}
\begin{figure}[!htb]
	
	\centering
		\includegraphics[width=5in]{./figures/BAC}
\caption{Binary asymmetric channel}
	\label{Fig:BAC}
\end{figure}

In Figure~\ref{Fig:BAC}, suppose the BAC channel has a probability $\sigma_0$ such that an input $0$ will be flipped into a $1$ and probability $\sigma_1$ such that an input $1$ will be flipped into a $0$. 

Unlike with the BSC~\cite{IT}, the conditional probability of the received bit value, given the send bit $P(y|x)$, can be expressed as the following piecewise function:
\begin{equation} \label{Equ:Corr_x_y}
f(x, y, \sigma_0, \sigma_1) = P(y|x, \sigma_0, \sigma_1) = \left\{
		\begin{array}{ll}
			\sigma_0^y(1-\sigma_0)^{(1-y)}&: x = 0 \\
			\sigma_1^{1 - y}(1-\sigma_1)^{y}&: x = 1
		\end{array}.
	\right.
\end{equation}
The above piecewise function can be unified as
\begin{equation}
f(x, y, \sigma_0, \sigma_1) = \sigma_x^{x \oplus y}\left(1-\sigma_x\right)^{1\oplus x \oplus y}. \label{equ:BAC_Contional}
\end{equation}

Suppose
\begin{equation} \label{Equ:Corr_x_y}
p(x) = \left\{
		\begin{array}{ll}
			\pi_0&: x = 0 \\
			\pi_1&: x = 1
		\end{array},
	\right.
\end{equation}
then the joint distribution $p(x,y)$ is given as
\begin{equation}
\begin{split}
p(x,y)&=p(x)p(y|x) \\
&=\pi_x\sigma_x^{x \oplus y}\left(1-\sigma_x\right)^{1\oplus x \oplus y}\\
&=\left(\pi_0\sigma_0^y(1-\sigma_0)^{(1-y)}\right)^{(1-x)} \left(\pi_1\sigma_1^{(1-y)}(1-\sigma_1)^{y}\right)^{x}.
\end{split}.\label{Equ:BACJointD}
\end{equation}


\begin{figure}[!htb]
	
	\centering
		\includegraphics[width=6in]{./figures/FactorGraphOfBAC.eps}
\caption{Factor graph of binary asymmetric channel (BAC) sources} 
	\label{BAC:FactorGraphBAC}
\end{figure}

In the following, the correlation between the sources $X$ and $Y$ is modeled as the BAC channel in Equation~\eqref{equ:BAC_Contional}, and an adaptive source coding method will be proposed. 

For simplicity, it is assumed that both $\pi_0$ and $\pi_1$ in~\eqref{Equ:BACJointD} are known in advance, and that both are equal to $0.5$. However, the proposed VB-based method will assume that $\pi_0$ and $\pi_1$ are unknown and will estimate them, and priors will be given for this algorithm for Bayesian inference.

\section{Proposed Methods} \label{section:bac_proposed_method}
To encode the source data $\mathbf{x}$ (vector form), the low-density parity-check code (LDPC) is utilized. LDPC is one of the linear block codes, which can be defined by a sparse parity-check matrix. Given an LDPC parity-check matrix  $\mathbf{H}$, the syndrome $\mathbf{s}$ (codeword) can be easily obtained by $\mathbf{s} = \mathbf{H}\times \mathbf{x}$.
 
To decode the syndromes, a factor graph (see section \ref{sec2:factor_graph}) is first proposed in Figure \ref{BAC:FactorGraphBAC}. As shown in Figure~\ref{BAC:FactorGraphBAC}, the factor graph can be divided into three regions: Region I is the low-density parity-check code (LDPC) area, which encodes the constraints between the sources $\{x_1, \dots, x_L\}$ and the compressed syndromes $s_1, \dots, s_K$. Region II is the correlation between the sources $\{x_1, \dots, x_L\}$ and the side information (SI) $\{y_1, \dots, y_L\}$. The BAC channel parameters $\sigma_0$ and $\sigma_1$ are estimated in the Region III. In the following sections, each of these three regions will be explained one by one.

\subsection{Region I - LDPC Codes} \label{BAC_subsection:LDPC}
In Region I of Figure \ref{BAC:FactorGraphBAC}, there are $L$-bit source and $M$-bit syndrome, which are connected by the check nodes $\{h_1,\dots,h_M\}$, and the sparse parity-check matrix $\mathbf{H}$ has the size of $L \times M$, i.e., $\mathbf{H} \in \mathbb{F}_2^{L \times M}$.

%For each entry $h_{i,j}$ of $H$, if the $j$-th bit of the source is connected with the $i$-th syndrome, then $h_{i,j}=1$; otherwise $h_{i,j}=0$.

To satisfy the LDPC constraint, each factor $h_a$ ($a \in \{1,\dots,M\}$) is defined as
\begin{equation}
h_a({\mathbf{x}}_a, s_a)=\sum{{\mathbf{x}}_a}+s_a=0. \label{Equ:h_function}
\end{equation}
All of the summations in Equation~\eqref{Equ:h_function} are implemented in $\mathbb{F}_2$, and ${\mathbf{x}}_a$ denotes all of the sources that connect with the factor $h_a$. The vectorized form of~\eqref{Equ:h_function} can be written as
\begin{equation}
\mathbf{H}\times \mathbf{x} + \mathbf{s} = \mathbf{0} \label{Equ:h_function_vectorized}
\end{equation}

The BP algorithm (see section~\ref{subsection:BP_algorithm}) will be applied to solve the LDPC decoding, and the detailed steps are listed in Table \ref{tab:algorithm_LDPC}. 

\begin{table*}
\caption{Summary of LDPC decoding using BP}
\vspace{0.1 in}
\centering
%\small
\begin{tabular}{p{5in}} \hline 
{\textbf{Input:}} the syndromes $s_1,\dots,s_K$.\\
{\textbf{Initialization:}} 
\begin{equation}
L_{ai} = 0, \forall \mbox{ check node } a
\end{equation}. \\
\textbf{Repeat until all messages $m_{a \rightarrow i}(x_i)$ and  $m_{i \rightarrow a}(x_i)$ are converged or reach maximum number of iterations:}
\begin{itemize}
\item update all the messages $m_{i \rightarrow a}(x_i)$ via Equation (\ref{Equ:LDPC_V2F});
\item update all the messages $m_{a \rightarrow i}(x_i)$ via Equation (\ref{Equ:LDPC_F2V});
\end{itemize}
{\textbf{Output:}}  the log-likelihood ratio for variable node $i$ can be updated according to Equation (\ref{Equ:LDPC_Belief}).
\\ \hline 
\end{tabular}
\label{tab:algorithm_LDPC}
\end{table*} 

Initially, no information is obtained from the check factor node. Therefore the message from those nodes should bear no information, i.e., $m_{ai}(1)=m_{ai}=0.5$ for any factor node $a$ to variable node $i$. 

After assigning the initial values to the messages from factor to variable, we can update the messages from variable to factor and the messages from factor to variable alternatively as described in section~\ref{subsection:BP_algorithm}. Since the unknown variables are binary, it is more convenient to represent the messages using likelihood or log-likelihood ratios. Define
\begin{equation} \label{equ:ratioform}
\begin{split}
l_{ai} \triangleq \frac{m_{ai}(0)}{m_{ai}(1)}\\
L_{ai} \triangleq \log l_{ai}
\end{split},
\end{equation}
and
\begin{equation}
\begin{split}
l_{ia}\triangleq \frac{m_{ia}(0)}{m_{ia}(1)}\\
L_{ia} \triangleq \log l_{ia}
\end{split}
\end{equation}
for any variable node $i$ and factor node $a$.

Therefore, the message from variable node $i$ to check node $a$ can be updated as
\begin{equation} \label{Equ:LDPC_V2F}
L_{ia} \leftarrow \sum_{b \in N(i)\backslash a} L_{bi}.
\end{equation}

For the message from the check node $a$ to the variable node $i$, according to Equation~\eqref{Equ:h_function}, a simplified update rule can be divided as follow
\begin{equation} \label{Equ:LDPC_F2V}
L_{ai} \leftarrow \left\{
		\begin{array}{ll}
			2\tanh^{-1}\left(\prod_{j \in N(a)\backslash i}\tanh\left(\frac{L_{ja}}{2}\right)\right)&: s_a = 0 \\
			2\tanh^{-1}\left(-\prod_{j \in N(a)\backslash i}\tanh\left(\frac{L_{ja}}{2}\right)\right)&: s_a = 1
		\end{array}.
	\right.
\end{equation}

The belief of variable node $i$ can be updated as 
\begin{equation} \label{Equ:LDPC_Belief}
\beta_i = \sum_{a \in N(i)} L_{ai}.
\end{equation}

Hard thresholding can be applied to the final belief to estimate the value of each bit. That is,
 \begin{equation} 
\hat{x}_i = \left\{
		\begin{array}{ll}
			0&: \mbox{if } \beta_i \ge 0 \\
			1&: \mbox{otherwise}
		\end{array}.
	\right.
\end{equation}
\subsection{Region II - Correlation between the sources and the side information}
\delete{The correlation between the sources and the SI is Equation~\eqref{equ:BAC_Contional}. }
For the messages between the Region \Rmnum{2} and Region \Rmnum{3}, the deterministic approximate inference methods - the EP and VB algorithms will be applied to update these messages. The details of the updating rules will be derived in the next subsection. In this subsection, only the message sent from the Region \Rmnum{2} to Region \Rmnum{1}, i.e. $m_{f_i \rightarrow X_i }(x_i)$, will be derived.

For both of the deterministic algorithms, each message sent from the Region \Rmnum{3} to Region \Rmnum{2} has the form of Beta distribution. For example, the messages $m_{\Sigma_0^j \rightarrow f_i}(\sigma_0)$ and $m_{\Sigma_1^j \rightarrow f_i}(\sigma_1)$ can be given as
\begin{equation}
\begin{split}
m_{\Sigma_0^j \rightarrow f_i}(\sigma_0) = \mbox{Beta}(\sigma_0|\alpha_{0,i}^j, \beta_{0,i}^j)= \frac{1}{\mbox{B}(\alpha_{0,i}^j, \beta_{0,i}^j)} \sigma_0^{\alpha_{0,i}^j - 1} (1- \sigma_0)^{\beta_{0,i}^j - 1}\\
m_{\Sigma_1^j \rightarrow f_i}(\sigma_1) = \mbox{Beta}(\sigma_1|\alpha_{1,i}^j, \beta_{1,i}^j)= \frac{1}{\mbox{B}(\alpha_{1,i}^j, \beta_{1,i}^j)} \sigma_1^{\alpha_{1,i}^j - 1} (1- \sigma_1)^{\beta_{1,i}^j - 1}
\end{split},
\end{equation}
where $\mbox{B}(\cdot)$ is defined in~\eqref{Equ:BetaFunction1}.

Then according to Equation~\eqref{BP:f2v}, we have
\begin{equation} \label{BACequ:msg_2_3}
\begin{split}
m_{f_i \rightarrow X_i }(x_i) &= \int m_{\Sigma_0^j \rightarrow f_i}(\sigma_0) m_{\Sigma_1^j \rightarrow f_i}(\sigma_1) f(x,y,\sigma_0,\sigma_1)d\sigma_0 d\sigma_1 \\
&=\int \mbox{Beta}(\sigma_0) \mbox{Beta}(\sigma_1) \sigma_x^{x \oplus y}\left(1-\sigma_x\right)^{1\oplus x \oplus y}d\sigma_0 d\sigma_1
\end{split}.
\end{equation}

By using the following properties
\begin{equation}
\begin{split}
\mbox{B}(\alpha,\beta) &= \int \mbox{Beta}(\sigma|\alpha, \beta) d \sigma \\
\mbox{B}(\alpha+1,\beta) &= \mbox{B}(\alpha,\beta)\frac{\alpha}{\alpha+\beta}\\
\mbox{B}(\alpha,\beta+1) &= \mbox{B}(\alpha,\beta)\frac{\beta}{\alpha+\beta}
\end{split},
\end{equation}

The ratio form (see~\eqref{equ:ratioform}) of $m_{f_i \rightarrow x_i }(x_i)$ in Equation~\eqref{BACequ:msg_2_3} can be simplified as
\begin{equation}
m_{f_i \rightarrow X_i }(x_i) = \left\{
		\begin{array}{ll}	\frac{\beta_{0,i}}{\alpha_{1,i}}\frac{\alpha_{1,i}+\beta_{1,i}}{\alpha_{0,i}+\beta_{0,i}}: y = 0 \\
\frac{\alpha_{0,i}}{\beta_{1,i}}\frac{\alpha_{1,i}+\beta_{1,i}}{\alpha_{0,i}+\beta_{0,i}}: y = 1
		\end{array}.
	\right.
\end{equation}

\subsection{Region III}
To make the algorithm more adaptive, in Figure \ref{BAC:FactorGraphBAC}, the channel parameters $\sigma_0$ and $\sigma_1$ are assumed to be unknown, and both of these parameters will be estimated. In addition, it is assumed that the parameters $\sigma_0$ and $\sigma_1$ are changing over time $j$, i.e., $\sigma_0^j$ and $\sigma_1^j$. Here, the number of factor nodes in Region II connecting to each variable node $\sigma_k^j$ is called the connection ratio $C$, where $k \in \{0, 1\}$. In Figure \ref{BAC:FactorGraphBAC}, it is seen that  the connection ratio is $3$.  

The Bayesian inference requires one to estimate  the posterior distribution $p(\sigma_k^j | {\mathbf{y}}_j)$, where ${\mathbf{y}}_j = \{y_i|i \in N^{\backslash g_k^j}(\Sigma_k^j)\}$, and $N^{\backslash g_k^j}(\Sigma_k^j)$ denotes the set of all neighboring indices for a variable node $\sigma_k^j$ excluding $g_k^j$. The posterior is given as

\delete{\begin{equation} \label{Equ:sigma_posterior}
\begin{split}
p(\sigma_k^j|{\mathbf{y}}_j) &= \frac{1}{Z_j} \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} p(\sigma_k^j) p(y_i|\sigma_k^j) \\
&= \frac{1}{Z_j} \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} \int p(\sigma_k^j) p(\sigma_{1-k}^j) p(x_i) p(y_i|x_i,\sigma_k^j, \sigma_{1-k}^j) dx_i\ d\sigma_{1-k}^j \\
\delete{&= \frac{1}{Z_j} \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} \int p(\sigma_k^j) p(\sigma_{1-k}^j) p(x_i) p(y_i|x_i,\sigma_k^j) dx_i\ d\sigma_{1-k}^j \\}
&= \frac{1}{Z_j} g(\sigma_k^j)\prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} \sum_{x_i}{m_{X_i \rightarrow f_i}(x_i)}\int p(y_i|x_i,\sigma_i^j, \sigma_{1-k}^j) p(\sigma_{1-k}^j) d\sigma_{1-k}^j,
\end{split} 
\end{equation}}

\begin{equation} \label{Equ:sigma_posterior}
\begin{split}
p(\sigma_k^j|{\mathbf{y}}_j) &= \frac{1}{Z_j} \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} p(\sigma_k^j) p(y_i|\sigma_k^j) \\
&= \frac{1}{Z_j} \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} \sum_{x_i} p(\sigma_k^j) p(x_i) p(y_i|x_i,\sigma_k^j)\\
&= \frac{1}{Z_j} g(\sigma_k^j)\prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} \sum_{x_i}{m_{X_i \rightarrow f_i}(x_i)}  p(y_i|x_i,\sigma_k^j),
\end{split} 
\end{equation}  
where $Z_j=\int \prod_{i \in \mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)} p(\sigma_k^j) p(y_i|\sigma_k^j)d\sigma_k^j$ is a normalization constant, $p(\sigma_k^j)=g(\sigma_k^j)$, and $p(x_i)$ is captured by the message $m_{X_i \rightarrow f_i}(x_i)$ with binary sources $x_i$ taking $0$ or $1$. The connection ratio $C$ equals to $|\mathcal{N}^{\backslash g_k^j}(\Sigma_k^j)|$.

According to Equation~\eqref{Equ:sigma_posterior}, the exact form of $p(\sigma_k^j | {\mathbf{y}}_j)$ contains $2^C$ terms, which are infeasible to estimate when $C$ is large, say $50$. Considering the trade off between complexity and accuracy, deterministic approximate methods are introduced. Two deterministic approximate methods are adopted, and they are  Expectation propagation (EP) (see Section~\ref{subsection:VB_algorithm}) and Variation Bayesian (VB) (see Section~\ref{subsection:EP_algorithm}). In the following, I will illustrate how to apply both of them to estimate the BAC parameters and derive each of the both updating equations for message passing.

\subsubsection{Estimating parameters by EP}
Only the estimate method for parameter $\sigma_0^j$ ($1 \le j \le J$) will be provided, and the estimate method for $\sigma_1^j$ ($1 \le j \le J$) can follow the same manners. In addition, all of the following steps are generalized for all $j$, hence the index $j$  will be omitted\delete{ in all of the symbols or notations}. For example, $\sigma_0^j$ is shorted as $\sigma_0$.

For $\sigma_0$, the Equation~\eqref{Equ:sigma_posterior} can be rewritten as the following
\begin{equation}
\begin{split}
p(\sigma_0|{\mathbf{y}})
&= \frac{1}{Z} m_{g_0 \rightarrow \Sigma_0}(\sigma_0) \prod_{i \in \mathcal{N}^{\backslash g_0}(\Sigma_0)} m_{f_i \rightarrow \Sigma_0}(\sigma_0).
\end{split} 
\end{equation}

The generalized EP algorithm is listed in Table \ref{Algorithm:EP}. For this specific problem, $q_i(\sigma_0)$ will be used to estimate each  $m_{f_i \rightarrow \Sigma_0}(\sigma_0)$, and confine all of $q_i(\sigma_0)$ into beta distribution. 

\delete{
The beta distribution is defined as
\begin{equation} \label{Equ:beta_distribution}
\begin{split}
p(x|\alpha, \beta) &= \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1}\\
&= \frac{1}{\mbox{B}\left(\alpha,\beta\right)}x^{\alpha-1}(1-x)^{\beta-1},
\end{split}
\end{equation} 
where $\Gamma(\cdot)$ is gamma function, $B(\cdot,\cdot)$ is beta function, $0\le x \le 1$ and $\alpha,\beta > 0$.} %The likelihood function (\ref{equ:BAC_Contional)} can be expressed in term of beta distribution with $\sigma_x$ as
%\begin{equation}
%p(y|x,\sigma_x)=\sigma_x^{x\}
%\end{equation}

The prior $g(\sigma_0)$ is given as the reference prior. Since the beta distribution is used to estimate the likelihood function (\ref{equ:BAC_Contional}), according to Equation~\eqref{Equ:Reference_Prior}, $g(\sigma_0)$ can be calculated as
\begin{equation} \label{Equ:Reference_prior_beta_distribution}
g(\sigma_0) = \frac{1}{B(\alpha_0,\beta_0)} \sigma_0^{\alpha_0-1}(1-\sigma_0)^{\beta_0-1},
\end{equation}
where $\alpha_0=\beta_0=0.5$. In addition, it can be seen that the reference prior (\ref{Equ:Reference_prior_beta_distribution}) has the same form as conjugate prior of beta distribution, which will lower the computational complexity.

In the process of estimating $\sigma_0$, the approximate distribution for $\sigma_1$ is fixed, and let's suppose
\begin{equation}
q(\sigma_1)=\frac{1}{\mbox{B}(\alpha_1,\beta_1)}\sigma_1^{\alpha_1-1}
\left(1-\sigma_1\right)^{\beta_1-1}.
\end{equation}

The following is the detail of the proposed EP algorithm for estimating $\sigma_0$. 

\begin{enumerate}
\item Initialize the prior $g_0(\sigma_0)$ according to Equation~\eqref{Equ:Reference_prior_beta_distribution}.
\item Initialize the likelihood messages as
\begin{equation}\label{Equ:Initialize_likelihood_messages}
\delete{\tilde{m}_{f_i\rightarrow \Sigma_0}\left(\sigma_0\right)}q_i(\sigma_0)= z_{0,i}\sigma_0^{\alpha_{0,i}-1}\left(1-\sigma_0\right)^{\beta_{0,i}-1}
\end{equation}
with $\alpha_{0,i}=1$, $\beta_{0,i}=1$ and $z_{0,i}=1$.\delete{, where the values selection method for the above parameters guarantees that $\sigma_0$ is equally likely before LDPC decoding.}
\item Initialize the posterior distribution as
\begin{equation}
q(\sigma_0)= g_0(\sigma_0).
\end{equation}
\item
For $1$ to the maximum number of iteration

For $i\in \mathcal{N}^{\backslash g_0} \left(\sigma_0\right)$
\begin{enumerate}
\item Remove \delete{$\tilde{m}_{f_i \rightarrow \sigma_0}\left(\sigma_0\right)$} $q_i(\sigma_0)$ from the posterior $q\left(\sigma_0\right)$
\begin{equation}
\begin{split}
\alpha_0^{'}=\alpha_0-\left(\alpha_{0,i}-1\right);\\
\beta_0^{'}=\beta_0-\left(\beta_{0,i}-1\right).
\end{split}
\end{equation}
\item Update $\alpha_0$ and $\beta_0$ according to moment matching method, which are 
\begin{equation} 
\alpha_0=\frac{m_1\left(m_1-m_2\right)}{m_2-m_1^2};
\end{equation}
and
\begin{equation} 
\beta_0=\alpha_0\left(\frac{1}{m_1}-1\right),
\end{equation}
where
\begin{equation}\label{Equ:sigma1_m1}
m_1=\frac{P_1}
{Q_1},
\end{equation}

\begin{equation}\label{Equ:sigma1_m2}
m_2=\frac{P_2}
{Q_2},
\end{equation}

\begin{equation}
\begin{split}
P_1 = \left(\alpha_0^{'}+y\right)\left(1+\left(\frac{\alpha_1}{\beta_1}\right)
^{2y-1}\right)Lr\left(x\right)+\left(\frac{\alpha_0^{'}}{\beta_0^{'}}\right)^
{1-y}\left(\alpha_0^{'}+\beta_0^{'}+1\right)
\end{split}
\end{equation}

\begin{equation}
Q_1 = \left(\alpha_0^{'}+\beta_0^{'}+1\right)\left(\left(1+\left(\frac{\alpha_1}
{\beta_1}\right)^{2y-1}\right)Lr\left(x\right)+1+\left(\frac{\beta_0^{'}}
{\alpha_0^{'}}\right)^{2y-1}\right),
\end{equation}

\begin{equation}
\begin{split} 
P_2= & \left(\alpha_0^{'}+y\right)\left(\alpha_0^{'}+y+1\right)\left(1+\left(\frac{\alpha_1}{\beta_1}\right)
^{2y-1}\right)Lr\left(x\right)+\left(\frac{\alpha_0^{'}}{\beta_0^{'}}\right)^
{1-y}\\
&\left(\alpha_0^{'}+1\right)\left(\alpha_0^{'}+\beta_0^{'}+2\right)\\,
\end{split}
\end{equation}

\begin{equation}
\begin{split} 
Q_2= & \left(\alpha_0^{'}+\beta_0^{'}+1\right)\left(\alpha_0^{'}+\beta_0^{'}+2\right)\\
&\left(\left(1+\left(\frac{\alpha_1}
{\beta_1}\right)^{2y-1}\right)Lr\left(x\right)
+1+ 
\left(\frac{\beta_0^{'}}
{\alpha_0^{'}}\right)^{2y-1}\right),\\
\end{split}
\end{equation}

and $Lr(x)=\frac{m_{X_i \rightarrow f_i(0)}}{m_{X_i \rightarrow f_i(1)}}$
\item Set approximated message
\begin{equation}
\begin{split}
\alpha_{0,i} &= \alpha_0-\left(\alpha^{'}-1\right),\\
\beta_{0,i} &= \beta_0-\left(\beta^{'}-1\right).\\
\end{split}
\end{equation}
\end{enumerate}
End For

End For

\item
Output $q(\sigma_0)$
\end{enumerate}

Then one can follow the similar above procedures to estimate $\sigma_1$, and all of the procedures are the same except the formulas of $P_1, Q_1$ in (\ref{Equ:sigma1_m1}), and $P_2, Q_2$ in (\ref{Equ:sigma1_m2}) should be modified as

\begin{equation}
P_1 = \left(\frac{\alpha_1^{'}}{\beta_1^{'}}\right)^y\left(\alpha_1^{'}+
\beta_1^{'}+1\right)Lr\left(x\right)+\left(\alpha_1^{'}+1-y\right)\left(1+
\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right),
\end{equation}

\begin{equation}
Q_1 = \left(\alpha_1^{'}+
\beta_1^{'}+1\right)\left(\left(1+\left(\frac{\alpha_1^{'}}{\beta_1^{'}}
\right)^{2y-1}\right)Lr(x)+1+\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right),
\end{equation}

\begin{equation}
\begin{split}
P_2 &= \left(\frac{\alpha_1^{'}}{\beta_1^{'}}\right)^y\left(\alpha_1^{'}+1\right)
\left(\alpha_1^{'}+
\beta_1^{'}+2\right)Lr\left(x\right)+\left(\alpha_1^{'}+1-y\right)\\
&\left(\alpha_1^{'}+2-y\right)\left(1+
\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)
\end{split}
\end{equation}
and
\begin{equation}
\begin{split}
Q_2 = &\left(\alpha_1^{'}+
\beta_1^{'}+1\right)\left(\alpha_1^{'}+
\beta_1^{'}+2\right)\\
&\left(\left(1+\left(\frac{\alpha_1^{'}}{\beta_1^{'}}
\right)^{2y-1}\right)Lr(x)+1+\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)
\end{split}.
\end{equation}

\delete{
\begin{equation}
m_1=\frac{\left(\frac{\alpha_1^{'}}{\beta_1^{'}}\right)^y\left(\alpha_1^{'}+
\beta_1^{'}+1\right)Lr\left(x\right)+\left(\alpha_1^{'}+1-y\right)\left(1+
\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)}{\left(\alpha_1^{'}+
\beta_1^{'}+1\right)\left(\left(1+\left(\frac{\alpha_1^{'}}{\beta_1^{'}}
\right)^{2y-1}\right)Lr(x)+1+\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)}
\end{equation}
 and
\begin{equation}
m_2=\frac{\left(\frac{\alpha_1^{'}}{\beta_1^{'}}\right)^y\left(\alpha_1^{'}+1\right)\left(\alpha_1^{'}+
\beta_1^{'}+2\right)Lr\left(x\right)+\left(\alpha_1^{'}+1-y\right)\left(\alpha_1^{'}+2-y\right)\left(1+
\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)}{\left(\alpha_1^{'}+
\beta_1^{'}+1\right)\left(\alpha_1^{'}+
\beta_1^{'}+2\right)\left(\left(1+\left(\frac{\alpha_1^{'}}{\beta_1^{'}}
\right)^{2y-1}\right)Lr(x)+1+\left(\frac{\beta_0}{\alpha_0}\right)^{2y-1}\right)}
\end{equation}
respectively.}

The $\sigma_0$ and $\sigma_1$ can be estimated alternately until both  $\sigma_0$ and $\sigma_1$ converge or reach the maximum number of iteration.

\subsubsection{Estimating parameters by VB}
The VB algorithm is also utilized to estimate $\sigma_0$ and $\sigma_1$. Before continuing, some notations should be classified:
\begin{itemize}
\item $N$ is the number of data (connection ratio in Region III of Figure \ref{BAC:FactorGraphBAC});
\item Define  $1$-of-$2$ random variables $\mathbf{x}$ as $\mathbf{x}=\left(1-x, x\right)^{T}$, where $x$ is the values of $X$ in Figure \ref{BAC:FactorGraphBAC}. It can be seen that only one element of $\mathbf{x}$ equals to  $1$ and the other one equals to $0$. $x_{nk}$ is the $k-th$ element of $\mathbf{x}$ for data $n$. Note $k=\{0,1\}$;
\item
let $\tilde{\sigma_0}=\sigma_0$, $\tilde{\sigma_1}=1-\sigma_1$, and $\tilde{\sigma} = [\tilde{\sigma_0}, \tilde{\sigma_1}]^T$. 
\end{itemize}

Considering the dependent relations among $\mathbf{x}$, $y$, $\pi$ and $\tilde{\sigma}$, their joint distribution can be factorized as
\begin{equation}\label{Equ:VB_PrimaryFunc}
p(\mathbf{x},y,\pi, \tilde{\sigma})=p(y|\mathbf{x},\tilde{\sigma}) p(\mathbf{x}|\pi)p(\pi)p(\tilde{\sigma}).
\end{equation}
Then one has
\begin{equation}\label{Equ:VB_PrimaryFunc1}
p(y|\mathbf{x}, \tilde{\sigma})=\prod_{n=1}^{N}\prod_{k=0}^1 \mbox{Bern}( y_n|\tilde{\sigma_{k}})^{x_{nk}},
\end{equation}
where $\mbox{Bern}$ is  Bernoulli distribution.

For the distribution $p(\mathbf{x}|\pi)$, $p(\pi)$ and $p(\sigma)$, they are assumed as
\begin{equation}\label{Equ:VB_PrimaryFunc2}
p(\mathbf{x}|\pi)=\prod_{n=1}^N\prod_{k=0}^1\pi_k^{x_{nk}},
\end{equation}
\begin{equation}\label{Equ:VB_PrimaryFunc3}
p(\pi)=\mbox{Dir}\left(\pi|\alpha_{0}, \alpha_{1}\right),
\end{equation}
\begin{equation}\label{Equ:VB_PrimaryFunc4}
p(\tilde{\sigma})=p(\tilde{\sigma}_0)p(\tilde{\sigma}_1) = \mbox{Beta}(\tilde{\sigma}_0|\beta_{00},\beta_{01}) \mbox{Beta}(\tilde{\sigma}_1|\beta_{10},\beta_{11}).
\end{equation}

The exact form of joint distribution (\ref{Equ:VB_PrimaryFunc}) can be obtained by multiplying (\ref{Equ:VB_PrimaryFunc1}), (\ref{Equ:VB_PrimaryFunc2}), (\ref{Equ:VB_PrimaryFunc3}) and (\ref{Equ:VB_PrimaryFunc4}). However, this multiplication form has much higher complexity, or it is infeasible. 

Under VB framework in Section \ref{subsection:VB_algorithm}, much more easier approximate distributions for $q({\mathbf{x}})$, $q(\pi)$ and $q(\sigma)$ will be given in order to make inference on ${\mathbf{x}}$, $\pi$ and $\sigma$ respectively. The detail of the algorithm is listed in the following

\begin{enumerate}
\item Initialize the parameters $\alpha_0$, $\alpha_1$, $\beta_{00}$, $\beta_{01}$, $\beta_{10}$ and $\beta_{11}$:

Suppose $\{x_1,\dots,x_N\}$ are current decoded words in Region I of Figure \ref{BAC:FactorGraphBAC}, then the parameters can be initialized as
\begin{equation}
\begin{split}
\alpha_0 &= \frac{\#\{x_n==0\}}{N}P_{\alpha} \\
\alpha_1 &= (1 - \alpha_0)P_{\alpha} \\
\beta_{00} &= \frac{\#\{x_n==0\&\&y_n==1\}}{\#\{x_n==0\}} P_{\beta_0}\\
\beta_{01} &= (1-\beta_{00})P_{\beta_0}\\
\beta_{10} &= \frac{\#\{x_n==1\&\&y_n==1\}}{\#\{x_n==1\}} P_{\beta_1}\\
\beta_{11} &= (1-\beta_{10})P_{\beta_1}\\
\end{split},
\end{equation}
where $P_{\alpha}$, $P_{\beta_0}$ and $P_{\beta_1}$ are constants to adjust the values of the parameters of the distributions. The bigger values the constants are, the more influences the corresponding distributions will take effect on the final estimations.
\item Repeat the following steps until all $q(\cdot)$ are converged or reach the maximum number of iterations.
\begin{enumerate}
\item Update $q(x)$:
\begin{equation}
q(x)=\prod_{n=1}^{N}\prod_{k=0}^1 r_{nk}^{x_{nk}},
\end{equation}
where
\begin{equation}
\begin{split}
r_{nk}&=\frac{\rho_{nk}}{\sum_{j=0}^1\rho_{nj}},\\
\ln \rho_{nk} &= y_n \mathbb{E}[\ln \tilde{\sigma_k}] + (1-y_n)\mathbb{E}[\ln(1-\tilde{\sigma_k})] + \mathbb{E}(\ln \pi_k),
\end{split}
\end{equation} 
where
\begin{equation}
\begin{split}
\mathbb{E}[\ln \tilde{\sigma_k}]&=\psi(\beta_{k0})-\psi(\beta_{k0} + \beta_{k1})\\
\mathbb{E}[\ln (1-\tilde{\sigma_k})]&=\psi(\beta_{k1})-\psi(\beta_{k0} + \beta_{k1})\\
\mathbb{E}(\ln \pi_k) &= \psi(\alpha_k)-\psi(\sum_{k=0}^1\alpha_k)
\end{split},
\end{equation}
where $\psi(\cdot)$ is digamma function.
\item Update $q(\pi)$:
\begin{equation}
q(\pi) = \mbox{Dir}(\pi|\alpha_0',\alpha_1'),
\end{equation}
where 
\begin{equation}
\begin{split}
\alpha_{k}{'}=\alpha_k+N_k\\
N_k = \sum_{n=1}^N r_{nk}
\end{split},
\end{equation}
\item Update $\tilde{\sigma}$:
\begin{equation}
q(\tilde{\sigma})=\prod_{k=0}^1 \mbox{Beta}(\tilde{\sigma_k}|\beta_{k0}{'}, \beta_{k1}{'}).
\end{equation}
where
\begin{equation}
\begin{split}
\beta_{k0}{'} &= S_k + \beta_{k0}\\
\beta_{k1}{'} &= T_k + \beta_{k1}
\end{split},
\end{equation}
where
\begin{equation}
\begin{split}
S_k &=\sum_{n=1}^N y_n r_{nk} \\
T_k &=\sum_{n=1}^N (1-y_n) r_{nk}
\end{split}.
\end{equation}
\end{enumerate}
\item Output $\sigma_0$ and $\sigma_1$ as
\begin{equation}
\begin{split}
\sigma_0 &= \tilde{\sigma_0}\\
\sigma_1 &= 1 - \tilde{\sigma_1}
\end{split}.
\end{equation}
\end{enumerate}
 
Both of the EP and VB algorithms can be applied to all of the pair nodes $\sigma_0^j$ and $\sigma_1^j$ so that provides the estimates for  $\sigma_0^j$ and $\sigma_1^j$ to Region II to continue message passing algorithm.

\section{Experimental Results}
The proposed algorithms will be tested in two experiments\delete{, where the correlations of the two sources can be modeled as the BSC channel and Z-channel}, and the details of setups are listed in Table \ref{TAB:BAC_SETUP}. All of the algorithms are implemented by Java and evaluated on an Intel 3.0 GHz CPU machine.

\begin{table}[ht]
\caption{Experimental setups}  \label{TAB:BAC_SETUP}
\begin{center}
\begin{tabular}{ |l|l|l| } 
\hline \hline
 & \textbf{Exp I} & \textbf{Exp II} \\ \hline\hline
$\sigma_0^t$ & $\mathcal{N}(0.15, 0.01^2)$ & $\mathcal{N}(0.2, 0.01^2)$\\
\hline
$\sigma_1^t$ & $\mathcal{N}(0.15, 0.01^2)$ & $\mathcal{N}(0.1, 0.01^2)$\\
\hline
$N$& \multicolumn{2}{|c|}{$10000$} \\ \hline
$C$ & \multicolumn{2}{|c|}{$500$} \\ \hline 
\# of iterations per each trial& \multicolumn{2}{|c|}{$60$} \\ \hline 
learning period & \multicolumn{2}{|c|}{every $40$ iterations}\\ \hline
\# of trials & \multicolumn{2}{|c|}{$100$}\\
\hline\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[!htb]
	\centering
		\includegraphics[width=6in]{./figures/BAC_EXP1}
\caption{Bit error rates versus code rates curves of Exp I} 
	\label{Fig:BAC_BSC}
\end{figure}
First, let's set both $\sigma_0^t$ and $\sigma_1^t$ as Gaussian distribution $\mathcal{N}(0.15, 0.01^2)$. The length of block $N$ is $10000$, and the connection ratio $C$ equals $500$. The crossover probabilities $\sigma_0^t$ and $\sigma_1^t$ for the proposed algorithm are initialized $0.1$ away from the mean of the true probabilities. The algorithm will run $60$ times for each testing, and the VB/EP only run every $40$ iterations. All the results are averaged over $100$ trials. The bit error rate is plotted in Figure \ref{Fig:BAC_BSC}. To make comparison, standard BP algorithm with no parameter estimation is also implemented. In Figure \ref{Fig:BAC_BSC}, It can be seen that the VB-based algorithm outperforms the standard BP algorithm and EP-based algorithm. The BP algorithm and EP-based algorithm achieve the same performance. The largest gap between the VB and the standard BP is more than $10$ dB. 
\begin{figure}[!htb]
	
	\centering
		\includegraphics[width=6in]{./figures/BAC_EXP2}
\caption{Bit error rates versus code rates curves of Exp II}
	\label{Fig:BAC_ZChannel}
\end{figure}

Second, $\sigma_0$ and $\sigma_1$ are set as Gaussian distributions $\mathcal{N}(0.2, 0.01^2)$ and $\mathcal{N}(0.1, 0.01^2)$, respectively. Other parameters are the same as Exp I. In addition, standard BP is also implemented for comparison. The bit error rate versus bit rate is plotted in Figure \ref{Fig:BAC_ZChannel}. In Figure \ref{Fig:BAC_ZChannel}, it  can be seen that the VB-based algorithm still achieves the best performance among all of the three algorithms. The EP algorithms even can not converge. 

Contrary to the results of \cite{cui2012adaptive}, the EP algorithm has the worst performance in above experiments, which can be explained as follow:
\begin{enumerate}
\item The algorithm in \cite{cui2012adaptive} assumes the channel type as BSC; however, none of above experiments is BSC channel. As a matter of fact, the EP algorithm has a better performance in Experiment I than Experiment II. That's because the channel type in Experiment I more approaches to BSC than that of Experimental II.
\item The property of EP algorithm are more suitable to learn the single mode model problem, however, the estimation of BAC problem is a multiple-mode estimation problem. 
\end{enumerate}
\delete{  
\section{Conclusion and Future Work}
In this chapter, we solve a BAC problem, whose parameters are changing with time. We use LDPC code to encode the source data. For the decoding, we establish a three-region factor graph, and they are LDPC region, correlation region and parameter estimation region. For the LDPC and correlation regions, we apply the the standard BP algorithm to implement message passing algorithm, and estimate the channel parameters in Region III by EP and VB algorithms. The proposed algorithm have been tested. The experimental results show that the VB-based algorithm achieve the best performance, and the performance EP-based algorithm are even worse than that of non-estimation method. We point the reason are the essence of EP algorithm, which are not suitable to estimate multi-mode problem.}

\chapter{COMPRESSION OF CORRELATED TEMPERATURE DATA OF SENSOR NETWORK}\label{chap:WZ}
\section{Introduction}\label{WZ:introduction}
In the real world, to monitor the temperatures of different locations in a neighborhood, we usually install sensors at each location in order to obtain their readings periodically. However, due to the limited memory of each sensor, the readings should be sent quickly after they are sampled. As the frequency of the samplings increases, the channel between each sensor and the destination may not transmit the data correctly due to its limited capacity. Since the sensors are installed in the same neighborhood, there may exist spatial correlations among the sensors. In addition, for each sensor, temporal correlations may also exist between the consecutive readings. In this section, both spatial and temporal correlations are utilized to reduce the transmission data rate between each sensor and the destination. To verify the validity of the proposed algorithm, a small wireless network is established in the laboratory for testing.

%Hence, let's first suppose there are two correlated sources $X$ and $Y$ in a network, where $Y$ is used as the side information for the compression of $X$. On the coding side: each of $X$ are firstly quantized into a given bits, and then encode all these bits by LDPC method and send the syndromes and side information to the decoding side. We use belief propagation (BP) algorithm to decode: the side information $Y$ offers the spatial information for the $X$ decoding. In addition, we also take the temporal correlations among $(X_1,\dots ,X_N)$ into account. We applied our method to a set of practical temperature network data, and the results illustrated that the proposed approach outperforms the method that only considers side information.

\section{Problem Formulation} \label{section2}
A Crossbow temperature network~\cite{websiteCrossbow} is installed in the laboratory. The network contains four sensors (senders) and one common destination (receiver). The four sensors are located in different locations and will send their readings to the destination every hour. 

In the network, one of the data source is used as the side information (SI), and the other three sources are compressed based on their correlations with the SI. This type of coding is called "asymmetric coding". Under this framework, this four-source DSC problem can be viewed as the three individual two-source DSC problems.

For simplicity, two sources, $X=[X_1,X_2,\dots,X_L ]$ and $Y=[Y_1,Y_2,\dots,Y_L]$, are assumed in the problem, where $X$ is the source to be compressed, $Y$ is the SI, and each of both sources has $L$-length block data.   

Suppose $X$ and $Y$ are jointly Gaussian, i.e.,
\begin{equation} \label{sec_wz:equ_spatial_correlation}
\displaystyle Y_i=X_i+\mathcal{N}(0,\sigma_s^2), 
\end{equation}
where $\mathcal{N}(0,\sigma_s^2)$ is a Gaussian distribution with $0$ mean and $\sigma_s^2$ variance. 
 
In addition, for each $X_i,1\leq i\leq L-1$, it is assumed that 
\begin{equation} \label{sec_wz:equ_temporal_correlation}
\displaystyle X_{i+1}=X_i+\mathcal{N}(0,\sigma_t^2 ),
\end{equation}
where $\mathcal{N}(0,\sigma_t^2 )$ is a Gaussian distribution with $0$ mean and $\sigma_t^2$ variance.

Data of the source $X$ are first quantized into bits (its binary representation), and then the bits are encoded by an LDPC matrix (see section \ref{section:bac_proposed_method}) in order to obtain the syndromes. These syndromes are then sent to the destination. Data from source $Y$ (SI) are sent to the destination directly.
   	
Once it receives the syndromes and the SI, the decoder will implement the message passing algorithm in the proposed factor graph to recover the data from the source $X$. In the next section, a detailed explanation of the proposed factor graph and how to implement the message algorithm in it will be given. 

\section{Proposed Method} \label{chapter_wz:sec_proposed_method}
\delete{
First, let's classify some notations. $N(a)$ denotes all the neighbors of factor node $a$. In addition, ${\mathbf{x}}_a$ denotes the set containing all variables connecting to $a$. That is,  ${\mathbf{x}}_a=\{x_i |i\in N(a)\}$ .}

\begin{figure}[!htb]
	\centering
		\includegraphics[width=5in]{./figures/DSCFactorGraph}
\caption{Factor Graph of Our Proposed Method} 
	\label{Fig:Kalman}
\end{figure}

Figure \ref{Fig:Kalman} is the factor graph of the proposed algorithm. This factor graph captures all of the aforementioned correlations and constraints. Based on the types of correlations and constraints, the factor graph can be divided into three regions. The rest of this chapter will elaborate on them.

\subsection{Region \Rmnum{1}}
Region \Rmnum{1} captures the constraints of the LDPC algorithm. Suppose each $x_i\ (1\leq i\leq L)$ are quantized into $K$ bits (see Figure \ref{Fig:Kalman}, where $K=3$), and $x_{i,k}\ (1\leq i\leq L,1\leq j\leq K)$ denotes the $k-th$ bit of $x_i$. And all of the $x_{i,k}$ are encoded together by the LDPC algorithm (see Section \ref{BAC_subsection:LDPC}). Suppose the LDPC compression ratio is $M/L$. Then $M$ syndromes $\mathbf{s}=\{s_1,s_2,\dots,s_M\}$ are generated for transmission to the decoder. At the joint decoder, with the constraint imposed by the received syndromes, the factor $f_m^\Rmnum{1}$ can be expressed as
\begin{equation} \label{sec_wz:equ_func1}
f_m^\Rmnum{1}(\textbf{X}_m,s_m) = \left\{\begin{array}{c l}
1 & ,s_m\oplus(\bigoplus \textbf{X}_m)\\
0 & ,otherwise\,
\end{array},\right.
\end{equation}
where $\mathbf{X}_m$ denotes the set of neighbors of the factor node $f_m^I$, $\oplus$ is the bit-xor operation, and $\bigoplus$ is bit-sum of the $\mathbf{X}_m$ set.

All of the messages in Region \Rmnum{1} can be updated according to the algorithm listed in Table \ref{tab:algorithm_LDPC}.

\subsection{Region \Rmnum{2}} 
In Region \Rmnum{2}, we use $f_l^{\Rmnum{2}}$ to denote the relation between the source $x_i$ and its quantized bits $\{x_{i1},x_{i2},\dots,x_{iL}\}$. Noticeably, all the variables on the left of $f_l^{\Rmnum{2}}$ are continuous, and due to the assumption of Gaussian distribution among them, the messages left to $f_l^{\Rmnum{2}}$ are also Gaussian distributions. Since Gaussian distributions can be fully specified by their mean and variance, only two parameters are required to express these messages. Let $Q[X]$ denote one of the candidates of the quantized source $\{x_{i1},x_{i2},\dots,x_{iL}\}$. The factor node $f_l^{\Rmnum{2}}$ can then be given as
\begin{equation} \label{sec_wz:equ_func2}
f_l^{\Rmnum{2}}(Q[X],\mu_l,\sigma_l^2)=
\int_{P(Q[X])-\Delta}^{P(Q[X])+\Delta} 
\frac{1}{\sqrt{2\pi}\,{\sigma_l}} \exp\left(-\frac{(z-\mu_{l})^2}{2\sigma_l^2}\right)dz,
\end{equation}
where $\mu_l$,$\sigma_l^2$ are the parameters of the message $m_{X_l\rightarrow f_l^{\Rmnum{2}}} (\mu_l,\sigma_l^2 )$, $P(Q[X])$ is the de-quantized value of $Q[X]$, and $\Delta$ is half of the step size of the quantization.

For message $m_{f_l^{\Rmnum{2}}\rightarrow X_{lk}}(x_{lk})$, $k=1,2,\dots,K$, since $x_{lk}$ is a binary variable, the BP algorithm in~\eqref{BP:f2v} can be implemented directly to update this type of messages. However, for message $m_{f_l^{\Rmnum{2}}\rightarrow X_{l}}(\mu_l,\sigma_l^2 )$, there are $2^K$ items, which make it difficult or even infeasible for further implementation of message passing in Region \Rmnum{3}. To cope this difficulties, we can use a deterministic approximation method-EP algorithm, and the details will be given in the following section.

\subsection{Region \Rmnum{3}}
According to~\eqref{sec_wz:equ_spatial_correlation} and~\eqref{sec_wz:equ_temporal_correlation}, the functions $f_l^{\Rmnum{3}}$ and $f_l^{\Rmnum{4}}$ in Region \Rmnum{3} can be given as
\begin{equation} \label{sec_wz:equ_func3}
f_l^{\Rmnum{3}}(x_l)=\frac{1}{\sqrt{2\pi}\,{\sigma_s}} \exp\left(\frac{-(x_l-y_l)^2}{2\sigma_s^2}\right)
\end{equation}
and
\begin{equation} \label{sec_wz:equ_func4}
f_l^{\Rmnum{4}}(x_l,x_{l+1})=\frac{1}{\sqrt{2\pi}\,{\sigma_t}} \exp\left(\frac{-(x_{l+1}-x_{l})^2}{2\sigma_t^2}\right)
\end{equation}
respectively.

For each of the variable nodes $X_l$, there are three types of messages that flow into this variable, and they are from factors $f_l^{\Rmnum{2}}$, $f_l^{\Rmnum{3}}$ and $f_l^{\Rmnum{4}}$. The message from the $f_l^{\Rmnum{2}}$ contains $2^K$ items. The message from $f_l^{\Rmnum{3}}$ can be viewed as the prior, which has a Gaussian form. The exact forms of the messages from $f_l^{\Rmnum{4}}$ are also Gaussian distributions. To continue implementing the message passing algorithm in region \Rmnum{3}, the EP algorithm in section~\ref{subsection:EP_algorithm} will be utilized to estimate all of the messages that flow out of the variable node $X_l$. The EP estimations of these messages are expressed as $\tilde{m}_{X_l \rightarrow f_l^{\Rmnum{2}}}(x_l)$, $\tilde{m}_{X_l \rightarrow f_l^{\Rmnum{4}}}(x_l)$, and all of these messages are restricted into Gaussian distribution. Then, the EP updating rules can be derived by using the following equations:
\begin{equation} \label{sec_wz:equ_multiplyGaussian}
\mathcal{N}(x;m_1,\nu_1 )\mathcal{N}(x;m_2,\nu_2 )=\mathcal{N}(m_1;m_2,\nu_1+\nu_1 )\mathcal{N}(x;m,\nu),
\end{equation}
where $\displaystyle v=\frac{\nu_1 \nu_2}{\nu_1+\nu_2}$, $\displaystyle m=v\left(\frac{m_1}{\nu_1}+\frac{m_2}{\nu_2}\right)$. The details of the algorithm can be found in~\cite{Minka:01}. 

The messages $\tilde{m}_{X_l \rightarrow f_l^{\Rmnum{2}}}(x_l)$ and $\tilde{m}_{X_l \rightarrow f_l^{\Rmnum{4}}}(x_l)$ can then be fed into $X_L$'s neighbor factors in order to update the messages that flows out of factors $f^{\Rmnum{2}}$ (see Region \Rmnum{2}) and $f^{\Rmnum{4}}$. The messages that flow out of the factor $ f_l^{\Rmnum{4}}$ can be updated according to the equation

\begin{equation}
\int \mathcal{N}(x_1;m_1,\nu_1 )\mathcal{N}(x_2;x_1,\nu_2)\ dx_1 \propto \mathcal{N}(x_2;m_1,\nu_1+\nu_2).
\end{equation}

\begin{remark}
Actually, region \Rmnum{3} is mainly composed by a Kalman filter. The SI $Y$ can be viewed as an observable value, and $X$ is the hidden variable. 

In addition, when comparing with the prior work~\cite{icip2011a}, there is no need to implement the procedure of de-quantization since the beliefs of the source $X$ in Region \Rmnum{3} are continuous.
\end{remark}

\section{Results} \label{section4}
In this section, the proposed algorithm will be applied to the network. For all of the $4$ sensors, since the data from one sensor will be utilized as the SI, there are $3$ combinations for testing. They can be called Group $1$, Group $2$, and Group $3$. In addition, $50$ readings are obtained from each sensor, i.e., the length of data block $L=50$. The numbers of quantized bits are set as $K = 3, 4, 5, 6, 7$, and the range of quantization is $[0,60]$. Hence, the compression rate $R$ can be calculated as $R=K*M/L$.

\begin{table}[ht]
\caption{The rate/bits of the three data groups given the number of bit planes $(K)$}
\label{sec_wz:tab_results}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
	\hline\hline
     & $K=3$ & $K=4$ & $K=5$ & $K=6$ & $K=7$ \\ \hline
   Group 1($\sigma_s=0.72$,$\sigma_t=1.1$) & $0.4$ & $0.6$ & $1.0$ & $2.3$ & $3.4$ \\
   Group 2($\sigma_s=0.68$,$\sigma_t=0.9$) & $0.2$ & $0.3$ & $1.2$ & $2.0$ & $3.2$ \\
   Group 3($\sigma_s=0.65$,$\sigma_t=1.0$) & $0.2$ & $0.4$ & $1.2$ & $2.3$ & $3.5$ \\ 
   \hline\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[!htb]
  \centering
  \subfigure[ ]{\label{fig:wz_Err1}\includegraphics[width=0.55\textwidth]{./figures/WZ-ErrGroup1.eps}}
  \subfigure[ ]{\label{fig:wz_Err2}\includegraphics[width=0.55\textwidth]{./figures/WZ-ErrGroup2.eps}}
\subfigure[ ]{\label{fig:wz_Err3}\includegraphics[width=0.55\textwidth]{./figures/WZ-ErrGroup3.eps}}
  \caption{Distortion versus the rate listed in Table \ref{sec_wz:tab_results}: (a) Group $1$; (b) Group 2; (c) Group 3}
  \label{fig:wz_Errs}
\end{figure}
 
Since the message passing algorithm is an iterative algorithm, the stop criterion should be given in advance. In all of these experiments, there are two stop criteria: The first one is that all of the beliefs of the decoding bits $x_{i,k}$ successfully satisfy the LDPC syndrome constraints in~\eqref{sec_wz:equ_func1}. The second one is when the implementation reaches the maximum number of iterations. The maximum number of iterations is set to be $60$. In addition, the EP algorithm will be implemented every $20$ iterations. 

For a given number of bit plane, the minimum rates $R_{min}$ of compressed syndromes that can be successfully decoded for the proposed algorithm are calculated in Table \ref{sec_wz:tab_results}. In addition, the spatial and temporal standard deviations of each data group are also listed.

Figure~\ref{fig:wz_Errs} gives the distortions versus the rates listed in Table~\ref{sec_wz:tab_results}. Here, the distortion is defined as the mean square error, and it is given in decibels (dB) unit. To make a comparison, the algorithm in~\cite{icip2011a}(by Wang et al.) is also implemented, and~\cite{icip2011a} only considers spatial correlations~\cite{icip2011a}. Unlike~\cite{icip2011a}, the proposed algorithm assumes the correlations are priors and estimates them offline. For a fair comparison, the same assumption is also applied to~\cite{icip2011a}. From the results in Figure~\ref{fig:wz_Errs}, the proposed algorithm outperforms the algorithm in~\cite{icip2011a}. In the Figure~\ref{fig:wz_Err1} and Figure~\ref{fig:wz_Err3},it is seen that the distortions of the algorithm in~\cite{icip2011a} increase as the rates of data increase. That's because, for a given bit plane and  LDPC setting, the algorithm in~\cite{icip2011a} cannot decode successfully. However, since the correlations among the data have been fully utilized in the proposed algorithm, the probability of decoding successfully will increase.

\chapter{STREAMLINED GENOME SEQUENCE COMPRESSION USING DISTRIBUTED SOURCE CODING}\label{chapter:genome}
\section{Introduction}

\delete{
\begin{wrapfigure}{r}{0.45\textwidth}
\includegraphics[width=3in]{./figures/MinION}
\vspace{-0.1in}
\caption{An example of a miniaturised sequencing instrument.}
\vspace{-0.2in}
\label{fig:MinION}
\end{wrapfigure}}

%The miniaturized sequencing device decked with tiny batteries, sensor chip and USB port is no longer just a vision. 
The vision of miniaturized sequencing devices is turning into reality with the emergence of MinION by Oxford Nanpore \cite{MinION2012_ONLINE}.
Such devices are promising in a variety of potential applications, ranging from studying of wildlife \com{and} clinical capture of sequenced genes, to food inspection for identifying pathogens.
However, such portable devices are commonly \com{subject} to the constraints in processing capabilities, power budget, and storage and communication \com{limitations}.
With these constraints, the traditional view of genome compression architecture as simple decoder and complex encoder needs to be changed.
It is urgent to develop novel techniques to satisfy the emerging reality challenges.
Data compression methods (for reducing the storage space with significantly lower computational complexity and memory requirements) become crucial for the efficient management of genomic data in portable devices.
%The biomedical community urgently calls for novel techniques to tackle the reality challenges in dealing with the ever growing genomic data. Data compression methods for reducing the storage space and saving the data transfer bandwidth become crucial for the efficient management of large genomic data.

%With the emergence of MinION by Oxford Nanpore, the miniaturized sequencing device decked with tiny batteries, sensor chip and USB port is no longer just a vision.
%Recent progress in DNA sequencing has led to the emergence of the portable sequencing devices, e.g. MinION by Oxford Nanopore, which enables the simple and free gene sequencing at a low cost, e.g. clinical, wildlife food inspectation.
%Such miniaturized sequencers with the size of a USB memory stick are adaptable for techniques such as DNA sequencing and protein sensing and compatible with the sequencing machine with the delivery of same workflows.
In both situations (with or without reference sequences), \com{traditional genome compression is computationally expensive at the encoder}. The complexity is dominated by matching (approximately) repeated patterns of nucleotides\ths{---namely Adenine (A), Cytosine (C), Guanine (G) and Thymine (T)---}between or within the DNA sequences.
These patterns are also accompanied by insertions, deletions and substitutions of single nucleotides.

%Due to the nature of DNA data with long, sparse matches and low contextual correlation, almost all proposed substitutional methods are superior to the solutions based only on Markov models.
To date, a number of specialized DNA sequence compression algorithms have been proposed.
In the spirit of Ziv and Lempel \cite{Ziv1978_TIT}, Grumbach and Tachi \cite{Grumbach1993_DCC} proposed the first DNA sequence compressor, \delete{known as} Biocompress, to compress the exact repeating patterns with a specially designed Fibonacci coder. \com{The algorithm} was then improved in Biocompress-2 \cite{Grumbach1994_JIPM} by introducing \com{a} Markov model for encoding the non-repeated regions. \com{Chen et al.} \cite{Chen2001_MEMB,Chen2002_BIO} extended the earlier approach to cover approximated repeats by further exploiting the nature of DNA sequences. Meanwhile, \com{the work in} \cite{Matsumoto2000_GI} introduced a combined CTW+LZ algorithm for searching approximate repeats and palindrome using hash and dynamic programming. \com{Behzadi and Fessant} \cite{Behzadi2005_CPM} proposed a dynamic programming approach for the optimal selection of approximate repeats with promising compression efficiency being witnessed. However, such methods are heuristic
\replace{as the underlying statistics of the sequence patterns are generally ignored}{since no information on theoretic principles can be obtained in the approximate matching}. \com{The authors in} \cite{Tabus2003_DCC, Korodi2007_MSP, Korodi2007_DCC} proposed to combine the matching and substitution of approximate repeats and \replace{a}{the} specific normalized maximum likelihood model\com{, obtaining} a much higher compression ratio. \com{Subsequently}, statistical modeling for predicting the generation of symbols and arithmetic coding for such symbols in DNA sequences were proposed for more \replace{efficient}{intensive} compression. \replace{Cao {\it et al.}}{XM} \cite{Cao2007_DCC} proposed to estimate the probability distribution of symbols with a panel of ``expert'' \replace{to tackle the approximate repeat problem.}{ with the repeat expert concerning the approximate repeats.} \com{Alternatively,} finite context models are proposed \delete{and compared} to \delete{rapidly} capture different aspects of statistical information along the sequence \cite{Pinho2009_ICASSP, Pratas2011_PACBB}\delete{. In spite of the evolutionary development of compression technique}, such reference free methods are \replace{plagued by}{subjected to} their low compression \com{rates} (not greater than 6:1) and prohibitive computational \replace{consumption}{consuming} for large DNA sets.

Recognizing reference-free architectures do not fully utilizing information, a series of algorithms are proposed to compress \replace{sequences}{a sequencing} by matching \delete{the} approximate repeats with \replace{a}{the} reference sequence. The RLZ algorithm proposed by Kuruppu \textit{et al.} \cite{Kuruppu2010_SPIR} performed relative Lempel-Ziv compression of DNA sequences with the collection of related sequences. Wang \textit{et al.} \cite{Wang2011_NAR} proposed the GRS compressor, that is able to compress a sequence with \com{using a} reference without any additional information. Applying the copy model into the matching of exact repeats in reference sequences, GReEn \cite{Pinho2012_NAR} achieved even larger gains when compared to \cite{Kuruppu2010_SPIR} and \cite{Wang2011_NAR}. \com{Recently, reference-based algorithms \cite{Kozanitis2011_JCB, Fritz2011_GR} achieved highly efficient compression performance for the fastq data format, by matching and comparing repeated subsequences in the reference sequences.} Although the reference-based architectures can achieve hundreds of folds compression, the requirement of reference sequences makes it impractical for miniaturized devices, which have very limited storage space and communication bandwidth.

In this section, a novel and pioneering architecture for the genome compression application in miniaturized devices with limited processing capabilities, power budget, storage space and communication bandwidth is proposed. The contribution of the proposed method is three-fold.

\begin{enumerate}
\item To the best of our knowledge, the proposed architecture is the first practical one to meet the demands of miniaturized devices. Motivated by the distributed source coding (DSC) for sensor networks \cite{Xiong2004_MSP}, \com{the proposed scheme} includes a simplified encoder \com{without having access to reference sequences or communicating with other encoders}, and a complex decoder that detects repeated subsequences in the stored reference sequences and decompress the received encoded bits with the specifically designed graphical model. \com{Hence, the proposed compression system can successfully meet the constraints} and requirements of the miniaturized sequencing devices.

\item A flexible encoding and decoding mechanism is proposed. \com{Using feedback from the decoder, the encoder} transmits either hashes conducting the detection of variable-size exact repeats in decoder or syndromes \com{obtained with low-complexity Slepian-Wolf coding \cite{Pradhan2003_TIT} of the non-repeated subsequences.} The proposed encoder and decoder perform efficiently by taking consideration of both exact repeats and approximately repeated subsequences (e.g. insertion, deletion and substitution).

\item With  the syndrome and reference sequences at the decoder~\cite{Pradhan2003_TIT}, a novel factor graph model is constructed to tackle the challenge in detecting insertion, deletion and substitution between the reference and original source. Experimental results show that the proposed architecture \com{can} achieve an efficient compression performance with significantly low encoding complexity when compared to the benchmark compressor GRS.

\end{enumerate}

The rest of this section is organized as follows. In Section~\ref{section:system}, the proposed DSC based genome compression system is introduced, and the system includes the implementation details of the proposed hash based (exactly) repeated sequence coding with adaptive length and an overview of syndrome based non-repeated sequence coding. Then, in Section~\ref{section:DSC}, a design of the syndromes based non-repeated sequence coding is proposed, and this design can handle the insertion, deletion, and substitution between sources and reference. The experimental results can be found in Sections~\ref{section:results}.



\section{System architecture}\label{section:system}

\begin{figure}[!t]
\centering
\includegraphics[width=5.0in]{./figures/workflow}
\vspace{-0.1in}
\caption{Workflow of genome compression based on DSC.}
\vspace{-0.2in}
\label{fig:workflow}
\end{figure}

The block diagram of the proposed Genome compression framework is depicted in Figure~\ref{fig:workflow}. Suppose that there are two correlated DNA sequences (i.e., source and reference sequences) available at the encoder and decoder, respectively, where the variations between two sequences are modeled by insertion, deletion and substitution. The alphabet of our studied DNA sequence is confined within the set $\{\mbox{\lq}A\mbox{\rq}, \mbox{\lq}C\mbox{\rq}, \mbox{\lq}G\mbox{\rq}, \mbox{\lq}T\mbox{\rq}, \mbox{\lq}N\mbox{\rq}\}$, where \com{$\mbox{\lq}N\mbox{\rq}$ denotes an unknown base} due to a low sequencing quality. Figure~\ref{fig:logicalFlow} shows the logical flow of the proposed framework.

\begin{figure}[!t]
\centering
\includegraphics[width=6.0in]{./figures/logicalFlow}
%\vspace{-0.1in}
\caption{Workflow of genome compression based on DSC.}
%\vspace{-0.1in}
\label{fig:logicalFlow}
\end{figure}

At the encoder (see the left hand side of Figure \ref{fig:logicalFlow}), \com{a} streaming DNA sequence obtained from the portable sequencer will be first stored in the incoming data buffer for further processing. Second, a sub-sequence $\mathbf{x}_i^L$, which starts with the $i$-th to be compressed base in the source sequence, is extracted from the incoming data buffer, where its length $L$ and the corresponding coding method are decided by the adaptive code length and types selection module. The compressed sequence can be either LDPC Accumulate (LDPCA\footnote{LDPCA is an extension of LDPC in section \ref{BAC_subsection:LDPC}, which has a feedback with the decoder so that the encoder can change the LDPC length of syndrome adaptively.}) syndromes $\mathbf{s}_{\mathbf{x}_i^L} = \mathbf{H}\mathbf{x}_i^L$ or hash bits $\mathbf{h}_{\mathbf{x}_i^L}$ depending on whether variations are presented between \com{the} source and \com{the} reference \com{sequence,} based on the decoder feedback, where $\mathbf{H}$ is the parity check matrix in LDPC codes (see section~\ref{BAC_subsection:LDPC}). Third, \com{the} encoded sequence will be temporally stored in the forward data buffer and send to the decoder.

\begin{figure}[!t]
\centering
\includegraphics[width=4.5in]{./figures/hash_match}
\vspace{-0.1in}
\caption{The diagram of hash-based coding}
\vspace{-0.1in}
\label{fig:slidingWindow}
\end{figure}

At the decoder (see the right hand side of Fig. \ref{fig:logicalFlow}), the received streaming data in the incoming data buffer will be processed by one of the following modules based on the corresponding data \com{compression mode} (i.e., either hash bits or syndromes).
\begin{enumerate}
\item For the received hash data $\mathbf{h}_{\mathbf{x}_i^L}$, it will be compared with the hashes generated from a bunch of sub-sequence candidates $\mathbf{y}_{j+U}^L, \cdots, \mathbf{y}_{j+V}^L$ within the reference sequence for \com{$V-U+1$ total} candidates, where $j$ is the current offset compensated start location, and $U$ and $V$ are predefined lower and upper bounds of the search region for start locations. Then, the comparison result can be further processed as follows.
\begin{enumerate}
\item \label{stepA} If a matched hash $\mathbf{h}_{\mathbf{y}_{k}^L}$ for $k = j+U, \cdots, j+V$ is detected (i.e., $\mathbf{h}_{\mathbf{y}_{k}^L} = \mathbf{h}_{\mathbf{x}_{i}^L}$), the next offset compensated start location of the sliding window can be updated as $j = k+L$ (see Fig. \ref{fig:slidingWindow}). Moreover, we claim that $\mathbf{y}_{k}^L$ will be identical to $\mathbf{x}_{i}^L$, if $\mathbf{h}_{\mathbf{y}_{k}^L}$ and $\mathbf{h}_{\mathbf{x}_{i}^L}$ are matched with each other, which is the fundamental assumption of our proposed system. Intuitively, the aforementioned assumption can be enforced by choosing a strong hash code with a small search region. The experimental results based on sequences~\cite{huala2001arabidopsis, rhee2003arabidopsis} with total more than 238 million bases demonstrate that a 16-bit cyclic redundancy check (CRC) hash code with a search region $U=-2$ and $V= 10$ provides a strong assertion of such assumption. In addition, the decoder will inform the success to the encoder and request a longer code length based on a predefined protocol as updating $L^{\mbox{\small current}} = bL_0$, where $L_0$ is a predefined initial length and the scaling factor $b$ is updated as $b = b  + d_b$, $d_b$ is an incremental constant, and $b$ is initialized as 0.  For example, at the beginning, $L^{\mbox{\small current}} = L_0$, if a matched hash is detected, the adaptive length  $L^{\mbox{\small current}}$ will be updated as  $L^{\mbox{\small current}} = d_bL_0$, as the scaling factor $b = 0+ d_b$. Similarly, if $n_h$ number of successively matched hashes are detected, the adaptive length and its corresponding scale factor will be $L^{\mbox{\small current}} = n_hd_bL_0$ and $b = n_hd_b$, respectively.
\item If no matched hash can be detected, the following two conditions will be checked.
\begin{enumerate}
\item if $L^{\mbox{\small current}} = L_0$, the decoder will inform the hash matching failure to the encoder and request syndromes from \com{the} encoder for further action.
\item Otherwise, the decoder also informs the hash matching failure to the encoder, but requests a shorter code length by setting $L^{\mbox{\small current}} = L_0$.
\end{enumerate} 
\end{enumerate}
\item For the received syndromes, the decoder will pass \com{the} syndrome to the proposed factor graph based LDPCA decoder with the capability of handling deletion, insertion and substitution between \com{the} source and \com{the} reference (see the next section for more implementation details). The following two conditions will be checked.
\begin{enumerate}
\item If the decoded source $\mathbf{\hat{x}}_i^L$ satisfies both the parity check constraint (i.e., $\mathbf{s}_{\mathbf{x}_i} = H\mathbf{\hat{x}}_i^L$ ) and the hash constraint (i.e., $\mathbf{h}_{\mathbf{x}_i^L} = \mathbf{h}_{\mathbf{\hat{x}}_i^L}$), \com{the} decoder will send \com{an} LDPCA success message back to the encoder and update \com{the} offset compensated start  location $j$ through the Smith-Waterman local alignment between the reference and the decoded source. Moreover, the encoder will send hash codes to the decoder for the next sub-sequence.
\item Otherwise, the decoder will request additional LDPCA syndromes from the encoder.
\end{enumerate}

\end{enumerate}


\section{Syndrome Based non-repeated sequence coding}\label{section:DSC}

\begin{figure}[!t]
\centering
\includegraphics[width=5in]{./figures/DCC_factorgraph}
\vspace{-0.1in}
\caption{Factor graph of genome compression based on DSC.}
\vspace{-0.2in}
\label{fig:factorGraph}
\end{figure}

As previously mentioned in our system architecture, if an exact repeat cannot be identified by hash coding, the decoder will request syndromes from \com{the encoder through a feedback channel}. In this section, \com{the} codec design of the proposed syndrome based non-repeated sequence coding is proposed.

\subsection{Syndrome based non-repeated sequence encoding}\label{sec:syndrom_encoding}

The first step of the proposed syndrome based non-repeat encoder is to convert DNA \com{data} into a binary source, such that they can be compressed \com{under a} binary LDPCA encoder. Suppose the following mapping rule for the letters within the alphabet, i.e., $\mbox{\lq}A\mbox{\rq} \rightarrow 000$, $\mbox{\lq}C\mbox{\rq} \rightarrow 001$, $\mbox{\lq}G\mbox{\rq} \rightarrow 010$, $\mbox{\lq}T\mbox{\rq} \rightarrow 011$, $\mbox{\lq}N\mbox{\rq} \rightarrow 100$, a DNA subsequence $\bf{x}$ can be represented by the corresponding binary vector $\bf{x}_b$. For instance, given a DNA subsequence $\bf{x} = $ [\lq A\rq \lq T\rq \lq G\rq \lq C\rq \lq T\rq \lq N\rq$]^T$ with length $N = 6$, its corresponding binary vector will be $\bf{x}_b = [000\ 011\ 010\ 001\ 011\ 100]^T$ with length $3N$. Thus, for LDPC based Slepian-Wolf (SW) coding (i.e., lossless DSC), the compressed syndromes will be generated through $\bf{s_x} = \bf H \bf{x}_b$, where $\bf {H}$ is a sparse parity check matrix with size $M \times 3N$ and $M < 3N$. Thus, the resulting code rate can be expressed as \replace{$R=M/N$ bits per base.}{$R = M/3N $ code rate.} It is worth mentioning that the computational complexity of the aforementioned encoder is ultra-low, since the only operation is the bit-wise multiplication between \com{the} sparse matrix $\bf H$ and the original source. Moreover, LDPCA codes \com{to implement} rate adaptive decoding is employed, where \com{the} decoder can incrementally request additional LDPCA syndromes from \com{the encoder through a feedback channel}, when facing decoding errors.

\subsection{Syndrome based non-repeated sequence decoding}
To perform syndrome based decoding for non-repeat DNA subsequence $\bf{x}$ with the reference sequence as side information $\bf{y}$, the key factor is to be able to explore the variations between \com{the} source subsequence $\bf{x}$ and \com{the} reference sequence $\bf{y}$, where the variations are modeled by the insertion, deletion, and substitution between the source and reference. Moreover, a substitution can be expressed as an insertion in the source sequence followed by a deletion in the corresponding location in the reference sequence. In this section, it is \com{demonstrated} that such variations can be effectively estimated through Bayesian inference on graphical models. The graphical model of the proposed syndrome based decoding with variation is depicted in Figure \ref{fig:factorGraph}. In Figure \ref{fig:factorGraph}, the variable nodes (usually depicted by a circle) denote variables such as source symbol, binary source bits, local offset introduced by variation, and syndromes. The detail of the proposed factor graph will be explained in the followings.

The parity check constraint imposed by the received syndromes is studied, where  $s_1, \cdots, s_M$, the realization of variable node $S_l$, $l=1, \cdots, M$, denotes the received syndromes in Figure~\ref{fig:factorGraph}. Similar as the LDPC codes in Section \ref{section:bac_proposed_method}, \com{the} factor nodes $c_l$, $l=1, \cdots, M$, take into account the parity check constraints, where the corresponding factor function can be expressed as

\begin{eqnarray}
\label{eqn:factor_compression_constraints}
 c_l(\mathbf{x}_{c_l}, {s}_l) = \left\{
\begin{matrix}
 1,  & \mbox{if } s_l \oplus \bigoplus \mathbf{x}_{c_l} = 0 \\
0, & \mbox{otherwise}
\end{matrix},
\right.
\end{eqnarray}
where $\mathbf{x}_{c_l}$ denotes the set of neighbors of \com{the} factor node $c_l$, and $\bigoplus \mathbf{x}_{c_l}$ denotes the binary sum of all elements of the set $\mathbf{x}_{c_l}$.

Moreover, $x_i^1, x_i^2, x_i^3$, the realization of variable node $X_i^r$ with $i =1, \cdots, N$, $r=1, 2,3$, are the binary representation for the $i$-th base $x_i$ in the DNA sequence according to the mapping rule introduced in Section \ref{sec:syndrom_encoding}, where the mapping rule is captured by the factor node $f_i$, $i=1, \cdots, N$ with corresponding factor function as follows


\begin{eqnarray}
f_i(x_i,x_i^1, x_i^2, x_i^3) = \left\{
\begin{matrix}
 1,  & \mbox{if } \mathbf{map} (x_i^1, x_i^2, x_i^3) = x_i, \\
0, & \mbox{otherwise}.
\end{matrix}
\right.
\end{eqnarray}
where $\mathbf{map}(\cdot)$ denotes the mapping from the binary bits \lq\lq $\cdot$ \rq\rq to a letter in the alphabet, e.g. the output of $\mathbf{map}(0,1,1)$ corresponds to the letter \lq\lq T\rq\rq.


Moreover, since the alphabet is not uniformly distributed in an arbitrary DNA sequence, the prior distribution for the alphabet is captured by the factor node $h_{x_i}$, where learning prior through training DNA sequences will be discussed shortly in the results section.

Now, an additional erasure variable node $M_i$ to capture the variation between reference $y_i$ and source $x_i$ is introduced, where the variable $m_i = 1$ indicates the presence of variations, $m_i = 0$ means the existence of matches $y_{i+t_i} = x_i$, where $t_i = -T, \cdots, T$ are all possible local offsets within the search region $[-T, T]$.  Moreover, the corresponding prior distribution of variable $m_i$ is captured by the factor node $h_{m_i}$ with factor function defined as
\begin{eqnarray}
h_{m_i}(m_i) = \left\{
\begin{matrix}
 1 - p_e & \mbox{if exist matches $y_{i+t_i} = x_i$}\\
p_e, & \mbox{otherwise}
\end{matrix},
\right.
\end{eqnarray}
\noindent where $p_e$ can be learnt through training DNA sequence.

For the existence of matches $y_{i+t_i} = x_i$, the local offset $t_i$ is captured by the variable node $T_i$ and its corresponding prior is represented by the factor node $h_{t_i}$ with $h_{t_i}(t_i) = p_{t_i}$, where $p_{t_i}$ can be learnt through training DNA sequences. Furthermore, as the local offsets between adjacent DNA bases \com{do} not vary significantly in this work, it is expected that adjacent \com{variables} $t_i$ will not differ much in value. Such characteristic is captured by the additional factor node $q_{t_i, t_{i+1}}$, where the corresponding factor function is defined as
 \begin{equation}
\label{eq:Laplacian}
q_{t_i, t_{i+1}}({t_i}, t_{i+1}, \alpha) = \frac{\alpha}{2}e^{-\alpha | t_i - t_{i+1} |},
\end{equation}
where $\alpha$ is the scale parameter of the Laplace distribution.

The factor node $g_i$ and its corresponding factor function $g_i(m_i, y_i, x_i, t_i)$ are introduced to combine the impact imposed by the side information $y_i$, erasure $m_i$ and local offset $t_i$. For $m_i = 0$, \com{the }factor function can be expressed as,

\begin{eqnarray}
\label{eq:correlation_Constraint}
g_i(m_i =0, y_i, x_i, t_i) = \left\{
\begin{matrix}
1  & \mbox{if}\ y_{i + t_i} = x_i, \\
0, & \mbox{otherwise}.
\end{matrix}
\right.
\end{eqnarray}

\noindent For $m_i = 1$, \com{the} variable nodes $T_i$ and $Y_i$ will be disconnected from the factor node $g_i$ . Therefore, the simplified factor function $g_i(m_i = 1, x_i) = 1$ can be used to \com{take} the impact of erasure into account.

\delete{Finally, in the context of Bayesian inference, estimating unknown variables corresponds to the evaluations of their marginal distribution, which can be efficiently achieved by performing message passing on a factor graph.} With the factor graph defined above in Figure~\ref{fig:factorGraph}, BP algorithm in Section~\ref{subsection:BP_algorithm} can be implemented so that obtain the posterior distribution (belief) of each variable. The original DNA sequence can be recovered \com{by} the posterior distribution of each source $x_i$.

\section{Results}\label{section:results}

\fc{Two genome sequence data - The Arabidopsis Information Resrouce (TAIR) \cite{websiteTAIR} and The Institute for Genomic Research (TIGR) \cite{websiteTIGR} are adopted for testing in this experiment. These two database are collected by professional groups or institutes, and have been widely used by research communities. }

\fc{The TAIR maintains a database of genetic and molecular biology data for the model higher plant Arabidopsis thaliana \cite{websiteTAIR}. In this experiment,} TAIR8 \cite{huala2001arabidopsis} dataset and TAIR9 \cite{rhee2003arabidopsis} dataset are tested, where each dataset contains five chromosomes \com{with over $238$ millions bases in total}. Moreover, the genome of TAIR9 is used for testing the compression performance with TAIR8 as reference only available at the decoder. \fc{For this experiment}, all the hyper-parameters are initialized as follows, the initial code length $L_0 = 528$, the incremental constant $b_d = 3$, the scale parameter of Laplace distribution $\alpha = 1$, the maximum local offset search region $T = 4$ and the erasure probability $p_e = 0.01$. The proposed codec is implemented in MATLAB and evaluated on an Intel 3.0GHz CPU machine.

\begin{figure}[!htb]
  \centering
  \subfigure[ ]{\label{fig:stat_DNABase}\includegraphics[width=0.4\textwidth]{./figures/stat_DNABase}}
  \subfigure[ ]{\label{fig:stat_offset}\includegraphics[width=0.4\textwidth]{./figures/stat_offset}}
  \caption{The empirical statistics of (a) \com{the} DNA bases \{\lq A\rq, \lq T\rq, \lq G\rq, \lq C\rq, \lq N\rq\} and these of (b) \com{the} local offsets with the range from $-4$ to $4$.}
  \label{fig:DNABase_offset}
\end{figure}

First, the empirical \com{marginal} statistics of \com{the} DNA bases \{\lq A\rq, \lq T\rq, \lq G\rq, \lq C\rq, \lq N\rq\} and these of \com{the} local offsets $t_i$ within the range from $-4$ to $4$ \com{are} shown in Figures~\ref{fig:stat_DNABase} and~\ref{fig:stat_offset}, respectively, which will be used as the priors in the syndrome based non-repeated sequence decoding. In Figure~\ref{fig:stat_DNABase}, it is \com{verified} the assumption that the alphabets of DNA sequences are usually non-uniformly distributed. Moreover, Figure~\ref{fig:stat_offset} depicts that the maximum local offset with $T= 4$ is sufficiently large for capturing shifts between \com{the reference and the source}.

\begin{figure}[!htb]
  \centering
  \subfigure[ ]{\label{fig:codeRateVSoffset}\includegraphics[width=0.4\textwidth]{./figures/DSC_codeRateVSoffset}}
  \subfigure[ ]{\label{fig:compressionSizeVSchrNum}\includegraphics[width=0.4\textwidth]{./figures/compressionSizeVSchrNum}}
  \caption{Compression performance of the proposed  codec \fc{on TAIR dataset}, (a) the average code rates vs. the different maximum local offsets in syndrome coding; (b) the overall compression performance (i.e., hash bits + syndromes) for all 5 chromosomes. }
  \label{fig:codeRateVSoffset_SizeVSchrNum}
\end{figure}

Figure \ref{fig:codeRateVSoffset} illustrates the relationship between the average code rates and the different maximum local offsets in syndrome coding based on all $5$ chromosomes.  In Figure~\ref{fig:codeRateVSoffset}, it shows that the code rates decrease as the maximum local offsets increase, \com{due to the fact that} a larger maximum local offset offers a wider search region for exploring the reference. However, a larger maximum local offset may also \com{result} in a higher decoding complexity.  Figure~\ref{fig:compressionSizeVSchrNum} shows the overall compression performance (i.e., hash bits + syndromes) for all 5 chromosomes in terms of compressed file size. Moreover, 
Figure~\ref{fig_compare} shows a side-by-side comparison of the compression rate and compression time. It shows \com{that} both the proposed method and GRS algorithm achieve significant file size reductions (i.e., up to $8252 \mbox{x}$ file size reduction). 


\begin{figure}[ht]
  \begin{center}
    \includegraphics[width=0.7\textwidth]{./figures/GRScompare}
  \end{center}
  \caption{Performance comparison between GRS and our proposed codec on \fc{TAIR dataset}. }
  \label{fig_compare}
\end{figure}

\fc{For the TIGR data, the chromosome $4$ ($35.8$MB) of the TIGR5 dataset is tested using the chromosome 4 of the TIGR6 as the reference by varying the LDPC code length (i.e., 528, 1056, 1584, 2112 and 2640) as shown in Table \ref{TIGRTable}. Moreover, GRS method is also implemented on this dataset, and the result was also listed in Table~\ref{TIGRGRS} as reference. It can be seen that the compression performance decreases as the LDPC code length increases. The proposed algorithm achieved better compression performance comparing with GRS in this case. It is because that the reference chromosome $4$ in TIGR6 includes a significant amount of insertions when comparing with the same chromosome in TIGR $5$, where the insertion information in the reference chromosome has no contribution to the size of DSC compressed data.}

The proposed encoder shows a significantly lower encoding complexity. It is worth mentioning that the proposed codec is implemented by MATLAB, where a potential performance boost is highly expected by using more efficient programming languages e.g., C/C++. To the best of my knowledge, this is the first study of DSC based genome compression. There is no doubt that it opens many possibilities for the portable miniaturized applications in which energy consumption and bandwidth usage are of paramount importance.

\begin{table}[t]
\begin{center}
\caption{Performances of our proposed method on Chromosome $4$ of TIGR$5$ ($35.8$ MB)} \label{TIGRTable}
\begin{tabular}{|c|c|c|c|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  LDPC Length & Compression Size (KB) & Enc Time (s) & Dec Time (s)
  \\
  \hline\hline
  528 & 3.68 & 0.04 & 298 \\
  1056 & 4.58 & 0.009 & 596 \\
  1584 & 5.67 & 0.01 & 787 \\
  2112 & 6.97 & 0.01 & 1102 \\
  2640 & 6.72 & 0.08 &1374 \\
  \hline\hline
\end{tabular}

\end{center}
\end{table}

\begin{table}[t]
\begin{center}
\caption{Performance of GRS on Chromosome $4$ of TIGR$5$ ($35.8$ MB)} \label{TIGRGRS}

\begin{tabular}{|c|c|c|}
  \hline\hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Compression Size (KB) & Encoding Time (Seconds) & Decoding Time (Seconds)
  \\
  \hline
  26.34 & 12 & 6 \\
  \hline\hline
\end{tabular}
\end{center}
\end{table}


% !TEX root = ../thesis_main.tex
 \chapter{CONCLUSIONS}\label{chap:conclusion}
DSC has attracted the interests of many researchers in recent years. The S-W and W-Z problems are two of the most important topics with regard to DSC. Although the theoretical compression limits of these two problems have been known for decades, a systematic method for solving both problems is still unavailable. Therefore, more and more researchers have focus on and tried to solve these problems. This dissertation has attempted to solve three main problems: adaptive Slepian-Wolf decoding for two binary sources (ASWDTBS) problem, the compression of correlated temperature data of sensor network (CCTDSN) problem and the streamlined genome sequence compression using distributed source coding (SGSCUDSC) problem. 

ASWDTBS problem is an S-W problem. LDPC code was used to encode the source data. For the decoding, a three-region factor graph was established. The three regions are the LDPC region, the source data and SI correlation region,  and the parameters estimation region. For the LDPC and correlation regions, the BP algorithm was applied to implement the message passing algorithm. The channel parameters in Region III were estimated by EP and VB algorithms. The proposed algorithm was tested with the synthetic data and compare with the standard BP algorithm with no parameter estimation. The experimental results show that the VB-based algorithm performed the best, but that the EP-based algorithm performed the worse. 

CCTDSN can be viewed as an extension of the W-Z problem.  A Distributed Source Coding algorithm, taking both spatial and temporal correlations into account, is proposed. To test this method, a wireless temperature networks was established in our laboratory, and each sensor sent its temperature readings every hour. The results showed that the algorithm achieved much better results than the algorithm of prior work on the temperature network. In the future, more research will be done on the quantization step of this W-Z problem, and an optimization between the number of bit planes and each bit plane's compression ratio for a given rate will be sought.

Finally, a DSC based genome compression architecture was presented. This area of research can be classified as a S-W problem. To the best of our knowledge, the proposed framework is the first study of its kind. It is especially targeted at the low complexity genome encoding for miniaturized devices, which have limited processing capabilities, power budgets, storage space and communication bandwidth. Compared to the traditional reference based DNA compression algorithm (e.g., GRS), the proposed framework offers ultra-low encoding complexity (non-repeated subsequences are encoded \replace{using}{through a} low complexity DSC encoding), while (exactly) repeated subsequences are compressed through adaptive length hash coding based on the decoder feedback. The customized factor graph based decoder tackles the challenges of detecting insertion, deletion and substitution between the reference and \com{the} original source, and it recovers the non-repeated subsequences based on received syndromes. Last but not least, our proposed genome compression framework incorporates LDPCA codes for rate adaptive decoding. Experimental results show that the proposed architecture could achieve an efficient compression performance with significantly lower encoding complexity when compared to the benchmark compressor GRS.


%
%\nocite{*} % Use to exclude specific citations from *.bib file

\bibliographystyle{ieeetr}
\singlespacing
\bibliography{./ref/ref2,./ref/referencesGenome}
%\bibliography{./ref/references0,./ref/references1,./ref/ref0,./ref/ref1}
%\appendix
%\input{appendix}
%\input{vita}
%       \pagestyle{empty}
% \begin{abstract}{Your Name}{Master of Science}{My Major} % Creates abstract
%    \input{abstract}
% \end{abstract}

\end{document}

%\floatstyle{ruled}
%\newfloat{algorithm}{htbp}{loa}
%\floatname{algorithm}{Algorithm}
%
%\begin{algorithm}
%\caption{Expectation Propagation}
%\label{alg:GS}
%
%\begin{tabbing}
%1. \=Initialize the compatibility approximation terms $\tilde{\psi}_i(x)$, e.g. $\tilde{\psi}_i(x) = 1$\\
%2. \=Compute the approximated distribution \\
%\begin{equation}
%q(x;\theta^i) = \frac{\prod_i\tilde{\psi}_i(x)}{\int_x \prod_i \tilde{\psi}_i(x)}
%\end{equation}\\
%
%\>2. For \=$i=1$ to $n$\\
%\>\>Set 
%\begin{math}
%x_i^{(k)} = 
%\frac{b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}
%         -\sum_{j=i+1}^{n}a_{ij}x_j^{(k-1)}}{a_{ii}}
%\end{math}
%\\
%\>3. If $\|\vec{x}^{(k)}-\vec{x}^{(k-1)}\| < \epsilon$, 
%where $\epsilon$ is a specified stopping criteria, stop.
%\end{tabbing}
%\end{algorithm}
