\section{Introduction}

With the fast advances in next-generation sequencing (NGS) techniques, the monetary cost for DNA sequencing has been declining at a rate that is outpacing Moore's law~\cite{dnacost}.
However, on the other side of the coin, this rapid throughput in sequencing means that data processing has become a more severe bottleneck.
For example, performing read mapping of the human genome on an Intel Xeon processor now takes more than 20$\times$ the sequencing time~\cite{goyal2017ultra}.

\emph{Read mapping}, a process in sequence alignment, maps a piece of query DNA (generated from sequencing) to matching locations of a reference DNA.
However, the reference and query DNAs do not exactly match due to sequencing errors and/or mutations.
Therefore, the problem falls into the category of approximate string matching.

One popular strategy to solve this problem is \emph{seed-and-extend}, where the seeding phase locates a few exact matches, and the extension phase performs approximate matching based on dynamic programming (DP) such as the Smith--Waterman algorithm~\cite{sw}.
Unfortunately, the Smith--Waterman algorithm has quadratic complexity with the input length, which continues to increase as NGS techniques evolve.

In this work, we focus on the extension, which consumes a significant portion of the execution time for read mapping~\cite{primer}.
Many proposals have been made for dedicated hardware accelerators with ASICs~\cite{fujiki2018genax,turakhia2018darwin,fujiki2020seedex} or FPGAs\cite{ham2020genesis, houtgast2015fpga,zhang2007implementation,ahmed2015heterogeneous,blastp} to cope with the computational bottleneck problem.
However, such accelerators are yet to be widely used in practice.
The ASICs are expensive to produce and would have difficulty appearing in the market unless mass production was guaranteed.
The FPGAs are less prone to this issue because many FPGA accelerator cards are already available on the market.
However, FPGA accelerator cards are not yet the dominant option because they are expensive and hard to program.

Under such circumstances, GPU-based accelerations can be viable options because they are cheaper, easier to find on the market, and require less effort to develop software.
Therefore, many GPU-based libraries have been developed for bioinformatics~\cite{cushaw1,gasal2,nvbio}.
With the fast growth in the computational capabilities of GPUs~\cite{ampere}, GPUs will likely continue to be the primary option for accelerating sequence alignments.

However, by analyzing the state-of-the-art GPU-accelerated seed extension, we identified several missing opportunities for performance improvement that have become critical, especially with long string queries.
Considering that the lengths of sequence reads are rapidly growing with third-generation sequencers~\cite{third}, the performance gap between the ideal case and the existing software will widen.
Specifically, we found that existing GPU kernels 1) inefficiently utilize memory and 2) suffer from load imbalance.

In this paper, we present SALoBa\xspace (\underline{\textbf{S}}equence \underline{\textbf{A}}lignment with Data \underline{\textbf{Lo}}cality and Workload  \underline{\textbf{Ba}}lance), which addresses the mentioned problems to achieve superior performance.
Taking lessons from several ASIC/FPGA accelerators~\cite{ahmed2015heterogeneous,zhang2007implementation, fujiki2020seedex}, SALoBa\xspace puts together a kernel that makes better use of the computational capability of modern GPUs.

SALoBa\xspace utilizes intra-query parallelism by allocating multiple CUDA threads to a query-reference pair.
Even though this option has been studied by some prior art~\cite{korpar2013sw,adept},
these suffer from inefficient memory access patterns and resource underutilization, being outperformed by the current state-of-the-art, which only uses inter-query parallelism.
SALoBa\xspace addresses those issues by introducing two novel techniques: lazy spilling to global memory and subwarp scheduling.
In \emph{lazy spilling}, the data are first accumulated to the CUDA shared memory and later spilled to the global memory in a coalesced manner.
By utilizing the double-buffered shared memory region in a rotating manner, the redundancy in global memory is greatly reduced without much overhead to the shared memory.
In addition, \emph{subwarp scheduling} divides a warp into multiple subwarps to mitigate the underutilization problem.
While this slightly increases the workload imbalance, the gain from better resource utilization often dominates the overhead from workload imbalance.

According to the experimental results, SALoBa\xspace performs significantly faster than the state-of-the-art GPU aligner in the seed extension kernel over several sequence lengths.
When tested on real-world datasets with various lengths, we obtained superior results due to better workload balancing.
Our contributions can be summarized as follows:
\begin{itemize}
    \item We identify several unexploited opportunities from the current libraries for performance improvements.
    \item We propose SALoBa\xspace, which provides significant speedup compared to the current state-of-the-art GPU-based alignment libraries.
    \item We perform an extensive amount of evaluation, including synthetic and real-world data from a popular read alignment software to demonstrate the efficiency of SALoBa\xspace.
\end{itemize}

\section{Background}

\subsection{Seed Extension}
The heart of the sequence alignment problem is approximate string matching. Because there could be multiple possible paths to take every time a mismatch occurs, its complexity quickly increases along with the length of the input pair.

The seed-and-extend strategy is an approach to perform alignment efficiently. Based on the observation that a good alignment usually contains many matches, the strategy first finds multiple exact matches, called \emph{seeds}.
Based on these, approximate matching scores are calculated by \emph{extending} to both directions from the found seeds.

\begin{figure}
    \centering
    \includegraphics[width=0.95\columnwidth]{figs/sw.pdf}
    \caption{\vspace{-4mm}An example seed extension algorithm.\vspace{-1mm}}
    \label{fig:sw}
\end{figure}

Seed extension is often performed using a DP algorithm that fills a two-dimensional DP table with the complexity of $O(N^2)$.
Popular algorithms are Smith--Waterman~\cite{sw} and Needleman--Wunsch~\cite{nw}.
The score of each cell is calculated based on predetermined scores for the type of differences (i.e., insert, delete, and mismatch).
With some adjustments for considerations for frequent long gaps, the affine gap function calculates the score of a DP cell $H(i, j)$ as follows:
\begin{align}
    H(i,j) &= max\begin{cases}0\\
    E(i,j)\\
    F(i,j)\\
    H(i-1,j-1) + S(i,j)
    \end{cases},\\
    E(i,j) &= max\begin{cases}    H(i,j-1) - \alpha \\
    E(i,j-1)-\beta \\
    \end{cases},\\
    F(i,j) &= max\begin{cases}    H(i-1,j) - \alpha \\
    F(i-1,j)-\beta \\
    \end{cases},
    \label{eq:sm}
\end{align}
where $S(i, j)$ is a score function that returns a positive value when the reference at $i$ and query at $j$ match and returns a negative value otherwise.
In addition, $E$ and $F$  are auxiliary variables that keep track of the continuing gaps.
Last, $\alpha$ and $\beta$ account for different gap penalties according to new gaps or continued gaps, respectively.

\fig{sw} presents an example DP table, where $Q$ represents a query string, and $R$ represents a reference string.
To calculate a cell, the cells from three adjacent cells, the top, left, and top left, are needed, represented by green arrows.
The trace of the highest score in the DP table represents the best match found by the algorithm, denoted with red arrows and numbers.

\subsection{Baseline GPU-based Seed Extension}
\label{sec:gasal}

There have been several attempts to accelerate seed extension using GPUs.
\JL{change this to intra-query and inter-query parallelism. Also talk about 2bit issues.}
In sequence read data, each base is represented by a character data type of eight bits.
However, only five bases exist within the sequences: A, C, G, T/U, and N, where T is used for DNA and U is used for RNA.
The N denotes unknown bases.
Having five bases indicates that at least three bits are needed to represent each.
Because three-bit representations are inefficient to deal with in modern architectures, a four-bit packed representation is often used~\cite{gasal2,nvbio} and some work utilizes eight-bit representation~\cite{korpar2013sw}.

