# Deletion Correcting Codes for Efficient DNA Synthesis

Johan Chrisnata,1 Han Mao Kiah,1 and Van Long Phuoc Pham1
<EMAIL_ADDRESS>
1 School of Physical and Mathematical Sciences, Nanyang Technological
University, Singapore

###### Abstract

The synthesis of DNA strands remains the most costly part of the DNA storage
system. Thus, to make DNA storage system more practical, the time and
materials used in the synthesis process have to be optimized. We consider the
most common type of synthesis process where multiple DNA strands are
synthesized in parallel from a common alternating supersequence, one
nucleotide at a time. The synthesis time or the number of synthesis cycles is
then determined by the length of this common supersequence. In this model, we
design quaternary codes that minimizes synthesis time that can correct
deletions or insertions, which are the most prevalent types of error in array-
based synthesis. We also propose polynomial-time algorithms that encode binary
strings into these codes and show that the rate is close to capacity.

##  I Introduction

DNA storage has been progressing rapidly in the last decade due to its
viability in storing information, namely its durability and high storage
density (see [1, 2] and the references therein). The process of storing
information in the form of synthetic DNA strands involves converting original
binary information into quaternary strings (nucleotide bases) and writing
those information into DNA strands using the synthesis machine. The
corresponding DNA strands are then stored in DNA pools. To read/retrieve the
stored data, the user employs a sequencing platform that creates multiple
possibly erroneous copies of the same strand which are clustered and decoded
to the original binary strings.

Although many works have been done on improving the feasibility of DNA
storage, they are still far from being practical in storing large data, mainly
because of its cost. Specifically, the DNA synthesis remains the most costly
part of the DNA storage model. In order to make the synthesis part efficient
and less costly, we need to understand how the synthesis process works.

To limit the error in the synthesis process, typically each DNA strand
contains no more than 250 nucleotides [3] (see also Table 1.1 in [2] for
experiments for DNA data storage). Thus, the user splits the converted
quaternary strings into multiple short sequences to synthesize and stores them
in a container in an unordered manner. The task of the synthesizing these
multiple strands is usually array-based, which means that they are done in
parallel [4]. To be exact, the synthesizer scans a fixed supersequence of DNA
nucleotides and appends a nucleotide at a time to the corresponding strands.
In other words, in each cycle of synthesis, the machine only appends one
nucleotide to a subset of the DNA strands that require that particular
nucleotide. As a result, all of the synthesized DNA strands must be
subsequences of the fixed supersequence. The length of the fixed supersequence
then determines the number of cycles required to generate all strands which
corresponds to the total synthesis time. Our goal now is thus to minimize the
total synthesis time which is consistent with decreasing the total cost of DNA
synthesis since each cycle requires reagents and chemicals [4].

Lenz et al.[5] showed that the supersequence that maximizes the information
rate is the alternating quaternary sequence. Therefore, for this paper, we
always use the alternating quaternary sequence as our supersequence. For
example, a naive scheme in which the binary information is encoded into DNA
strings is using the rule $(00)\xrightarrow{}A$, $(01)\xrightarrow{}C$,
$(10)\xrightarrow{}G$, $(11)\xrightarrow{}T$. In this case, an alternating
sequence which repeats the substring ACGT is scanned and the synthesis cycle
begins. This scheme achieves an information rate of 0.5 bits/cycle since four
cycles are required to synthesize one nucleotide which is equivalent to two
bits of information. They then showed that by introducing redundancy to the
synthesized strands using constrained coding, they can improve the information
rate from 0.5 bits/cycle to 0.947 bits/cycle. However, the encoding scheme
does not correct errors that happen during synthesis.

In this paper, we design codes for efficient DNA synthesis in the presence of
errors. The most prevalent types of errors that stem from array-based
synthesis are deletion and insertion [4]. We use $\mathcal{B}$ to denote an
error-function ball that maps a quaternary sequence $\mathbfsl{x}$ to the set
of its possible errors, $\mathcal{B}(\mathbfsl{x})$. To this end, we define an
$(n,T,\mathcal{B})$-synthesis code to be a quaternary code $\mathtt{C}$ of
length $n$ which is a subset of the set of all subsequences of the alternating
quaternary supersequence of length $T$, such that for any two codewords
$\mathbfsl{x}$ and $\mathbfsl{y}$ in $\mathtt{C}$,
$\mathcal{B}(x)\cap\mathcal{B}(y)=\emptyset$. Observe that there is a trade
off between the synthesis time $T$ or the length of the fixed supersequence
with the rate of the code $\mathtt{C}$. Thus, we aim to construct optimal
codes with high information rate per symbol with different synthesis time
constraints that also corrects deletion or insertion errors. For this paper,
we fix $\mathcal{B}$ to be a single insertion or deletion error ball (indel),
which means that $\mathcal{B}(\mathbfsl{x})$ is the set of all words that can
be obtained by at most a single deletion or at most a single insertion from
the word $\mathbfsl{x}$. In addition, we also design efficient encoders for
our codes.

To correct a single indel, we have the famous binary Varshamov-Tenengolts (VT)
codes [9]. This class of codes was later extended to the nonbinary case by
Tenengolts [10], who also provided a linear-time systematic encoder. However,
in the worst case, a length-$n$ codeword of a quaternary VT-code may have
synthesis time $4n$. Hence, in this paper, we fix a synthesis time $n<T\leq
4n$ and study ways to encode into quaternary VT-codewords whose synthesis time
is at most $T$.

Before we formulate our problem, we mention some works that adapt VT codes for
DNA based data storage. In [11, 12], the authors provided efficient encoders
that map messages into VT codewords that obey certain GC-content and
homopolymer constraints. [13, 14] modified VT codes to correct adjacent
transpositions, chromosome-inversions, while [15, 16] adapted VT codes to
correct deletions using more reads.

Figure 1: Plot for information rate for $W(n,\gamma n)$

##  II Preliminaries

$(n,T^{\prime})$-synthesis code$(n,T,\mathcal{B})$-synthesis code${\rm
VT}_{n}(a,b,T)$Binary messages Block encoder ${\rm ENC}{A}$ and Special
encoder ${\rm ENC}{B}$ synthesis-constrained encoder systematic single-indel
correcting encoder Direct encoder ${\rm ENC}{C}$ Figure 2: Single-indel
correcting synthesis-constrained encoder

Let $\Sigma=\\{1,2,3,4\\}$ be an alphabet of size four to represent the four
nucleotides of DNA. Throughout this paper, we use the shifted modulo operator
$(a\bmod 4)\in\\{1,2,3,4\\}$. Let $\Sigma^{n}$ denote the set of all
quaternary words of length $n$ and
$\Sigma^{n*}=\Sigma^{n-1}\cup\Sigma^{n}\cup\Sigma^{n+1}$. For a positive
integer $T$, let $\mathbf{a}_{T}=12341234\cdots$ denote the alternating
quaternary sequence of length $T$. For a quaternary word
$\mathbfsl{x}=x_{1}x_{2}\cdots x_{n}\in\Sigma^{n}$, we define the differential
word of $\mathbfsl{x}$ to be $D(\mathbfsl{x})\triangleq(x_{1}\bmod
4,x_{2}-x_{1}\bmod 4,x_{3}-x_{2}\bmod 4,\ldots,x_{n}-x_{n-1}\bmod 4)$. Note
that the map $D:\Sigma^{n}\xrightarrow{}\Sigma^{n}$ is a one to one
correspondence and hence is invertible. For a quaternary word
$\mathbfsl{x}\in\Sigma^{n}$, we define
$||\mathbfsl{x}||\triangleq\sum_{i=1}^{n}{x_{i}}$ to be the $L_{1}$-norm of
$\mathbfsl{x}$. We define the synthesis time of $\mathbfsl{x}$, denoted
$S(\mathbfsl{x})$, with respect to the alternating quaternary sequence to be
the smallest integer $T$ such that $\mathbfsl{x}$ is a subsequence of
$\mathbf{a}_{T}$. Furthermore, we can show that $S(x)=||D(\mathbfsl{x})||$.

######  Example 1.

Let $\mathbfsl{x}=23144$. Note that $\mathbfsl{x}$ is a subsequence of
$\mathbf{a}_{12}=123412341234$, and thus $S(\mathbfsl{x})=12$. On the other
hand, $D(\mathbfsl{x})=21234$, and $||D(\mathbfsl{x})||=12=S(\mathbfsl{x})$.

Let $\mathcal{C}\subseteq\Sigma^{n}$. For $\alpha\in\Sigma$, we denote
$\mathcal{C}\circ\alpha=\\{\mathbfsl{x}\circ\alpha:\mathbfsl{x}\in\mathcal{C}\\}$,
where $\mathbfsl{x}\circ\alpha$ is the word obtained by appending $\alpha$ to
$\mathbfsl{x}$. We say that $\mathcal{C}$ ia an $(n,T)$-synthesis code if for
every $\mathbfsl{x}\in\mathcal{C}$, we have $S(\mathbfsl{x})\leq T$.
Furthermore, $\mathcal{C}$ is an $(n,T,\mathcal{B})$-synthesis code if
$\mathcal{C}$ is an $(n,T)$-synthesis code and for every
$\mathbfsl{x},\mathbfsl{y}\in\mathcal{C}$, we have
$\mathcal{B}(\mathbfsl{x})\cap\mathcal{B}(\mathbfsl{y})=\emptyset$. We define
its redundancy to be $2n-\log_{2}|\mathcal{C}|$ bits. Let $W(n,T)$ denote the
largest $(n,T)$-synthesis code and let $A(n,T)$ denote its size. Likewise, we
define $A(n,T,\mathcal{B})$ to be the maximum size of an
$(n,T,\mathcal{B})$-synthesis code. For $n<T\leq 4n$, we define the
information rate of an $(n,T)$-synthesis code $\mathcal{C}$ to be
$R(\mathcal{C},T)=\frac{\log_{2}|\mathcal{C}|}{T}$ bits per cycle. Thus, for
some $1<\gamma\leq 4$, we define the corresponding capacity to be
$C(\gamma)=\limsup_{n\to\infty}{\frac{\log_{2}A(n,\gamma n)}{\gamma n}}$, and
the asymptotic rate
$C(\gamma,\mathcal{B})=\limsup_{n\to\infty}{\frac{\log_{2}A(n,\gamma
n,\mathcal{B})}{\gamma n}}$.

Although there is no explicit closed formula for $A(n,T)$ for any value of $n$
and $T$, we can find its value using recursion (see Section III-A for
details). Nevertheless, the generating fuction for $A(n,T)$ was provided by
Lenz et al. [5].

######  Proposition 1 (Lenz et al.[5]).

We have that

$\sum_{T\geq
0}A(n,T)z^{T}=\frac{1}{1-z}(z+z^{2}+z^{3}+z^{4})^{n}=\frac{z^{n}(1-z^{4})^{n}}{(1-z)^{n+1}}\,.$

The graph for the values of information rate of $W(n,T)$ for various values of
$T=\gamma n$ is given in Figure 1.

Let $\mathcal{B}$ be a single deletion or a single insertion error (indel
error), which means that $\mathcal{B}(\mathbfsl{x})$ is the set of all words
that can be obtained from $\mathbfsl{x}$ by at most one deletion or at most
one insertion. Our interest is now to estimate the value of
$A(n,T,\mathcal{B})$, namely to give upper bound and lower bound of the value.
In the next section, we present a construction of an
$(n,T,\mathcal{B})$-synthesis code for any $T$, where $n<T\leq 4n$. We also
provide encoding algorithm/encoder to construct such codes. To be more formal,
we define the following notion of encoder.

######  Definition 2.

The map ${\rm ENC}:\\{0,1\\}^{m}\xrightarrow{}\Sigma^{n}$ is a single-indel
correcting synthesis-constrained encoder if there exists a decoder map ${\rm
DEC}:\Sigma^{n*}\xrightarrow{}\\{0,1\\}^{m}$ such that

  1. (i)

For every $\mathbfsl{x}\in\\{0,1\\}^{m}$, if $\mathbf{c}={\rm
ENC}(\mathbfsl{x})$ and $\mathbfsl{y}\in\mathcal{B}(\mathbf{c})$, then
$\text{DEC}(\mathbfsl{y})=\mathbfsl{x}$,

  2. (ii)

For every $\mathbfsl{x}\in\\{0,1\\}^{m}$, $S({\rm ENC}(\mathbfsl{x}))\leq T$.

This means that the code $\mathcal{C}=\\{{\rm
ENC}(\mathbfsl{x}):\mathbfsl{x}\in\\{0,1\\}^{m}\\}$ is an
$(n,T,\mathcal{B})$-synthesis code.

In this paper, we provide several single-indel correcting synthesis-
constrained encoders that map binary strings to an
$(n,T,\mathcal{B})$-synthesis code. Specifically, for $T\geq 2.5n$, the
encoder in Algorithm 7 runs in linear time, while for general $T$, the
encoders in Algorithm 5 and Algorithm 11 run in $O(n^{3})$ and $O(n^{5})$ time
respectively. Before we describe these encoders, we describe the main
ingredient of the construction – the VT-codes.

###  II-A Varshamov-Tenengolts (${\rm VT}$) codes

######  Definition 3.

Let $\mathbfsl{x}$ be a binary word of length $n$. Then the VT-syndrome of
$\mathbfsl{x}$ is defined as ${\rm VT}(\mathbfsl{x})=\sum_{i=1}^{n}{ix_{i}}$.

######  Definition 4.

Let $\mathbfsl{x}\in\Sigma^{n}$ be a quaternary word of length $n$. We define
the auxiliary binary sequence of $\mathbfsl{x}$ to be the binary sequence
$\tilde{\mathbfsl{x}}$ of length $n-1$ such that
$\tilde{x}_{i}=\begin{cases}1&\quad\text{if }x_{i+1}\geq x_{i}\\\
0&\quad\text{if }x_{i+1}<x_{i},\\\ \end{cases}$ for all $1\leq i\leq n-1$.

######  Construction A (Tenengolts[10]).

