text:
\begin{document}

\maketitle

\begin{abstract}

LT codes are modern and efficient \emph{rateless forward error correction} (FEC) codes with close to channel capacity performance. Nevertheless, in \emph{intermediate range} where the number of received encoded symbols is less than the number of source symbols, LT codes have very low recovery rates.

In this paper, we propose a novel algorithm which significantly increases the intermediate recovery rate of LT codes, while it preserves the codes' close to channel capacity performance. To increase the intermediate recovery rate, our proposed algorithm \emph{rearranges} the transmission order of the encoded symbols exploiting their structure, their transmission history, and an estimate of the channel's erasure rate. We implement our algorithm for conventional LT codes, and numerically evaluate its performance.
\end{abstract}

\section{Introduction}

\emph{LT codes} are modern, efficient, and flexible \emph{rateless forward error correction} (FEC) codes \cite{Luby02lt}. An LT encoder at a source $S$ can potentially generate a limitless number of encoded symbols $c_i,i\in\{0,1,\ldots\}$ from $k$ source symbols $\underline{x} = \{x_1, x_2,\ldots,x_k\}$. The decoder at a destination $D$ can successfully recover $\underline{x}$ from any $k\gamma$ of received encoded symbols, where $\gamma$ is the received \emph{coding overhead}. Let $\gamma_{succ}$ denote the required overhead for a successful decoding. Conventional LT codes can achieve $z\approx1$ with high probability at $\gamma_{succ}$'s slightly larger than one, where $z$ denotes the ratio of the number of recovered source symbols to $k$ at $D$. Note that, for LT codes, $\gamma_{succ}$ is fixed and known \cite{Luby02lt}.

Each LT code is identified by a \emph{degree distribution}. In LT encoding, first an encoded symbol degree $d$ is chosen from a degree distribution, $\{\Omega_{1},\Omega_{2},\ldots,\Omega_{k}\}$, where $\Omega_{i}$ is the probability that $d=i$. This degree distribution is also shown by its generator polynomial $\Omega(y)=\sum_{i=1}^{k}\Omega_{i}y^{i}$. Next, $d$ $x_j$'s are chosen uniformly at random from $\underline{x}$, and are \emph{XOR}ed to generate an encoded symbol $c_i$. Let $\varepsilon\in [0,1)$ denote the channel erasure rate. Therefore, $S$ needs to generate on average $m=\frac{k\gamma_{succ}}{1-\varepsilon}$ encoded symbols so that $D$ can collect $k\gamma_{succ}$ of them for a successful decoding.

In LT decoding at $D$, if a newly delivered $c_i, i \in\{1,2,\ldots,m\}$ has degree-one, it can decode one $x_j,j\in\{1,2,\ldots,k\}$. If the degree of $c_i$ is larger than one then previously recovered $x_j$'s, if any, are removed from $c_i$ to decrease its degree. If the degree of $c_i$ is reduced to one, this symbol is similar to a new degree-one $c_i$. On the other hand, if the remaining degree is larger than one, the symbol is \emph{buffered}. If any new $x_j$ is recovered in the previous step, it is removed from the buffered $c_i$'s. Hence, some buffered $c_i$'s may acquire a degree equal to one, resulting in recovery of further $x_j$'s. This procedure is repeated iteratively until no more degree-one $c_i$ emerges.