Existing methods for GPU-based seed extension can be roughly categorized by the parallelism they utilize: \emph{intra-query parallelism} and \emph{inter-query parallelism}.
Intra-query parallelism refers to processing multiple cells in a DP table concurrently.
As shown in \fig{sw}, parallelism exists in an anti-diagonal form in the DP table. For example, after the cell (0,0) has been processed, cells (0,1) and (1,0) can be processed in parallel.
Afterwards, cells (0,2), (1,1), and (2,0) can be processed in parallel.
Many ASIC/FPGA accelerators successfully take advantage of this, commonly using one-dimensional systolic array structures~\cite{zhang2007implementation,fujiki2020seedex,ahmed2015heterogeneous}.
Some GPU-based approaches utilize this type of parallelism~\cite{korpar2013sw,adept}.
However, they often fail to achieve sufficient speedup due to resource under-utilization and inefficient memory access~\cite{adept}.
On the other hand, inter-query parallelism refers to simply processing multiple query-reference pairs in parallel.
Because the latter is easier to optimize resource utilization, many successful libraries adopt this approach~\cite{gasal2,nvbio,cushaw2gpu,soap3}.

Among the libraries using inter-query parallelism, GASAL2~\cite{gasal2} is known to show the current state-of-the-art performance.
In this section, we describe the details of its strategy.
The GPU registers are 32-bit wide; thus, eight bases from each query and reference are fetched in a single step.
Therefore, it is natural to process 8\texttimes8 cells at once that correspond to a single-word reference and a single-word query.
After processing the block of 8\texttimes8 cells, the thread advances to the right.
For this, a column of eight rightmost cells must be stored for dependency. The dependency is fulfilled by keeping the values of these DP cells within the registers.
When a thread reaches the rightmost end of the table, it moves to the bottom part of the DP table.
Then, the cells from the top must be accessed for dependency, and the thread stores the bottom eight cells of each 8\texttimes8 block in the global memory.
Thus, considering the matrix size of $N\times N$, the amount of global memory access becomes $N\times N/4$ for reading and writing dependent cell values.

\input{eval/histo250}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{figs/parallel.pdf}
    \caption{Intra-query parallelism. (a) shows parallelization strategy using a warp, (b) illustrates the terms used in this figure. (c) demonstrates how the dependency between cells are handled and (d) depicts the use of anti-diagonal structure.}
    \label{fig:parallel}
\end{figure*}

\section{Motivation - Diagnosis}
\label{sec:moti}
In this section, we diagnose the state-of-the-art aligner GASAL2~\cite{gasal2} and discuss potential opportunities for further performance improvement.
First, we reveal a significant load imbalance between individual queries.
Second, we identify that the use of global memory is not coalesced, which requires a significant amount of redundant memory access.

\subsection{Load Imbalance}
As discussed in Section~\ref{sec:gasal}, the main parallelization scheme of GASAL2 is inter-query parallelism.
In the strategy of letting each CUDA thread handle a single query, the main advantage is that there is no complicated inter-thread synchronization or costly communication.

However, one drawback of this approach is that it ignores the variance in the workload for the seed extension.
The seeding step provides the query and reference sizes as input to the seed extension.
Because of this, the lengths vary by individual inputs, and a significant imbalance occurs.

\fig{histo250} plots the distribution of the query and reference sequences from two types of reads (see \sect{eval} for detailed settings) from the popular alignment software BWA-MEM~\cite{bwamem}.
As illustrated, both distributions range from zero to several hundred or thousand and are not well clustered, implying a substantial amount of warp divergence within GPUs.
As depicted in the figure, the difference of length between the shortest and longest strings can be up to $10\times$ for both the query and reference string.
Provided that the computational complexity is proportional to the multiple of the two string lengths, the workload imbalance could be very large in practice.

\subsection{Memory Inefficiency}
With four-bit sequence packing, it is rational to process $8\times8$ DP cells at a time.
When moving to the next cell,\footnote{Without loss of generality, we assume that cells are processed left to right and top to bottom.} the current cell content can be captured within the registers.
However, the dependency structure of the seed extension also requires long-term storage of the cell information to process the next row of $8\times8$ cells.
This intermediate data well exceed the size of the register file and are  usually stored in the global memory (i.e., DRAM).

However, this scheme incurs two inefficiencies.
First, the intermediate data need not remain in the global memory at the end of the kernel and are considered overhead.
Therefore, reducing this overhead contributes to better performance.
Second, the minimum access size of the GPU’s global memory is 128 B (or 32 B from Volta~\cite{volta} architecture, as identified by~\cite{khairy2018exploring}), whereas the individual cell data size is only 4 B.
If not captured by the L2 cache, this will incur a number of redundant access instances for the same intermediate data, leading to an inefficient kernel.

\begin{table}[t]
\normalsize
    \centering
    \caption{Amount of Data Stored and Accessed for the Existing GPU Aligner}
    \begin{tabular}{lc}
    \toprule
        Data & Quantity \\
        \midrule
        Necessary & $2N$ \\
        Stored & $2N + N^2/4 $ \\
        Accessed (Until Pascal~\cite{pascal}) &$128N + 16N^2$  \\
        Accessed (After Volta~\cite{volta}) & $32N + 4N^2$  \\
        \bottomrule
    \end{tabular}
    \label{tab:moti}
\end{table}

\tab{moti} lists  the data volume, input data necessary for the extension (Necessary), data stored for intermediate data (Stored), and data accessed due to the access granularity (Accessed).
The table reveals that the inefficiency is multifold even with modern architectures, which provides another opportunity for improvements.
In this work, we demonstrate that SALoBa\xspace can remove much of this access, leading to superior performance.

\section{SALoBa\xspace Design}
\subsection{Intra-query Parallelism}
\label{sec:parallel}
Because each cell in the DP table has dependencies from the top, left, and top-left cells, it is known that intra-query parallelism exists in an anti-diagonal form, as depicted in \fig{parallel}.
\JL{chk the sentence here}
Adopting the idea, we use multiple CUDA threads to concurrently handle anti-diagonal elements.
A slight difference is that instead of individual threads computing a single cell at a time, the cells are assigned to threads in $8\times8$ \emph{blocks} because of the 4-bit packing on the reference and query sequences.

In the first version of the library, we decided that 32 threads at most---a warp---should collaborate in processing a query.
Because intra-warp synchronization is relatively cheap or free (for GPUs before introduction of independent thread scheduling~\cite{volta}), this becomes an attractive design choice.
While using more than 32 threads is theoretically possible, it would require threadblock level synchronizations (i.e., \texttt{\small\_\_syncthreads()}) that cause non-negligible overhead.

For communication between the threads, we used the CUDA shared memory.
Using shared memory as a communication channel and reusing it every iteration, only a fixed amount of shared memory is needed, and all access to the shared memory is conflict-free (see \sect{shm}).

\fig{parallel} (a) and (b) illustrate the procedure of SALoBa\xspace.
The largest units of computation are called \emph{chunks},  which are horizontal partitions of the DP table.
The chunks have the height of 32 blocks and the width of the entire query sequence.