Let $n$ be a positive integer. We define ${\rm
VT}_{n}(a,b)=\\{\mathbfsl{x}\in\Sigma^{n}:{{\rm
VT}}(\tilde{\mathbfsl{x}})=a\bmod n,$ and $\sum_{i=1}^{n}{x_{i}}=b\bmod 4\\}.$

######  Theorem 2 (Tenengolts[10]).

For any $a\in\mathbb{Z}_{n}$ and $b\in\mathbb{Z}_{4}$, the set ${\rm
VT}_{n}(a,b)$ is an $(n,4n,\mathcal{B})$-synthesis code. Moreover, there exist
$a\in\mathbb{Z}_{n}$ and $b\in\mathbb{Z}_{4}$ such that the redundancy of
${\rm VT}_{n}(a,b)$ is at most ${\left\lceil\log_{2}n\right\rceil}+2$.
Furthermore, there exists a systematic quaternary encoder that maps binary
sequence to a quaternary code that can correct single insertion or single
deletion error with redundancy at most $\log_{2}n+12$.

In the same paper, Tenengolts also provided a linear-time decoder that
corrects a single indel.

######  Theorem 3 (Tenengolts[10]).

Let $a\in\mathbb{Z}_{n}$ and $b\in\mathbb{Z}_{4}$. Given $\mathbfsl{x}\in{\rm
VT}_{n}(a,b)$, for any $\mathbfsl{y}\in\mathcal{B}(\mathbfsl{x})$, there
exists a linear-time decoder ${\rm
DEC}_{a,b}:\Sigma^{n*}\xrightarrow{}\Sigma^{n}$, such that ${\rm
DEC}_{a,b}(\mathbfsl{y})=\mathbfsl{x}$.

Next, we give lower bound on the size of an $(n,T,\mathcal{B})$-synthesis
code.

######  Definition 5.

We define the set ${\rm
VT}_{n}(a,b,T)\triangleq\\{\mathbfsl{x}\in\Sigma^{n}:S(\mathbfsl{x})\leq
T,{\rm VT}(\tilde{\mathbfsl{x}})=a\bmod n,$ and $\sum_{i=1}^{n}{x_{i}}=b\bmod
4\\}.$

######  Theorem 4.

For $n<T\leq 4n$, we have that ${\rm VT}_{n}(a,b,T)$ is an
$(n,T,\mathcal{B})$-synthesis code with size at least $A(n,T)/4n$.

###### Proof.

Consider the set of all quaternary words of length $n$, whose synthesis time
is at most $T$. We know that there are $A(n,T)$ such words. Note that ${\rm
VT}_{n}(a,b,T)$ is a subset of ${\rm VT}_{n}(a,b)$, therefore
$B(\mathbfsl{x})\cap\mathcal{B}(\mathbfsl{y})=\emptyset$ for any
$\mathbfsl{x},\mathbfsl{y}\in{\rm VT}_{n}(a,b,T)$. Thus ${\rm VT}_{n}(a,b,T)$
is an $(n,T,B)$-synthesis code. There are $4n$ choices of parameters ($4$
choices of $b$ and $n$ choices of $a$), so by the pigeon hole principle, there
exist a choice of $a$ and $b$ such that $VT_{n}(a,b,T)$ has size at least
$A(n,T)/4n$. ∎

######  Corollary 5.

For $2.5n\leq T\leq 4n$, there exists an $(n,T,\mathcal{B})$-synthesis code
with size at least $\frac{4^{n-1}}{2n}$. Hence the redundancy of the code is
at most $3+\log_{2}n$ bits.

Now, Theorem 4 and Corollary 5 are existence results. In other words, a greedy
method that constructs codes meeting the lower bound takes exponential time.
Therefore, in this work, we provide explicit construction of codes whose rates
are close to the ones promised in Theorem 4 and Corollary 5. Specifically, we
present several ways to encode binary messages to an
$(n,T,\mathcal{B})$-synthesis code. The encoders are summarized in Figure 2.

  1. (A)

Block encoder ${\rm ENC}{A}$. The first idea is to encode binary messages to
an $(n,T^{\prime})$-synthesis code, where $T=T^{\prime}+O(\log_{2}n)$. We do
this by splitting the binary messages into $\ell$ blocks and encode each block
using ranking/unranking to the set of all quaternary strings of length $k$
with synthesis time at most $T^{\prime}/\ell$. Finally we concatenate all the
blocks and encode it to an $(n,T,\mathcal{B})$-synthesis code using a
systematic single-indel correcting code, ${\rm ENC}{H}$. The whole process is
presented as a single-indel correcting synthesis-constrained encoder, ${\rm
ENC}{A}$, in Algorithm 5 that runs in $O(n^{3})$ time. This is discussed in
Subsection III-A.

  2. (B)

Special encoder ${\rm ENC}{B}$. For the special case where $T\geq 2.5n$, we
borrow ideas from [5] to provide another encoder, ${\rm ENC}{B}$, that incurs
smaller redundancy (as compared to ${\rm ENC}^{\mathbfsl{A}}$). This encoder
is given in Algorithm 7 and it runs in $O(n)$ time.

  3. (C)

Direct encoder ${\rm ENC}{C}$. Finally we provide a method that encodes binary
messages directly to a ${\rm VT}_{n}(a,b,T)$ code using ranking/unranking
algorithm for any general $T$. This encoder ${\rm ENC}{C}$ is given in
Algorithm 11 and has the best redundancy among other encoders. However, it
runs in $O(n^{5})$ time.

We compare the information rates of the three encoders in Table I.

$\gamma$, where $T=\gamma n$ | ${\rm ENC}^{\mathbfsl{A}}$ | ${\rm ENC}^{\mathbfsl{B}}$ | ${\rm ENC}^{\mathbfsl{C}}$
---|---|---|---
1.1 | 0.000 |  | 0.364
1.5 | 0.315 |  | 0.841
1.9 | 0.646 |  | 0.890
2.3 | 0.698 |  | 0.821
2.7 | 0.665 | 0.691 | 0.714
3.1 | 0.594 | 0.602 | 0.622
3.5 | 0.526 | 0.533 | 0.551
3.9 | 0.472 | 0.479 | 0.495
4.0 | 0.461 | 0.467 | 0.482
TABLE I: Information rates for different encoding methods when $n=127$.

##  III Polynomial-Time Encoders

###  III-A Block encoder ${\rm ENC}{A}$

In 1984, Tenengolts [10] presented a systematic encoder that maps quaternary
sequences to an $(n,4n,\mathcal{B})$-synthesis code with at most
${\left\lceil\log_{2}n\right\rceil}+12$ bits of redundancy. With a little
modification of the encoder from Tenengolts, we provide a systematic single
indel-correcting encoder that runs in linear time with only
${\left\lceil\log_{2}n\right\rceil}+6$ bits of redundancy. We remark that a
non-systematic linear-time encoder with ${\left\lceil\log_{2}n\right\rceil}+2$
redundant bits were given in [11]. However, we need a systematic encoder for
the purpose of this paper.

######  Definition 6.

The map ${\rm ENC}:\Sigma^{m}\xrightarrow{}\Sigma^{n}$ is a systematic single-
indel correcting encoder if there exists a decoder map ${\rm
DEC}:\Sigma^{n*}\xrightarrow{}\\{0,1\\}^{m}$ such that

  1. (i)

For every $\mathbfsl{x}\in\Sigma^{m}$, if $\mathbf{c}={\rm ENC}(\mathbfsl{x})$
and $\mathbfsl{y}\in\mathcal{B}(\mathbf{c})$, then ${\rm
DEC}(\mathbfsl{y})=\mathbfsl{x}$,

  2. (ii)

There exists an injective map
$\omega:\\{1,2,\ldots,m\\}\xrightarrow{}\\{1,2,\ldots,n\\}$, such that for
every $\mathbfsl{x}\in\Sigma^{m}$, if ${\rm ENC}(\mathbfsl{x})=\mathbf{c}$,
then $x_{i}=c_{\omega(i)}$ for all $1\leq i\leq m$.

This means that the code $\mathcal{C}=\\{{\rm
ENC}(\mathbfsl{x}):\mathbfsl{x}\in\Sigma^{m}\\}$ is a systematic single-indel
correcting code.

In order to have an encoder from binary sequences to quaternary sequences, we
use a natural bijection $\phi$ from quaternary alphabet $\Sigma=\\{1,2,3,4\\}$
to the following two-bit binary sequences:

$1\longleftrightarrow 00,\quad 2\longleftrightarrow 01,\quad
3\longleftrightarrow 10,\quad 4\longleftrightarrow 11.$

Then, we extend the mapping $\phi$ to map a quaternary sequence
$\mathbfsl{x}\in\Sigma^{n}$ to a binary sequence in $\\{0,1\\}^{2n}$. For
example, $\phi(2314)=01100011$.

First, we present our systematic single-indel correcting encoder in Algorithm
1. The input of this encoder is a quaternary sequence of length
$m=n-{\left\lceil\log_{4}n\right\rceil}-3$, where $n$ is the length of the
output quaternary sequence. Thus the redundancy of the encoder in bits is
${\left\lceil\log_{2}n\right\rceil}+6$ bits.

Algorithm 1 Systematic single-indel correcting encoder ${\rm ENC}{H}$

Input $\quad\mathbfsl{x}=x_{1}x_{2}\cdots x_{m}\in\Sigma^{m}$
Output $\mathbf{c}=c_{1}c_{2}\cdots c_{n}={\rm
ENC}{H}(\mathbfsl{x})\in(n,4n,\mathcal{B})$-synthesis code

1: Set $c_{1}\cdots c_{m}=\mathbfsl{x}$.

2: Set $c_{m+1}=c_{m+2}=c_{m}+2\bmod 4$.

3: Set $c_{m+3}=\sum_{i=1}^{m}{c_{i}}\bmod 4$.

4: Set $c_{m+4}\cdots c_{n}$ to be the quaternary representation of ${\rm
VT}(\tilde{\mathbfsl{x}})\bmod n$.

5: return $\mathbf{c}$.

For completeness, we also present our decoder map ${\rm
DEC}{H}:\Sigma^{n*}\xrightarrow{}\Sigma^{m}$ in Algorithm 2.

######  Theorem 6.

${\rm ENC}{H}$ is a systematic single-indel correcting encoder that maps
quaternary sequences of length $m=n-{\left\lceil\log_{4}n\right\rceil}-3$ to a
quaternary single-indel correcting correcting code of length $n$ with
${\left\lceil\log_{2}n\right\rceil}+6$ bits of redundancy. Furthermore, the
encoder runs in $O(n)$ time.

Algorithm 2 Systematic single-indel correcting decoder ${\rm DEC}{H}$

Input $\quad\mathbf{c}\in\Sigma^{n*}$
Output $\mathbfsl{x}\in\Sigma^{m}={\rm DEC}{H}(\mathbf{c})$

1: Let $\ell=|\mathbf{c}|$ denote the length of the input $\mathbf{c}$.

2: If $\ell=n$, then we return $\mathbfsl{x}=c_{1}\cdots c_{m}$.

3: If $\ell=n+1$, then a single insertion happened. If $c_{m+1}=c_{m+2}$, then
there is no error in the data part, hence we output $\mathbfsl{x}=c_{1}\cdots
c_{m}$. Otherwise, set $a=c_{m+4}$ and $b$ to be decimal representation of the
quaternary word $c_{m+5}\cdots c_{n+1}$. Using Theorem 3, we return
$\mathbfsl{x}={\rm DEC}_{a,b}(c_{1}\cdots c_{m+1})$

4: If $\ell=n-1$, then a single deletion happened. If $c_{m}=c_{m+1}$, then
set $a=c_{m+2}$ and $b$ to be decimal representation of the quaternary word
$c_{m+3}\cdots c_{n-1}$. Using Theorem 3, we output $\mathbfsl{x}={\rm
DEC}_{a,b}(c_{1}\cdots c_{m-1})$. Otherwise, we return
$\mathbfsl{x}=c_{1}\cdots c_{m}$.

We can verify easily that ${\rm ENC}{H}$ satisfies the definition of
systematic single-indel correcting encoder in Definition 6.

Now, we present a synthesis-constrained encoder which encodes binary messages
to an $(n,T^{\prime})$-synthesis code using ranking/unranking algorithm. Let
$\mathcal{S}$ be a finite set with cardinality $N$. A ranking function for a
set $\mathcal{S}$ is a bijection
$\text{rank}:\mathcal{S}\xrightarrow{}\\{1,2,\ldots,N\\}$. There is a unique
unranking function associated with the function rank, which is
$\text{unrank}:\\{1,2,\ldots,N\\}\xrightarrow{}\mathcal{S}$, so that
$\text{rank}(s)=i$ if and only if $\text{unrank}(i)=s$ for all
$s\in\mathcal{S}$ and $i\in\\{1,2,\ldots,N\\}$. Recall that $W(n,T)$ denote
the set of all quaternary words of length $n$ whose synthesis time is at most
$T$ and $A(n,T)$ is its size. We define the map

$\sigma^{n,k}_{T}:W(n-1,T-k)\xrightarrow{}W(n,T),$

such that

$\sigma^{n,k}_{T}((x_{1},x_{2},\ldots,x_{n-1}))=(x_{1},\ldots,x_{n-1},(x_{n-1}+k)\bmod
4),$

if $(x_{1},x_{2},\ldots,x_{n-1})\in W(n-1,T-k)$.

######  Theorem 7.

The map $\sigma^{n,k}_{T}$ is injective and

$W(n,T)=\bigsqcup_{k=1}^{4}{\sigma^{n,k}_{T}(W(n-1,T-k))},$ (1)

where $\bigsqcup$ denotes the disjoint union of sets.

###### Proof.

First, we prove injectivity. If $\mathbfsl{x},\mathbfsl{x}^{\prime}\in
W(n-1,T-k)$ where $\mathbfsl{x}\neq\mathbfsl{x}^{\prime}$ then
$\sigma^{n,k}_{T}(\mathbfsl{x})\neq\sigma^{n,k}_{T}(\mathbfsl{x}^{\prime})$
since the first $n-1$ positions of $\sigma^{n}_{T}(\mathbfsl{x})$ and
$\sigma^{n}_{T}(\mathbfsl{x}^{\prime})$ are different.