\emph{Intermediate range} of LT codes refers to the case where the transmission is still in progress, i.e., $0\leq\gamma<1$ \cite{Sanghavi07intermeidate}. Since conventional LT codes \cite{Luby02lt} are designed to have an almost complete recovery of $\underline{x}$ ($z\approx1$) for certain $\gamma_{succ}\text{'s}>1$, they have low intermediate recovery rates ($z \approx 0$) at $0\leq\gamma<1$. However, in many applications where \emph{partial} recovery of the source symbols from the incomplete received encoded symbols is still beneficial, intermediate recovery rate becomes important. For instance, in multimedia transmission the receiver can play a lower quality of the multimedia contents from the incomplete recovered data. This motivates the design of an algorithm to improve the intermediate recovery rate of existing LT codes.

In order to obtain a high $z$ (as close as possible to $\gamma$) in $0\leq\gamma<1$, each delivered $c_i$ should decode (on average) one $x_j$ instead of being buffered. We propose an algorithm to \emph{rearrange} the transmission order of $c_i$'s based on their transmission history, structure of each $c_i$, and an estimate of $\varepsilon$, such that each delivered $c_i$ can decode one $x_j$ with high probability. By employing our proposed algorithm, the intermediate performance of a given LT code improves significantly, while the recovery rate of $\underline{x}$ at $\gamma_{succ}$ remains intact. In other words, the code remains capacity achieving similar to its original setup.

This paper is organized as follows. Section \ref{related_work} reviews the exiting work on intermediate performance of LT codes. In Section \ref{algorithm}, we propose our novel sorting algorithm. Next, Section \ref{evaluation} reports the performance evaluation of our proposed algorithm. Finally, Section \ref{conclusion} concludes the paper.

\section{Related Work} \label{related_work}

In \cite{Sanghavi07intermeidate}, Sanghavi has shown that $0\leq\gamma<1$ and equivalently $0\leq z<1$ can be divided into three regions of $z \in [0,\frac{1}{2}]$, $z \in [\frac{1}{2},\frac{2}{3}]$, and $z \in (\frac{2}{3},1)$. In each region, the upper bound on $z$ for all rateless codes versus $\gamma$ is formulated, and the optimum $\Omega(y)$'s to gain these upper bounds are provided. The $\Omega(y)$'s provided for each region perform optimally in that specific region only, thus they have a low $z$ compared to $\gamma$ in other two regions. Further, provided $\Omega(y)$'s are not capacity achieving, and they are designed for asymptotic case (infinite $k$) and may not be employed in practice where $k$ is finite.

Authors in \cite{Kamra06growthcodes} have proposed \emph{Growth codes}, which are designed to increase the number of recovered $x_j$'s in intermediate range in wireless sensor networks. In Growth coding, $S$ gradually increases the degree of $c_i$'s on-the-fly based on the value of $z$ (which is known to $S$ by feedbacks received from $D$), such that the instantaneous decoding probability of each delivered $c_i$ is maximized. Growth codes only consider the instantaneous recovery probability of each $c_i$, thus they do not have a close to capacity performance. More importantly, Growth codes require a lot of feedbacks from the receiver.

Authors in \cite{Kim09improved} have proposed to employ multiple feedbacks to transmit $c_i$'s in the order of their \emph{degree} to increase the intermediate recovery rate. Since the decoding of $c_i$'s with lower degrees depends on the recovery of a smaller subset of $x_j$'s, they have a higher probability of decoding an $x_j$ in $D$ at the beginning of transmission. Consequently, an improvement is observed in the intermediate recovery rate of LT codes. However, the algorithm proposed in \cite{Kim09improved} cannot outperform the code of \cite{Kamra06growthcodes, Sanghavi07intermeidate} in intermediate range. Besides, we show that our proposed algorithm always surpasses the algorithm in \cite{Kim09improved}.

In our recent work \cite{talari09intermediate}, we have designed several $\Omega(y)$'s for LT codes to obtain optimum intermediate performance \emph{throughout} the intermediate range rather than a single $\gamma$ \cite{Sanghavi07intermeidate}. The LT codes designed in \cite{talari09intermediate} do not require channel information or feedbacks in contrast to \cite{Kamra06growthcodes, Kim09improved}. However, similar to \cite{Sanghavi07intermeidate} the codes designed in \cite{talari09intermediate} cannot achieve channel capacity.

\section{The Proposed Algorithm} \label{algorithm}

\subsection{Discussion and Idea}

In the previous section, we observed a \emph{trade-off} between being channel capacity achieving and having a high intermediate recovery rate. Some LT codes have optimum intermediate performance but they cannot achieve channel capacity \cite{talari09intermediate, Sanghavi07intermeidate}. On the other hand, the algorithms that are proposed to improve the intermediate recovery rate of LT codes with close to channel capacity such as \cite{Kim09improved} cannot outperform codes of the first group in intermediate range.

The reason for this trade-off is that the codes in the first group have $\Omega(y)$'s that result in generation of a large fraction of low-degree $c_i$'s. Low-degree $c_i$'s have a higher probability of decoding a source symbol when $z$ is small. However as $z$ grows, some of the received $c_i$'s become redundant and cannot recover any $x_j$ due to earlier recovery of all their adjacent $x_j$'s. Therefore, these codes cannot have a close to channel capacity performance.

On the other hand, $\Omega(y)$'s of codes with close to channel capacity performance result in generation of $c_i$'s with much higher degrees. Decoding of high-degree $c_i$'s depends on the recovery of many $x_i$'s at $D$, thus in intermediate range these $c_i$'s are mostly buffered for a later decoding. Therefore, for these codes $z$ does not grow considerably with $\gamma$. However, the buffered high-degree $c_i$'s are simultaneously decoded together and recover $\underline{x}$ at $\gamma_{succ}$, which makes the code capacity achieving.

We can see that in LT codes with capacity achieving performance, each $c_i$ eventually decodes close to one $x_j$ on average, since all $k$ $x_j$'s are decoded from $k \gamma_{succ}$ $c_i$'s, which are slightly more than $k$ symbols. Consequently, if $c_i$'s are transmitted in the \emph{order} that they are \emph{decoded}, we can significantly improve the intermediate recovery rate of \emph{capacity achieving} LT codes. We propose our algorithm to transmit $c_i$'s in this \emph{correct order} for two cases of constant and varying $\varepsilon$.

\subsection{Algorithm for Constant $\varepsilon$}

Similar to \cite{Gummadi08relaying}, we assume that an estimate of the channel erasure rate, $\varepsilon$, is available at $S$. Since we assume that $D$ generates no feedbacks, our algorithm is designed to be implemented on the encoder side. Therefore, our algorithm can exploit the information available in $S$ only, and the decoder remains intact.

In conventional LT coded symbol transmission, $S$ generates $m=\frac{k\gamma_{succ}}{1-\varepsilon}$ random $c_i$'s from a capacity achieving degree distribution $\Omega(.)$ such as \emph{Robust-Soliton} degree distribution \cite{Luby02lt}. After transmission, $D$ receives $k\gamma_{succ}$ encoded symbols, which results in a successful decoding of $\underline{x}$. This method has a poor intermediate recovery rate as we later show.

In our proposed scheme, after generating $m$ $c_i$'s, $S$ performs as follows. $S$ maintains a probability vector $\underline{\rho}=[\rho(1),\rho(2),\ldots,\rho(k)]$, in which $\rho(j)$ represents the probability that $x_j$ is still \emph{not recovered at} $D$. Clearly, $S$ sets $\underline{\rho}$ to an all-one vector when the transmission has not started yet since no $x_j$ is recovered at $D$. At each transmission and based on $\underline{\rho}$, $S$ finds a $c_i$ that has the highest probability of recovering an $x_i$ at $D$ (as we later describe in Algorithm \ref{sorting_algth}). Next, $S$ transmits $c_i$ and updates $\rho(j), j\in \mathcal{N}(c_i)$, where $\mathcal{N}(c_i)$ represents the set of indices of $x_i$'s \emph{XOR}ed together to generate $c_i$. $S$ continues until all $m$ $c_i$'s are transmitted.

An encoded symbol $c_i$ with degree $d$, i.e., $|\mathcal{N}(c_i)|=d$, where $|.|$ represents the cardinality of a set, can recover a source symbol $x_j$ iff all $x_k,k \in \{\mathcal{N}(c_i)-j\}$ has already been recovered. Let $\underline{p_{dec}}=\{p_{dec}(1), p_{dec}(2),\ldots,p_{dec}(m)\}$, where $p_{dec}(i)$ is the probability that $c_i$ can recover a source symbol at $D$, and $p_{dec}(i)=0$ if $c_i$ has been previously transmitted. Since at the beginning of transmission no source symbol is still recovered, we have $p_{dec}(i)=0$ if $|\mathcal{N}(c_i)|>1$, i.e., $c_i$'s with degrees larger than one cannot decode any $x_j$ at $D$. Besides, for $|\mathcal{N}(c_i)|=1$ we have $p_{dec}(i) = (1-\varepsilon)$, i.e. only degree-one encoded symbols that are not erased by the channel (with probability $(1-\varepsilon)$) can recover a source symbol.

We can see that at the beginning of transmission degree-one $c_i$'s have the highest probabilities of decoding an $x_j$ at $D$. Therefore, $S$ transmits degree-one $c_i$'s with $\mathcal{N}(c_i)=\{j\}$, and updates $\rho(j)= \varepsilon \rho_{old}(j)$, where $\rho_{old}(j)$ is the value of $\rho(j)$ before $c_i$ was transmitted.

Next, we consider a degree-two $c_i, \mathcal{N}(c_i)=\{j,k\}$. $c_i$ can recover $x_j$ with probability $(1-\varepsilon) (1-\rho(j))\rho(k)$, which is the probability that $c_i$ is not dropped on channel, $x_j$ has not been recovered previously, and the $x_k$ has already been recovered. Similarly, $c_i$ can recover $x_k$ with probability $(1-\varepsilon) (1-\rho(k))\rho(j)$. Consequently, $p_{dec}(i) = (1-\varepsilon) [(1-\rho(k))\rho(j)+(1-\rho(j))\rho(k)]$. Assume $\forall l\not =i, p_{dec}(i) > p_{dec}(l)$, i.e. $c_i$ has the highest probability of decoding an $x_j$ at $D$. Therefore, $S$ transmits $c_i$ next and sets $\rho(j)= \rho_{old}(j) (1-(1-\varepsilon) (1-\rho_{old}(k)))$ and $\rho(k)= \rho_{old}(k) (1-(1-\varepsilon) (1-\rho_{old}(j)))$.

Further, we consider $c_i, |\mathcal{N}(c_i)|=d$. If $c_i$ is successfully delivered to $D$, it can possibly decode $x_j, j \in \mathcal{N}(c_i)$. Similar to low degree $c_i$'s, $x_j$ can be decoded with probability $(1-\varepsilon) \rho(j) \prod \limits_{v \in \mathcal{N}(c_i), v \not = j}(1-\rho(v))$. Therefore, $p_{dec}(i) = (1-\varepsilon) \sum \limits_{l \in \mathcal{N}(c_i)} [\rho(l) \prod \limits_{v \in \mathcal{N}(c_i),v \not = l}(1-\rho(v))]$. If $p_{dec}(i)>p_{dec}(l),\forall l\not =i$, $S$ transmits $c_i$ and updates $\rho(j)=\rho_{old}(j)[1-(1-\varepsilon)\prod \limits_{v \in \mathcal{N}(c_i),v \not = j}(1-\rho_{old}(v))], j \in \mathcal{N}(c_i)$.

We summarize our proposed sorting scheme in Algorithm \ref{sorting_algth}. The output of Algorithm \ref{sorting_algth} is the suitable rearranged transmission order $\underline{\pi}$ of $c_i$'s that substantially improves $z$ in $0\leq\gamma<1$. In this algorithm, $argmax(\underline{p_{dec}})$ is a function that returns $i$ where $\forall j\not = i, p_{dec}(i) > p_{dec}(j)$. Further, if $c_i$ and $c_l$ both have the highest probability of decoding of an $x_j$, i.e., $p_{dec}(l) = p_{dec}(i)$, then $argmax(\underline{p_{dec}})$ returns the index of $c_i$ or $c_l$, whichever has the \emph{lowest} degree. Further, if $c_i$ and $c_l$ have equal degrees (similar to degree one $c_i$'s at the beginning of transmission), $argmax(\underline{p_{dec}})$ randomly returns one of the indices. Clearly, we are assigning earlier transmission priority to lower degree symbols, which we later show is important for $\varepsilon\rightarrow 1$.

\begin{algorithm}
\caption{The proposed symbol sorting algorithm}
\label{sorting_algth}
\begin{algorithmic}
\STATE Initialize: $\underline{\pi}=\emptyset$, $\underline{\rho}=\{1\}_{1\times k}$
\WHILE{$|\underline{\pi}|<m$}
    \FOR{$j = 1$ to $m, j \not\in \underline{\pi}$}
        \STATE $p_{dec}(j)=(1-\varepsilon)\sum \limits_{l \in \mathcal{N}(c_i)} [\rho(l) \prod \limits_{v \in \mathcal{N}(c_i),v \not = l}(1-\rho(v))]$
    \ENDFOR
    \STATE $i^* = argmax(\underline{p_{dec}})$
    \STATE  $\underline{\pi} =[\underline{\pi} , i^*]$
    \FOR{$j \in \mathcal{N}(c_{i^*})$}
        \STATE $\rho(j)=\rho_{old}(j)[1-(1-\varepsilon)\prod \limits_{v \in \mathcal{N}(c_{i^*}),v \not = j}(1-\rho_{old}(v))]$
    \ENDFOR
\ENDWHILE
\end{algorithmic}
\end{algorithm}