\begin{figure*}[t]
    \includegraphics[width=0.96\textwidth]{figs/shm}
    \caption{Lazy spill optimization. (Left): Naive global memory access. (Right): Proposed optimized access. Each small rectangle represents 8\texttimes1 cells stored for dependency. A colored rectangle labeled `t\#' denotes that it is being accessed by thread\# and a gray rectangle labeled `block \#' represent that it is from block id \# and is ready to be spilled to the global memory.}
    \label{fig:shm}
\end{figure*}

As displayed in the figure, a thread processes a block in a \emph{strip} (i.e., a row of blocks) of the chunk per step.
In the first step, the first thread starts the processing by fetching a 32-bit word from each of the query and reference sequences, which is enough to process an  $8\times$8 \emph{block} of cells.

The number of ready-to-process cells increases by one per step in the upcoming steps.
Therefore, the number of threads participating in the computation increases until the 32nd step, forming a 31-step long \emph{prologue}.
In the main loop stage,\footnote{This pattern is also called the `kernel,’ but we decided to call it the main loop stage because it would be confusing with the term GPU kernel.} all threads in the warp simultaneously process blocks in an anti-diagonal manner, advancing one block to the right per step.
This process  continues for $Q - 31$ steps, where $Q$ is the number of the blocks in the query sequence (i.e., $\lceil query\_length/8\rceil$).
When the first thread reaches the end of the row, the \emph{epilogue} starts, and the number of participating threads decreases by one per step until the last thread reaches the end of the row.
At the end of the epilogue, the total number of steps taken for processing the entire chunk is $Q + 31$.

Each time the threads advance to the next step, they access the dependency data from the cells on the top, left, and top-left (\eq{sm}).
\fig{parallel} (c) demonstrates how this is done.
When a block is calculated, the data from the eight cells to the left are what the current thread computed in the previous step.
Therefore, the data can be stored in the register.
The data from the eight cells at the top were computed by the adjacent thread in the previous step.
At the end of processing a block, each thread writes the bottom 8\texttimes1 cells in the shared memory.
In the next step, each thread reads the data in the next position of the shared memory so that it can receive the 8\texttimes1 cells from the above block. \JL{this paragraph has been edited a lot by the editor. Check for its meaning again}
Last, one cell at the top left has dependency on the current block, which is what the thread received from the shared memory in the previous step.
Fortunately, it can also be passed using the register.
Thus, the number of cells stored in the register becomes nine instead of eight.

This scheme not only reduces the warp divergence but also improves memory access.
In the previous technique, all bottom cells in each strip must be stored and read back to/from the global memory.
In contrast, with intra-query parallelism, only the bottom cells of each chunk (not of the strip) are stored and read to/from the global memory.
For a 32-thread warp, this reduces the amount of intermediate data access to 1/32.

\subsection{Lazy Spill to Global Memory Access}
\label{sec:shm}
For processing a chunk, the data from the bottom-most cells must be stored for processing the top-most cells of the next chunk.
Usually, this easily exceeds the shared memory capacity and must be stored in the global memory.
In a naive scheme, as in \fig{shm} (left), the last thread ($t31$) stores the bottom-most cells of a block in the global memory, and the first thread ($t0$) reads them from the global memory as the next chunk is processed.
However, as diagnosed in \sect{moti}, this results in an abundant amount of non-coalesced access, becoming a reason for inefficiency.

To address this problem, we used \emph{lazy spilling} implemented using double-buffered shared memory for each warp to reduce the amount of global memory access.
\fig{shm} (right) indicates the shared memory location where the threads write at each step.
For each warp, a shared memory region is allocated with the size of $2\cdot dim(block)\cdot \#threads$. \JL{is it a cell? or 8 cells?}
As explained in \sect{parallel}, the first thread ($t0$) starts writing on the first location of the shared memory, which is passed to the second thread before the next step.
In each step, the threads read data from the shared memory, process a block, and overwrite the bottom cells in the block to the shared memory.
Then, all threads shift one block location for the next step.

\begin{figure*}[t]
    \includegraphics[width=0.96\textwidth]{figs/vwc_rev}
    \caption{Subwarp scheduling.\vspace{-2mm}}
    \label{fig:vwc}
\end{figure*}

After the prologue, the threads leave a data trail, written by the last ($t31$) thread.
Luckily, these are the data to be stored on the global memory.
When the first thread ($t0$) reaches the end of the buffer at step 64, the shared memory has a consecutive region of data at the left side which are large enough to be processed by a warp.
The threads stop processing and gather together to spill the data to the global memory in a coalesced write.
Afterward, the threads wrap around in the shared memory, and the trail of block data starts forming on the right half of the shared memory, which is later spilled to the global memory.

The spilled data from the previous chunk must be read to process the next chunk  so that the first thread can use the data to calculate the cell values.
This process is also performed in the same way using the double-buffered shared memory in the opposite direction of the writes.

\subsection{Subwarp Scheduling: Trade-off Between Parallelism and Workload Balancing}
The proposed strategy---a warp per query---can eliminate the workload imbalance between threads in a warp because they work together.
However, as demonstrated in \sect{parallel}, a problem exists in the prologue and epilogue.
The average thread utilization over these phases is 50\% (because it has a triangular form), and their size is proportional to the warp size.
The DP table must be significantly larger than the number of threads in a warp to amortize this overhead.
However, according to \fig{histo250}, the size of the table dimension is only around several hundred, which is only a few times larger than the size of 32, and sometimes even smaller.

We chose to split the warp into multiple \emph{subwarps}  and let them process a single query, as in \fig{vwc}, similar to the idea of VWC~\cite{vwc} used in graph processing.
Its effect is essentially a trade-off between workload balancing and resource utilization.

When many subwarps are used (i.e., the subwarp size is smaller), the sizes of the prologue and epilogue become smaller; therefore, the average utilization increases.
However, multiple subwarps make the scheme suffer more often from warp divergence due to workload imbalance.
In contrast, when fewer subwarps are used (i.e., the subwarp size is larger), it is less likely to suffer from imbalance at the cost of lower utilization from increased portions of the prologue and epilogue.
Empirically, we found that using two or four subwarps yields the best results in our setting (see \sect{eval}).

One potential drawback is that using subwarps complicates the memory access optimization.
If the same double-buffered technique is used, the number of threads accessing the consecutive memory address becomes less than 32.
Subwarps larger than eight threads do not incur a considerable problem for architectures later than NVIDIA Volt
~\cite{khairy2018exploring}.
For older architectures, the problem can be solved by allocating $N+32$ slots of shared memory instead of $2N$.
By making the active region rotate around and making all threads in the entire warp
write to 32 slots of data to the global memory together, a full coalescing can be achieved at the expense of a slightly higher shared memory capacity requirement.

\section{Evaluation}
\label{sec:eval}
\subsection{Experimental Setup}
\label{sec:env}
We evaluated SALoBa\xspace on two platforms.
First, as an `affordable’ system, we used a GTX1650 GPU card based on the Turing~\cite{turing} architecture.
The server has a six-core Intel i5-9600K CPU with 16 GB of RAM.
We also conducted experiments on a `high-end’ system with an RTX3090 GPU card based on the Ampere~\cite{ampere} architecture.
The server has a single socket, 12-core AMD EPYC 7272 CPU and 128 GB of RAM.
Both machines run on Ubuntu 18.04 with CUDA version 11.2 and NVIDIA driver 460.27.04.

\input{eval/kernel_all}

For comparison, we used the seed extension kernels from the following libraries as the baseline listed in \tab{baselines}.
\textbf{SOAP3-dp~\cite{soap3dp}} is an early algorithm for GPU-assisted short read alignment, which utilizes inter-query parallelism.
\textbf{CUSHAW} is a family of GPU-assisted short read alignment algorithms.
While most of its members \cite{cushaw1,cushaw2,cushaw3} use GPU only for seeding phase, \textbf{CUSHAW2-GPU}~\cite{cushaw2gpu} accelerates the seed extension part with a GPU using a similar strategy to SOAP3-dp.
Better performance is obtained by compacting the global memory storage format and using CUDA texture memory.
\textbf{NVBIO~\cite{nvbio}} is a library of reusable components designed by NVIDIA to accelerate bioinformatics applications using CUDA.
\textbf{SW\#~\cite{korpar2013sw}} is an algorithm that utilizes intra-query parallelism.
It splits the DP table into multiple anti-diagonal partitions and processes one partition per each kernel launch.
Because it launches the kernel multiple times per query, it targets very long sequences.
\textbf{ADEPT~\cite{adept}} is the most recent work among algorithms that use intra-query parallelism.
It uses shuffle instruction along with binary masking to make use of the anti-diagonal parallelism.
Lastly, \textbf{GASAL2~\cite{gasal2}} is the state-of-the-art library.
It achieves superior performance by further executing the sequence packing on GPU devices.

\begin{table}[b]

\small
    \centering
    \caption{Baseline Kernels Under Comparison}
    \begin{tabular}{lccc}
    \toprule
        Kernel & Parallelism & Bitwidth & Mapping \\
        \midrule
        SOAP3-dp~\cite{soap3} & inter-query & 2 bits & one-to-one \\
        CUSHAW2-GPU~\cite{cushaw2gpu} & inter-query & 2 bits & one-to-many \\
        NVBIO~\cite{nvbio} & inter-query & 2,4,8 bits & one-to-many \\
        GASAL2~\cite{gasal2} & inter-query & 4 bits & one-to-one \\
        SW\#~\cite{korpar2013sw} & INTRA-query & 8 bits & one-to-many \\
        ADEPT~\cite{adept} & INTRA-query & 8 bits & one-to-one \\
        \bottomrule
        \multicolumn{4}{r}{\scriptsize *(All kernels have been modified to have at least four-bit  } \\
        \multicolumn{4}{r}{\scriptsize GPU-assisted packing and to support one-to-one mapping.)}
    \end{tabular}
    \label{tab:baselines}

\end{table}

For a fair comparison, we have put our best efforts into optimizing the existing methods under the same environment.
Most importantly, we assume on-GPU sequence packing for all methods.
GASAL2 is known to achieve state-of-the-art performance using a custom sequence packing kernel executed on GPUs.
While the strategy is successful, the packing is orthogonal to the seed extension itself and the same method can be applied to all other kernels.

We have modified SOAP3-dp, CUSHAW2-GPU, and ADEPT kernels to support five possible types of literals by unifying the kernels to work on four-bit packing.
We left SW\# to use its original eight-bit packing, because modifying it to support packing would
complicate the memory access behavior.
In addition, we have modified some kernels (CUSHAW2-GPU, SW\#, and NVBIO) that only have one-to-many alignment modes to support one-to-one alignment.

\subsection{Kernel Performance Measurements}
In this subsection, we compare the performance of the seed extension kernels without load imbalances using input reads of equal lengths.
To measure the performance under various input lengths, we used an in-house sequence read simulator similar to Wgsim~\cite{wgsim} to generate synthetic reads for each length in the range of 64 to 4096.
Each kernel processed 5,000 reads per call 200 times, and these results were averaged to output the results in \fig{kernelall}.
The execution times were measured with the $cudaEventElapsedTime()$ API.

\fig{kernelall} (a) and (b) present the performance comparison of the methods on the GTX1650 card, and (c) and (d) present that on the RTX3090 card.
For lengths equal to or longer than 128 bp, SALoBa\xspace outperforms all other methods.
For a very short length of 64 bp, the execution time of NVBIO is slightly shorter than SALoBa\xspace (0.42 ms vs 0.51 ms in GTX 1650 and 0.21 ms vs 0.24 ms in RTX 3090).
This is reasonable because with intra-query parallelism, the effect of resource under-utilization from prologue and epilogue is relatively significant.
However, the break-even point is found at 128 bp, where the reduced global memory access of SALoBa\xspace becomes dominant.
The methods other than SALoBa\xspace that utilize intra-query parallelism (SW\# and ADEPT) perform poorly compared to the others.
SW\# is especially slow because it divides a single DP table into multiple kernel calls where one kernel call processes an anti-diagonal group, resulting in very low resource utilization.

A rather surprising observation is that GASAL2 does not always show superior performance over other baselines.
While GASAL2 reports multi-fold speedup compared to its previous work, its speedup is mainly from the on-GPU packing, not the extension strategy.
Because we combined CUSHAW2-GPU with the on-GPU packing module taken from GASAL2, it shows comparable results on GTX1650 for all lengths and slightly better performance on RTX3090 for long sequence lengths.

\begin{comment}

\JL{how does adept do here?}
The gain fluctuates over the sequence length in the range of shorter reads.
At this range, the portions of the prologue and epilogue sharply decrease and provide more benefit.
However, the portion of reduced global memory access from using more threads also rapidly decreases, compensating the gain.
The best acceleration is obtained at 64 bp with 0.39 ms vs. 0.21 ms for 74.3\% acceleration.
The lowest acceleration  is at 256 bp with 13.8\% acceleration.
\end{comment}

The best speedup obtained at short lengths by SALoBa\xspace is observed at 512 bp.
The speedup is 27.7\% on GTX1650 and 43.6\% on RTX3090 against GASAL2.
At longer sequence lengths, some baseline kernels fail to run due to structural limitation (ADEPT) or bounded device memory (NVBIO and SOAP3-dp). \JL{who fails for what}
The performance trend becomes more consistent for longer lengths because the portion of the overhead of global memory access in execution time diminishes.
Another cause for this trend would be the decrease of the relative size of prologue/epilogue compared to the sequence length.
For inputs equal to or longer than 1,024 bp, the speedup of SALoBa\xspace against GASAL2 is consistently around 30\% on GTX 1650 and 50\% on RTX 3090.
Against CUSHAW2-GPU, the speedup is around 40\% on GTX1650 and 20\% on RTX 3090.

\input{eval/ablation}

\subsection{Ablation Study}

\fig{ablation} breaks down the contributions of the three techniques composing SALoBa\xspace.
The speedups are normalized against GASAL2, which performs reasonably well along all sequence lengths we target.
The most effective technique for shorter lengths ($\leq$1024) is subwarp scheduling, which is expected because at shorter inputs, the portions of the prologue and the epilogue are very large.
Therefore, switching to intra-query parallelism yields performance degradation.
Using subwarp scheduling directly reduces the overhead and provides a substantial speedup.

A seemingly large speedup is observed at 64 bp for both GTX1650 and RTX3090, which is counter-intuitive because there is less intra-query parallelism that SALoBa\xspace can exploit.
However, further analyzing \fig{kernelall} (a) and (c) reveals that it is the inefficiency of GASAL2, not speedup of SALoBa\xspace.
GASAL2 has a relatively large memory initialization cost at the beginning, and using a small sequence length fails to amortize it.
Compared to NVBIO instead of GASAL2 at 64 bp, it correctly reflects the intuition that inter-query parallelism is better suited for very short queries.

As the input becomes longer, the gain from using subwarps becomes marginal, and the effects of the other two techniques become significant.
Intra-query parallelism and lazy spilling both contribute to less global memory access.
The former reduces the amount of data stored in the global memory, whereas the latter reduces redundant access by coalescing access better.
Therefore, these become the main driver for performance gains in longer input ranges.
Interestingly, in  RTX3090, intra-query parallelism has more effect, whereas lazy spilling is more effective for GTX1650.
We believe this is explainable by the fact that RTX3090 has a higher computation/memory bandwidth ratio.
RTX3090 has a peak performance of 35.58 TFlops, and its memory bandwidth is 936.2 GB/s using GDDR6X.
On the other hand, GTX has a peak performance of 2.98 TFlops with 128.1 GB/s memory bandwidth.
The ratio is 38.91 Flops/B in RTX3090 and 23.82 Flops/B in GTX1650.
This means that  RTX3090 is generally more bottlenecked at its memory.
The major drawback of intra-query parallelism is the low CUDA thread utilization.
However, as RTX3090 is more bottlenecked in memory bandwidth, it is likely that there are more computational resources to process the data even when some are being idle.
This partially offsets the utilization problem in RTX3090.
\JL{check this again}

\subsection{Real-world Data Experiments}
To test SALoBa\xspace under workload imbalance with the realistic distribution of workload sizes,
we generate seeds from BWA-MEM~\cite{bwamem} using real-world datasets.
It takes a whole reference genome sequence and many
sequence reads to generate seeds that are composed of multiple pairs that can be processed by extension kernels.

For the reference genome sequence, we used the latest release of the human genome assembly project, Build 38 patch release 13 (GRCh38.p13)~\cite{refsequence} that has 3.1G base pairs (bp).  \JL{SY: how long is it?}
For the input sequence reads, we used two datasets downloaded from the sequence read archive~\cite{sequencearxiv}.
The \textbf{dataset A} (SRR835433) is from a 2nd generation sequencer Illumina MiSeq and represents short reads
where each sequence has a length of 250 bp.
The dataset comprises 8.3M sequences which are randomly read genome parts from a human.
Each sequence is read twice, resulting in the total number of base pairs to be 250$\times$8.3M$\times$2=4.1 Gbp.
The \textbf{dataset B} (SRP091981) is from a 3rd generation sequencer PacBio RS and represents long reads.
It also contains randomly read genome parts from a human where there are 82K sequences with variable lengths that averages around 2,000 bp to a total of 182.4 Mbp.
The distributions of the seeds processed by BWA-MEM~\cite{bwamem} are presented in \fig{histo250}.

\input{eval/realkernel}

The results are plotted in \fig{realkernel}.
\fig{realkernel} (a) shows the performance for short read dataset A.
The best speedup of SALoBa\xspace over the baseline GASAL2 is 32.5\%
for GTX1650 and 20.2\% for RTX3090.
SOAP3-dp could not complete the workload on GTX1650 as some of the inputs exceeded the length it could process.
The speedup values observed from SALoBa\xspace are slightly larger than that of \fig{kernelall} (a) and (c) due to the fact that SALoBa\xspace suffers less from workload imbalance.
While the performance of SOAP3-dp and NVBIO are inferior to GASAL2, CUSHAW2-GPU exhibits some speedup over GASAL2 (with the help of GASAL2's on-GPU packing kernel).
However, its speedup is smaller than that of SALoBa\xspace.
ADEPT achieves a similar speedup compared to SALoBa\xspace only on RTX3090, but the speedup comes from placing all the intermediate values in the shared memory (no global memory access), which fundamentally limits ADEPT up to 1024 bp.

An interesting difference between GTX1650 and RTX3090 is that the optimal subwarp size was 16 for GTX1650 and eight for RTX3090 (see \fig{realkernel} (c)).
This outcome is due to the different effectiveness of each technique.
According to \fig{ablation}, the benefit of applying subwarp scheduling for shorter sequences was 2.26$\times$ for GTX1650 and 2.85$\times$ for RTX3090 in geomean.
In GTX1650, the benefit of using subwarps is partially offset by the overhead from workload imbalance, shifting the optimal configuration to 16 threads per subwarp.

\fig{realkernel} (b) reveals the performance for dataset B with longer reads.
On this dataset, SOAP3-dp, ADEPT, and NVBIO fail to run due to their input length limitations.
The speedup of SALoBa\xspace greatly improves compared to that of \fig{kernelall} (b) and (d) because the amount of workload imbalance increases which works in favor of SALoBa\xspace.
The best speedup is around 2.1$\times$ for both GTX1650 and RTX3090.
Subwarps with 16 threads per subwarps gained the best performance, which is the sweet spot between exploiting intra-query parallelism and reducing the workload imbalance.

\section{Related Work}
\subsection{GPU-based Smith--Waterman Algorithm Implementations}

 Accelerating Smith--Waterman algorithm~\cite{sw} with GPUs has been studied for many years.
Some earlier work based on OpenGL library~\cite{liu2006gpu, liu2006bio} addressed issues for mapping the algorithm on graphics pipeline.

A popular CUDA implementation can be found from the CUDASW++ family~\cite{cudasw1,cudasw2,cudasw3}.
They suggest utilizing both the intra- and inter-query parallelism based on a pre-determined sequence length threshold~\cite{cudasw1}.
Similar approaches can be found from GSWABE~\cite{gswabe} and gpu-pairalign~\cite{gpu-pairalign}.
In addition, focused on all-to-all patterns for protein DB alignments, query profiling optimization~\cite{cudasw2} and CPU-GPU hybrid parallelism~\cite{cudasw3} have been suggested.
However, such optimizations are often too specific to protein alignments, and cannot be easily generalized to other domains such as DNAs.

In such regard, the CUDAlign family~\cite{cudalign1, cudalign2, cudalign3, cudalign4} focus on a general Smith--Waterman algorithm using intra-query parallelism.
CUDAlign~\cite{cudalign1} splits the DP table into multiple blocks and distributes them to threadblocks.
Then, the communication between the threadblocks are done using a dedicated region in the global memory called the `horizontal bus' and `vertical bus'.
The algorithm is later extended to support linear space algorithm to reduce the memory size~\cite{cudalign2}.
Its recent versions support multi-GPU alignment with execution time prediction~\cite{cudalign3} and traceback~\cite{cudalign4}.
When the DP table is too large, approaches such as MultiBP~\cite{multiBP}, MASA~\cite{masa}, and SW\#~\cite{korpar2013sw} suggest block pruning~\cite{cudalign21} that removes some blocks that can never achieve the optimal score.

However, these algorithms are optimized for huge sequences that sometimes do not even fit into a single GPU card~\cite{khajeh2010acceleration}.
For DNA alignment softwares with seed-and-extend strategies, the input sequence length for the extension ranges around several hundreds even with the long DNA reads.
Because of this, libraries that rely on inter-query parallelism~\cite{nvbio,gasal2} often perform much better for DNA alignments as demonstrated in \cite{adept}.
On the other hand, SALoBa\xspace outperforms the previous approaches on DNA alignment scenarios by adopting several careful optimizations.

\subsection{GPU-accelerated Sequence Alignment Softwares}

Some approaches design an end-to-end DNA alignment software that are friendly to GPU architectures.
SARUMAN~\cite{saruman} and GPU-RMAP~\cite{gpu-rmap} are some early methods that
implement hashtable lookups on GPUs to perform short read mappings.
With the increased use of the Burrows--Wheeler transform (BWT)~\cite{bwt} based on suffix-trees, BarraCUDA~\cite{barracuda}, SOAP3~\cite{soap3} and CUSHAW1~\cite{cushaw1} were introduced with GPU implementations of BWT.
SOAP3-dp~\cite{soap3dp}, CUSHAW2-GPU~\cite{cushaw2gpu}, and LOGAN~\cite{zeni2020logan} are later methods that adopt seed-and-extend strategy,
and Arioc~\cite{arioc} further expands the strategy to multiple GPUs.
However, these methods often sacrifice alignment quality for the sake of throughput.
For example, SOAP3 and LOGAN only supports limited number of mismatches, and CUSHAW2-GPU packs the sequences to two bits by converting the fifth `N' bases to random bases.

Although the approaches above show good throughput with only a small amount of quality degradation,
the industry has grown to value quality over speed.
As high-quality algorithms such as BWA-MEM~\cite{bwamem} became a de facto standard, the usability of the GPU-aware alignment softwares were limited.
Some approaches~\cite{zaid1,zaid2,zaid3,adept} tackle this problem and design a seed extension kernel general enough to be used for BWA-MEM using intra-query performance.
However, later approaches based on inter-query parallelism outperformed these kernels, which is the strategy adopted by the current state-of-the-art methods such as NVBIO~\cite{nvbio} or GASAL2~\cite{gasal2}.

Finally, it is worth mentioning cuBLASTP~\cite{cublastp}, a GPU accelerated protein search algorithm that uses a variant of seed-and-extend strategy.
However, cuBLASTP only accelerates the seeding part using GPU, and leaves the extension part to CPU for utilizing CPU-GPU hybrid parallelism.

\begin{comment}

\subsection{Algorithmic Approaches}
As the environment around both the sequencing and processing systems develops, it is accompanied by advances in processing algorithms.
\cite{vasimuddin2019efficient} provided some algorithmic optimizations for BWA-MEM to perform better on modern architectures.
While BWA-MEM~\cite{bwamem} is used as an almost de facto standard for aligning short reads, it is not as efficient for long reads, leading to the development of algorithms designed for long reads that often use the seed-chain-fill strategy~\cite{blasr, minimap2}.

Aside from standard techniques, there have been some pioneering proposals.
AnySeq~\cite{muller2020anyseq} is an abstraction for high-performance sequencing, targeting various hardware platforms, and \cite{besta2020communication} proposed using Jaccard similarity as a metric to perform genome comparisons.
In addition, GenieHD~\cite{kim2020geniehd} used hyperdimensional computing to perform DNA pattern matching and proposed an accelerator for the algorithm.

Another stream is to scale the system for distributed environments.
PaKman~\cite{ghosh2019pakman} is a scalable algorithmic solution for K-mer and contig generation that forms a genomics pipeline.
Moreover, \cite{guidi2021parallel} designed a distributed  memory algorithm for de novo assembly.
\cite{nisa2021distributed} suggested GPU-based distributed memory K-mer counting.
\cite{dash2021scaling} used a large-scale cluster of several thousand GPUs to identify combinations of genes with mutations.

\end{comment}

\subsection{FPGA/ASIC Acceleration}
As interest in the genomics pipeline has increased, more researchers have considered designing dedicated accelerators using FPGAs or ASICs.
\cite{zhang2007implementation} provided a implementation of the Smith--Waterman algorithm on an Altera FPGA, followed by a few other studies using banded algorithms~\cite{chen2013hardware, blastp} or flexible systolic arrays~\cite{houtgast2015fpga,ahmed2015heterogeneous}.
Moreover, DRAGEN~\cite{dragen} is an FPGA-based platform from Illumina that implements a full end-to-end analysis into an FPGA.
 Darwin~\cite{turakhia2018darwin} is an ASIC accelerator that speeds up the whole genome sequencing through the co-design of both the software and hardware, powered by gapped filtering.
GenAx~\cite{fujiki2018genax} is an automata-based accelerator for both seeding and extension.

Although these approaches provide significant speedup,
FPGAs or ASICs are much harder to design and it requires a long time to reach the market.
Compared to them, GPU-based alignment software could be a quick solution easily adopted by any system that has GPUs attached.

\section{Discussion}

\subsection{CUDA Shuffle Instructions}
The CUDA shuffle instruction is an alternative to
shared memory by allowing a direct register-to-register data exchange for inter--thread communication.

However, using shuffle-based communication did not add any additional speedup on top of SALoBa\xspace.
This outcome aligns with the CUDA specification that the throughput of the shuffle instructions is almost the same as the nonconflicting shared memory reads~\cite{warpperf}.

One potential gain from shuffle instructions is reducing shared memory usage.
However, the portion of data cached in the shared memory was negligible in the current scenario and did not provide noticeable speedup.

\subsection{Banded Algorithms}
Many seed extension methods use the banded algorithm to reduce the computational burden.
Taking advantage of the fact that the best matching path of the DP table is usually near the diagonal, processing cells within some predetermined width (the \emph{band}) often yields solutions of sufficient quality.
However, most GPU-based extension libraries do not exploit this.

One problem is that the band sizes for each query are often different, which worsens load balancing.
However, we envision that banded algorithms would have more benefits and could be adopted for better performance with longer reads.

\subsection{Multiple GPUs}
It is often beneficial to use multiple GPUs, especially when equipped on a single machine.
We believe that extending this work by splitting the queries into equal numbers and assigning them to multiple GPUs would be straightforward.
One possible drawback of such a strategy would be the load imbalance between the GPUs.
However, the penalty would be small compared to the thread-level imbalance problem.
If this becomes significant, it could be solved by dynamic assignment or preprocessing with approximate sorting.

\section{Conclusion}
We proposed SALoBa\xspace, a new GPU-based seed extension library, exhibiting  state-of-the-art performance over the existing work. The experiments reveal that the software performs well on modern devices, and the performance gain increases with recent GPUs. We believe this software will be useful in fields where sequence alignment times are critical for diagnosing fatal diseases.

\section*{Acknowledgements}
\small{
This work has been supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2022R1C1C1008131,
2022R1C1C1011307)
and
Institute of Information \& communications Technology Planning \& Evaluation
(IITP) grant funded by the Korea government (MSIT)
(2020-0-01361,
Artificial Intelligence Graduate School Program (Yonsei University),
2021-0-00853, Developing Software Platform for Programming of PI
).
}

\JL{Fix the capitalizations in ref}

\bibliographystyle{IEEEtranS}
\balance

\subsection{FPGA/ASIC Acceleration}

Due to its flexibility, FPGAs often become the target platform for emerging applications requiring fast deployment time.
\cite{zhang2007implementation} provided  a detailed implementation of the Smith--Waterman algorithm on an Altera FPGA, followed by a few other studies using banded algorithms~\cite{chen2013hardware, blastp} or flexible systolic arrays~\cite{houtgast2015fpga,ahmed2015heterogeneous}.
Moreover, DRAGEN~\cite{dragen} is an FPGA-based platform from Illumina that provides speedup by implementing the full end-to-end analysis into an FPGA.
These approaches adopt similar intra-query parallelization strategies with the GPU-based ones.
However, FPGAs are still hard to find from the market as accelerators.

Some researchers have considered dedicated ASIC accelerators~\cite{primer}.
Compared to FPGAs, ASICs provide better performance and energy efficiency at the expense of flexibility.  Darwin~\cite{turakhia2018darwin,turakhia2019darwin} is an accelerator that speeds up the whole genome sequencing through the co design of both the software and hardware, powered by gapped filtering.
GenAx~\cite{fujiki2018genax} is an automata-based accelerator for both seeding and extension that supports configurable scoring and traceback functionality.
Moreover, GenASM~\cite{genasm} is an accelerator framework for bitvector-based algorithms.

\begin{comment}
The EXMA~\cite{jiang2021exma} accelerator designs a row-buffer friendly algorithm for the FM  index to build a dedicated accelerator.
Further, SeedEx~\cite{fujiki2020seedex} is an extension accelerator that uses a banded algorithm. By opportunistically using a narrower band than required and providing a function that checks the optimality, it greatly reduces the width of the band used for extension algorithms.

\end{comment}

Although these researchers have claimed that they can achieve multi fold speedup, sometimes by a few hundred, ASICs fall into a solution that requires at least a few years to reach the market.
Compared to them, GPU-based alignment software could be a quick solution easily adopted by any system that has GPUs attached.

}


\section{Introduction}
Optical heterodyne detection is a commonly used method for measurement of amplitude and phase of spectral components of an optical beam \cite{kuri2003optical, delange1968optical, chtcherbakov2007optical, hall1981optical, jacobs1988optical}. Despite its frequent usage, complete and concise mathematical descriptions of the heterodyne measurement process are difficult to come by, with many sources simply stating a non-generalised final result \cite{keyes2013optical, menzies2005laser}. The lack of a readily available comprehensive description can make it difficult to deeply understand the method, including underlying assumptions, which increases the chance of experimental or analytical errors.

Here, we present a complete and concise derivation of optical heterodyne measurement of a beam consisting of regularly spaced frequency components. Such beams are commonly generated using electro-optic modulators, and ultrafast laser techniques, and where the number components is large, the optical field is referred to as a frequency comb. We attempt to use consistant notation, and where appropriate, give a physical interpretation to the mathematical results.

At a given position any electric field strength can be written as
\begin{equation}
	E(t) = \frac{1}{2}\mathcal{E}(t)e^{-{\rm{i}}\omega_0t} + \frac{1}{2}\mathcal{E}^*(t)e^{{\rm{i}}\omega_0t}
	\label{a}
\end{equation}
where $\mathcal{E}(t)$ is the complex envelope function, $\mathcal{E}^*$ is its complex conjugate, and $\omega_0$ is the reference angular frequency. The choice of reference frequency is completely arbitrary, and can be positive or negative. However, the reference is \emph{usually} chosen so that the variation in time of $\mathcal{E}(t)$ is as small as possible. In practice, the reference frequency is often taken as the carrier frequency of a beam, so the complex envelope fully represents only the modulations to this carrier (both in phase, amplitude, and any additional frequency components).

The intensity of a given electric field is given by
\begin{equation}
	I(t) = c\varepsilon \frac{1}{T}\int_{-T/2}^{T/2} E^*(t)E(t){\rm{d}}t,
	\label{ba}
\end{equation}
where $T$ is an integer number of periods of the optical cycle, with minimum one period.

Subbing equation \ref{a} into this expression yields:
\begin{align}
	I(t) &= \frac{c\varepsilon}{4T}\int_{-T/2}^{T/2} \left[\mathcal{E}^*(t)\mathcal{E}(t) + \mathcal{E}(t)\mathcal{E}^*(t) + \mathcal{E}^2(t)e^{-{\rm{i}}2\omega_0t} + \mathcal{E}^{*2}(t)e^{{\rm{i}}2\omega_0t}  \right]{\rm{d}}t \label{b0}\\
	     & = \frac{c\varepsilon}{2T}\int_{-T/2}^{T/2} \mathcal{E}^*(t)\mathcal{E}(t){\rm{d}}t + \frac{c\varepsilon}{4T}\int_{-T/2}^{T/2} \left[ \mathcal{E}^2(t)e^{-{\rm{i}}2\omega_0t} + \mathcal{E}^{*2}(t)e^{{\rm{i}}2\omega_0t}  \right]{\rm{d}}t
	\label{b}
\end{align}
If we have chosen $\omega_0$ appropriately such that $\mathcal{E}$ is slowly varying in time, the expression in the first integral can be pulled out of the integral, and the integral evaluates to $T$. The two terms on the right oscillate at twice the optical frequency, and so will integrate to zero. Note this definition of intensity breaks down where the amplitudes vary significantly over the time scale of an optical cycle, such as with ultra-short laser pulses.

Averaging over the fast oscillations gives a convenient expression for intensity written in terms of the complex amplitudes:
\begin{align}
	I(t) &= \frac{c\varepsilon}{2}\mathcal{E}^*(t)\mathcal{E}(t)\label{c0}\\
	     &=\frac{c\varepsilon}{2}|\mathcal{E}(t)|^2.
	\label{c}
\end{align}

Now, when we interfere two beams, we will have the signal $E_{\rm{sig}}(t)$ and the local oscillator $E_{\rm{LO}}(t)$ fields.
The signal can be written as
\begin{equation}
	E_{\rm{sig}}(t) = \frac{1}{2}\mathcal{E}_{\rm{sig}}(t)e^{-{\rm{i}}\omega_0t} + \frac{1}{2}\mathcal{E}_{\rm{sig}}^*(t)e^{{\rm{i}}\omega_0t}.
	\label{d}
\end{equation}

A frequency comb with tooth spacing $\omega_{\rm{SB}}$ (`SB' is signifies `side band'), is written as a sum of fields each with a constant amplitude and a frequency that differs by an integer multiple of $\omega_{\rm{SB}}$. For compactness, a reference frequency can be pulled out of each term, and the total signal field can be written as equation \ref{d}, where the complex envelope is given by:
\begin{equation}
	\mathcal{E}_{\rm{sig}}(t) = \sum_{n}a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}
	\label{e}