Next, we prove that
$W(n,T)=\bigsqcup_{k=1}^{4}{\sigma^{n,k}_{T}(W(n-1,T-k))}$. If
$\mathbfsl{x}\in W(n-1,T-k)$ and $\mathbfsl{x}^{\prime}\in
W(n-1,T-k^{\prime})$, where $k\neq k^{\prime}$ but $\mathbfsl{x}$ and
$\mathbfsl{x}^{\prime}$ are not necessarily distinct, then
$\mathbfsl{y}=\sigma^{n,k}_{T}(\mathbfsl{x})\neq\sigma^{n,k^{\prime}}_{T}(\mathbfsl{x}^{\prime})=\mathbfsl{y}^{\prime}$
since $k=y_{n}-y_{n-1}\bmod 4\neq y^{\prime}_{n}-y^{\prime}_{n-1}\bmod
4=k^{\prime}$. Thus, the four sets are disjoint. Furthermore, note that for
any $\mathbfsl{x}\in W(n-1,T-k)$, we have
$S(\sigma^{n,k}_{T}(\mathbfsl{x}))=||D(\sigma^{n,k}_{T}(\mathbfsl{x}))||=||D(\mathbfsl{x})||+k$.
Thus for any $\mathbfsl{x}\in W(n-1,T-k)$, we have
$\sigma^{n,k}_{T}(\mathbfsl{x})\in W(n,T)$ and hence
$\sigma^{n,k}_{T}(W(n-1,T-k))\subset W(n,T)$. Next, suppose that
$\mathbfsl{y}=(y_{1},y_{2},\ldots,y_{n})\in W(n,T)$. Let $k=y_{n}-y_{n-1}\bmod
4$. Let $\mathbfsl{x}=(y_{1},y_{2},\ldots,y_{n-1})$. Note that
$\mathbfsl{x}\in W(n-1,T-k)$, and
$\sigma^{n,k}_{T}(\mathbfsl{x})=\mathbfsl{y}$. Thus,
$W(n,T)\subset\bigsqcup_{k=1}^{4}{\sigma^{n,k}_{T}(W(n-1,T-k))}$. ∎

######  Corollary 8.

$A(n,T)=\sum_{k=1}^{4}{A(n-1,T-k)}.$

######  Proposition 9.

The base cases of the recursion in (1) are

  * •

$W(1,T)=\\{1,\ldots,T\\}$, if $1\leq T\leq 4$,

  * •

$W(1,T)=\\{1,2,3,4\\}$, if $T>4$,

  * •

$W(n,T)=\emptyset$, if $T\leq 0$ and $n\geq 1$.

Thus from (1) and Proposition 9, we can build $W(n,T)$ recursively from the
base cases.

######  Example 7.

We can construct $W(3,6)$ from $W(2,5),W(2,4),W(2,3)$ and $W(2,2)$ using the
map $\sigma^{3,k}_{6}$. Then $W(2,5)$ can in turn be constructed from
$W(1,4)=\\{1,2,3,4\\}$, $W(1,3)=\\{1,2,3\\}$, $W(1,2)=\\{1,2\\}$, and
$W(1,1)=\\{1\\}$ using the map $\sigma^{2,k}_{5}$. Similarly for
$W(2,4),W(2,3)$ and $W(2,2)$.

Thus, we obtain

$\displaystyle W(2,5)$
$\displaystyle=\sigma^{2,1}_{5}(W(1,4))\cup\sigma^{2,2}_{5}(W(1,3))$
$\displaystyle\cup\sigma^{2,3}_{5}(W(1,2))\cup\sigma^{2,4}_{5}(W(1,1))$
$\displaystyle=\\{12,23,34,41\\}\cup\\{13,24,31\\}\cup\\{14,21\\}\cup\\{11\\}.$
$\displaystyle W(2,4)$
$\displaystyle=\sigma^{2,1}_{4}(W(1,3))\cup\sigma^{2,2}_{4}(W(1,2))\cup\sigma^{2,3}_{4}(W(1,1))$
$\displaystyle=\\{12,23,34\\}\cup\\{13,24\\}\cup\\{14\\}.$ $\displaystyle
W(2,3)$ $\displaystyle=\sigma^{2,1}_{3}(W(1,2))\cup\sigma^{2,2}_{3}(W(1,1))$
$\displaystyle=\\{12,23\\}\cup\\{13\\}.$ $\displaystyle W(2,2)$
$\displaystyle=\sigma^{2,1}_{2}(W(1,1))$ $\displaystyle=\\{12\\}.$

And finally, we have

$\displaystyle W(3,6)$
$\displaystyle=\sigma^{3,1}_{6}(W(2,5))\cup\sigma^{3,2}_{6}(W(2,4))$
$\displaystyle\cup\sigma^{3,3}_{6}(W(2,3))\cup\sigma^{3,4}_{6}(W(2,2))$
$\displaystyle=\\{123,234,341,412,134,241,312,141,212,112\\}$
$\displaystyle\cup\\{124,231,342,131,242,142\\}$
$\displaystyle\cup\\{121,232,132\\}\ \cup\\{122\\}.$ (2)

In fact, they give rise to a polynomial time algorithm, via dynamic
programming, to build $W(n,T)$ and to calculate its size. In this paper, we
focus on explaining the ideas behind the design of our ranking/unranking
algorithms and not the optimization of the detailed running time analysis of
these algorithms. The idea of ranking/unranking algorithm is the unfolding of
the recursion given in (1). We can follow the natural total ordering of the
elements in $W(n,T)$ where the elements in $\sigma^{n,k}_{T}(W(n-1,T-k))$ are
before after the elements in $\sigma^{n,k^{\prime}}_{T}(W(n-1,T-k^{\prime}))$
if $k<k^{\prime}$. Our goal is to propose a polynomial time ranking/unranking
algorithm for the set $W(n,T)$. The formal unranking algorithm is presented in
Algorithm 3.

Algorithm 3 $\text{unrank}(j,n,T)$

Input Integers $j;n,T$, where $1\leq j\leq A(n,T).$
Output The $j$-th element of $W(n,T)$.

1: If $n=1$ and $1\leq j\leq T\leq 4$, then return $(j)$.

2: If $n=1$, $T>4$ and $1\leq j\leq 4$, then return $(j)$.

3: Let $k\leq 4$ be the smallest positive integer such that
$j\leq\sum_{i=1}^{k}{A(n-1,T-k)}$.

4: Let $S=\sum_{i=1}^{k-1}{A(n-1,T-k)}$ if $k>1$, and let $S=0$ if $k=1$.

5: return $\sigma^{n,k}_{T}(\text{unrank}(j-S;n-1,T-k))$

The corresponding ranking algorithm can be defined analogously in Algorithm 4.

Algorithm 4 $\text{rank}(\mathbfsl{x},n,T)$

Input Integers $n,T$ and string $\mathbfsl{x}\in W(n,T)$.
Output Integer $j$ s.t. $\mathbfsl{x}$ is the $j$-th element of $W(n,T)$.

1: If $n=1$ then return $\mathbfsl{x}$.

2: Let $k=x_{n}-x_{n-1}\bmod 4$.

3: Let $S=\sum_{i=1}^{k-1}{A(n-1,T-i)}$ if $k>1$, and let $S=0$ if $k=1$.

4: return $S+\text{rank}((x_{1},x_{2},\ldots,x_{n-1}),n-1,T-k)$.

Now, we are ready to present our single-indel correcting synthesis-constrained
encoder for general $T$. We assume that the input binary sequence
$\mathbfsl{x}$ consists of $\ell$ blocks, each of length $m$, where the $i$-th
block is denoted by $\mathbfsl{x}^{(i)}$, for $1\leq i\leq\ell$. Each block of
$\mathbfsl{x}^{(i)}$ is converted to its decimal representation $c^{(i)}$,
which is then mapped to a quaternary string $\mathbfsl{y}^{(i)}$ of length $k$
in $(k,T^{\prime}/\ell)$-synthesis code using unrank function. Note that we
must have $2^{m}\leq A(k,T^{\prime}/\ell)$, and therefore we set
$m={\left\lfloor\log_{2}A(k,T^{\prime}/\ell)\right\rfloor}$. Finally, we
concatenate all the blocks $\mathbfsl{y}^{(i)}$’s and encode it to a
quaternary string in $(n,T,\mathcal{B})$-synthesis code using ${\rm ENC}{H}$.
Here $T=T^{\prime}+\epsilon$, where
$\epsilon=4({\left\lceil\log_{4}n\right\rceil}+\ell+3)$, where the additional
$4\ell$ might come from the concatenated parts. Therefore, if the length of
the output quaternary sequence is $n$, then $\ell k=n-\log_{4}n-3$. Thus the
redundancy of the encoder in bits is $2n-\ell
m=2n-\ell{\left\lfloor\log_{2}A(\frac{n-\log_{4}n-3}{\ell},\frac{T-\epsilon}{\ell})\right\rfloor}$
bits. The formal encoder is described in Algorithm 5 and the result is
summarized in Theorem 10.

Algorithm 5 Single-indel correcting synthesis-constrained encoder ${\rm
ENC}{A}(\ell,m,\mathbfsl{x};n,T)$

Input $\mathbfsl{x}=x_{1}x_{2}\cdots x_{\ell m}\in\\{0,1\\}^{\ell m}$ such
that $2^{m}\leq A(\frac{n-\log_{4}n-3}{\ell},\frac{T-\epsilon}{\ell})$, where
$\epsilon=4({\left\lceil\log_{4}n\right\rceil}+\ell+3)$.
Output $\mathbf{c}=c_{1}c_{2}\cdots c_{n}={\rm
ENC}{A}(\ell,m,\mathbfsl{x};n,T)\in(n,T,\mathcal{B})$-synthesis code

1: Set $\mathbfsl{x}^{(t)}=x_{m(t-1)+1},x_{m(t-1)+2},\ldots,x_{m(t-1)+m}$, for
all $1\leq t\leq\ell$ which represents the $t$-th block of $\mathbfsl{x}$.

2: Let $c^{(t)}=\sum_{i=1}^{m}{\mathbfsl{x}^{(t)}_{i}2^{\ell-i}}$ be the
decimal representation of the binary string $\mathbfsl{x}^{(t)}$.

3: Let
$\mathbfsl{y}^{(t)}=\text{unrank}(c^{(t)};\frac{n-\log_{4}n-3}{\ell},\frac{T-\epsilon}{\ell})$.

4: Set
$\mathbfsl{y}=\mathbfsl{y}^{(1)}\mathbfsl{y}^{(2)}\cdots\mathbfsl{y}^{(\ell)}\in\Sigma^{n-\log_{4}n-3}$.

5: Set $\mathbf{c}={\rm ENC}{H}(\mathbfsl{y})\in\Sigma^{n}$ and return
$\mathbf{c}$.

######  Theorem 10.

The encoder ${\rm ENC}{A}$ is a single-indel correcting synthesis-constrained
encoder that maps a binary string $\mathbfsl{x}$ of length $\ell m$ to a
quaternary string of length $n$ with synthesis time at most $T$. This encoder
runs in $O(n^{3})$ time with
$2n-\ell{\left\lfloor\log_{2}A(\frac{n-\log_{4}n-3}{\ell},\frac{T-\epsilon}{\ell})\right\rfloor}$
bits of redundancy, where
$\epsilon=4({\left\lceil\log_{4}n\right\rceil}+\ell+3)$. Thus the code
$\mathcal{C}=\\{{\rm
ENC}{A}(\ell,m,\mathbfsl{x};n,T):\mathbfsl{x}\in\\{0,1\\}^{\ell m\\}}$ is an
$(n,T,\mathcal{B})$-synthesis code for a fixed $n<T\leq 4n$.

###### Proof.

We want to show that the synthesis time of the output $\mathbf{c}$ of encoder
${\rm ENC}{A}$ is indeed at most $T$. Note that the synthesis time of each
$\mathbfsl{y}^{(t)}$ is at most $\frac{T-\epsilon}{\ell}$. Thus the synthesis
time of $\mathbfsl{y}$ is at most $T-\epsilon+4\ell$, where the additional
$4\ell$ synthesis time may come from the concatenated parts between each
block. Recall that $\mathbf{c}={\rm ENC}{H}(\mathbfsl{y})$ and ${\rm ENC}{H}$
is a systematic encoder. Furthermore, the encoder appends a quaternary
substring of length ${\left\lceil\log_{4}n\right\rceil}+3$, which increases
the synthesis time by at most $4({\left\lceil\log_{4}n\right\rceil}+3)$. Thus
$S(\mathbf{c})\leq
T-\epsilon+4\ell+4({\left\lceil\log_{4}n\right\rceil}+3)=T$. ∎

The comparison of the information rates for the encoder ${\rm ENC}{A}$ for
different values of block lengths $k$ with fixed output length $n=127$ is
given in Figure 3.

For completeness, we also provide our decoder map ${\rm
DEC}{A}:\Sigma^{n*}\xrightarrow{}\Sigma^{m}$ in Algorithm 6.

Algorithm 6 Single-indel correcting synthesis-constrained decoder ${\rm
DEC}{A}(\ell,m;n,T,\mathbf{c})$

Input $\mathbf{c}\in\Sigma^{n*}.$
Output $\mathbfsl{x}=x_{1}x_{2}\cdots x_{\ell m}={\rm
DEC}{A}(\ell,m;n,T,\mathbf{c})$.

1: Set $\mathbfsl{y}={\rm DEC}{H}(\mathbf{c})\in\Sigma^{n-\log_{4}n-3}$.

2: Split $\mathbfsl{y}$ into $\ell$ blocks of equal length, such that
$\mathbfsl{y}=\mathbfsl{y}^{(1)}\mathbfsl{y}^{(2)}\cdots\mathbfsl{y}^{(\ell)}$.

3: Let
$c^{(t)}=\text{rank}(\mathbfsl{y}^{(t)},\frac{n-\log_{4}n-3}{\ell},\frac{T-\epsilon}{\ell})$,
where $\epsilon=4({\left\lceil\log_{4}n\right\rceil}+\ell+3)$, for all $1\leq
t\leq\ell$.