Our proposed algorithm increases the coding complexity of LT codes from $O(k\log k)$ \cite{Luby02lt} to $O(k^2)$, while it does not deteriorate the decoding complexity.

Further, in our proposed algorithm all $c_i$'s need to be generated and sorted before the transmission can start in contrast to the conventional LT coding where each $c_i$ can be independently transmitted upon generation. Therefore, some delays may be introduced.

However, this delay can be easily eliminated with the following procedure. Clearly, the performance of our proposed scheme is independent of $\underline{x}$'s contents and only depends on $\mathcal{N}(c_i),i\in\{1,2,\ldots,m\}$ and $\underline{\pi}$. Therefore, before the transmission starts, $S$ generates $c_i$'s from a dummy $\underline{x}$, and obtains an off-line version of $\underline{\pi}_{\text{off-line}}$. $S$ saves $\mathcal{N}_{\text{off-line}}(c_i)$, and $\underline{\pi}_{\text{off-line}}$ for a later use. When the actual encoding starts, $\underline{x}$ of interest replaces the dummy $\underline{x}$, and $S$ generates $m$ $c_i$'s \emph{in the order} dictated by  $\underline{\pi}_{\text{off-line}}$, \emph{XOR}ing $x_j , j\in\ \mathcal{N}_{\text{off-line}}(c_i)$. In this way, each $c_i$ can be transmitted upon generation and the delay is completely eliminated at the cost of some data storage.

In the next section, we extend Algorithm \ref{sorting_algth} to the case where $\varepsilon$ varies.

\subsection{Algorithm for Varying $\varepsilon$}

Assume that $S$ has generated $m$ $c_i$'s considering $\varepsilon$. Assume that the erasure rate of the channel changes to $\varepsilon_{new}$ when $\frac{k\gamma_{c}}{1-\varepsilon}$ symbols has already been transmitted so that $\frac{k(\gamma_{succ}-\gamma_{c})}{1-\varepsilon}$ $c_i$'s are still remaining to be transmitted.

If  $\varepsilon_{new} > \varepsilon$, less than $k\gamma_{succ}$ $c_i$'s would be collected by $D$, making the full decoding impossible. In this case, $S$ generates $t=(\frac{1}{1-\varepsilon_{new}}-\frac{1}{1-\varepsilon})k(\gamma_{succ}-\gamma_{c})$ \emph{new} $c_i$'s, and adds them to the queue of $c_i$'s to be transmitted to ensure the delivery of $k\gamma_{succ}$ $c_i$'s to $D$. Next, $S$ \emph{rearranges} all $c_i$'s in the queue according to Algorithm \ref{sorting_algth} and continues the transmission.

In the second case for $\varepsilon_{new} < \varepsilon$, $S$ randomly drops $1-\frac{1-\varepsilon}{1-\varepsilon_{new}}$ fraction of remaining $c_i$'s from the transmission queue. This limits the number of delivered $c_i$'s to $k\gamma_{succ}$, hence the code maintains its close to channel capacity performance.

If the erasure rate of the channel varies several times, the same procedures are followed after each change. We assume that $\varepsilon_{new}$ is known to $S$, which can be realized by receiving few feedbacks from the receiver.

Note that the symbol dropping procedure described above is similar to \emph{puncturing LDPC} \cite{Pishro-Nik05Nonuniform} and \emph{turbo} codes \cite{Douglas97Rate-compatible} to achieve a certain higher coding rate for these codes.

\vspace{-1mm}

\section{Evaluation of the Proposed Algorithm} \label{evaluation}

\vspace{-1mm}

We emphasize that our proposed algorithm can be applied to an LT code with any degree distribution to increase its intermediate recovery rate when $\varepsilon$ is available at $S$. The advantage of this algorithm is that if the code is capacity achieving, it remains capacity achieving after employing our algorithm.

To evaluate the performance of our proposed algorithm, we implement it for two well-known LT codes. The first code we employ is the LT code used in \emph{Raptor codes} \cite{Shokrollahi06raptor} with degree distribution $\Omega(y)$ given below and $k=1000$. Since low error floors cannot be achieved in intermediate range, the precoding phase of Raptor codes can be skipped.

\vspace{-3mm}

\begin{equation*}
\footnotesize{
\begin{split}
\Omega(y)=0.00797x+0.49357x^{2}+0.16622x^{3}+0.07265x^{4}+0.08256x^{5}\\
             +0.05606x^{8}+0.03723x^{9}+0.05559x^{19} + 0.02502x^{65}+0.00314x^{66}.
\end{split}
}
\end{equation*}

The second code is an LT code with $k=1000$ and Robust-Soliton degree distribution \cite{Luby02lt} with parameters $c=0.05$ and $\delta=0.01$.

Figures \ref{result_shok} and \ref{result_rsd} illustrate the improvement made in $z$ for aforementioned Raptor and LT codes employing our proposed sorting algorithm.