\end{equation}
where $n$ can take any integer value (positive, negative or zero), $a_n$ is the real amplitude of the $n^{\rm{th}}$ comb tooth, and $\phi_{n}$ is the corresponding phase of that frequency component (cf. \ref{apxc1}). Note that this is simply a general expression that represents a frequency comb, it is not a recipe for how to generate one. It is well known that modulating either the phase or amplitude of a carrier periodically puts sidebands on that carrier, but here we are simply concerned with how to represent (and then measure) a field that \emph{does} have evenly spaced sidebands.

The full signal electric field is therefore written as:
\begin{equation}
	E_{\rm{sig}}(t) = \frac{1}{2}\sum_{n}a_{n}\left[e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}e^{-{\rm{i}}\omega_0t} +e^{{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}e^{{\rm{i}}\omega_0t}\right].
	\label{f}
\end{equation}

In a heterodyne measurement, the signal beam is combined with a local oscillator beam, which is shifted in angular frequency relative to the signal carrier by an amount $\Delta\omega$. The complex envelope can therefore be written as:
\begin{equation}
	\mathcal{E}_{\rm{LO}}(t) = be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})},
	\label{g}
\end{equation}
where $b$ is the local oscillator real amplitude.

The full local oscillator electric field can then be written as:
\begin{equation}
	E_{\rm{LO}}(t) = \frac{b}{2}\left[e^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})}e^{-{\rm{i}}\omega_0t} + e^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})}e^{{\rm{i}}\omega_0t}\right].
	\label{h}