4: Let $\mathbfsl{x}^{(t)}\in\\{0,1\\}^{m}$ be the binary representation of
length $m$ of $c^{(t)}$.

5: Set $\mathbfsl{x}$ to be the concatenation of
$\mathbfsl{x}^{(1)}\mathbfsl{x}^{(2)}\cdots\mathbfsl{x}^{(\ell)}\in\\{0,1\\}^{\ell
m}$ and return $\mathbfsl{x}$.

Figure 3: Information rate of block encoder ${\rm ENC}{A}$ for $n=127$ with
multiple block size values

###  III-B Special encoder ${\rm ENC}{B}$ for $T\geq 2.5n$.

For the special case of $T\geq 2.5n$, we only need two bits of redundancy when
we encode the binary messages to an $(n,T^{\prime})$-synthesis code, where
$T=T^{\prime}+O(\log n)$. And then similar to the general encoder presented in
the previous subsection, we use systematic single-indel correcting encoder,
where the extra redundant bits are again encoded to minimize the synthesis
time. This special encoder ${\rm ENC}{B}$ is presented in Algorithm 7. For
simplicity, we assume that the length of the input sequence is even. The input
of this encoder is a binary sequence of length
$2m=2n-2{\left\lceil\log_{4}n\right\rceil}-10$, where $n$ is the length of the
output quaternary sequence. Thus the redundancy of the encoder ${\rm ENC}{B}$
is ${\left\lceil\log_{2}n\right\rceil}+10$ bits.

Algorithm 7 Single-indel correcting synthesis-constrained encoder ${\rm
ENC}{B}(m,\mathbfsl{x};n)$ for $T\geq 2.5n$

Input $\quad\mathbfsl{x}=x_{1}x_{2}\ldots x_{2m}\in\\{0,1\\}^{2m}$
Output $\mathbf{c}=c_{1}c_{2}\ldots c_{n}={\rm
ENC}{B}(m,\mathbfsl{x};n)\in(n,T,\mathcal{B})$-synthesis code

1: Set $\mathbfsl{y}=\phi^{-1}(\mathbfsl{x})$ to be a quaternary word of
length $m$.

2: Let $\mathbfsl{y}^{\prime}\in\Sigma^{m+1}$ such that $y^{\prime}_{1}\cdots
y^{\prime}_{m}=D(\mathbfsl{y})$ is the differential word of $\mathbfsl{y}$,
and $y^{\prime}_{m+1}=1$.

3: If $||\mathbfsl{y}^{\prime}||>2.5(m+1)$, then for all $1\leq i\leq m+1$,
replace $y_{i}$ with $5-y_{i}$.

4: Set $\mathbfsl{z}={\rm ENC}{H}(D^{-1}(\mathbfsl{y}^{\prime}))$, where
$z_{1}\cdots z_{m+1}=D^{-1}(\mathbfsl{y}^{\prime})$ is the systematic part of
the encoding and $z_{m+2}\cdots z_{n-1}$ is the suffix of $\mathbfsl{z}$.

5: Let $\mathbfsl{z}^{\prime}$ be a sequence of length $n$ where
$z^{\prime}_{1}\cdots z^{\prime}_{n-1}=D(\mathbfsl{z})$, and
$z^{\prime}_{n}=1$.

6: If $||z^{\prime}_{m+2}\cdots z^{\prime}_{n}||>2.5(n-m-1)$, then replace
$z^{\prime}_{i}$ with $5-z^{\prime}_{i}$ for all $m+2\leq i\leq n$.

7: Let $\mathbf{c}=D^{-1}(\mathbfsl{z}^{\prime})$ and return $\mathbf{c}$.

For completeness, we also provide the nontrivial decoder ${\rm DEC}{B}$ for
${\rm ENC}{B}$ in Algorithm 8.

Algorithm 8 Single-indel correcting synthesis-constrained decoder ${\rm
DEC}{B}(m;n,\mathbf{c})$ for $T\geq 2.5n$

Input $\quad\mathbf{c}\in\Sigma^{n*}$.
Output $\mathbfsl{x}\in\\{0,1\\}^{2m}$

1: Let $\ell=|\mathbf{c}|$ denote the length of the input $\mathbf{c}$.

2: If $\ell=n+1$ and $c_{m+2}\neq c_{m+3}$, then the insertion happened at the
first $m+3$ positions. Set $\mathbf{c}^{\prime}=D(\mathbf{c})$. If
$c^{\prime}_{n+1}=4$, replace $c^{\prime}_{i}$ with $5-c^{\prime}_{i}$ for all
$m+5\leq i\leq n+1$ and set
$\mathbfsl{z}=D^{-1}(c^{\prime}_{1}c^{\prime}_{2}\cdots c^{\prime}_{n})$. Set
$a=z_{m+5}$ and $b$ to be decimal representation of the quaternary word
$z_{m+6}\cdots z_{n}$. Using Theorem 3, we compute $\mathbfsl{y}={\rm
DEC}_{a,b}(z_{1}\cdots z_{m+2})$.

3: If $\ell=n+1$ and $c_{m+2}=c_{m+3}$, then there is no error in the first
$m+1$ positions. Set $\mathbfsl{y}=c_{1}c_{2}\cdots c_{m+1}$.

4: If $\ell=n-1$ and $c_{m+1}=c_{m+2}$, then the deletion happened at the
first $m+1$ positions. Set $\mathbf{c}^{\prime}=D(\mathbf{c})$. If
$c^{\prime}_{n-1}=4$, replace $c^{\prime}_{i}$ with $5-c^{\prime}_{i}$ for all
$m+1\leq i\leq n-1$ and set
$\mathbfsl{z}=D^{-1}(c^{\prime}_{1}c^{\prime}_{2}\cdots c^{\prime}_{n-2})$.
Set $a=z_{m+3}$ and $b$ to be decimal representation of the quaternary word
$z_{m+4}\cdots z_{n-2}$. We compute $\mathbfsl{y}={\rm DEC}_{a,b}(z_{1}\cdots
z_{m})$.

5: If $\ell=n-1$ and $c_{m+1}\neq c_{m+2}$, then there is no error in the
first $m+1$ positions. Set $\mathbfsl{y}=c_{1}c_{2}\cdots c_{m+1}$.

6: If $\ell=n$, then set $\mathbfsl{y}=c_{1}c_{2}\cdots c_{m+1}$.

7: Set $\mathbfsl{y}^{\prime}=D(\mathbfsl{y})$. If $y^{\prime}_{m+1}=4$,
replace $y^{\prime}_{i}$ with $5-y^{\prime}_{i}$ for all $1\leq i\leq m$. Set
$\mathbfsl{x}=\phi(D^{-1}(y^{\prime}_{1}y^{\prime}_{2}\cdots y^{\prime}_{m}))$
and return $\mathbfsl{x}$.

The result is formalized in Theorem 11.

######  Theorem 11.

The encoder ${\rm ENC}{B}$ is a single-indel correcting synthesis-constrained
encoder that maps a binary string $\mathbfsl{x}$ of length $2m$ to a
quaternary string of length $n$ with synthesis time at most $2.5n$. This
encoder runs in linear time with $10+{\left\lceil\log_{2}n\right\rceil}$ bits
of redundancy. Thus the code $\mathcal{C}=\\{{\rm
ENC}{B}(m,\mathbfsl{x};n):\mathbfsl{x}\in\\{0,1\\}^{2m}\\}$ is an
$(n,T,\mathcal{B})$-synthesis code for $T\geq 2.5n$.

###### Proof.

Let $\mathbfsl{x}\in\\{0,1\\}^{2m}$ be the input sequence. Note that
$\mathbfsl{y}=\phi^{-1}(\mathbfsl{x})$ is the quaternary representation of the
input binary string. Observe that Step 3 guarantees that the synthesis time of
$D^{-1}(\mathbfsl{y}^{\prime})$ is
$||D(D^{-1}(\mathbfsl{y}^{\prime}))||=||\mathbfsl{y}^{\prime}||\leq 2.5(m+1)$.

Furthermore, Step 6 guarantees that the suffix of $\mathbfsl{z}$ after the
systematic part of $z$ in Step 4 also has synthesis time at most $2.5(n-m-1)$.
Thus the total synthesis time for the systematic part and the suffix of
$\mathbfsl{z}$ after Step 6 is at most $2.5n$. Finally, we can verify that
${\rm ENC}{B}$ satisfies both conditions in Definition 2. ∎

###  III-C Direct encoder ${\rm ENC}{C}$

Finally, in this subsection, we provide ranking/unranking algorithm that
encodes binary messages directly into a ${\rm VT}_{n}(a,b,T)$ code for any
$n<T\leq 4n$. This encoder should give the best redundancy among other
encoders, however the drawback is that this encoder takes longer time to run
compared to the previous encoders. Similar to Subsection III-A, we need to
introduce the following recursion.

Let $W_{n}(\ell,T,a,b,\alpha)$ denote the set of all quaternary words
$\mathbfsl{x}=x_{1}x_{2}\cdots x_{\ell}$ of length $\ell$, such that
$x_{\ell}=\alpha$, $S(\mathbfsl{x})\leq T$ and ${\rm
VT}(\tilde{\mathbfsl{x}})=a\bmod n,\text{ and
}\sum_{i=1}^{\ell}{\mathbfsl{x}_{i}}=b\bmod 4$. Note that $b$ and $\alpha$ are
always in $\\{1,2,3,4\\}$, namely the shifted modulo $4$, and $a$ is always
mod $n$, where $n$ is fixed in the beginning. We have the following recursion.

$\displaystyle W_{n}(\ell,T,a,b,\alpha)$
$\displaystyle=\bigsqcup_{i=1}^{4}W_{n}(\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq
i}(\ell-1),b-\alpha,i)\circ\alpha,$ (3)

where $c=\alpha-i\bmod 4$ and $\ell\geq 2$. For example, we have
$W_{n}(\ell,T,a,b,3)=W_{n}(\ell-1,T-2,a-\ell+1,b-3,1)\circ 3\sqcup
W_{n}(\ell-1,T-1,a-\ell+1,b-3,2)\circ 3\sqcup
W_{n}(\ell-1,T-4,a-\ell+1,b-3,3)\circ 3\sqcup W_{n}(\ell-1,T-3,a,b-3,4)\circ
3$.

######  Proposition 12.

The base cases of the recursion in (3) are

  * •

$W_{n}(1,T,a,b,\alpha)=\emptyset$, if $a\neq 0$ or $\alpha\neq b$ or
$\alpha>T$,

  * •

$W_{n}(1,T,0,\alpha,\alpha)=\\{\alpha\\}$, if $\alpha\leq T$

Similar to Subsection III-A, using ranking/unranking algorithm, from (3) and
Proposition 12, we can build
$W_{n}(n,T,a,b,*)\triangleq\bigsqcup_{\alpha=1}^{4}{W_{n}(n,T,a,b,\alpha)}$
recursively from the base cases. The formal unranking algorithm is presented
in Algorithm 9, in which our goal is to find the $j$-th element of
$W_{n}(n,T,a,b,*)$.

Algorithm 9 $\text{unrank}(j;n,\ell,T,a,b,\alpha)$

Input Integers $j,n,\ell,T,a,b,\alpha$, where
$\alpha\in\mathbb{Z}_{4}\cup\\{*\\},b\in\mathbb{Z}_{4},1\leq\ell\leq
n,a\in\mathbb{Z}_{n}$ and $1\leq j\leq|W_{n}(\ell,T,a,b,\alpha)|.$
Output The $j$-th element of the set $W_{n}(\ell,T,a,b,\alpha)$.

1: If $\ell=j=1$, $b=\alpha\leq T$, and $a=0$ then return $(\alpha)$.

2: If $\ell=n$ and $\alpha=*$, then let $\beta\leq 4$ be the smallest positive
integer such that $j\leq\sum_{i=1}^{\beta}{|W_{n}(n,T,a,b,i)|}$. If $\beta>1$,
then set $S=\sum_{i=1}^{\beta-1}{|W_{n}(n,T,a,b,i)|}$ , otherwise, if
$\beta=1$ then set $S=0$. return unrank$(j-S;n,n,T,a,b,\beta)$.

3: If $\alpha\neq*$, then let $k\leq 4$ be the smallest positive integer such
that $j\leq\sum_{i=1}^{k}{|W_{n}(\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq
i}(\ell-1),b-\alpha,i)|}$, where $c=\alpha-i\bmod 4$.

4: If $k>1$, then let
$S=\sum_{i=1}^{k-1}{|W_{n}(\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq
i}(\ell-1),b-\alpha,i)|}$, where $c=\alpha-i\bmod 4$. Otherwise, if $k=1$ set
$S=0$.

5: return unrank$(j-S;n,\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq
j}(\ell-1),b-\alpha,j)\circ\alpha$, where $c=\alpha-k\bmod 4$.

The corresponding ranking algorithm can be defined analogously in Algorithm
10.

Algorithm 10 $\text{rank}(\mathbfsl{x},n,\ell,T,a,b,\alpha)$

Input $\mathbfsl{x}\in W_{n}(\ell,T,a,b,\alpha)$ and integers
$n,\ell,T,a,b,\alpha$, where
$\alpha\in\mathbb{Z}_{4}\cup\\{*\\},b\in\mathbb{Z}_{4},1\leq\ell\leq
n,a\in\mathbb{Z}_{n}$
Output Integer $j$ s.t. $\mathbfsl{x}$ is the $j$-th element of
$W_{n}(\ell,T,a,b,\alpha)$.

1: If $\ell=1$ and $\alpha=*$ then return $\alpha$, else if $\ell=1$ and
$\alpha\neq*$, then ${\bf return}$ $1$.

2: If $\alpha=*$, then let $\beta=x_{\ell}$. If $\beta>1$, then let
$S=\sum_{i=1}^{\beta-1}{|W_{n}(\ell,T,a,b,i)|}$, otherwise, if $\beta=1$, then
$S=0$. return  $S+\text{rank}(\mathbfsl{x},n,\ell,T,a,b,\beta)$