\begin{figure}[h]
\includegraphics[width=3.4in,angle=0]{result_shok.eps}
\caption{Improvement made in the intermediate recovery rate of a Raptor code proposed in \cite{Shokrollahi06raptor} employing proposed algorithms for various $\varepsilon$'s versus $\gamma$.}
\label{result_shok}
\end{figure}

\begin{figure}[h]
\includegraphics[width=3.4in,angle=0]{result_rsd.eps}
\caption{Improvement made in the intermediate recovery rate of an LT code \cite{Luby02lt} employing proposed algorithms for various $\varepsilon$'s versus $\gamma$.}
\label{result_rsd}
\end{figure}

From Figures \ref{result_shok} and \ref{result_rsd}, we can see a significant improvement in $z$ for both codes. For instance at $\gamma=1$, for the first and the second code we can see $382\%$ and $843\%$ improvement in $z$, respectively. We can observe that the amount improvement depends on the value of $\varepsilon$. The rationale behind this is that when $\varepsilon$ decreases our algorithm can estimate the recovery probability of $x_j$'s more accurately, which results in a more efficient reordering of $c_i$'s. As $\varepsilon$ becomes larger, the ordering of $c_i$'s becomes less accurate. It is worth noting that the conventional transmission of LT codes results in the same curve of $z$ regardless of $\varepsilon$'s value.

\subsection{Upper and Lower Bounds on Algorithm's Performance}

As $\varepsilon\rightarrow 1$, $S$ cannot make a meaningful estimation about the recovery of $x_j$'s at $D$, and $\underline{\rho}$ always remains an all-one vector. Since in our proposed algorithm, $c_i$'s with lower degrees have higher priority of transmission, for $\varepsilon\rightarrow 1$ our algorithm approaches to the case where $c_i$'s are only sorted based on their degrees in $\underline{\pi}$. Consequently, in this case our proposed algorithm boils down to the algorithm proposed in \cite{Kim09improved}. As a result, the improvement made by our algorithm is lower bounded by the results of \cite{Kim09improved}. Moreover, for $\varepsilon\rightarrow 0$, $S$ can estimate which $x_j$'s are being decoded with a high accuracy, thus more exact $\underline{\pi}$ can be acquired, and the intermediate performance approaches the ideal upper bound, i.e., $z=\gamma$.

The upper and the lower bounds on our proposed scheme are depicted in Figures \ref{result_shok_bounds} and \ref{result_rsd_bounds} for distribution of Raptor and LT codes, respectively. For comparison, we have also provided the recovery rate curves of Growth codes \cite{Kamra06growthcodes} and conventional LT transmission for the same codes. Further, the lower bound illustrates the performance of the scheme proposed in \cite{Kim09improved}.

\begin{figure}[h]
\includegraphics[width=3.4in,angle=0]{result_shok_bounds.eps}
\caption{The upper and the lower bounds on the improvement made by our proposed scheme for Raptor codes \cite{Shokrollahi06raptor}, compared to the intermediate recovery rate of Growth codes and conventional LT codes.}
\label{result_shok_bounds}
\end{figure}

\begin{figure}[h]
\includegraphics[width=3.4in,angle=0]{result_rsd_bounds.eps}
\caption{The upper and the lower bounds on the improvement made by our proposed scheme for LT codes \cite{Luby02lt}, compared to the intermediate recovery rate of Growth codes and conventional LT codes.}
\label{result_rsd_bounds}
\end{figure}

Figures \ref{result_shok_bounds} and \ref{result_rsd_bounds} show considerable improvement for intermediate recovery rate of both employed LT codes for any value of $\varepsilon$. We can see that Growth codes outperform our propose algorithm only for a small region of $\gamma$, while they cannot have a close to channel capacity performance.

\subsection{Algorithm's Performance for Varying $\varepsilon$}

As described earlier, if $\varepsilon$ increases, $S$ needs to generate some new $c_i$'s. Assume that $S$ is transmitting $m$ Raptor encoded $c_i$'s (without precoding) generated for $\varepsilon=0.3$. Also assume that $\varepsilon$ increases to $\varepsilon_{new}=0.5$ at $\gamma_c = 0.5$. Based on our proposed algorithm, $S$ adds $t = \lceil 0.5714 k (\gamma_{succ}-\gamma_c)\rceil$ new $c_i$'s to the queue, and updates $\underline{\pi}$ accordingly. We have depicted $z$ versus $\gamma$ for this case in Figure \ref{varying}. For comparison, we have also depicted $z$ for constant $\varepsilon\in\{0.3,0.5\}$.

\begin{figure}[h]
\includegraphics[width=3.2in,angle=0]{varying.eps}
\caption{The intermediate performance of a Raptor code employing our proposed algorithm for the case where $\varepsilon$ increases from $0.3$ to $0.5$ with $\gamma_c=0.5$.}
\label{varying}
\end{figure}

Figure \ref{varying} shows that fluctuation in the $\varepsilon$ is well compensated by our algorithm. Since the newly generated $c_i$'s at $\gamma_c$ disturb the original sorting order, a slight degradation in $z$ can also be observed. These new $c_i$'s might have been transmitted earlier than $\gamma_c$ if they were present in the queue from the beginning. In spite of all the fluctuation in intermediate performance, the code remains capacity achieving since $D$ can collect $k\gamma_{succ}$ encoded symbols in total.

\subsection{Comparison with Fixed-Rate Codes}

Since we assumed that an estimate of the channel loss rate, $\varepsilon$, is available at the source, $S$ can employ existing \emph{fixed-rate} codes instead of LT codes to encoded $\underline{x}$. Therefore, we need to compare the performance of our proposed scheme with the intermediate recovery rate of existing fixed-rate codes.

\emph{Systematic irregular repeat-accumulate (Systematic IRA)} codes \cite{jin00irregular} are capacity achieving fixed-rate codes on erasure channels, which can provide a high intermediate recovery rate compared to other existing fixed-rate codes. Figure \ref{IRA_comparison} shows the intermediate recovery rate of a systematic IRA code with rate $R=0.5$ versus our proposed algorithm employing Raptor codes without precoding for two cases of $\varepsilon=0.44$ and $\varepsilon=0.46$ with $k=10000$.

\vspace{-3mm}

\begin{figure}[h]
\includegraphics[width=3.2in,angle=0]{IRA.eps}
\caption{The intermediate recovery rate of systematic IRA codes \cite{jin00irregular} compared to our proposed algorithm employing Raptor codes \cite{Shokrollahi06raptor} for $k=10000$.}
\label{IRA_comparison}
\end{figure}
\vspace{-2mm}

From Figure \ref{IRA_comparison}, we can see that the employed systematic IRA code outperforms our proposed scheme for small values of $\gamma$. This high intermediate recovery is due to the systematic part of these codes which results in transmission of uncoded $\underline{x}$. When the systematic part of the code is being transmitted, each delivered $c_i$ is itself a source symbol, hence $z$ and $\gamma$ are equal. However, we can observe that when the transmission of systematic part ends, $z$ does not increase linearly with $\gamma$ anymore and our proposed scheme outperforms this code for larger values of $\gamma$. Further, as $\varepsilon\rightarrow 1$ the systematic part cannot be delivered and the gain from the systematic part is eliminated. In this case, our proposed algorithm always outperforms systematic IRA codes.

Furthermore, fixed-rate codes seriously suffer from their fixed rates since they cannot compensate slight variations in $\varepsilon$. For instance, we can see that when $\varepsilon$ is increased from $0.44$ to $0.46$, $z$ decreases from $1$ to $0.6150$ at the end of transmission. However, as observed in Figure \ref{varying} our algorithm employed along with a Raptor code exhibit good performance in spite of very large variations in $\varepsilon$.

\vspace{-1mm}

\section{Conclusion}\label{conclusion}

In this paper, we proposed an algorithm to increase the intermediate recovery rate of capacity achieving LT codes. In our proposed algorithm, the transmitter exploits the structure of LT encoded symbols, a history of previously transmitted encoded symbols, and an estimate of channel's erasure rate to sort the transmission order of the encoded symbols to gain a high intermediate recovery rate.

Using numerical results, we showed that our algorithm can increase the intermediate recovery rate of LT codes to a great extent while the code remains capacity achieving. We also showed that this algorithm performs well for fluctuating channel erasure rates.

To extend this work, we intend to implement the proposed algorithm for real multimedia transmission and observe the improvement made in the quality of the received stream.

\vspace{-1mm}

\bibliographystyle{ieeetr}
\bibliography{ref}

\end{document}

config: arxiv


TITLE: I need to prove a simple proposition about polynomials
QUESTION [0 upvotes]: I am working on systems theory and I have reduced one problem that I have into a simple polynomial proposition.
All (non-constant) polynomials (with integer degrees) in which the sum of positive coefficients is exactly 1 and the sum of negative coefficients is in [-1,0] have maximum one solution in the range [0-1].
Can somebody either help me or find one polynomial that disproves the proposition?
Thanks a lot in advance...
SB7

REPLY [0 votes]: $\frac{1}{2}x^4-\frac{1}{2}x^2-\frac{1}{2}x+\frac{1}{2}$ has two roots in $[0,1]$ as well, there are many examples. (Unless I don't understand what you mean by integer degree).


TITLE: GR time dilation cares only about local curvature, right?
QUESTION [2 upvotes]: Regarding this remark on Worldbuilding SE and the discussion leading up to it:
Can someone properly knowledgeable on the subject please explain whether the time dilation due to being in a gravitation field is in fact due only to the curvature of spacetime in the region where the clock is found, or if a flat "plateau" somehow knows that it is just a shelf along a larger slope?
In particular, would using compensating masses to form a nearly flat spacetime in a small region (as described by Robert L. Forward) cause a clock to experience only "flat spacetime" or would the clock be slowed anyway because this prepared region is itself in a gravitation well?
E.g. a prepared flat region in the inner solar system vs. one out in interstellar space: is there any difference in time rate when the local curvature of spacetime surrounding the clock is the same?
Edit for those who found it unclear: look at previous paragraph for a specific physics question.

REPLY [1 votes]: A simple way of illustrating gravitational time dilation is using the gravitational redshift. If I shine light of frequency $\nu_{me}$ from the Earth's surface up to you on the International Space Station then when you received it the light will have reduced in frequency to $\nu_{you}$. We explain this by saying that my time is running more slowly than yours, so I measure the time between peaks in the light waveform to be smaller than you do. The relative time dilation between me and you is then simply:
$$ \frac{\tau_{me}}{\tau_{you}} = \frac{\nu_{you}}{\nu_{me}} $$
This should make it obvious that time dilation is not a local phenomenon because it depends on where I am, where you are and what is in between us.
In most cases we are likely to encounter gravity is relatively weak, and we can describe the spacetime geometry using the weak field approximation. In this case the relative time dilation for any two observers is given by:
$$ \frac{\Delta t_1}{\Delta t_2} = \sqrt{1 - \frac{2\Delta\Phi}{c^2}} $$
where $\Delta\Phi$ is the difference in the (Newtonian) gravitational potential difference between the two observers.


TITLE: Does $f_n(x) = \sum_{k=0}^n x^k$ converge uniformly where f maps from $[0, 3/4]$ to reals
QUESTION [0 upvotes]: Does  $f_n(x) = \sum_{k=0}^n x^k$ converge uniformly where f maps from $[0, 3/4]$ to reals?
I found the pointwise limit to be $\frac{1}{1-x} $for $ x \in (0, 3/4]$ and $0$ for $x=0$. I do not know how to proceed

REPLY [1 votes]: Observe by the geometric series summation formula that
$$
f_n(x)=\frac{1-x^{n+1}}{1-x},
$$
and also let $f(x)=1/(1-x)$ denote the limit function. Then
$$
|f(x)-f_n(x)|=\frac{|x|^{n+1}}{|1-x|}.
$$
To establish uniform convergence, we must show that for every $\epsilon>0$ there exists an $N$ such that for all $n>N$, we have $|f(x)-f_n(x)|<\epsilon$ - where the bounds do not depend on the value of $x$. In this case, observe from $1-x\geq \tfrac14$ and $x\leq \tfrac34$ that we obtain the uniform bound
$$
|f(x)-f_n(x)|\leq 4 \left(\frac{3}{4}\right)^{n+1},
$$
allowing you to choose the appropriate $N$ given $\epsilon$, as required.


TITLE: computation of a line integral
QUESTION [1 upvotes]: in an excersise I have to compute a few line integrals but one of them I can't solve. It is not even written as a line integral but the others are. I am talking about:
$$\int_0^{2\pi}e^{it+e^{it}}dt.$$
Is there any possibility to write that integral in the form $\int_{\gamma}f(z)dz$ with a line $\gamma:[0,2\pi]\rightarrow\mathbb{C}$? Maybe then I can use some results of Cauchy.
Thanks and regards
N.Sch

REPLY [1 votes]: Writing $z = e^{it}$, you can recognize this as
$$\int_{\gamma} e^z \, \frac{dz}{i}$$
where $\gamma$ is the unit circle oriented counterclockwise.

REPLY [0 votes]: Let $f = \exp$ be the exponential function, and consider the integral of this function over the unit circle $S^1$. This circle can be parameterized by the curve
$$ \gamma(t) = e^{it} $$
It then follows that
$$ \int_{S^1} f(z) dz = \int_0^{2 \pi} (f \circ \gamma)(t) \gamma'(t) dt = i \int_0^{2\pi} e^{e^{it}} e^{it} dt $$
So
$$ \int_0^{2\pi} e^{e^{it}} e^{it} dt = -i \int_{S^1} f(z) dz $$


TITLE: Choice of Axioms
QUESTION [2 upvotes]: I was studying the Foundations of Mathematics and a question was bugging me for quite some while. Does our choice of axioms depend upon the system we are studying? Could it be that if some statement is an axiom in one system it is a theorem in another? For example, I was learning Peano Axioms where we prove theorems establishing commutativity, associativity, distributivity etc. of number systems (natural, integers, rational, real, etc.) and then when I began Abstract Algebra, it treated these propositions as axioms. Even in Elementary Algebra these properties of natural or real numbers are taken to be demonstrable axioms. So is our choice of axioms subjective? It would be great if someone could shed some light on this.

REPLY [0 votes]: The field you are talking about is call Reverse Mathematics, and there is a number of related questions under the tag reverse-math. A typical procedure would be, instead of using the usual axioms to prove, say, the extremal value theorem (EVT), postulating the EVT to see what former "axioms" can be deduced from it or are equivalent to it. Here a certain bare minimum is assumed which is usually denoted RCA$_0$; see reverse mathematics.


TITLE: Probability of at least one k-length run of heads in n coin flips
QUESTION [0 upvotes]: Question says it all - the flips are uniform, independent events that can be Heads or Tails. My strategy has been to find the probability of having no k-length runs of heads in n flips, and take 1 minus that. I think the solution has to do with recursion/induction, but I’m not sure exactly how it works...

REPLY [0 votes]: The probability of any particular outcome of heads and tails in n flips is $0.5^n$. The number of outcomes is $\frac{1}{0.5^n}$. Obviously $0.5^n\cdot \frac{1}{0.5^n}=1$ so the probability of not getting a single k length run of heads is:
$$1 - (n-(k-1))0.5^n$$
Where $(n-(k-1))$ is the number of ways to get a k length run of heads in n flips.
In considering more than one run of k heads where k is considerably smaller than n, we add in extra terms to subtract from 1 until a maximum possible number of k runs is reached. The following is the probability for one and two runs of k. Another series of terms is required for the probability of more than two runs of k.
$$1 - (n-(k-1))0.5^n - ((n-(k+1))-(k-1))0.5^n - ((n-(k+2))-(k-1))0.5^n..........(1)0.5^n$$
Additionally, runs of heads greater than k are not considered.


TITLE: Comparing two explanatory variables explanatory variables
QUESTION [1 upvotes]: I am new to stats, I am working on a problem and I need some direction.
I have a binary dataset and I have three explanatory variables x,y, and z that I want to analyze. There is no guarantee of independence between any of them. I would like to perform a test to compare the explanatory power of two of them e.g x,y. I think that if for example, I do logistic regression with x,z, and another y,z and compare the goodness of fit or I compare coefficients in one logistic regression with all variables, it wouldn't work because none of the variables are independent.
So I am not sure how to proceed.
Please let me know if the post needs more information/clarity.

REPLY [0 votes]: You can add an interaction term. So if your response is w, your formula would be $w \sim x + y + x*y$ if you just wanted to use x and y as predictors, or $w \sim  x + y + z + x*y + y*z + x*z + x*y*z$. Usually three factors interactions aren’t common so you could omit the last term.


TITLE: Are all metric translations isometries
QUESTION [1 upvotes]: Let $(M, d)$ be a metric space. I define a translation of $M$ to be a function $f$ from $M$ to $M$ such that $d(x, f(x)) = d(y, f(y))$ for all $x$ and $y$ in $M$. My conjecture is that every translation on $M$ is an isometry under the same metric.
Can anyone prove this, or give me a counterexample?

REPLY [2 votes]: If we don't require continuity of $f$, we have a counterexample $f\colon \mathbb{R}\to\mathbb{R}$ with
$$f(x) = \begin{cases}x + 1 &, x \in \mathbb{Q}\\ x-1 &, x \notin\mathbb{Q}.\end{cases}$$
For continuous $f$, we get a counterexample by considering a discrete metric space $(M,d)$, with $d(x,y) = 1$ for $x\neq y$, that contains at least three points. Let $x_0,x_1$ be two distinct points, and
$$f(x) = \begin{cases}x_0 &, x \neq x_0\\ x_1 &, x = x_0. \end{cases}$$
Then $d(x,f(x)) = 1$ for all $x$, but $d(f(x),f(y)) = 0$ for all $x,y \in M\setminus \{x_0\}$.


TITLE: Evaluate $(1-\frac1{2^2})(1-\frac1{3^2})\ldots(1-\frac1{2015^2})$
QUESTION [0 upvotes]: Evaluate
$$\prod_{k=2}^{2015} \left(1-\frac1{k^2}\right) = \left(1-\frac{1}{2^2}\right)\left(1-\frac{1}{3^2}\right)\ldots\left(1-\frac{1}{2014^2}\right)\left(1-\frac{1}{2015^2}\right)$$

REPLY [0 votes]: I saw really no effort put into the question.
Even if you are not really good at maths, you can evaluate the first few products. Always helps, just calculating a few results so you get an idea how the problem works.
The first product is 3/4.
Multiplied by 8/9 = 2/3.
Multiplied by 15/16 = 5/8.
Multiplied by 24/25 = 3/5.
Multiplied by 35/36 = 7/12.
Multiplied by 48/49 = 4/7.
Now do you see a pattern? If you multiply numerator and denominator each by 2 in every second line, you get the pattern
3/4, 4/6, 5/8, 6/10, 7/12, 8/14
Now you can write down a hypothesis, and prove it using induction. Really easy. Even for a beginner.
What you added in a comment to your question is really nonsense. Look at the product. Each of the numbers multiplied is very, very close to 1. All the numbers that you multiply in your comment are tiny. It should have been obvious that this cannot be right.


TITLE: A subset that is closed under multiplication but not addition?
QUESTION [3 upvotes]: I can't get my head around subspaces despite having studied on them quite a lot. Here goes:
The problem statement, all given variables and data
Give an example of a non-empty subset U of R^2 such that U is closed under scalar multiplication but is not a subspace of R^2.
Attempt at a solution
So a set such as [x1, x2] | x1 >= 0 is NOT closed under scalar multiplication but is closed under addition. I get this but I cannot find an answer in reverse.
A wild guess is x1 + x2 + 2 = 0. I could multiply it and 2x1 + 2x2 + 4 = 0 would still hold true. However, if I wish to perform an addition on it, I don't exactly understand what I add to what equation exactly.. it's all so bogged up. Help would be appreciated. :-)

REPLY [3 votes]: Since you only need a subset for this example- we can "cheat" a bit. We shall define our subset $U$ to be the union of $x$ axis, and the $y$ axis, i.e.
$U$= $\{(x,0):x \in \mathbb{R}\} \cup \{(0,y):y \in \mathbb{R}\}$


TITLE: Wave Function Question
QUESTION [1 upvotes]: For part a) the value of a=60° and k=2, therefore the answer is 2sin(x+60)°.
I am struggling in part b) on how to find an answer, the answer is := min. of 2 @ x=30°.
Can someone help on how to solve part b).

REPLY [1 votes]: Hint: To minimize $f$ you want to maximize the denominator.  $\sin \theta$ is maximized at $\theta = 90$.


TITLE: Exercise 6, section 2.7 of Hammack's Book of Proof
QUESTION [0 upvotes]: I have two questions about exercise 6 of section 2.7 of  Book of Proof from Hammack.
The exercise asks the reader to write the following statement in English and then to say if the statement is true or false.
$$ \exists \ n \in \mathbb{N}, \forall X \in \mathcal{P}(\mathbb{N}), \vert X \vert < n $$
My first question is whether or not I am reading this correctly.
I interpret this statement as saying "there exists a positive integer greater than the cardinality of all subsets of $\mathbb{N}$". Is that correct?
My second question, assuming my understanding of the statement is correct, is whether the statement is true or not.
My fist line of thought was that the statement is false because for any positive integer $n$ there will always be a subset of $\mathbb{N}$ with cardinality equal or larger than $n$ since the set of natural numbers is infinitely large. However, the converse reasoning seems as valid: for any subset of $\mathbb{N}$, however large it is, there will always exist a positive integer larger than this subset's cardinality. It seems to me now that we cannot really say if the statement is true or false: the statement itself is sort of paradoxical. Would this assessment be the correct answer? Am I missing something here?

REPLY [1 votes]: My fist line of thought was that the statement is false because for any positive integer $n$ there will always be a subset of $\mathbb{N}$ with cardinality equal or larger than $n$ since the set of natural numbers is infinitely large.

Yes, pretty much ... but you want to be a little more explicit/concrete about exactly what subset you can always pick such that its cardinality is not smaller than $n$ for whatever number $n$ you may have picked initially. But it is easy to see what works: you can pick any infinite-sized subset, such as all even numbers ... or even just simply $\mathbb{N}$ itself ... which is also an element of $\mathcal{P}(\mathbb{N})$

However, the converse reasoning seems as valid: for any subset of $\mathbb{N}$, however large it is, there will always exist a positive integer larger than this subset's cardinality.

First of all, you are now considering the truth of the statement:
$$ \forall X \in \mathcal{P}(\mathbb{N}), \exists \ n \in \mathbb{N}, \vert X \vert < n $$
which is a different statement. But, like the original one, this one is also false. Again, you can use any infinite-sized subset as a counterexample.


TITLE: Indiscernibles and colorations.
QUESTION [1 upvotes]: Let $L$ be a language, and $(X,\leq)$ be a total order contained in an $L$-structure $\frak{A}$. Now if we denote by $[X]^n$ the set of $n$-sized sequences in $X$ and consider a set $\Gamma$ of $L$-formulas in $n$ free variables then is it possible to give a coloration $c:[X]^n\longrightarrow 2^{|L|}$ such that:
$$
c(\bar{a})=c(\bar{b}) \Longleftrightarrow \frak{A} \models \varphi[\bar{a}]\longleftrightarrow \varphi[\bar{b}] \tag{i}
$$
For all $\varphi \in \Gamma$
How to give a coloration such that indiscernible sets are of homogeneous color?. Is it possible to do so without actually stating $(i)$ in the definition?. Is there any way to do this under some other set of hypothesis maybe?

REPLY [3 votes]: It's easier to define $c$ to map $[X]^n$ into the set ${}^\Gamma\{\top,\bot\}$ of functions from $\Gamma$ to truth values, rather than into the cardinal $2^{|L|}$, and it doesn't matter for the result since there's a bijection between ${}^\Gamma\{\top,\bot\}$ and a subset of $2^{|L|}$.  So let $c(\overline a)$ be the function assigning to each $\phi\in\Gamma$ the truth value of $\phi(\overline a)$ in $\mathfrak A$.


TITLE: What is a physically intuitive description of the "dot product" and "cross product"?
QUESTION [1 upvotes]: I just learned about the definition of the vector dot product and cross product.
Almost every source says that the dot product is the:

Product of the $\cos(\theta)$ component of a vector along another vector times the magnitude of the other vector.

While cross product is the:

Area of the parallelogram formed by the two vectors to be crossed.

There seems to be no meaning of what is done and why it is done so. Why do we need to find the area of the parallelogram formed when two vectors are cross produced?
It may seem as a weird question for people who have had a long bond with vectors but this question would have surely popped in the minds of everyone who is beginning the journey of vectors!
What I would like to know is the physical and a more intuitive meaning of the dot and cross products.

REPLY [2 votes]: There are many ways of answering this question - here is one which I find interesting.
If you think of vectors as columns of 3 numbers, ${\bf a}=(a_1,a_2,a_3)$ and ${\bf b} =(b_1,b_2,b_3)$, and you think about how you might multiply them together in some sensible way, you could write $\tilde{\bf a} {\bf b}$
or ${\bf a} \tilde {\bf b}$ (where the tilde means taking the transpose).
The first way gives the scalar product, and it's easy to see (by writing out the components) that $\tilde {\bf a} {\bf b}=\tilde{\bf b}{\bf a}$.
The second way - it's called the Cartesian product -  gives a 3 by 3 matrix: $M_{ij}=a_ib_j$. It is not symmetric:${\bf a} \tilde {\bf b}\neq {\bf b} \tilde {\bf a}$. You can subtract the two and get
$$\left(\matrix{ 0 &a_1b_2-b_1a_2&a_1b_3-b_1a_3\\a_2b_1-b_2a_1&0&a_2b_3-b_2a_3\\a_3b_1-b_3a_1&a_3b_2-b_3a_2&0}\right)$$
This is antisymmetric, so there are only 3 independent components. You can assign these to the three elements of a vector - the (1,2) matrix element goes in the (3) vector component, the (2,3) goes in (1) and (1,3) goes in (2).  (Technically you do this using the Levi-Civita symbol).
The interesting point is that this only works in 3 dimensions. In 2D there is only one non-zero element, and in 4D there are 6. Scalar products are very very general, but vector products are peculiar to our 3D space.
So in terms of components the scalar product is the sum of the components multiplied together. The vector product is the difference of the two possible Cartesian products.
This happens to give the area of the parallelogram, but that's just a useful spin-off.


text:
\subsection{Inferring aquifer transmissivity}

Now we demonstrate the usefulness of LA-MCMC in a large-scale example: inferring the spatially heterogeneous transmissivity of an unconfined aquifer. In this problem, the likelihood contains a set of coupled partial differential equations that model transport of a nonreactive tracer through the aquifer; the tracer concentration is then observed, with noise, at selected locations in the domain. Each likelihood evaluation is thus  computationally intensive.

More specifically, the tracer is advected and diffused through the unconfined aquifer (a groundwater resource whose top boundary is not capped by an impermeable layer of rock/soil) by a steady velocity field. The aquifer is divided into ten irregularly shaped regions with a river running along the right-hand side (Figure \ref{fig:tracer-transport-domain}). A well is located at $(x_w, y_w) = (0.2, 0.75)$. Soil properties are parameterized by the transmissivity $\kappa(x,y)$, which we take to be piecewise constant in each region. Given a value of $\kappa$, we first compute the steady state hydraulic head $h$ by solving
\begin{subequations}
\begin{equation}
    - \nabla \cdot (\kappa h \nabla h) = f_h
\end{equation}
on the unit square $[0,1]^2 \ni (x,y)$, with Dirichlet boundary conditions $h(0, y) = (y + 5)/2$, $h(1, y) = (y + 3)/2$, and homogeneous Neumann boundary conditions on the $y=0$ and $y=1$ boundaries. We set the source term above to:
\begin{equation*}
    f_h(x,y) = 10 \exp{(-((x-x_w)^2+(y-y_w)^2)/0.005)},
\end{equation*}
\label{eq:hydraulic-head}
\end{subequations}
which models hydraulic forcing due to well pumping. The velocity is then
\begin{equation}
    u = -\kappa h \nabla h.
    \label{eq:tracer-velocity}
\end{equation}
The hydraulic head and velocity resulting from the ``true'' parameter values in Table \ref{tab:tracer-log-kappa} are illustrated in Figure \ref{fig:tracer-transport-hydraulic-head}. Given the steady-state velocity field $u$ and an initial tracer concentration,
\begin{equation*}
    c(x, y, 0) = \exp{(-((x-x_w)^2+(y-y_w)^2)/0.005)},
\end{equation*}
advection and diffusion of the tracer throughout the domain is modeled by a time-dependent concentration field $c(x,y,t)$ that obeys the following transport equation:
\begin{equation}
    \frac{\partial c}{\partial t} + \nabla \cdot ((d_m \mathbf{I} + d_{\ell} u u^T) \nabla c) - u^T \nabla c = -f_t
    \label{eq:tracer-concentration}
\end{equation}
with $d_m=0.01$, $d_{\ell}=0.01$, and $f_t(x,y) = f_h(x,y)$. The initial conditions $c_0(x,y)$ and tracer forcing $f_t(x,y)$ model an initial tracer release at the well, followed by continued leakage into the domain. The concentration field at time $t = 1$ is shown in Figure \ref{fig:tracer-transport-hydraulic-head}. We numerically implement this tracer-transport model using the finite element method, via the software package FEniCS \citep{FENICS}. We evolve the tracer concentration in time using an adaptive time integration scheme implemented in Sundials \citep{sundials}.

\begin{figure}[h!]
\centering
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_Domain.png}
    \caption{Spatial domain for the unconfined aquifer model. Each of the $10$ colored areas corresponds to a region of constant transmissivity. The white diamond in the top left denotes a well location, and there is a river along the right side of the domain.}
    \label{fig:tracer-transport-domain}