\end{equation}

Now, when the signal and local oscillator beams are combines on a beamsplitter, the output of a single port is just the sum of their fields. We will only be interested in a single port, so the phase shift introduced by reflection of a mirror is irrelevant to us and will be ignored (as will the halving of their amplitudes). The result is:
\begin{align}
	E_{\rm{tot}}(t) &= E_{\rm{sig}}(t) + E_{\rm{LO}}(t)\\
	     &= \frac{1}{2}\left[\sum_{n} a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right]e^{-{\rm{i}}\omega_0t}
	        +\frac{1}{2}\left[\sum_{n} a_{n}e^{{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right]e^{{\rm{i}}\omega_0t}
	\label{i}
\end{align}
where the terms have been grouped by pulling out their common exponential factor.

We can see that equation \ref{i} is in the form of equation \ref{a}, where the complex envelope function is given by the term in the square brackets. We can therefore immediately write down the resulting time varying intensity for the combined field by substituting the envelope and its conjugate into equation \ref{c0}. Note that when calculation the intensity the summations over the signal teeth will multiply and hence the summations must be performed separately. The index label of one summation is therefore changed to $m$.

\begin{align}
	I_{\rm{tot}}(t) &= 	\frac{c\varepsilon}{2} \left[\sum_{n} a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right] \times
					\left[\sum_{m}a_{m}e^{{\rm{i}}(m\omega_{\rm{SB}}t + \phi_{m})} + be^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right] \label{j0}\\
	     &=  \frac{c\varepsilon}{2}\left\{
	     b^2 + \sum_{n, m} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}
	      + \sum_{m}a_{m}be^{{\rm{i}}\left[(m\omega_{\rm{SB}}-\Delta\omega)t + \phi_{m} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \right\}\\
	     &=  \frac{c\varepsilon}{2}\left\{
	     b^2 + \sum_{n}a_{n}^2
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}
	      + \sum_{n}a_{n}be^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \right\}
	      \label{j1}
\end{align}

So, equation \ref{j1} tells us exactly what the \emph{intensity} of the combined field will be at a given point in time. Clearly the intensity will vary periodically in time so this is commonly referred to as the \emph{beat note}, although in the case of a signal beam which is made up of many comb teeth, there will be multiple frequency components present.

We are often primarily interested in knowing something about these different frequency components, so a common way to analyse the heterodyne signal is using a spectrum analyser to display at the power spectrum of the photocurrent from a photodiode. The photocurrent is directly proportional to intensity, so the generated power spectrum is directly proportional to the power spectrum of the intensity. A spectrum analyser simply Fourier transforms the signal, and then displays the amplitude (or power, which is proportional to the square of the amplitude), so it is useful to see what the analytic form is of the Fourier transform of the time varying intensity, $\breve{I}(\omega)$:
\begin{equation}
	\breve{I}(\omega) = \frac{1}{\sqrt{2\pi}}\int I(t)e^{{\rm{-i}}\omega t}{\rm{d}}t,
	\label{k}
\end{equation}
where the integral is understood to be over $-\infty$ to $\infty$.

Substituting equation \ref{j1} into the Fourier transform yields:

\begin{align}
	\breve{I}_{\rm tot}(\omega) &=    \frac{c\varepsilon}{2\sqrt{2\pi}}  \int
	    \Big\{
	     b^2 + \sum_{n}a_{n}^2
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]} \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}be^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \Big\}
	     e^{{\rm{-i}}\omega t}{\rm{d}}t  \label{l0}\\
	     &=  \frac{c\varepsilon}{2\sqrt{2\pi}}
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\int e^{{\rm{-i}}\omega t}{\rm{d}}t
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}\int e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b\int e^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t
	      + \sum_{n}a_{n}b\int e^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t
	     \Big\}
	       \label{l2}\\
	     &=  \frac{c\varepsilon}{2\sqrt{2\pi}}
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\int e^{{\rm{-i}}\omega t}{\rm{d}}t
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\int e^{{\rm{i}}(m-n)\omega_{\rm{SB}}t}e^{{\rm{-i}}\omega t}{\rm{d}}t \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}be^{{\rm{i}} (\phi_{n} - \phi_{\rm LO})}\int e^{{\rm{i}}(n\omega_{\rm{SB}}-\Delta\omega)t }e^{{\rm{-i}}\omega t}{\rm{d}}t
	      + \sum_{n}a_{n}be^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\int e^{{\rm{-i}}(n\omega_{\rm{SB}}-\Delta\omega)t}e^{{\rm{-i}}\omega t}{\rm{d}}t
	     \Big\}
	       \label{l3}