3: If $\alpha\neq*$, then let $\beta=x_{\ell-1}$. If $\beta>1$ then let
$S=\sum_{i=1}^{\beta-1}{|W_{n}(\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq
i}(\ell-1),b-\alpha,i)|}$, where $c=\alpha-i\bmod 4$, otherwise if $\beta=1$,
then let $S=0$. return  $S+\text{rank}(x_{1}x_{2}\cdots
x_{\ell-1},n,\ell-1,T-c,a-\mathbbm{1}_{\alpha\geq\beta}(\ell-1),b-\alpha,\beta)$,
where $c=\alpha-\beta\bmod 4$.

Now, we are ready to present the direct encoder ${\rm ENC}{C}$ in Algorithm
11. The input of this encoder is a binary string of length $m\leq|{\rm
VT}_{n}(a,b,T)|$ and the output is a quaternary sequence in ${\rm
VT}_{n}(a,b,T)$.

Algorithm 11 Direct encoder ${\rm ENC}{C}(m,\mathbfsl{x};n,T,a,b)$

Input $\mathbfsl{x}=x_{1}x_{2}\cdots x_{m}\in\\{0,1\\}^{m}$ such that
$2^{m}\leq|{\rm VT}_{n}(a,b,T)|$
Output $\mathbf{c}=c_{1}c_{2}\cdots c_{n}={\rm
ENC}{C}(m,\mathbfsl{x};n,T,a,b)\in\text{VT}_{n}(a,b,T)$

1: Let $j=\sum_{i=1}^{m}{x_{i}2^{m-i}}$ be the decimal representation of the
binary string $\mathbfsl{x}$.

2: return $\text{unrank}(j;n,n,T,a,b,*)$

For completeness, we also provide our decoder map ${\rm DEC}{C}$ for ${\rm
ENC}{C}$ in Algorithm 12.

Algorithm 12 Direct decoder ${\rm DEC}{C}(m;n,T,a,b,\mathbf{c})$

Input $\mathbf{c}\in\Sigma^{n*}$
Output $\mathbfsl{x}=x_{1}x_{2}\cdots x_{m}={\rm
DEC}{C}(m;n,T,a,b,\mathbf{c}))$

1: Let $\mathbfsl{y}={\rm DEC}_{a,b}(\mathbf{c})\in\Sigma^{n}$

2: Let $j=\text{rank}(\mathbfsl{y},n,n,T,a,b,*)$

3: Let $\mathbfsl{x}$ be the binary representation of length $m$ of $j$, and
return $\mathbfsl{x}$.

The result is formalized in the Theorem 13.

######  Theorem 13.

The encoder ${\rm ENC}{C}$ is a single-indel correcting synthesis-constrained
encoder that maps a binary string $\mathbfsl{x}$ of length $m$ to a quaternary
codeword in ${\rm VT}_{n}(a,b,T)$. Furthermore, ${\rm ENC}{C}$ runs in
$O(n^{5})$ time and there exist $a\in\mathbb{Z}_{n}$ and $b\in\mathbb{Z}_{4}$
such that this encoder has at most
$2n-{\left\lfloor\log_{2}A(n,T)\right\rfloor}+\log_{2}n+2$ bits of redundancy.
Thus the code $\mathcal{C}=\\{{\rm
ENC}{C}(m,\mathbfsl{x};n,T,a,b):\mathbfsl{x}\in\\{0,1\\}^{m}\\}$ is an
$(n,T,\mathcal{B})$-synthesis code for a fixed $n<T\leq 4n$.

## Acknowledgement

This research of Han Mao Kiah is supported by the Ministry of Education,
Singapore, under its MOE AcRF Tier 2 Award MOE-T2EP20121-0007.

## References

  * [1] S. Yazdi, H. M. Kiah, E. R. Garcia, J. Ma, H. Zhao, and O. Milenkovic, “DNA-based storage: Trends and methods.” IEEE Trans. Molecular, Biological, Multi-Scale Commun., vol. 1, no. 3, pp. 230-248, 2015.
  * [2] Shomorony, I., and Heckel, R. (2022). “Information-Theoretic Foundations of DNA Data Storage,” Foundations and Trends®in Communications and Information Theory, 19(1), 1-106.
  * [3] D. Carmean, L. Ceze, G. Seelig, K. Stewart, K. Strauss, and M. Willsey, “DNA data storage and hybrid molecular-electronic computing,” Proc. IEEE, vol. 107, no. 1, pp. 63–72, Jan. 2019.
  * [4] L. Ceze, J. Nivala, and K. Strauss, ”Molecular digital data storage using DNA,” Nat. Rev. Genet., vol. 20, no. 8, pp. 456–466, Aug. 2019.
  * [5] A. Lenz, Y. Liu, C. Rashtchian, P. H. Siegel, A. Wachter-Zeh and E. Yaakobi, ”Coding for Efficient DNA Synthesis,” 2020 IEEE International Symposium on Information Theory (ISIT), 2020, pp. 2885-2890, doi: 10.1109/ISIT44484.2020.9174272.
  * [6] L. Organick, S. D. Ang, Y.-J. Chen, R. Lopez, S. Yekhanin, K. Makarychev, M. Z. Racz, G. Kamath, P. Gopalan, B. Nguyen, C. N. Takahashi, S. Newman, H.-Y. Parker, C. Rashtchian, K. Stewart, G. Gupta, R. Carlson, J. Mulligan, D. Carmean, G. Seelig, L. Ceze, and K. Strauss, “Random access in large-scale DNA data storage,” Nat. Biotechnol., vol. 36, no. 3, pp. 242–248, Mar. 2018.
  * [7] S. M. H. T. Yazdi, R. Gabrys, and O. Milenkovic, “Portable and errorfree DNA-based data storage,” Sci. Rep., vol. 7, no. 1, p. 5011, Dec. 2017.
  * [8] M. H. Caruthers, “The chemical synthesis of DNA/RNA: Our gift to science,” J. Biol. Chem., vol. 288, no. 2, pp. 1420–1427, Jan. 2013.
  * [9] Levenshtein, V. I. (1966, February). Binary codes capable of correcting deletions, insertions, and reversals. In Soviet physics doklady (Vol. 10, No. 8, pp. 707-710).
  * [10] G. Tenengolts, “Nonbinary codes, correcting single deletion or insertion,” IEEE Trans. Inf. Theory, vol. IT-30, no. 5, pp. 766–769, Sep. 1984.
  * [11] K. Cai, Y. M. Chee, R. Gabrys, H. M. Kiah and T. T. Nguyen, ”Correcting a Single Indel/Edit for DNA-Based Data Storage: Linear-Time Encoders and Order-Optimality,” in IEEE Transactions on Information Theory, vol. 67, no. 6, pp. 3438-3451, June 2021, doi: 10.1109/TIT.2021.3049627.
  * [12] Nguyen, T. T., Cai, K., Immink, K. A. S., and Kiah, H. M. (2021). Capacity-approaching constrained codes with error correction for DNA-based data storage. IEEE Transactions on Information Theory, 67(8), 5602-5613.
  * [13] Gabrys, Ryan, Eitan Yaakobi, and Olgica Milenkovic. ”Codes in the Damerau distance for deletion and adjacent transposition correction.” IEEE Transactions on Information Theory 64, no. 4 (2017): 2550-2570.
  * [14] Nguyen, T. T., Cai, K., Song, W., and Immink, K. A. S. (2022, June). Optimal Single Chromosome-Inversion Correcting Codes for Data Storage in Live DNA. In 2022 IEEE International Symposium on Information Theory (ISIT) (pp. 1791-1796). IEEE.
  * [15] K. Cai, H. M. Kiah, T. T. Nguyen and E. Yaakobi, ”Coding for Sequence Reconstruction for Single Edits,” _IEEE Transactions on Information Theory_ , vol. 68, no. 1, pp. 66-79, Jan. 2022.
  * [16] Chrisnata, J., Kiah, H. M., and Yaakobi, E. (2022). Correcting deletions with multiple reads. IEEE Transactions on Information Theory, 68(11), 7141-7158.
  * [17] R. N. Grass, R. Heckel, M. Puddu, D. Paunescu, and W. J. Stark, “Robust chemical preservation of digital information on DNA in silica with errorcorrecting codes,” Angew. Chem. Int. Ed., vol. 54, no. 8, pp. 2552–2555, Feb. 2015.
  * [18] A. Lenz, P. H. Siegel, A. Wachter-Zeh, and E. Yaakobi, “Coding over sets for DNA storage,” IEEE Trans. Inf. Theory, 2019.
  * [19] S. Newman, A. P. Stephenson, M. Willsey, B. H. Nguyen, C. N. Takahashi, K. Strauss, and L. Ceze, “High density DNA data storage library via dehydration with digital microfluidic retrieval,” Nat. Commun., vol. 10, no. 1, p. 1706, Dec. 2019.
  * [20] C. Rashtchian, K. Makarychev, M. Racz, S. Ang, D. Jevdjic, S. Yekhanin, L. Ceze, and K. Strauss, “Clustering billions of reads for DNA data storage,” in Proc. Conf. Neural Inf. Process. Sys., Long Beach, CA, Dec. 2017, p. 12.
  * [21] M. Willsey, K. Strauss, L. Ceze, A. P. Stephenson, C. Takahashi, P. Vaid, B. H. Nguyen, M. Piszczek, C. Betts, S. Newman, and S. Joshi, “Puddle: A dynamic, error-correcting, full-stack microfluidics platform,” in Proc. Int. Conf. Architectural Support for Programming Languages and Operating Systems. Providence, RI, USA: ACM Press, 2019, pp. 183–197.


Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

KDH@IJCAI’23: Knowledge Discovery From Healthcare Data, August 20, 2023,
Macao, S.A.R

<EMAIL_ADDRESS>] [1]

<EMAIL_ADDRESS>] [1]<EMAIL_ADDRESS>] [1]

[1]Corresponding author. [1]These authors contributed equally.

# Anderson Acceleration For Bioinformatics-Based Machine Learning

Sarwan Ali  Georgia State University, Atlanta, USA  Prakash Chourasia  Murray
Patterson

(2023)

###### Abstract

Anderson acceleration (AA) is a well-known method for accelerating the
convergence of iterative algorithms with applications in various fields,
including deep learning and optimization. Despite its popularity in these
areas, the effectiveness of AA in classical machine learning classifiers has
not been thoroughly studied. Tabular data, in particular, presents a unique
challenge for deep learning models, and classical machine learning models are
known to perform better in these scenarios. However, the convergence analysis
of these models has received limited attention. To address this gap in
research, we implement a support vector machine (SVM) classifier variant
incorporating AA to speed up convergence. We evaluate the performance of our
SVM with and without Anderson acceleration on several datasets from the
biology domain and demonstrate that the use of AA significantly improves
convergence and reduces the training loss as the number of iterations
increases. Our findings provide a promising perspective on the potential of
Anderson acceleration in training simple machine learning classifiers and
underscore the importance of further research in this area. By showing the
effectiveness of AA in this setting, we aim to inspire more studies that
explore the applications of AA in classical machine learning.

###### keywords:

Anderson Acceleration SVM Sequence Analysis

##  1 Introduction

Anderson acceleration is a method that can be used to enhance the convergence
of gradient descent algorithms. Based on the difference between the current
and prior weight vectors, a correction term is added to the weight vector
updates at each iteration. When the gradients are changing quickly, or the
optimization landscape is very non-convex, this correction term can aid in
reducing oscillations and speeding convergence. Consider the optimization
issue as a trajectory in the weight space, where the weight vector reflects
the position at each iteration, to appreciate this concept better. Without
Anderson acceleration, the gradients at each location alone control the
trajectory of the optimization process. While Anderson acceleration can smooth
out the trajectory and minimize oscillations, the trajectory is also affected
by the difference between the current and prior weight vectors.

Solving the convex problem in finding gradient descent is a typical problem in
optimization. Newton’s methods use the inverse Hessian matrix [1] to
accelerate gradient descent, and they are successful in achieving a faster
rate of convergence compared to gradient descent or accelerated gradient
descent, but it is very expensive. By utilizing knowledge of the curvature of
the loss function landscape and quasi-Newton algorithms [2] that compute a
low-rank approximation of the Hessian, it is feasible to accelerate the
training of machine learning(ML) models. The approximate replacement matrix
for the Hessian Inverse can be found using quasi-Newton methods, described in
detail by authors in [3]. The alternative option is to use the existing data
that is already available for the fixed point approach, such as natural
gradient descent (NGD) [4] and extrapolation methods for vector sequences [5].
Approximation-based approaches such as popular techniques, including momentum
acceleration methods [6], and stochastic Newton methods [2], can be used to
calculate the approximate replacement matrix. However, as dimensionality and
complexity rise, the benefits of these methods are outweighed by the
computational expense. Due to the increased dimensionality and complexity of
the data, evaluating the training loss function is getting increasingly
expensive; consequently, fixed-point [7] approaches are more suited for
hastening the training of ML models.

The machine learning (ML) community has become interested in extrapolation
techniques like Anderson acceleration (AA) [5]. These techniques have been
demonstrated to speed up the training of sequential deep learning models (DL)
[8]. Since the optimization landscape for deep learning models can be more
complicated to improve than simpler models, Anderson acceleration is typically
more successful when used in optimization methods for deep learning [9].
Although it might not offer as big of a gain in terms of convergence speed
compared to more complex models, Anderson acceleration may still be effective
for optimizing simpler classical ML models. Ultimately, the precise
characteristics of the optimization problem being solved will determine how
effective Anderson acceleration is in each given scenario.

In this work, we propose a robust approach to perform Anderson acceleration
(AA) to speed up the training of SVM classifier models for multi-dataset
training from the domain of biological sequencing. We regularize AA by
including it in the loss optimization of simple linear classifier models
(SVMs) and classical ML training, in contrast to previous work in complex deep
learning models. We numerically demonstrate the effectiveness of the proposed
acceleration by comparing the training loss with an increasing number of
iterations on different sets of biological sequences. The results show that
using AA significantly improves convergence and efficiently accelerates the
training of traditional ML models.

##  2 Related Work