\end{figure}
\begin{table}[h!]
    \centering

    \begin{tabular}{c|c|c}

        Parameter & ``Truth'' & Prior mean $\bar{\theta}$

        \\ \hline \hline

        $\theta_1 = \log{(\kappa_1)}$ & $-1$ & $-1.1$ \\
        $\theta_2 = \log{(\kappa_2)}$ & $-0.25$ & $-0.1$ \\
        $\theta_3 = \log{(\kappa_3)}$ & $-0.5$ & $-0.6$ \\
        $\theta_4 = \log{(\kappa_4)}$ & $-2$ & $-1.8$ \\
        $\theta_5 = \log{(\kappa_5)}$ & $-0.5$ & $-0.4$ \\
        $\theta_6 = \log{(\kappa_6)}$ & $-3$ & $-2.5$  \\
        $\theta_7 = \log{(\kappa_7)}$ & $-1$ & $-1.2$  \\
        $\theta_8 = \log{(\kappa_8)}$ & $-1$ & $-0.9$  \\
        $\theta_9 = \log{(\kappa_9)}$ & $-1.5$ & $-1.3$  \\
        $\theta_{10} = \log{(\kappa_{10})}$ & $-1.25$ & $-1.25$
    \end{tabular}
    \caption{Parameter values used to generate the data, along with the prior mean for each of the ten log-transmissivity parameters. Each parameter corresponds to a numbered region of Figure \ref{fig:tracer-transport-domain}.}
    \label{tab:tracer-log-kappa}