\end{align}

The Fourier integrals can be evaluated by noting that:
\begin{equation}
	\frac{1}{\sqrt{2\pi}}\int e^{{\rm i}\alpha t} e^{-{\rm i}\omega t} {\rm d}t= \sqrt{2\pi}\,\delta(\omega - \alpha),
	\label{m}
\end{equation}
where $\delta(\omega)$ is the Dirac delta function.

Substituting this in equation \ref{l2}, taking the corresponding value of $\alpha$ in each integrand gives:
\begin{align}
	\breve{I}_{\rm tot}(\omega) &=   \frac{c\varepsilon\sqrt{2\pi}}{2}
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\delta(\omega) \notag\\
	     &\qquad\qquad+ \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\delta(\omega-[m-n]\omega_{\rm{SB}}) \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})}\delta(\omega-[n\omega_{\rm{SB}}-\Delta\omega]) \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\delta(\omega+[n\omega_{\rm{SB}}-\Delta\omega])
	     \Big\}
	       \label{n}
\end{align}

Equation \ref{n} describes three combs of teeth the in frequency domain, each with spacing $\omega_{SB}$, and a peak also at zero frequency representing the DC level. The two combs that sum only over $n$ represent actual optical frequencies present in the signal beam. There are two of these combs because each actual unique signal comb tooth has a positive and negative frequency component when the original electric field is represented in a form such as \ref{apxd}.