Iterative optimization methods like gradient descent and its variants are
widely used for training ML models, but convergence can be slow, especially
for high-dimensional problems. Anderson acceleration (AA) is a technique for
speeding up the convergence of these methods by exploiting the geometry of the
search space. It was first introduced by Anderson [10] as a way to accelerate
the convergence of the conjugate gradient method and has since been applied to
a variety of optimization techniques, such as Newton’s method [11], stochastic
gradient descent [12], and the Nelder-Mead simplex algorithm [13].

In recent years, there has been a growing interest in using Anderson
acceleration for training deep neural networks, where it has been applied to a
variety of tasks, such as image classification [14], natural language
processing and reinforcement learning [15]. Anderson acceleration is
particularly well-suited for deep learning problems, where it has been shown
to improve convergence and generalization performance [16, 17]. A method to
estimate a sparse generalized linear model with convex or non-convex separable
penalties using Anderson acceleration is also proposed in [17]. In these
approaches, Anderson acceleration has been shown to improve convergence and
generalization performance compared to traditional optimization methods, such
as gradient descent. In addition, it has also been applied to logistic
regression [18] and other ML models.

However, despite these advances, Anderson acceleration has not been widely
applied to classical machine learning classifiers, such as support vector
machines (SVM), despite the potential for improved convergence rates. This is
particularly relevant for tabular data, where classical machine learning
classifiers are widely used. The limited exploration of Anderson acceleration
in classical machine learning classifiers is surprising, given that the
technique is effective in improving convergence in other optimization problems
[10, 11, 12, 13].

Another area where Anderson acceleration has shown promise is in training
sparse models, such as sparse coding and dictionary learning [19]. In these
applications, Anderson acceleration effectively improves convergence and
achieves sparsity, an essential consideration in many machine-learning models.

In recent years, researchers have also explored the use of Anderson
acceleration in the training of generative adversarial networks (GANs) [20].
In these applications, Anderson acceleration has been shown to improve
convergence and stability and to produce high-quality synthesized images.

Finally, it’s worth noting that Anderson acceleration has also been applied to
the training of robust models that are robust to outlier examples and to
adversarial attacks [21]. In these applications, Anderson acceleration
effectively improves the robustness of machine learning models and defends
against adversarial attacks.

##  3 Proposed Approach

This section first discusses the algorithm we use for the proposed method.
Later, we discuss the theoretical understanding of Anderson Acceleration and
the assumptions considered.

Anderson Acceleration (AA) attempts to make greater use of previous data than
the fixed-point iteration, which only takes the most recent iteration to
produce a new estimate, $y_{k+1}=g(y_{k})$. The proposed method’s algorithmic
pseudocode is provided in Algorithm 1, and the model training flow chart is
shown in Figure 1. For model training, given the feature embedding X (or
$\phi$) made from SARS-CoV-2 sequences and its lineage (variants) as labels Y,
the first step involves the embedding generation using the methods discussed
in Section 3.2, the feature vector generated and the labels for the sequences
are then supplied to the algorithm. In the algorithm, firstly, the weight
vector is initialized with random values (1 $\times$ length of sequence). We
then initialize the $iterLoss$ and $gradient$ values (lines 2-3 in Algorithm
1) for each iteration. Afterward, for each input sample X and its label y, we
predict using weight vector $\vec{w}$ (line 5 in Algorithm 1), the predicted
value is normalized, and the gradient is updated (lines 5 and 6 in Algorithm
1). Sample loss is updated, and the iteration loss list $iterLoss$ is
maintained (lines 8 to 9 in Algorithm 1). After every sample is processed, the
gradient is averaged out, and weight history is maintained for the iteration
(lines 11 and 12 in Algorithm 1), also shown in Figure 1-e. Anderson
acceleration is used to update the weight vector from the third iteration
since we need at least two weight histories. The difference between the last
two weight histories is computed and is multiplied with Anderson factor
$\alpha$ as shown in lines 14 and 15 in Algorithm 1, also shown in Figure
Figure 1-ii. The loss and accuracy for the iteration are saved, and the next
iteration is performed to do the same steps. Finally, after all iterations,
the loss list is returned for the given input feature vectors. The loss for
each Iteration is captured and argued to be the better option for faster
convergence using Anderson Acceleration.

###  3.1 Anderson Acceleration

One way to formally prove the convergence of Anderson acceleration is to use
the concept of “linear convergence”, which refers to the rate at which the
optimization process approaches the optimal solution. Specifically, we can
show that under certain conditions, the Anderson acceleration optimization
process converges linearly, meaning that the error decreases by a constant
factor at each iteration. This contrasts standard gradient descent, which
converges at a sublinear rate (e.g., the error decreases by a factor less than
1 at each iteration).

To prove this result, we can start by considering the optimization problem in
the form of a series of updates to the weight vector, where the update at each
iteration is given by:

$w_{t+1}=w_{t}-\alpha\nabla f(w_{t})$ (1)

where $w_{t}$ is weight vector at iteration $t$, $\alpha$ is the learning
rate, and $\nabla f(w_{t})$ is the gradient of the objective function $f$ at
$w_{t}$. Now, we can add the Anderson acceleration term to the update,
resulting in:

$w_{t+1}=w_{t}+\alpha(w_{t}-w_{t-1})-\alpha\nabla f(w_{t})$ (2)

Next, we can define the error at each iteration as:

$e_{t}=w_{t}-w^{*}$ (3)

where $w^{*}$ is the optimal weight vector. Now, we can substitute the
expression for the update into the expression for the error and rearrange it
to get:

$e_{t+1}=(1-\alpha)e_{t}+\alpha e_{t-1}-\alpha\nabla f(w_{t})$ (4)

where we have used the fact that $w^{*}=w_{t}-\nabla f(w_{t})$. Now, we can
define the “damping factor” as:

$\delta=1-\alpha$ (5)

and rewrite the expression for the error as:

$e_{t+1}=\delta e_{t}+(1-\delta)e_{t-1}$ (6)

This expression has the form of a weighted average, where the weight of the
current error is given by $\delta$, and the weight of the previous error is
given by $1-\delta$. Now, we can make the following assumptions:

  1. 1.

The objective function $f$ is continuously differentiable and has a Lipschitz
continuous gradient, i.e., there exists a constant $L$ such that

$|\nabla f(x)-\nabla f(y)|\leq L|x-y|$ (7)

for all $x,y$.

  2. 2.

The objective function $f$ is bounded below, i.e., there exists a constant
$f_{\min}$ such that $f(x)\geq f_{\min}$ for all $x$.

  3. 3.

The optimization algorithm is using a fixed step size $\alpha$, and the
sequence of points ${x_{k}}$ generated by the algorithm satisfies

$|x_{k+1}-x_{k}|\leq R$ (8)

for some constant $R$ and all $k$.

These assumptions are typically made in the analysis of gradient descent
algorithms. They allow us to establish certain convergence properties of the
algorithm. Specifically, under these assumptions, it can be shown that the
sequence of points generated by gradient descent with Anderson acceleration
converges to a stationary point (a point where the gradient is zero) of the
objective function $f$ at a rate of $O(1/k)$, where $k$ is the iteration
number. This convergence rate is faster than $O(1/k^{2})$ achieved by plain
gradient descent without Anderson acceleration.

Intuitively, Anderson acceleration can be thought of as a way to incorporate
information from past iterations into the current iteration to improve the
convergence rate of the optimization algorithm. This is achieved using a
weighted combination of the current gradient and the difference between the
current and previous iterates. The weights are chosen such that the resulting
update direction better approximates the true gradient at the current iterate,
leading to faster convergence.

To compute the loss, we use “cross-entropy loss" using the following
expression:

$\text{Cross Entropy Loss}=-sum(y\times log(yPred+1e-10))$ (9)

where y is the true label, and yPred is the predicted label. The $1e-10$ is
added to the yPred to avoid the log of zero, which will cause an infinity
error. The negative sign ensures the optimization problem is formulated as a
minimization problem (hence, our loss can be negative).

The flowchart for training with Anderson Acceleration(AA) is shown in Figure
1. We provide the Feature vectors X (or $\phi$) as input along with the labels
Y. Few parameter initializations are required, such as the number of
iterations and the weight vector initialized with random values. Anderson
acceleration factor $\alpha$, for which we tried several values to study its
impact and select the best value. An empty list for loss is also shown in
Figure 1-b. For a given number of iterations, we process the samples to
compute the gradient and loss for the sample. The gradient is averaged out,
and we update the weight using Anderson Acceleration for that iteration, also
shown in Figure 1-ii. The process is repeated for the given number of
iterations.

Figure 1: Flow chart for Anderson Acceleration training. Given the input data
(X,y), we First, predict with initial weight w; this w is updated using SGD
and modified with AA in consecutive steps to compute the loss and accuracy for
the prediction.  Algorithm 1 SVM With Anderson Acceleration