\end{table}

\begin{figure}[h!]
\centering
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_HydraulicHead.pdf}
    \caption{The steady-state hydraulic head and velocity field computed by solving \eqref{eq:hydraulic-head} and \eqref{eq:tracer-velocity} given the ``true'' parameter values in Table \ref{tab:tracer-log-kappa}. The white diamond represents the location of a well in the aquifer.}
    \label{fig:tracer-transport-hydraulic-head}
\end{figure}

\begin{figure}[h!]
\centering
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_ConcentrationStep5.pdf}
    \caption{Tracer concentration computed at time $t=1$ by solving \eqref{eq:tracer-concentration} given the ``true'' parameter values in Table \ref{tab:tracer-log-kappa} and the steady state velocity field shown in Figure \ref{fig:tracer-transport-hydraulic-head}. The white diamond represents the location of a well in the aquifer, and the circles are sensor locations used in the Bayesian inference problem.}
    \label{fig:tracer-transport-concentration}
\end{figure}

Our goal is to infer the log-transmissivity parameters $\theta_i = \log{(\kappa_i)}$ for $i \in \{1,\ldots,10\}$, given observations of the tracer concentration. We endow these parameters with a Gaussian prior distribution, $\theta \sim N(\bar{\theta}, 0.05 {I})$, where $\bar{\theta}$ is the prior mean given in Table \ref{tab:tracer-log-kappa} and ${I}$ is the identity matrix of appropriate dimension. Data are collected at multiple sensor locations: at the well $(x_w, y_w) = (0.2, 0.75)$, and at a $9 \times 9$ array of points evenly spaced in $x \in [0.1, 0.9]$ and $y \in [0.1, 0.9]$. These locations are marked in Figure \ref{fig:tracer-transport-concentration}. We make observations at $10$ evenly spaced times $t \in [0, 1]$. The observations are modeled as,
\begin{equation}
    y = f(\theta) + \varepsilon,