The term that sums over $n$ and $m$ represents the beatnotes formed when the signal beam comb teeth beat against \emph{each other}. Different pairs of frequency components with equal frequency spacing will beat at the same frequency, and may add or subtract to that beatnote depending on their phase. For a set of sidebands produced by a phase only electro-optic modulator, the different phases between all frequency components add in a way so this ``cross comb" does not show up at all.

Depending on the number of teeth, the sideband spacing, and the local oscillator frequency offset, it is possible for teeth belonging to all three combs to overlap. It is generally undesirable for the teeth to overlap, since this prevents measurement of the amplitude and phase of individual of the individual frequency components that make up the signal beam frequency comb. Figures \ref{fig1}, \ref{fig2}, \ref{fig3} and \ref{fig4} show how combs can begin to intersect as the number of teeth get larger. Note that $\Delta\omega$ has deliberately been chosen to be an integer plus one third times the sideband spacing, so that even though all three combs become interspersed, the individual comb teeth do not lie on top of one another. The amplitude of the signal comb teeth was was set to be a Gaussian simply as a visual aid.
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_1_teeth.pdf}
    \caption{Beatnote from comb with single tooth.}
    \label{fig1}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_3_teeth.pdf}
    \caption{Beatnote from comb with 3 teeth.}
    \label{fig2}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_10_teeth.pdf}
    \caption{Beatnote from comb with 10 teeth.}
    \label{fig3}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_30_teeth.pdf}
    \caption{Beatnote from comb with 30 teeth.}
    \label{fig4}