1:Input: Embedded training data $X$, target labels $y$, no. of iterations
$iter$, Anderson acceleration factor $\alpha$, initial weight vector
$\vec{w}\leftarrow\text{randn}(len(X[0])$, and loss list $loss$

2:Initialization:

3:$\vec{w}\leftarrow\text{randn}(X.\text{shape}[1])$ $\triangleright$ Initial
random weight

4:$loss\leftarrow[]$ $\triangleright$ Empty list

5:$wHistory\leftarrow\vec{[w]}$ $\triangleright$ Weight vector history

6:for $i\in\\{1,2,\ldots,iter\\}$ do

7: $iterLoss\leftarrow 0$ $\triangleright$ Initial loss

8: $grad\leftarrow\textbf{0}$ $\triangleright$ Gradient

9: for $(x,y)\in(X,y)$ do

10: $yPred\leftarrow w\cdot x$ $\triangleright$ Predict the class

11: $yPred\leftarrow\frac{yPred}{\text{sum}(yPred)}$ $\triangleright$
Normalize the prediction

12: $grad\leftarrow grad+y-yPred$

13: $sampleLoss-=\text{sum}(y\times\text{log}(yPred))$

14: $iterLoss+=sampleLoss$

15: end for

16: $grad\leftarrow\frac{grad}{\text{len}(X)}$ $\triangleright$ Avg the grad
over total samples

17: $wHistory+=\vec{[w]}$ $\triangleright$ Add updated $\vec{w}$ to history

18: if $\text{len}(wHistory)>2$ then

19: $diff\leftarrow wHistory[i-1]-wHistory[i-2]$

20: $\vec{w}\leftarrow\vec{w}+\alpha\times diff+grad$ $\triangleright$ Update
$\vec{w}$ using AA

21: else

22: $\vec{w}\leftarrow\vec{w}+grad$ $\triangleright$ Use gradient to update
$\vec{w}$

23: end if

24: $loss.append([\frac{iterLoss}{\text{len}(X)}])$ $\triangleright$ Append
loss

25:end for

26:return $loss$

###  3.2 Embedding Methods

We employ the three representation learning techniques below to convert the
biological sequences into low-dimensional embeddings.

####  3.2.1 Spike2Vec [22]

This technique offers numerical embedding of the supplied input spike
sequences to facilitate the use of ML models. Initially, it produces $k$-mers
of the supplied spike sequence because $k$-mers are known to maintain the
sequence’s ordering information. For a sequence of length $N$, the total
number of $k$-mers produced is $N-k+1$. For every particular sequence,
$k$-mers is a collection of (contiguous) amino acids (also known as mers) of
length $k$. (also called nGram in the NLP domain). To convert the $k$-mers
alphabetical data into a numerical representation, the Spike2Vec computes the
frequency vector based on $k$-mers. This vector comprises the counts of each
$k$-mer in the sequence. A fixed-length feature vector is then made using the
generated $k$-mers and their frequencies in a sequence. The character alphabet
$\Sigma$ and the length of the $k$-mers are used to calculate the length of
this feature vector, which is $|\Sigma|^{k}$.

####  3.2.2 Minimizer [23]

A minimizer, also known as an $m$-mer, is the substring of consecutive letters
(amino acids) of length $m$ from a given $k$-mer that is lexicographically the
smallest in both forward and backward order of the $k$-mer, where $m<k$ and is
fixed. The repetition of $k$-mers in a long sequence, which increases
computing and storage costs, is one of the fundamental issues with $k$-mers.
Minimizers could be used to get rid of this redundant information. The
lexicographically lowest $m$-mer from both the forward and backward $k$-mers
is then used to calculate a minimizer. Two parameters ($k,$m) are provided to
the minimizer. The $k$-mer’s length is $k$, where $k=9$ in our example, and
the $m$-mer’s size is $m$, where $m=3$. It extracts $k$-mers of a sequence
given to it. Then, it calculates a corresponding $m$-mer for each $k$-mer
(minimizer). A fixed-length feature vector is then made using the generated
minimizers $m$-mers and their frequencies in a sequence. The character
alphabet $\Sigma$ and the length of the $m$-mers are used to calculate the
length of this feature vector, which is $|\Sigma|^{m}$.

####  3.2.3 Spaced k-mer [24]

The performance of sequence classification is significantly impacted by the
size and sparsity of feature vectors for sequences based on $k$-mers
frequencies. The idea of employing non-contiguous length $k$ sub-sequences
($g$-mers), proposed by spaced $k$-mers, to create compact feature vectors
with reduced sparsity and size. It first computed $g$-mers using a spike
sequence as input. We calculate $k$-mers, where $kg$, from those $g$-mers. To
conduct the trials, we used $k=4$ and $g=9$. The gap’s dimensions are
determined by $g-k$. However, this approach still involves bin scanning, which
is computationally expensive and generates very high dimensional feature
representation. We took 500 Principle components by applying PCA [25] for high
dimensional embeddings (feature vector length $>1000$).

##  4 Experimental Evaluation

To perform evaluation, we use datasets including Genome and Host. The details
are as follows:

###  4.1 Dataset Statistics

####  4.1.1 Genome Dataset

Using the well-known and widely used database of SARS-CoV-2, GISAID [26], we
retrieve the full-length nucleotide sequences of the coronavirus. Our dataset
includes the COVID-19 variant information and $8220$ nucleotide sequences. In
our sample, there are 41 different Lineages altogether. The goal is to
classify the sequences and predict the Lineage it belongs to.

####  4.1.2 Host Dataset

The National Institute of Allergy and Infectious Disease (NIAID) Virus
Pathogen Database, Investigation Resource (ViPR) [27], and GISAID was used to
retrieve the Spike protein sequences from a collection of spike sequences from
several clades of the Coronaviridae family, along with details about the hosts
that each spike sequence has infected. The hostname is used as the class label
in our classification tasks for this dataset. It displays the distribution of
the dataset across the various host types (grouped by family).

###  4.2 Evaluation Metrics

For performance evaluation of SVM without and with Anderson acceleration, we
use cross-entropy loss. The cross-entropy loss, also known as the negative
log-likelihood loss, is commonly used in supervised learning problems with
categorical targets. The cross-entropy loss for a single sample can be
expressed mathematically as follows:
$L=-\log\left(\frac{e^{f_{y_{i}}}}{\sum_{j=1}^{K}e^{f_{j}}}\right)$, where
$f_{y_{i}}$ is the predicted score for the correct class and $K$ is the number
of classes. The cross-entropy loss is averaged over the entire training set to
obtain the final objective function optimized during training.

The cross-entropy loss penalizes the predicted scores for the incorrect
classes and rewards the predicted score for the correct class. During
training, the goal is to minimize the cross-entropy loss so that the predicted
scores for the correct class are as high as possible compared to those for
incorrect classes.

##  5 Results And Discussion

In this section, we report results comparison without and with Anderson
acceleration using cross-entropy loss for different biological sequence
datasets.

###  5.1 Results For Genome Data

The results for genome data using all embedding methods are reported in Figure
2 for the best value of Anderson Acceleration (AA) factor $\alpha$. We use
cross-validation to get the best value for $\alpha$ ranging from (0, 0.1, 0.2,
$\cdots$, 1.0) for respective embeddings, where 0 implies no AA and 1.0 shows
maximum AA. For Spike2Vec embedding, we can observe that although cross-
entropy loss without Anderson acceleration is smaller with fewer iterations,
as we increase the iterations, the loss increases too. On the other hand, the
loss does not increase significantly while using Anderson acceleration in SVM.
Moreover, with AA, the loss started to converge after $300$ iterations, which
is almost half compared to the loss convergence without AA (i.e., $\approx$
600 iterations). For Minimizer-based embedding, although we can observe more
fluctuation in loss compared to Spike2Vec, the loss (and convergence) is less
when SVM is used along with AA. Similarly, the behavior of spaced $k$-mers-
based embedding differs from both Spike2Vec and Minimizer-based embedding.
Although we can see an overall increasing trend in loss with an increasing
number of iterations, the SVM with AA loss is lower than without AA when the
number of iterations increases. Overall, it is evident from all three
embedding results that the loss with AA is less than the loss without AA for
different embedding methods as we increase the number of iterations, showing
the significance of using AA for the training of SVM.

(a) Spike2Vec

(b) Minimizer

(c) Spaced Kmer

Figure 2: Loss plots with AA (green line) and without AA (red line) for Genome
data. The x-axis shows the increasing number of iterations, while the y-axis
shows the cross entropy loss. The figure is best seen in color.

###  5.2 Results For Host Data

The results for host data using all embedding methods are reported in Figure 3
for the best value of Anderson Acceleration (AA) factor $\alpha$. We use
cross-validation to get the best value for $\alpha$ ranging from (0, 0.1, 0.2,
$\cdots$, 1.0) for respective embeddings, where 0 implies no AA and 1.0 shows
maximum AA. For Spike2Vec-based embedding, the behavior is not different from
the same embedding in the case of Genome data. Although SVM without and with
Anderson acceleration converges very fast (i.e., in $<100$ iterations), the
cross entropy loss with AA is smaller than SVM without AA. We observed some
improvement in the SVM without AA in the Minimizer and Spaced $k$-mers-based
embedding methods. However, when the number of iterations is smaller, we can
observe some fluctuation in the cross-entropy loss for SVM without AA,
compared to the smooth loss curve for SVM with AA, showing its significance in
efficient training of the SVM classifier.

(a) Spike2Vec

(b) Minimizer

(c) Spaced Kmer

Figure 3: Loss plots with AA (green line) and without AA (red line) for Host
data. The x-axis shows the increasing number of iterations, while the y-axis
shows the cross entropy loss. The figure is best seen in color.

Overall, we can observe from the results of all two datasets and all three
embeddings that SVM with Anderson acceleration converges faster and shows less
cross-entropy loss than SVM without Anderson acceleration. This behavior shows
the effectiveness of SVM when used with Anderson acceleration in terms of
training time and performance when applied in a real-world setting (as we show
its performance on a $3$ real-world set of biological sequences). We believe
this study can open more opportunities for researchers to explore further the
applications of Anderson acceleration for other datasets (e.g., from other
domains such as finance and images) when applied to other machine learning
models such as perceptron. These future studies can help us further understand
and improve the predictive performance of classical ML models when applied in
a real-world (biological) setting.

##  6 Conclusion

In summary, this study provides a novel use of Anderson acceleration in the
bioinformatics area of support vector machine (SVM) classifier. Our
experiments on several sequence-based bioinformatics datasets show that
Anderson acceleration results in a considerable decrease in training loss and
improved convergence compared to the standard SVM. In the future, we will
investigate more traditional linear classifier models, such as the Perceptron,
and bigger biological data to assess their scalability and resilience.
Moreover, evaluating the robustness and generalizability of the proposed
Anderson acceleration method is also an interesting future extension.

## References

  * Goldfarb et al. [2020] D. Goldfarb, Y. Ren, A. Bahamou,  Practical quasi-newton methods for training deep neural networks,  Advances in Neural Information Processing Systems 33 (2020) 2386–2396.
  * Tong et al. [2021] Q. Tong, G. Liang, X. Cai, C. Zhu, J. Bi,  Asynchronous parallel stochastic quasi-newton methods,  Parallel computing 101 (2021) 102721.
  * Shanno [1970] D. F. Shanno,  Conditioning of quasi-newton methods for function minimization,  Mathematics of computation 24 (1970) 647–656.
  * Naganuma and Yokota [2019] H. Naganuma, R. Yokota,  A performance improvement approach for second-order optimization in large mini-batch training,  in: 2019 19th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGRID), IEEE, 2019, pp. 696–703.
  * Brezinski et al. [2022] C. Brezinski, S. Cipolla, M. Redivo-Zaglia, Y. Saad,  Shanks and anderson-type acceleration techniques for systems of nonlinear equations,  IMA Journal of Numerical Analysis 42 (2022) 3058–3093.
  * Qian [1999] N. Qian,  On the momentum term in gradient descent learning algorithms,  Neural networks 12 (1999) 145–151.
  * Walker and Ni [2011] H. F. Walker, P. Ni,  Anderson acceleration for fixed-point iterations,  SIAM Journal on Numerical Analysis 49 (2011) 1715–1735.
  * Shi et al. [2019] W. Shi, S. Song, H. Wu, Y.-C. Hsu, C. Wu, G. Huang,  Regularized anderson acceleration for off-policy deep reinforcement learning,  Advances in Neural Information Processing Systems 32 (2019).
  * Keskar et al. [2016] N. Keskar, et al.,  On large-batch training for deep learning: Generalization gap and sharp minima,  arXiv preprint arXiv:1609.04836 (2016).
  * Anderson [1965] D. G. Anderson,  Iterative procedures for nonlinear integral equations,  Journal of the ACM (JACM) 12 (1965) 547–560.
  * An et al. [2017] H. An, X. Jia, H. F. Walker,  Anderson acceleration and application to the three-temperature energy equations,  Journal of Computational Physics 347 (2017) 1–19.
  * Mai and Johansson [2020] V. Mai, M. Johansson,  Anderson acceleration of proximal gradient methods,  in: International Conference on Machine Learning, 2020, pp. 6620–6629.
  * Barton and Ivey Jr [1991] R. R. Barton, J. S. Ivey Jr, Modifications of the Nelder-Mead simplex method for stochastic simulation response optimization, Technical Report, 1991.
  * Pasini et al. [2022] M. L. Pasini, J. Yin, V. Reshniak, M. K. Stoyanov,  Anderson acceleration for distributed training of deep learning models,  in: SoutheastCon 2022, 2022, pp. 289–295.
  * Zuo et al. [2022] G. Zuo, S. Huang, J. Li, D. Gong,  Offline reinforcement learning with anderson acceleration for robotic tasks,  Applied Intelligence (2022) 1–14.
  * Pasini et al. [2021] M. L. Pasini, J. Yin, V. Reshniak, M. Stoyanov,  Stable anderson acceleration for deep learning,  arXiv preprint arXiv:2110.14813 (2021).
  * Bertrand et al. [2022] Q. Bertrand, Q. Klopfenstein, P.-A. Bannier, G. Gidel, M. Massias,  Beyond l1: Faster and better sparse models with skglm,  arXiv preprint arXiv:2204.07826 (2022).
  * Bertrand and Massias [2021] Q. Bertrand, M. Massias,  Anderson acceleration of coordinate descent,  in: International Conference on Artificial Intelligence and Statistics, 2021, pp. 1288–1296.
  * Rodriguez [2021] P. Rodriguez,  Computational assessment of the anderson and nesterov acceleration methods for large scale proximal gradient problems,  in: Symposium on Image, Signal Processing and Artificial Vision (STSIVA), 2021, pp. 1–5.
  * He et al. [2022] H. He, S. Zhao, Y. Xi, C. J. Ho, Y. Saad,  Solve minimax optimization by anderson acceleration,  in: International Conference on Learning Representations, 2022.
  * Garstka et al. [2022] M. Garstka, M. Cannon, P. Goulart,  Safeguarded anderson acceleration for parametric nonexpansive operators,  in: European Control Conference (ECC), 2022, pp. 435–440.
  * Ali and Patterson [2021] S. Ali, M. Patterson,  Spike2vec: An efficient and scalable embedding approach for covid-19 spike sequences,  in: IEEE International Conference on Big Data (Big Data), 2021, pp. 1533–1540.
  * Roberts et al. [2004] M. Roberts, W. Hayes, B. Hunt, S. Mount, J. Yorke,  Reducing storage req for biological sequence comparison,  Bioinformatics 20 (2004) 3363–3369.
  * Singh et al. [2017] R. Singh, A. Sekhon, et al.,  Gakco: a fast gapped k-mer string kernel using counting,  in: Joint ECML and Knowledge Discovery in Databases, 2017, pp. 356–373.
  * Wold et al. [1987] S. Wold, K. Esbensen, P. Geladi,  Principal component analysis,  Chemometrics and intelligent laboratory systems 2 (1987) 37–52.
  * GISAID Website [2021] GISAID Website, https://www.gisaid.org/, 2021. [Online; accessed 17-October-2022].
  * Pickett et al. [2012] B. E. Pickett, E. L. Sadat, Y. Zhang, J. M. Noronha, R. B. Squires, V. Hunt, M. Liu, S. Kumar, S. Zaremba, Z. Gu, et al.,  Vipr: an open bioinformatics database and analysis resource for virology research,  Nucleic acids research 40 (2012) D593–D598.


# Improved Speech Enhancement with the Wave-U-Net

Craig Macartney
City, University of London
<EMAIL_ADDRESS>
Tillman Weyde
City, University of London
<EMAIL_ADDRESS>
###### Abstract

We study the use of the Wave-U-Net architecture for speech enhancement, a
model introduced by Stoller et al for the separation of music vocals and
accompaniment. This end-to-end learning method for audio source separation
operates directly in the time domain, permitting the integrated modelling of
phase information and being able to take large temporal contexts into account.
Our experiments show that the proposed method improves several metrics, namely
PESQ, CSIG, CBAK, COVL and SSNR, over the state-of-the-art with respect to the
speech enhancement task on the Voice Bank corpus (VCTK) dataset. We find that
a reduced number of hidden layers is sufficient for speech enhancement in
comparison to the original system designed for singing voice separation in
music. We see this initial result as an encouraging signal to further explore
speech enhancement in the time-domain, both as an end in itself and as a pre-
processing step to speech recognition systems.

##  1 Introduction

Audio source separation refers to the problem of extracting one or more target
sources while suppressing interfering sources and noise [18]. Two related
tasks are those of speech enhancement and singing voice separation, both of
which involve extracting the human voice as a target source. The former
involves attempting to improve speech intelligibility and quality when
obscured by additive noise [7, 9, 18]; whilst the latter’s focus is on
separating music vocals from accompaniment [12].

Most audio source separation methods operate not directly in the time-domain,
but with time-frequency representations as input and output (front-end). Since
2017, the U-Net architecture on magnitude spectrograms has achieved new state
of the art results in audio source separation for music [6] and speech
dereverbration [2]. Also, neural network architectures operating in the time
domain have recently been proposed for speech enhancement [9, 11]. These
approaches have been combined in the Wave-U-Net [12] and applied to singing
voice separation. In this paper we apply the Wave-U-Net to speech enhancement
and show that it produces results that are better than the current state of
the art.

The remainder of this paper is structured as follows. In section 2, we briefly
review related work from the literature. In section 3, we introduce briefly
the Wave-U-Net architecture and its application to speech. Section 4 presents
the experiments we conducted and their results including comparison to other
methods. Section 5 concludes this article with a final summary and
perspectives for future work.

##  2 Related work

Source separation of audio has seen great improvement in recent years through
deep learning models [5, 8]. These methods, as well as more traditional ones,
mostly operate in the time-frequency domain, from deep recurrent architectures
predicting soft masks, such as [4], to convolutional encoder-decoder
architectures like that of [1]. Recently, the U-Net architecture on magnitude
spectrograms has achieved new state of the art results in audio source
separation for music [6] and speech dereverbration [2].

Also recently, models operating in the time domain have been developed. The
development of Wavenet [17] inspired other developments, including [9, 11].
The SEGAN [9] architecture was developed for the purpose of speech enhancement
and denoising. It employs a neural network in the time-domain with an encoder
and decoder pathway that successively halves and doubles the resolution of
feature maps in each layer, respectively, and features skip connections
between encoder and decoder layers. It offers state-of-the-art results on the
Voice Bank (VCTK) dataset [14].