\end{equation}
where $f(\theta)$ is the \textit{forward model} induced by the partial differential equations above, mapping the log-transmissivity parameters $\theta$ to a prediction of the tracer concentration at the chosen sensor locations/times. The sensor noise is $\varepsilon \sim N(0, 10^{-2} {I})$, yielding the conditional distribution

\begin{equation}
y \vert \theta \sim  N(f(\theta), 10^{-2} I).
\end{equation}
The posterior density then follows from Bayes' rule:
\begin{eqnarray}
    \pi(\theta \vert y) &\propto& \pi(y \vert \theta) \pi(\theta)

    \label{eq:tracer-posterior}
\end{eqnarray}
To avoid an ``inverse crime,'' we generate the data $y$ by solving the forward model with a well-refined $100 \times 100$ numerical discretization using the ``true'' parameters in Table \ref{tab:tracer-log-kappa}, and then perform the inference using a $50 \times 50$ numerical discretization of the forward model.

Evaluating the forward model at every MCMC step is computationally prohibitive. The run time of the model is approximately $5$--$15$ seconds (on a 2.6 GHz ten-core Xeon processor), depending on the values of log-transmissivity parameters $\theta$. Generating $10^5$ samples would thus take roughly $10$ days of computation time. LA-MCMC is essential to making this Bayesian inference problem computationally feasible. Here, since prior density evaluations are trivial, we apply the approximation capability of LA-MCMC to the log-likelihood function, $\log p(y \vert \theta)$.

Recall that the local polynomial surrogate \eqref{eq:local-polynomial-estimate} requires computing $q = \text{dim}(\mathcal{P})$ coefficients, and in the total-degree setting of our theory, we have $q = {{n+p}\choose{p}}$. This quantity grows rapidly with the dimension $d$ of the parameters for $p>1$.
Retaining a total-degree construction, we could set $p=0$ and thus include only a constant term in $\mathcal{P}$. This results in a constant surrogate model within each ball; now the number of terms in the local polynomial approximation is always one, independent of the parameter dimension. Similarly, setting $p=1$ lets $q$ grow only linearly with $d$. Yet in the one-dimensional example of Section~\ref{sec:1d-example}, we showed that including higher-degree terms in the approximation can reduce the number of expensive target density evaluations required to achieve a given accuracy (see Figure \ref{fig:1d-example-order}).

In the present higher-dimensional example ($d=10$), we will implement these total-degree local polynomial approximations but also compare them empirically with other choices of polynomial approximation space---in particular, truncations that retain higher-order terms more selectively. A common practice in high-dimensional  approximation is to employ \emph{sparse truncations} of the relevant index set \citep{BlatmanSudret2011}. Let $\psi_s(\theta_i)$ be a univariate polynomial of degree $s$ in $\theta_i$. (Typically, we choose an orthogonal polynomial family, but this is choice is immaterial.) Each $d$-variate polynomial basis function $\phi_{\bm{\alpha}}(\theta) \coloneqq \prod_{i=1}^{d} \psi_{\alpha_i}(\theta_i)$ is defined by a multi-index $\bm{\alpha} \in \mathbb{N}_0^d$, where the integers $\alpha_i \geq 0$ are components of $\bm{\alpha}$. Now define the $\ell^\nu$ norm for $\nu > 0$ (a quasi-norm for $0 < \nu < 1$):

\begin{equation}
    \|\bm{\alpha}\|_{\nu} = \left(\sum_{i=1}^{d} \alpha_i^{\nu}\right)^{1/\nu}.
\end{equation}
Given a maximum degree $p$, we can build our \emph{local} polynomial approximations in the polynomial space $\mathcal{P}_\nu^{p,d}$ spanned by basis functions with multi-indices $\|\bm{\alpha}\|_{\nu} \leq p$. When $\nu=1$, this recovers a total-degree expansion, which is consistent with the error bounds in Section~\ref{sec:erroranalysis}. However, the dimensionality of $\mathcal{P}_\nu^{p,d}$ decreases as $\nu \rightarrow 0$. We define the special case of $\nu=0$ to be a polynomial space with no cross terms: each multi-index has one nonzero element, and that element is at most $p$.