\end{figure}

Equation \ref{n} gives the Fourier amplitude at a given frequency, but in reality, any measurement we make of the Fourier amplitudes is over some small range of frequencies. In a spectrum analyser, this range is set by the resolution bandwidth RBW, and in a numerical discrete Fourier transform of a time domain signal, it is set by the spacing of the frequency array (which is the inverse of the total length of time sampled). So experimentally, to extract the actual value for the Fourier amplitude at some frequency $\bar{\omega}$, we must integrate over some small frequency range : $\omega=\bar{\omega}\pm\epsilon$, and the resulting quantity will be denoted $\bar{\breve{I}}_{\rm tot}(\bar{\omega})$. Integration over the Dirac deltas means that the resulting values of $\bar{\breve{I}}_{\rm tot}(\bar{\omega})$ will be finite.

\begin{align}
	\bar{\breve{I}}_{\rm tot}(\bar{\omega}) &= \int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}  \breve{I}_{\rm tot}(\omega)  {\rm d}\omega  \label{o1} \\
	    &=\frac{c\varepsilon\sqrt{2\pi}}{2}
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)
	     \int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega){\rm d}\omega \notag\\
	     &\qquad\qquad+ \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega-[m-n]\omega_{\rm{SB}}) {\rm d}\omega \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega-[n\omega_{\rm{SB}}-\Delta\omega]) {\rm d}\omega \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega+[n\omega_{\rm{SB}}-\Delta\omega]) {\rm d}\omega
	     \Big\}
	       \label{o2}
\end{align}

Equation \ref{o2} provides us with a simple method to determine the amplitude and phase of the optical signal comb teeth by measuring the radio frequency beat note. If the time domain radio signal is measured on an oscilloscope, it can be Fourier transformed numerically to show the radio-frequency combs. Assuming the teeth of the different combs are not overlapped then the Fourier amplitude of at the frequency of the $n$th tooth has a contribution from only a single term of equation \ref{o2}. Under these conditions, the optical amplitude and phase of the $n$th tooth is related to the Fourier transform of the detected intensity simply by:
\begin{equation}
	\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n) = a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})},
	\label{p}
\end{equation}
so that the optical amplitude of the tooth is:
\begin{equation}
	a_{n} = \frac{|\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n)|}{b},
	\label{q}
\end{equation}
and the optical phase of the tooth is:
\begin{equation}
	\phi_{n} = {\rm Arg}[\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n)]  + \phi_{\rm LO}.
	\label{r}
\end{equation}

It is worth noting that the optical phase of the local oscillator $\phi_{\rm LO}$ will not usually be fixed in a heterodyne measurement. While it is certainly possible to set up some phase locking circuit by either feeding back to the laser or to a piezo controlled mirror, generally the $\phi_{\rm LO}$ will drift over time, as is seen in any unlocked interferometer. Therefore, the recovered value of $\phi_{n}$ will also drift over time. However, so long as $\phi_{\rm LO}$ stays constant over the time scale of a single measurement, then the \emph{difference} between all the $\phi_{n}$s will always be measured to be the same, even though $\phi_{\rm LO}$ may drift.

A final word about the beat note formed by the comb teeth beating against each: just because there is a comb, \emph{does not} mean you will actually observe any modulation of the intensity. Looking at the second term in equation \ref{n} (the term responsible for teeth beating with each other), it can be seen that pairs of $n,m$ with the same difference, do not necessarily have the same phase. If the phase of different pairs is not the same, then it is possible (indeed likely) that subsequent pairs will not all add in a way that increases the magnitude, and may totally cancel each other out. This is the case when a beam is phase modulated using an EOM. If the beam intensity is subsequently measured with only a photodiode, no intensity modulation is seen at all. However beating with a frequency shifted local oscillator reveals that there may be many other frequency components present.