The Wavenet for Speech Denoising [11], another architecture to operate
directly in the time domain, takes its inspiration from [17]. It has a non-
causal conditional input and a parallel output of samples for each prediction
and is based on the repeated application of dilated convolutions with
exponentially increasing dilation factors to factor in context information.

##  3 Wave-U-Net for Speech Enhancement

Figure 1: The Wave-U-Net architecture following [12].

The Wave-U-Net architecture of [12] combines elements of both of the
abovementioned architectures with the U-Net. The overall architecture is a
one-dimensional U-Net with down and upsampling blocks.

As per the spectrogram-based U-Net architectures (e.g. [6]), the Wave-U-Net
uses a series of downsampling and upsampling blocks to make its predictions,
whilst at each level of the network, the time resolution is halved.

In applying the Wave-U-Net architecture to the application of speech
enhancement, our objective is to separate a mixture waveform
$m\in[-1,1]^{L\times C}$ into K source waveforms $S^{1},...,S^{K}$ with
$S^{k}\in[-1,1]^{L\times C}$ for all $k\in{1,...,K},C$ as the number of audio
channels and $L$ as the numbers of audio samples. In our case of monaural
speech enhancement we have $K=2$ and $C=1$.

In doing so, for K sources to be estimated, a 1D convolution, zero-padded
before convolving, of filter size 1 with K · C filters is applied to convert
the stack of features at each audio sample into a source prediction for each
sample. This is followed by a tanh nonlinearity to obtain a source signal
estimate with values in the interval $(-1,1)$. All convolutions except the
final ones are followed by LeakyReLU non-linearities.

##  4 Experiments

###  4.1 Datasets

We use the same VCTK dataset [14] as the SEGAN [9], which is available
publicly, encouraging comparisons with future speech enhancement methods.

The dataset includes clean and noisy audio data at 48kHz sampling frequency.
However, like the SEGAN, we downsample to 16kHz for training and testing. The
clean data are recordings of sentences, sourced from various text passages,
uttered by 30 English-speakers, male and female, with various accents – 28
intended for training and 2 reserved for testing [16]. The noisy data were
generated by mixing the clean data with various noise datasets, as per the
instructions provided in [9, 14, 15].

With respect to the training set, 40 different noise conditions are considered
[9, 16]. These are composed of 10 types of noise (2 of which are artificially-
generated111available here: http://data.cstr.ed.ac.uk/cvbotinh/SE/data/ [14]
and 8 sourced from the DEMAND database [13], each mixed with clean speech at
one of 4 signal-to-noise ratios (SNR) (15, 10, 5, and 0 dB). In total, this
yields $11,572$ training samples, with approximately 10 different sentences in
each condition per training speaker.

Testing conditions are mismatched from those of the training. The speakers,
noise types and SNRs are all different. The separate test set with 2 speakers,
unseen during training, consists of a total of 20 different noise conditions:
5 types of noise sourced from the DEMAND database at one of 4 SNRs each (17.5,
12.5, 7.5, and 2.5 dB) [14, 15]. This yields 824 test items, with
approximately 20 different sentences in each condition per test speaker [14,
15].

###  4.2 Experimental setup

As per [12], our baseline model trains on randomly-sampled audio excerpts,
using the ADAM optimization algorithm, a learning rate of 0.0001, decay rates
$\beta 1$ = 0.9 and $\beta 2$ = 0.999 and a batch size of 16. We specify an
initial network layer size of 12, like in [12], although this is varied across
experiments, as described in the Results section below. 16 extra filters per
layer are also specified, with downsampling block filters of size 15 and
upsampling block filters of size 5 like in [12]. We train for 2,000 iterations
with mean squared error (MSE) over all source output samples in a batch as
loss and apply early stopping if there is no improvement on the validation set
for 20 epochs. We use a fixed validation set of 10 randomly selected tracks.
Then, the best model is fine-tuned with the batch size doubled and the
learning rate lowered to 0.00001, again until 20 epochs have passed without
improved validation loss.

###  4.3 Results

Table 1: Objective evaluation - comparing the mean results of the untreated noisy signal, the Wiener-, SEGAN- and Wave-U-Net-enhanced signals. Higher scores are better for all metrics. Metric | Noisy | Wiener | SEGAN | Wave-U-Net
---|---|---|---|---
PESQ | 1.97 |  2.22 |  2.16 |  2.40
CSIG | 3.35 |  3.23 |  3.48 |  3.52
CBAK | 2.44 |  2.68 |  2.94 |  3.24
COVL | 2.63 |  2.67 |  2.80 |  2.96
SSNR | 1.68 |  5.07 |  7.73 |  9.97
Table 2: Objective evaluation - mean results, comparing variations of the Wave-U-Net model with different numbers of layers, without fine-tuning applied. Metric | 12-layer | 11-layer | 10-layer | 9-layer | 8-layer
---|---|---|---|---|---
PESQ |  2.40 |  2.38 |  2.41 |  2.41 |  2.39
CSIG |  3.49 |  3.47 |  3.43 |  3.54 |  3.51
CBAK |  3.23 |  3.22 |  3.24 |  3.23 |  3.18
COVL |  2.95 |  2.92 |  2.92 |  2.97 |  2.95
SSNR |  9.79 |  9.95 |  9.98 |  9.87 |  9.30

To evaluate and compare the quality of the enhanced speech yielded by the
Wave-U-Net, we mirror the objective measures provided in [9]. Each measurement
compares the enhanced signal with the clean reference of each of the 824 test
set items. They have been calculated using the implementation provided in
[7]222available here: https://ecs.utdallas.edu/loizou/speech/software.htm. The
first metric is that of the Perceptual Evaluation of Speech Quality (PESQ) -
more specifically the wide-band version recommended in ITU-T P.862.2 (from
–0.5 to 4.5) [7, 9]. Secondly, composite measures of metrics that aim to
computationally approximate the Mean Opinion Score (MOS) that would be
produced from human perceptual trials are computed [11]. These are: CSIG, a
prediction of the signal distortion attending only to the speech signal [3]
(from 1 to 5); CBAK, a prediction of the intrusiveness of background noise [3]
(from 1 to 5); and COVL, a prediction of the overall effect [3] (from 1 to 5).
Last is the Segmental Signal-to-Noise Ratio (SSNR) [10] (from 0 to $\infty$).

Table 1 shows the results of these metrics for comparison across different
speech enhancement architectures versus the overall best-performing Wave-U-Net
model, that of the 10-layer with fine-tuning applied. As a comparative
reference, it also shows the results of these metrics when applied: directly
to the noisy signals; to signals filtered using the Wiener method, based on a
priori SNR estimation; and to the SEGAN-enhanced signal, as provided in [9].
The results indicate that the Wave-U-Net is the most effective model for
speech enhancement.

Table 2 shows the performance differences between different variations of the
Wave-U-Net, with different numbers of layers. No fine-tuning was performed to
obtain the results shown here, which explains the difference between the
10-layer Wave-U-Net in Tables 1 and 2. The results suggest that fine-tuning
does not make a meaningful difference, except on the CSIG measurement, and
that performance reaches a peak around the 10- and 9-layer models, smaller
than the best performing equivalent model for music vocals source separation
in [12]. This is likely due to the size of the receptive field, where for
speech the optimal size is probably smaller than for music.

##  5 Conclusions

###  5.1 Summary

The Wave-U-Net combines the advantages of several of the most recent
successful architectures for music and speech source separation and our
results show that it is particularly effective at speech enhancement. The
results improve over the state of the art by a good margin even without
significant adaptation or parameter tuning. This indicates that there is great
potential for this approach in speech enhancement.

###  5.2 Future work

In comparison to the SEGAN architecture, it is possible that the advantage
stems from the upsampling that avoids aliasing, which should be further
investigated. The results indicate that there is room for increasing
effectiveness and efficiency by further adapting the model size and other
parameters, e.g. filter sizes, to the task and expanding to multi-channel
audio and multi-source-separation.

## References

  * Chandna et al. [2017] Pritish Chandna, M. Miron, Jordi Janer, and Emilia Gómez.  Monoaural Audio Source Separation Using Deep Convolutional Neural Networks.  _13th International Conference on Latent Variable Analysis and Signal Separation (LVA ICA2017)_ , (2017), 2017.  URL http://mtg.upf.edu/node/3680.
  * Ernst et al. [2018] Ori Ernst, Shlomo E. Chazan, Sharon Gannot, and Jacob Goldberger.  Speech dereverberation using fully convolutional networks.  _CoRR_ , abs/1803.08243, 2018.  URL http://arxiv.org/abs/1803.08243.
  * Hu et al. [2008] Yi Hu, Philipos C Loizou, and Senior Member.  Evaluation of Objective Quality Measures for Speech Enhancement.  _IEEE transactions on audio, speech and language processing_ , 16(1):229, 2008.  doi: 10.1109/TASL.2007.911054.  URL http://www.utdallas.edu/loizou/speech/noizeus/.
  * [4] Po-Sen Huang, Minje Kim, Mark Hasegawa-Johnson, and Paris Smaragdis.  Singing-Voice Separation from Monaural Recordings Using Deep Recurrent Neural Networks.  Technical report.  URL http://www.isle.illinois.edu/sst/pubs/2014/huang-ismir2014.pdf.
  * Huang et al. [2015] Po-Sen Huang, Minje Kim, Mark Hasegawa-Johnson, and Paris Smaragdis.  Joint optimization of masks and deep recurrent neural networks for monaural source separation.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 23(12):2136–2147, 2015.
  * Jansson et al. [2017] Andreas Jansson, Eric J. Humphrey, Nicola Montecchio, Rachel M. Bittner, Aparna Kumar, and Tillman Weyde.  Singing voice separation with deep u-net convolutional networks.  In _Proceedings of the 18th International Society for Music Information Retrieval Conference, ISMIR 2017, Suzhou, China, October 23-27, 2017_ , pages 745–751, 2017.  URL https://ismir2017.smcnus.org/wp-content/uploads/2017/10/171_Paper.pdf.
  * Loizou [2013] Philipos C Loizou.  _Speech Enhancement: Theory and Practice_.  CRC Press, Inc., Boca Raton, FL, USA, 2nd edition, 2013.  ISBN 1466504218, 9781466504219.
  * Nugraha et al. [2016] Aditya Arie Nugraha, Antoine Liutkus, and Emmanuel Vincent.  Multichannel audio source separation with deep neural networks.  _IEEE/ACM Trans. Audio, Speech & Language Processing_, 24(9):1652–1664, 2016.
  * Pascual et al. [2017] Santiago Pascual, Antonio Bonafonte, and Joan Serrà.  Segan: Speech enhancement generative adversarial network.  In _Proc. Interspeech 2017_ , pages 3642–3646, 2017.  doi: 10.21437/Interspeech.2017-1428.  URL http://dx.doi.org/10.21437/Interspeech.2017-1428.
  * Quackenbush et al. [1988] Schuyler R. Quackenbush, T. P. (Thomas Pinkney) Barnwell, and Mark A. Clements.  _Objective measures of speech quality_.  Prentice Hall, 1988.  ISBN 0136290566.
  * [11] Dario Rethage, Jordi Pons, and Xavier Serra.  A Wavenet for Speech Denoising.  Technical report.  URL https://arxiv.org/pdf/1706.07162.pdf.
  * Stoller et al. [2018] Daniel Stoller, Sebastian Ewert, and Simon Dixon.  Wave-U-Net: A Multi-Scale Neural Network for End-to-End Audio Source Separation.  2018\.  URL https://arxiv.org/abs/1806.03185.
  * Thiemann et al. [2013] Joachim Thiemann, Nobutaka Ito, and Emmanuel Vincent.  The Diverse Environments Multi-channel Acoustic Noise Database (DEMAND): A database of multichannel environmental noise recordings The Diverse Environments Multi-channel Acoustic Noise Database (DEMAND): A database of multichannel environmental noise record-ings.  page 10, 2013.  doi: 10.5281/zen.  URL https://hal.inria.fr/hal-00796707.
  * Valentini-Botinhao [2017] Cassia Valentini-Botinhao.  Noisy speech database for training speech enhancement algorithms and TTS models, 2016 [sound].  _University of Edinburgh. School of Informatics. Centre for Speech Technology Research (CSTR)_ , 2017.  URL http://dx.doi.org/10.7488/ds/2117.
  * Valentini-Botinhao et al. [2016a] Cassia Valentini-Botinhao, Xin Wang, Shinji Takaki, and Junichi Yamagishi.  Investigating rnn-based speech enhancement methods for noise-robust text-to-speech.  In _9th ISCA Speech Synthesis Workshop_ , pages 146–152, 2016a.  doi: 10.21437/SSW.2016-24.  URL http://dx.doi.org/10.21437/SSW.2016-24.
  * Valentini-Botinhao et al. [2016b] Cassia Valentini-Botinhao, Xin Wang, Shinji Takaki, and Junichi Yamagishi.  Speech Enhancement for a Noise-Robust Text-to-Speech Synthesis System using Deep Recurrent Neural Networks.  2016b.  doi: 10.21437/Interspeech.2016-159.  URL http://dx.doi.org/10.21437/Interspeech.2016-159.
  * van den Oord et al. [2016] Aäron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew W. Senior, and Koray Kavukcuoglu.  Wavenet: A generative model for raw audio.  In _The 9th ISCA Speech Synthesis Workshop, Sunnyvale, CA, USA, 13-15 September 2016_ , page 125, 2016.  URL http://www.isca-speech.org/archive/SSW_2016/abstracts/ssw9_DS-4_van_den_Oord.html.
  * Vincent et al. [2018] Emmanuel Vincent, Tuomas Virtanen, and Sharon Gannot, editors.  _Audio Source Separation and Speech Enhancement_.  John Wiley & Sons Ltd, Chichester, UK, 9 2018.  ISBN 9781119279860.  doi: 10.1002/9781119279860.  URL http://doi.wiley.com/10.1002/9781119279860.