\begin{figure}[h!]
\centering
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_PosteriorMarginalDensities.pdf}
    \caption{One- and two-dimensional marginal distributions of the posterior distribution \eqref{eq:tracer-posterior} estimated using LA-MCMC with $\nu=1$, $\gamma_0=500$, $\gamma_1=0.5$, $\bar{\Lambda}=\infty$, $\tau_0 = 1$, $\eta = 1$, $k=75$, $p=2$, and $V(x) = \exp{(10^{-4} \|\theta-\theta_{\text{MAP}}\|)}$. Here, $\theta_{\text{MAP}} = \argmax_{\theta}{\pi(\theta \vert y)}$ is the posterior mode. In the diagonal subplots, the red and blue lines show the prior and posterior marginals, respectively. Vertical lines denote the ``true'' parameter values.}
    \label{fig:tracer-transport-marginals}
\end{figure}

We now apply the LA-MCMC algorithm to the posterior distribution of the tracer transport problem. We tested five different configurations of LA-MCMC, parameters of which are given in Table~\ref{tab:tracer-mcmc-parameters}. The Lyapunov function $V$ used in the tail correction is defined using the posterior mode, $\bar{\theta} = \argmax_{\theta}{\pi(\theta \vert y)}$. The first two rows of Table~\ref{tab:tracer-mcmc-parameters}, with $\nu =1$, correspond to using total-degree linear and quadratic local polynomial approximations (consistent with our theory). The lower three rows, with $\nu < 1$, are sparse quadratic and cubic approximations. The required number of nearest neighbors $k$ is chosen to be slightly more than the number of points $q = \text{dim}(\mathcal{P}_\nu^{p,d}) $ required to interpolate.

Reducing $\nu$ reduces $q$ and $k$. Chains produced by all these algorithmic configurations yield essentially identical posterior estimates. A trace plot of selected states from the LA-MCMC chain corresponding to the $p=1$ configuration is shown in  Figure \ref{fig:tracer-mixing-linear}; traces for LA-MCMC with higher-degree local polynomial approximations are qualitatively similar, and thus omitted here for brevity. One- and two-dimensional marginals of the posterior distribution are shown in Figure \ref{fig:tracer-transport-marginals}. While this figure is generated with $\nu=1$ and $p=2$, the results are essentially identical for the other algorithmic configurations.

\begin{table}[h!]
    \centering
    \begin{tabular}{c|c|c|c|c|c|c|c}
        degree $p$ & $\nu$ & $k$ & \# to interpolate ($q$) & $\gamma_0$ & $\gamma_1$ & $\eta$ & $\log{V(\theta)}$ \\ \hline \hline
        1 & $1$ & $15$ & $11$ & $3500$ & $1$ & $0.025$ & $10^{-2} \|\theta-\bar{\theta}\|$ \\
        2 & $1$ & $75$ & $66$ & $10^{4}$ & $1$ & $0.1$ & $10^{-2} \|\theta-\bar{\theta}\|$ \\
        2 & $0.5$ & $30$ & $21$ & $7500$ & $1$ & $0.01$ & $10^{-1} \|\theta-\bar{\theta}\|$ \\
        3 & $0.75$ & $80$ & $76$ & $7500$ & $1$ & $0.01$ & $0.05 \|\theta-\bar{\theta}\|$ \\
        3 & $0.5$ & $35$ & $31$ & $7500$ & $1$ & $0.01$ & $10^{-3} \|\theta-\bar{\theta}\|$
    \end{tabular}
    \caption{Parameter configurations for Algorithm \ref{alg:la-mcmc}, used to generate samples from the posterior distribution \eqref{eq:tracer-posterior} of the PDE/tracer transport problem. In all cases $\bar{\Lambda} = \infty$.}
    \label{tab:tracer-mcmc-parameters}
\end{table}

\begin{figure}
  \centering
  \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_MixingCubicDiagonal.pdf}
  \caption{PDE/tracer transport problem: trace plot of two states of an LA-MCMC chain using a locally cubic ($p=3$) surrogate model with $\nu = 0.5$, which excludes any cross terms. The additional parameters defined on the last row of Table \ref{tab:tracer-mcmc-parameters}.}
  \label{fig:tracer-mixing-linear}
\end{figure}

\begin{figure}
\centering
  \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_Refinements.pdf}
  \caption{PDE/tracer transport problem: number of likelihood evaluations as a function of MCMC steps, for LA-MCMC chains with parameters defined in Table \ref{tab:tracer-mcmc-parameters}. In all cases, we see a significant improvements (of three to four orders of magnitude) over the $10^6$ likelihood evaluations required by exact MCMC.}
  \label{fig:tracer-refinements}
\end{figure}

Though all these chains are successful in characterizing the posterior distribution, the required number of likelihood evaluations---and hence the overall computational cost---depend on a non-trivial relationship between the accuracy of the local approximations and the number of nearest neighbors required to define them. Figure \ref{fig:tracer-refinements} shows the number of likelihood evaluations $n$ as a function of the number of MCMC steps $t$ for each experiment in Table \ref{tab:tracer-mcmc-parameters}. We see that total-degree quadratic approximations (red line) require more likelihood evaluations for a given $t$ than linear approximations (purple line). Sparse ($\nu = 0.5$) quadratic models require fewer likelihood evaluations at a given $t$, however. And while cubic models with $\nu = 0.75$ require more likelihood evaluations than the other cases, reducing $\nu$ to 0.5 in the cubic case makes the number of  evaluations commensurate with those of the total-degree quadratic approximations. In all cases, we see significant improvements---between three and four orders of magnitude---over the $10^6$ likelihood evaluations that would be required by exact MCMC (at $t =10^6$, i.e., the right boundary of the horizontal axis).

It is important to note that our approach to controlling the bias-variance tradeoff, and thus the error thresholds for triggering refinement via \eqref{eq:local-error-bound} and \eqref{eq:error-threshold}, are derived for expansions of total degree $p$ and thus in principle applicable only for $\nu=1$. We use the same rules here for $\nu<1$, but this approach cannot guarantee that bias and variance decay at the same rate. Indeed, a full analysis of the local sparse approximations may depend on understanding the magnitudes of mixed derivatives of the target function $g$; this can become quite problem-specific, and we defer such an investigation to future work. The empirical results for $\nu<1$ here are intended to be practical and exploratory.

The choice of local approximation not only affects the overall computational cost of each chain, but also its mixing. Figure \ref{fig:tracer-autocorrelation-ess}(a) shows the decay of autocorrelation of the first variable $\theta_0$ along the chain for each of the five tested LA-MCMC configurations. Figure~\ref{fig:tracer-autocorrelation-ess}(b) shows the number of effectively independent samples (i.e., the effective sample size \citep{Wolffetal2004}) produced by each configuration, as a function of the number of MCMC steps $t$. Chains employing local linear approximations ($p=1$) tend to mix more poorly than the others, while the local cubic approximation appears to produce the best mixing on this problem.

Figure \ref{fig:tracer-autocorrelation-ess}(c), however, shows that increasing the complexity of the surrogate model does not necessarily make the algorithm more efficient \textit{as a function of the number of model evaluations} $n$. Under the ESS per model evaluation metric, we find that the most efficient scheme is a local quadratic model with no cross terms ($\nu=0.5$), followed by a total degree quadratic model. Intuitively, one would expect that higher-degree polynomials become less important as the local ball size decreases and, therefore, the extra model evaluations required to build a higher-degree surrogate yield diminishing returns. In theory, even a locally constant model will yield convergence, but in practice we find that local approximations of moderate degree (e.g., linear and quadratic) are helpful for efficiently exploring the target distribution in a reasonable number of MCMC steps.

\begin{figure}
  \centering

  \begin{tabular}{@{}c@{}}
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_Autocorrelation.pdf} \\[\abovecaptionskip]
    \small (a) Autocorrelation of $\theta_0$
  \end{tabular}

  \begin{tabular}{@{}c@{}}
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_ESS.pdf} \\[\abovecaptionskip]
    \small (b) ESS per MCMC step
  \end{tabular}

  \begin{tabular}{@{}c@{}}
    \includegraphics[width=0.475\textwidth]{Examples/Figures/Tracer/fig_ESS_refinements.pdf} \\[\abovecaptionskip]
    \small (c) ESS per model evaluation
  \end{tabular}

  \caption{PDE/tracer transport problem, computational efficiency for LA-MCMC chains with parameters defined in Table \ref{tab:tracer-mcmc-parameters}. (a) autocorrelation as a function of lag between samples; (b) effective sample size (ESS) as a function of the number of MCMC steps; (c) ESS as a function of the number of likelihood evaluations (refinements).}
  \label{fig:tracer-autocorrelation-ess}
\end{figure}