Further author information: E-mail<EMAIL_ADDRESS>Telephone: 1
(202) 767-9532

# How to find real-world applications for compressive sensing

Leslie N. Smith U.S. Naval Research Laboratory  Washington  D.C. 20375 USA

###### Abstract

The potential of compressive sensing (CS) has spurred great interest in the
research community and is a fast growing area of research. However, research
translating CS theory into practical hardware and demonstrating clear and
significant benefits with this hardware over current, conventional imaging
techniques has been limited. This article helps researchers to find those
niche applications where the CS approach provides substantial gain over
conventional approaches by articulating lessons learned in finding one such
application; sea skimming missile detection. As a proof of concept, it is
demonstrated that a simplified CS missile detection architecture and algorithm
provides comparable results to the conventional imaging approach but using a
smaller focal plane array. The primary message is that all of the excitement
surrounding CS is necessary and appropriate for encouraging our creativity but
we all must also take off our ”rose colored glasses” and critically judge our
ideas, methods and results relative to conventional imaging approaches.

###### keywords:

Compressed sensing, target detection

##  1 INTRODUCTION

The data deluge of the past few decades sparked a great deal of research in
image compression. In image compression, such as JPEG, most of the acquired
data is accurately approximated by much fewer coefficients of a basis such as
wavelets. An interesting research question was posed by Donoho[1] “why go to
so much effort to acquire all the data when most of what we get will be thrown
away? Can we not just directly measure the part that will not end up being
thrown away?” Compressive sensing (CS) bypasses this wasteful acquisition
process but rather only measures the relatively fewer values necessary to
reconstruct the signal or image. This amazing ability has been the topic of a
fast growing area of research as is illustrated in a paper by Elad[2] (see his
Figure 2).

However, translating the fundamentals of CS into practical hardware with
significant benefits over conventional imaging has been limited. As stated in
a JASON report for 2012[3], “The CS literature often deals with idealized
situations”. Despite thousands of paper published yearly on this subject, few
head to head comparisons between CS and conventional imaging techniques appear
in the literature. The advantages of CS are derived from a savings in terms of
size, weight, power, and costs compared to a conventional sensor. The
disadvantages come from the need to post-process the measurements and a
reduction in peak signal to noise ratio (PSNR). This tradeoff indicates
situations where implementing a CS architecture might be beneficial.

For example, our efforts to directly compare CS image construction and super-
resolution to conventional imaging with super-resolution for aerial imagery
proved conventional imaging methods to be preferable. More specifically, we
compared using the same low resolution focal plane array (FPA) in a CS camera
mode and in a conventional sensor with the goal of attaining a 4 times greater
resolution image reconstruction (double the width and height resolution of the
FPA). For the conventional camera approach we simply used bicubic
interpolation to obtain the higher resolution image. Interpolation is
computationally fast and only created some blurring in the final result. On
the other hand, a variety of compressive sensing approaches in the literature
were attempted, such as the CS imager described by Willett, et.al. [4] and the
algorithm provided by Romberg [5]. In some cases the results contained
artifacts and in all cases the CS approach was more computationally expensive.
These experiments indicated that whenever conventional techniques are easily
implemented, conventional imaging techniques are preferable to CS techniques.
These tests indicate that CS is most valuable in order to overcome limitations
of conventional techniques; that is, when conventional imaging is hard to do.

The aim of this article is to encourage and assist researchers in finding
those niche applications where the CS approach provide substantial gain
compared to conventional approaches. The main contributions of this paper are
twofold: first, to present lessons learned to help find practical, real-world
applications derived from the author’s recent experience in evaluating CS
applications; and second, these guidelines are illustrated by showing the
benefits of CS versus conventional imaging of a new application - sea skimming
missile detection.

Some of the past successful applications for CS provides guidance in finding
other successful applications. Two such early wins are SparseMRI[6] and radio
interferometry[3]. Analysis of the methods behind these two very different
applications show many similarities. From this analysis the following
guidelines for successful application of CS can be derived:

  1. 1.

CS most useful when it can overcome hardware, financial, bandwidth, or
battlefield limitations,

  2. 2.

A priori knowledge of what you will be imaging,

  3. 3.

The signal is sparse in the pixel domain or there exist a transformation to
create a sparse representation and additionally the sparsity is orders of
magnitude smaller than the number of pixels,

  4. 4.

The measurement matrix is (optimally) incoherent with the sparsifying
transform,

  5. 5.

There exists physically realizable hardware architecture corresponding to the
theoretical design.

The first two conditions come out of practical considerations and eliminate
many potential applications for CS. This implies that it is necessary to find
niche applications where CS has unusual strengths such that it is better to
use CS for those strengths than use a conventional sensor. The next two
conditions derive from CS theory and are prevalent throughout the literature.
The last condition grounds the architecture in a realizable CS camera.

In the next section we will provide a brief description of compressive sensing
prior to describing a new, practical application for CS.

###  1.1 Compressive Sensing

Compressive sensing (CS) grew out of the work of initial work by Emmanuel J.
Candès, Justin Romberg and Terence Tao[7], and David Donoho[1]. At it’s heart,
compressive sensing is about breaking the Nyquist barrier; that is the ability
to image when collecting much fewer measurements than twice the highest
frequency dictated by the Nyquist theorem. This is possible because there is
much less useful information than indicated by the frequency. The key to
compressive sensing’s success relies on the signal being sparse or
compressible in either its natural state or sparse in a known basis. In this
situation, solving an underdetermined, linear system of equations can uncover
the essential information within the few measurements. Transform coding[8, 9]
relies on finding a basis where the signal is sparse or compressible.

Here we provide a description of key terms:

  1. 1.

Sparsity - a signal is k sparse if it has at most k terms or if it can be
exactly represented by at most k atoms from a transform or dictionary.

  2. 2.

Compressible - the coefficients of a representation of a signal decrease
exponentially, implying that the signal can be accurately represented by it
largest k coefficients.

  3. 3.

Coherence/incoherence - the metric used for a measurement transform that is
the largest absolute product between the a sparse or compressible signal and
the measurement matrix. Intuitively, a small coherence means the signals are
spread out in the measurement domain, which leave ”clues” to their existence.
Fourier transforms is a good example of incoherence since a Dirac function or
spike in the signal domain is spread out throughout the Fourier/frequency
domain.

  4. 4.

Partial Fourier (PF) coefficients - undersampling the Fourier coefficients and
retaining only a fraction of the coefficients.

More formally, one seeks to recover a $k$-sparse signal $x\in\mathbf{R}^{N}$
based on the measurements $y\in\mathbf{R}^{M}$, each of which is a different
linear combination of the $x$. Here one expects $k<M\ll N$, where a lower
bound of $M\approx k\times log(N)$ is given in one of the original papers by
Candès, Romberg, and Tao[7] but more recent research claim even fewer
measurements are sufficient. In the event that $x$ is not naturally sparse, we
assume the existence of a linear transform $\Psi$ such that $x=\Psi\theta$ and
$\theta\in\mathbf{R}^{N}$ is $k$-sparse in the coefficient vector. Given this
notation, assume the measurement process is described by the equation

$y=A\Psi\theta+n$ (1)

where the $M\times N$ matrix $A$ models measurement process, and
$n\in\mathbf{R}^{M}$ represents noise in the measurements. CS theory has shown
that it is possible to accurately reconstruct $\theta$ given $y$ by solving
the optimization problem [7]

$\operatorname*{arg\,min}_{\theta}\parallel\theta\parallel_{0}~{}~{}subject~{}to\parallel
y-A\Psi\theta\parallel_{2}^{2}<\epsilon$ (2)

where $\epsilon$ is a small scalar related to the noise variance. If there
isn’t any noise in the measurements ($n=0$) this situation is called the exact
model and equation 2 reduces to

$\operatorname*{arg\,min}_{\theta}\parallel\theta\parallel_{0}~{}~{}subject~{}to~{}~{}y=A\Psi\theta$

In the past few years a significant amount of research was directed to solving
Equation 2 and reconstructing the original image. It is NP hard to solve this
equation and it is typically solved approximately either by greedy methods or
$L_{1}$ relaxation methods. Greedy methods include Orthogonal Matching
Pursuits (OMP)[10, 11], CoSaMP[12], Subspace Pursuits[13], thresholding[9],
and others. Since the objective function $\parallel.\parallel_{0}$ in Equation
2 in non-convex, it is common in the literature to “relax” the objective
function with a convex approximation $\parallel.\parallel_{1}.$ Hence one
instead solves

$\operatorname*{arg\,min}_{\theta}\parallel\theta\parallel_{1}~{}~{}subject~{}to\parallel
y-A\Psi\theta\parallel_{2}^{2}<\epsilon$ (3)

which can be solved by linear programming methods[14] for which efficient
solvers exist. An equivalent solution can be obtained by solving a modified
form of Equation 3

$\operatorname*{arg\,min}_{\theta}\left(\frac{1}{2}\parallel
y-A\Psi\theta\parallel_{2}^{2}+\lambda\parallel\theta\parallel_{1}\right)$ (4)

where $\lambda$ is a parameter to balance the solution between the first term
for integrity to the data and the second term to promote sparsity.

In addition to the sparsity requirement, the above solutions require that the
sensing matrix $A$ be “incoherent” with the signal model $\Psi$. Incoherence
can be understood heuristically in the following way; the incoherence property
suggests that the measurement system $A$ takes a sparse signal $\theta$ and
spreads it out over many observations $y$. If this property holds, information
sufficient to recover the vector $\theta$ are contained in the $M$
measurements $y$. For example, if the sparse signal consists of $k$ spikes we
might choose $\Psi={\bf I}$ (the identity matrix). In this case the Fourier
transform matrix $A$ is optimally incoherent and will spread out information
of the spikes throughout the frequency domain. One could also design a CS
imaging system so that $A$ is comprised of random entries (e.g., each entry
equiprobable “0” or “1”). This essentially guarantees that a vector that is
sparse in the basis $\Psi$ will be non-sparse in $A$.

###  1.2 Related Work

A great many papers in the compressive sensing literature discuss how to
implement practical compressed sensing systems (e.g., see [4, 15, 16]) but few
show CS imaging to be better than conventional imaging when all aspects of the
two approaches are weighed and compared. Most practical CS hardware
architectures can be found described on Igor Carron’s website at
https://sites.google.com/site/igorcarron2/compressedsensinghardware and the
latest information within the CS field may be discovered in his daily blog at
http://nuit-blanche.blogspot.com/.

As discussed above, two credible applications for CS are in its use for MRI[6]
and radio interferometry[17]. We used these wins as prototypes for
establishing the guidelines for successful application of CS. There are a
number of papers in the literature applying CS to detection of stars, such as
the paper by Gupta, et. al.[18]. Missile detection considered later in this
paper bears similarities to star point source detections algorithms. This also
bears some resemblance to the application of compressive sampling to radio
interferometry[17] where the use of partial Fourier reconstruction is
described.

The JASON report on ”Compressive sensing for DoD sensor systems”[3] is closely
related to the topic of this article as the report addresses the question
where CS can be beneficially utilized for DoD radar applications. Some of the
guidelines for successful application of CS can be found in this report but
are not clearly delineated as in this paper. Ultimately the report finds that
CS should be of interest to DoD but more research is required, such as this
paper, to determine useful applications.

The paper “Compressive sampling vs. conventional imaging” by Haupt and
Nowark[19] addresses the same question that is the focus of this paper.
However, that paper does not provide general guidelines nor an example of a
new CS application as is done in this paper. Several recent papers, such as
one by Adcock, et. al.[20], point out the gap between CS theory and use in
real world applications and show how to bridge it.

There are in the CS literature papers related to using compressive sampling
for target detection and tracking. The paper by Kashter, Levi, and Stern[21]
describe application of CS to motion detection. However, their method relies
on frame differencing to achieve sparsity; that is, only the moving objects
appear in the frame differenced image. Similarly the paper by Poon, et.
al.[22] describes a realizable hardware application that they applied to a
simplified target tracking problem. Again their method depends on frame
differencing, hence the objects moving to different pixels from frame to
frame. In the missile detection problem, the sub-pixel missile signal remains
in the same pixel for tens of frames and frame differencing eliminates all
possibility of detecting the already dim target. A paper by Duarte, et.
al.[23] describes the signal detection problem using CS principles without
reconstructing the image. The missile detection application described here
similarly follows this concept; that is, the goal is to detect the missile
regardless of the quality of the image reconstruction.

##  2 A NEW APPLICATION OF COMPRESSIVE SENSING

In this section we describe the conventional and new CS approaches to the
missile detection problem. Specifically we are able to demonstrate better
Receiver Operating Characteristic (ROC) performance with a lower resolution
focal plane array than a previously used, conventional approach. Rather than
try and recover a pristine (low PSNR) image, we show that by modifying the CS
recovery algorithm we can effect greater probability of detection for low
false alarm rates. This is, of course, the main goal of a missile warning
system and hence is a clear demonstration of the potential of CS
architectures.

Section 2.1 briefly describes the sensor system and field data collected to
support the system development. Section 2.2 describes how to adapt the sensor
system and algorithms for a CS approach.

Figure 1: Full 45x 3.7 degree sensor image at 0800.

###  2.1 Missile Detection

Several years ago the Naval Research Laboratory (NRL) developed a MWIR sensor
system for detection, tracking, and declaration of sub- and supersonic sea-
skimming, anti-ship cruise missiles. The objective was to find bright, sub-
pixel sized objects just above the horizon. As part of this effort, both the
sensor and detection algorithms were developed, and subsequently tested with
field data. This section describes the sensor, algorithms, and the collected
field data.

The sensor used in this study is a passive staring midwave IR sensor [24]. It
consists of an f/2 80 mm diameter optical system with a prism anamorphic
element that results in a 3.6 x 48 degree field of view. The focal plane for
the sensor is a 25 micron pitch, 512 x 2560, 5 micron cutoff HgCdTe array
cooled to 80K and operating at a 30 Hz frame rate. The sensor instantaneous
field of view is 130 (vertical) x 300 (horizontal) microradians. Either the
3.7 - 4.25 or 4.6 - 4.85 micron spectral band is chosen depending on viewing
conditions using a warm filter wheel that is part of the optical assembly. The
large focal plane array (FPA) size is required for detecting dim, sub-pixel
targets.

The data processing functions were broken down into two primary subdivisions,
the front-end and back-end processing. This division was required because of
the large amount of data from the sensor that must be handled in real time.
The primary goal of the front end processing is image processing and the
production of candidate detections, which are called exceedances (referred to
as ”xcds” for brevity in this paper). The primary goal of the back end
processing is the generation and maintenance of tracks, and for the sea
skimming missile at the horizon detection mode, correct declaration of threats
with a minimum of false alarms. This paper focuses on a realizable CS camera
and the front end algorithms while the back end processing remains identical.

The detection performance of the system against a surrogate missile target was
tested in a series of experiments in September 2005 [24]. The sensor was
mounted at a height of 65 feet above sea level on the roof of a building 100
meters from the shoreline. This location provided a clear view to the east
over the Atlantic Ocean for the full 45 degree field of view of the sensor.
Data was collected between 0800 and 1500 hours. The target consisted of a
radiometrically calibrated blackbody infrared source constructed using an
internally heated metal plate integrated on a towed decoy. It provided an in-
band signature of approximately 4 watts/ster, which is approximately the
expected signal from a subsonic missile. The target was towed from a Lear jet
flying at an altitude of approximately 2500 feet and speed of 270 knots. The
target altitude was controlled from the tow aircraft and varied between 50 and
150 feet above the sea surface depending on wind conditions. The target can be
seen to be above the geometric horizon at ranges exceeding 15 nautical miles
for target heights of greater than 50 feet Representative sensor imagery is
illustrated in Figure 1. Along with the towed target there were a variety of
small surface vessels, both commercial fishing and pleasure craft, that
presented point like sources to the sensor.

Perhaps the main challenge in the missile detection problem is the large
number of false alarms (FAs) that tend to arise due to clutter in the scene
such as sunlight reflecting off the sea. The front end processing must
minimize these point like sources without losing sight of the target, which is
is relatively dim compared to these bright FAs. To accomplish this goal,
background normalized values are computed by the front end. Specifically, this
two-dimensional spatial image processing consists of:

  1. 1.

Read in an image frame,

  2. 2.

Single frame spatial demeaning (i.e., convolve with a 1-D, 21 element, zero
mean filter),

  3. 3.

Four parallel, optical point spread function (PSF) matched filtering,

  4. 4.

Local neighborhood calculation of the spatial variance,

  5. 5.

Computing background normalized values as the filtered results divided by
spatial variance, and

  6. 6.

Thresholding of the background normalized values in order to identify
exceedances (xcds).

A summary of the original spatial detection algorithm is shown in Figure 2.
For simplicity, this algorithm was slightly modified for the results in this
paper by using only one PSF for the matched filtering. In this paper, applying
this algorithm on the collected field data constitutes the conventional
approach to missile detection. It is important to note that each of the above
listed steps constitutes a linear operation on the data. This will become
important later in the image recovery step of the CS approach.

Figure 2: Block diagram for single frame based spatial processing sea skimming
missile detection algorithm.

Figure 3: Sample of mask for the 128x640 FPA where 15 out of 16 pixels are
blocked (shown in black) and only 1 mask pixel is open (shown in white).

Figure 4: Direct image reconstruction from partial Fourier coefficients where
15 out of 16 pixels are blocked and only 1 mask pixel is open.

###  2.2 CS Missile Detection

The goal of a CS missile detection system is to maintain a high level of
correct true target detections while obtaining only a fraction of the number
measurements relative to the number of pixels in the original system.

Does the sea skimming missile problem fulfill the conditions for finding a
practical application of CS? This particular problem does fit nicely within
the CS requirements. The first requirement for a practical CS application is
when applying CS can overcome hardware, financial, or battlefield limitations.
In the missile detection scenario the trend is towards more resolution and
smaller pixels in order to better detect dim, sub-pixel missile. However, a
four times increase in number of pixels often translates into an order of
magnitude increase in cost. Hence, CS can combat this long term trend for more
expensive FPAs.

Next, one needs a priori knowledge of what you will be imaging in order to
develop the sparsifying and incoherent transforms. In this situation the
objects of interest are only the point sources within the scene. Since Fourier
transforms are incoherent with delta functions, which are point sources in the
discrete realm, they are ideal for this application. Third is the requirement
that the sparsity be orders of magnitude smaller than the number of pixels.
This is true in this case because there is only one missile and only hundreds
of point like pixels out of over a million pixels; hence a reduction of four
orders of magnitude. The fourth requirement is incoherence between measurement
matrix and the sparsifying transform. Here the point like objects are sparse
in the image domain and as mentioned, the Fourier transforms are ideally
incoherent.

The fifth guideline asks “Is the hardware for the CS sensor realizable?” Yes;
for the envisioned hardware one can use Fourier optics and a mask to measure
some fraction of the Fourier coefficients. For proof of concept and
prototyping purposes, we used a simple, random binary mask where each entry is
the result of a Bernoulli trial. As clearly explained in the paper by Lustig,
Donoho, and Paula[6],“Randomness is too important to be left to chance”. The
choice of the mask is critical to obtaining good results because the mask
directly affects the aliasing artifacts. Hence, several masks of this simple
type were tested with our imagery and the best one was used in our experiments
in Section 2.3. Two metrics were used to compare masks in order to pick the
best; one metric is to count the number of FAs and pick the mask that produced
the fewest FAs; the other metric is to compute the Mean Square Error (MSE) of
the image reconstructed with the partial Fourier coefficients and pick the
mask that produced the smallest error. The ROC results for both of these masks
are displayed in Figure 5, where the ”1st mask” is the result of the first
metric and the ”2nd mask” is the result of the second metric above.

Still, this type of mask is likely quite inefficient and there exists
abundance of numerical evidence demonstrating the advantage of variable
density sampling strategies [6, 25], which was not tested in this
demonstration. Furthermore, Romberg has shown [15] that greater efficiencies
are possible with random modulation masks.

For this simulation the scenario consists of a smaller FPA; that is, increase
the size of the pixel by 2 in each direction. However, the mask retains
original resolution/pixel size, which in this case is 512x2560. By creating a
mask that blocks out 3 of 4 pixels, each FPA ”fat” pixel measures the Fourier
coefficient of only the one open mask pixel. The mask is designed to randomly
pick which 3 are blocked and which is open. Similarly we tested reducing the
FPA further - to only 128x640 and 64x320 pixels and the 512x2560 pixel mask
blocks 15 out of 16 or 63 out of 64, respectively. A sample of a section of
the 1 out of 16 mask is shown in Figure 3.

An approach to solving the missile detection problem is to zero fill the
missing Fourier coefficient and reconstruct the image with an inverse Fourier
transform. Figure 4 shows the resulting image from this direct approach of
reconstruction from partial Fourier (PF) coefficients where 15 out of 16
pixels are blocked and only 1 mask pixel is open. Strong artifacts are clearly
visible in the image, which interfere with the processing described above in
section 2.1. However, as shown in the Results section, thresholding the
background normalized image for point like objects works reasonably well when
the SNR of the target is large. This implies that the PF coefficients retain
information on the point like objects in the frame.

Instead, previous successful examples of compressive sensing in sparse MRI [6]
and radio interferometry [17] point to more successful CS approaches for
solving the CS missile detection problem. Following the CLEAN approach[26]
used in radio interferometry and Donoho, Tsaig, and Starck’s stagewise
orthogonal matching pursuit [27], one can implement a greedy method by
iteratively thresholding the background normalized image and removing the
brightest points. In addition, this is equivalent to the nonlinear
thresholding scheme described by Lustig, Donoho, and Pauly[6]. On the other
hand, following the techniques from this same paper on sparse MRI [6] and a
paper describing radio interferometry [17] one can analogously solve the
following optimization problem as an $L_{1}$ relaxation approach:

$\operatorname*{arg\,min}_{x}\left(\parallel T(x)\parallel_{1}+0.5\parallel
y-MFx\parallel_{2}^{2}\right)$ (5)

where $x$ is the unknown image, $y$ is the observed undersampled Fourier
coefficients, $F$ is the Fourier transform, $M$ is a mask for under-sampling
the coefficients, and $T(x)$ is a function that represents the missile
detection algorithm described above in section 2.1. In other words, the
problem constraint has the detection algorithm “built in”. By minimizing
$\parallel T(x)\parallel_{1}$ we are minimizing the number of declared
targets. Given proper thresholding parameters in $T$, the false alarms will be
significantly reduced and only true positives survive the minimization. While
this is not at all appropriate for pristine image reconstruction, that is not
our goal.

There exists many algorithms for solving the $L_{1}$ optimization problem. In
the Results shown in Section 2.3 we used the conjugate gradient (CG) method
with a total variation prior described in the MRI [6] paper and the robust
TwIST method[28].

Table 1: Compares results the original algorithm on the original image to CS in various configurations. The fewer exceedences above the target pixel intensity, the better the algorithm. Image processing | Number of exceedences above target
---|---
Original 512x2560 image algorithm | 90
CS using 256x1280 PF Coeffs via TwIST or thresholding | 10
CS using 128x640 PF Coeffs via TwIST or thresholding | 23
CS using 128x640 PF Coeffs; Remove groups of 5 | 14
CS using 64x320 PF Coeffs via TwIST or thresholding | 32
CS using 64x320 PF Coeffs; Remove groups of 5 | 8

###  2.3 Results

This section presents two tests of this CS missile detection system; results
from a feasibility test to detect the target in a single frame and a more
comprehensive test of missile detection over 500 frames of video.

Initially we performed a first test to determine the feasibility of detecting
the surrogate missile target using CS methodology. This test used a frame
where the missile SNR is relatively large. Table 1 compares results the
original algorithm on the original image (512x2560) to CS in various
configurations. The original algorithm on the original image is able to find
missile target in the background normalized image but, unsurprisingly, there
are 90 false alarms (FAs) that are of stronger intensity than the true target
pixel. In this test, the fewer FAs above the target pixel intensity, the
better the likelihood of target detection and tracking.

If one doubles the width and height of the pixels, one gets 256x1280, or a
quarter of the pixels, which is the case using CS in the next row of Table 1.
In our experiments, using either TwIST for L1 optimization solution or simple
thresholding, all produced the same results, which was 10 FAs that are of
stronger intensity than the true target pixel. That the CS method improved
over the original algorithm with the full image was initially surprising but
later tests showed that this was related to the nature of the image frame used
in this experiment. However, this result does provide the proof we were
seeking that the concept of CS missile detection is valid. The next row of
Table 1 gives the results when the resolution of the FPA is reduced once again
(to 128x640) and here too the result of 23 brighter xcds is better than the
results from the original algorithm. Test of a non-linear thresholding, where
the brightest 5 xcds were iteratively removed as described in Section 2.2, is
shown in the next row. This demonstrates that the non-linear methodology can
improve on simple thresholding. Finally, the last two rows of the Table show
that the CS methods can work for this image with even a lower FPA resolution
of 64x320.

Given a satisfactory result from our proof of concept, a more comprehensive
test of the CS missile detection system was performed by comparing ROC curves
of the target detection over 500 frames of video. The target signal varies
significantly and randomly from frame to frame, offering a challenging test of
the CS technique. The goal is to detect the true target in all 500 frames with
a minimum number of FAs.

The conjugate gradient (CG) optimization [6] weights were both set to 0.0001
to include both the sparsity and total variation (TV) constraints. Figure 5
shows the ROC curve for the original algorithm and another curve for the CS
solution by CG with a total variation (TV) prior at 256x1280 resolution. Each
ROC is computed by processing each frame by its algorithm, applying a set of
thresholds, and recording at each threshold the number of FAs and if the true
target is detected. Afterward, the number of FAs and target detections is
summed over all of the frames and is used to plot the ROC line.

Figure 5 illustrates that our CS method using CG a TV prior at 256x1280
resolution provides a competitive result to the original algorithm. Of
particular interest is the initial behavior where up until 60% of the target
detections the CS method obtains the target with fewer FAs than the original
method. Of note is that the initial test described above used one of these
frames, explaining why in Table 1 the CS method obtained fewer FAs than the
original algorithm.

Finally we compare computation time for the CS method versus the original
algorithm. Running the original algorithm over 500 frames of video using a
Matlab function takes 166 seconds of execution time, which is an order of
magnitude slower than real time (at 30 Hz, 500 frames last 16.7 seconds). This
algorithm does run in real time when implemented on an FPGA. CS using simple
thresholding adds a negligible amount of execution time but the performance of
the thresholding was unsatisfactory. Running the CS method using CG with a TV
prior at 256x1280 resolution executes in 2300 seconds, which is much greater
than real time. Future work might show that utilizing GPU’s and parallel
processing techniques will permit running these codes within the real time
requirements of these types of systems.

Figure 5: ROC comparing the original missile detection algorithm to a CS
algorithm over 500 frames of video containing the surrogate missile target.
The CS method utilized a conjugate gradient method with sparse and TV priors.

##  3 CONCLUSIONS

The purpose of this article is to encourage and assist researchers in finding
those niche applications where the CS approach provides substantial gain over
conventional approaches. Now that the theory has advanced so significantly, it
is time to focus on finding applications where CS will win in real-world
comparisons to current conventional state-of-the-art techniques, which include
techniques such as super-resolution and image deblurring. This paper presented
lessons learned to help researchers find practical, real-world applications,
then used these guidelines to find and show the potential benefit of a new CS
application - sea skimming missile detection.

The results shown here for CS missile detection are only the beginning of the
required effort to demonstrate practical target detection. No effort was taken
to optimize the architecture or mask for the missile detection problem, yet
competitive results were obtained with this simple CS method using a lower
resolution focal plane array. This is proof of the potential of the CS method
for this application and additional research is required on masks, solvers and
architectures that certainly will improve on these initial results.

Finally, the primary message is that all of the excitement surrounding CS is
necessary and appropriate for encouraging our creativity but we all must take
off our ”rose colored glasses” and critically judge our ideas, methods and
results relative to conventional imaging approaches.

## Acknowledgments

The author gratefully thanks Jonathan M. Nichols and James R. Waterman for
their assistance and the Naval Research Laboratory for supporting this work.

## References

  * [1] Donoho, D. L., “Compressed sensing,” IEEE Transactions on Information Theory 52(4), 1289–1306 (2006).
  * [2] Elad, M., “Sparse and redundant representation modeling: What next?,” IEEE Signal Processing Letters 19(12) (2012).
  * [3] Group, J., “Compressive sensing for DoD sensor systems,” http://www.fas.org/irp/agency/dod/jason/compress.pdf (2012).
  * [4] Willett, R. M., Marcia, R. F., and Nichols, J. M., “Compressed sensing for practical optical imaging systems: a tutorial,” Optical Engineering 50(7), 72601–1–13 (2011).
  * [5] Romberg, J., “Imaging via compressive sampling,” Signal Processing Magazine, IEEE 25(2), 14–20 (2008).
  * [6] Lustig, M., Donoho, D., and Pauly, J. M., “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine 58(6), 1182–1195 (2007).
  * [7] Candès, E. J., Romberg, J., and Tao, T., “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory 52(2), 489–509 (2006).
  * [8] Eldar, Y. C. and Kutyniok, G., [Compressed Sensing: Theory and applications ], Cambridge University Press, Cambridge (2012).
  * [9] Elad, M., [Sparse and Redundant Representations: From theory to applications in signal and image processing ], Springer, New York (2010).
  * [10] Mallat, S. G. and Zhang, Z., “Matching pursuits with time-frequency dictionaries,” IEEE Transactions on Signal Processing 41(12), 3397–3415 (1993).
  * [11] A.Tropp, J. and Gilbert, A. C., “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on Information Theory 53(12), 4655–4666 (2007).
  * [12] Needell, D. and Tropp, J. A., “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis 26(3), 301–321 (2009).
  * [13] Dai, W. and Milenkovic, O., “Subspace pursuit for compressive sensing signal reconstruction,” IEEE Transactions on Information Theory 55(5), 2230–2249 (2009).
  * [14] Chen, S. S., Donoho, D. L., and Saunders, M. A., “Atomic decomposition by basis pursuit,” SIAM journal on scientific computing 20(1), 33–61 (1998).
  * [15] Romberg, J., “Compressive sensing by random convolution,” SIAM Journal on Imaging Sciences 2(4), 1098–1128 (2009).
  * [16] Wagadarikar, A., John, R., Willett, R., and Brady, D., “Single disperser design for coded aperture snapshot spectral imaging,” Applied Optics 47(10), B44 –B51 (2008).
  * [17] Li, F., Cornwell, T. J., and de Hoog, F., “The application of compressive sampling to radio astronomy,” Astronomy and Astrophysics 528 (2011).
  * [18] Gupta, R., Indyk, P., Price, E., and Rachlin, Y., “Compressive sensing with local geometric features,” in [Proceedings of the 27th annual ACM symposium on Computational geometry ], 87–96, ACM (2011).
  * [19] Haupt, J. and Nowak, R., “Compressive sampling vs. conventional imaging,” in [Image Processing, 2006 IEEE International Conference on ], 1269–1272, IEEE (2006).
  * [20] Adcock, B., Hansen, A. C., Poon, C., and Roman, B., “Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing,” arXiv preprint arXiv:1302.0561 (2013).
  * [21] Kashter, Y., Levi, O., and Stern, A., “Optical compressive change and motion detection,” Applied Optics 51(13), 2491–2496 (2012).
  * [22] Poon, P. K., Townsend, D., Wehrwein, S., Vera, E. M., Stenner, M. D., and Gehm, M. E., “A static architecture for compressive target tracking,” in [Proceedings of 2012 SPIE Defense Security and Sensing Symposium ], 8360, 83600G–1 (2012).
  * [23] Duarte, M. F., Davenport, M. A., Wakin, M. B., and Baraniuk, R. G., “Sparse signal detection from incoherent projections,” in [Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on ], 3, III–III, IEEE (2006).
  * [24] Smith, L., Waterman, J., and Smith, R., “Demonstration of novel detection algorithms for the navy das irst system with 2005 shore-based test data,” in [Military Sensing Symposium 2006 ], (2006).
  * [25] Puy, G., Vandergheynst, P., and Wiaux, Y., “On variable density compressive sampling,” IEEE Signal Processing Letters 18(10), 595–598 (2011).
  * [26] Hògbom, J. A., “Aperture synthesis with a non-regular distribution of interferometer baselines,” Astronomy and Astrophysics Supplement 15, 417 (1974).
  * [27] Donoho, D. L., Tsaig, Y., Drori, I., and Starck, J.-L., “Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit,” IEEE Transactions on Information Theory 58(2), 1094–1121 (2012).
  * [28] Bioucas-Dias, J. M. and Figueiredo, M. A., “A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Transactions on Image Processing 16(2), 2992–3004 (2007).


# Enumerating consistent subgraphs of directed acyclic graphs: an insight into
biomedical ontologies

Yisu Peng111Contributed equally to this work. Department of Computer Science,
Indiana University, Bloomington, Indiana, U.S.A.  Yuxiang Jiang11footnotemark:
1 Department of Computer Science,
Indiana University, Bloomington, Indiana, U.S.A.  Predrag
Radivojac222Corresponding author, email<EMAIL_ADDRESS>Department of
Computer Science,
Indiana University, Bloomington, Indiana, U.S.A.

###### Abstract

Modern problems of concept annotation associate an object of interest (gene,
individual, text document) with a set of interrelated textual descriptors
(functions, diseases, topics), often organized in concept hierarchies or
ontologies. Most ontologies can be seen as directed acyclic graphs, where
nodes represent concepts and edges represent relational ties between these
concepts. Given an ontology graph, each object can only be annotated by a
consistent subgraph; that is, a subgraph such that if an object is annotated
by a particular concept, it must also be annotated by all other concepts that
generalize it. Ontologies therefore provide a compact representation of a
large space of possible consistent subgraphs; however, until now we have not
been aware of a practical algorithm that can enumerate such annotation spaces
for a given ontology. In this work we propose an algorithm for enumerating
consistent subgraphs of directed acyclic graphs. The algorithm recursively
partitions the graph into strictly smaller graphs until the resulting graph
becomes a rooted tree (forest), for which a linear-time solution is computed.
It then combines the tallies from graphs created in the recursion to obtain
the final count. We prove the correctness of this algorithm and then apply it
to characterize four major biomedical ontologies. We believe this work
provides valuable insights into concept annotation spaces and predictability
of ontological annotation.

##  1 Introduction

Ontologies have become a common means of concept annotation in computational
biology and related fields [18]. A protein’s molecular function [1], an effect
of a genetic variant [28], or a patient’s diagnosis [19] are typical examples
in which biomedical entities such as macromolecules, mutations, or individuals
are associated with sets of mutually dependent descriptors. The dependencies
between these descriptors are often hierarchical, leading to the use of
directed acyclic graphs as concept space representations.

A directed acyclic graph is a pair $(V,E)$, where $V$ is a set of vertices
(nodes) and $E$ is a set of directed edges (links) between vertices such that
no cycles can be formed. Each vertex in the graph is associated with a unique
concept (term, description) and each edge is associated with a particular type
of relational tie. For example, when annotating proteins as biomedical
entities using the Gene Ontology graph [1], the terms “nucleic acid binding”
and “DNA binding” are linked by edges of the type $\mathsf{is}$-$\mathsf{a}$
asserting that DNA binding is a more specific form of nucleic acid binding.
Other types of relational ties include $\mathsf{part}$-$\mathsf{of}$,
$\mathsf{regulates}$, and so on.

A typical biomedical entity is associated with a set of terms determined
through experiment such as a molecular assay or a diagnostic procedure. A
protein, for example, may be assigned terms “DNA binding” and “RNA binding”,
neither of which is a generalization of the other. To avoid annotation
inconsistencies, this protein must also be annotated by the terms such as
“nucleic acid binding” and all others that generalize either of the
experimentally determined terms. More broadly, this implies that a biomedical
object can only be annotated by a set of terms that respect the hierarchy – a
_consistent subgraph_ of the ontology. Unfortunately, (manual) experimental
annotation is resource-demanding and often incomplete [16], giving rise to an
entire field of computational prediction [17, 11].

The development of computational prediction methods presents its own
challenges. Although it can be performed by building a separate binary
classifier for each concept in the ontology, this approach is currently
competitive only for specialized ranking tasks; e.g., disease-gene
prioritization [14], since it does not exploit relationships between the
terms. On the other hand, a more complete characterization is via learning
structured outputs [25] in which a method takes an object (e.g., a protein)
and is asked to provide the totality of concepts with which this object might
be associated (i.e., a consistent subgraph). However, the structured-output
formulation generally falls under the extreme classification umbrella because
the size of the output space is often exceedingly large. This poses problems
in measuring similarity between annotations, evaluating accuracy of
classification models, and optimization when solving the “argmax problem” [4,
5, 12].

We identify now what we believe is an open problem in computational biology
and computer science; that is, efficiently determining the exact number of
consistent subgraphs in a given ontology. This problem has a linear-time
solution for rooted trees [22], but to our knowledge no such algorithm exists
for directed acyclic graphs. This paper therefore proposes a practical
solution to this enumeration problem, proves its correctness, analyzes run-
time complexity, and introduces various computational speedups. Using this new
approach, we analyze four often-used ontologies from the biomedical domain and
explore the space of possible annotations. We believe that the algorithms,
software, and analysis carried out in this work will lead to better insights
into concept annotation spaces and facilitate ontology quality assurance.

##  2 A Motivating Example

A growing number of concept annotation problems are formulated as the manual
or computational assignment of a set of mutually related textual descriptors
to some objects of interest. One of such problems is the computational
prediction of protein function [5], which can be broadly operationalized as
follows:

> _Given_ : (1) an amino acid sequence with auxiliary data such as structure,
> expression, interactions, etc. of a protein $p$ with unknown or incomplete
> function; (2) training data that includes sequences, structures, or systems
> data corresponding to a (large) set of proteins, some of which have their
> true biological functions available; (3) a Gene Ontology (GO); i.e., a
> concept hierarchy used to represent biological functions of proteins in a
> structured and easy-to-compute-on form.
>
> _Objective_ : provide a set of GO terms that are most likely to be the true
> (experimental) annotation of $p$.

The objects of interest here are proteins and the set of textual descriptors
of protein function is given by GO – an ontology with a directed acyclic graph
structure where each node represents a textual descriptor and each edge
represents a particular type of a relational tie between two descriptors [1].

Figure 1: The functional annotation of the friend leukemia integration 1
transcription factor isoform 1 [_FLI1_ ; Homo sapiens] (RefSeq ID:
NP_002008.2) as a consistent subgraph of the molecular function ontology. The
arrows in this graph indicate an is-a relationship and are drawn in the
reverse direction.

An example of such an annotation is shown in Figure 1, where 8 terms from the
molecular function domain have been assigned to this protein. Due to the
hierarchical organization of GO, both the set of experimentally determined
terms and the set of computationally predicted terms must respect this
hierarchy. As shown in this example, the annotation of the term “DNA binding”,
implies the annotation of all the other GO terms that conceptually generalize
it; e.g., “nucleic acid binding”, “binding”, etc. Typically, the ontology used
to represent the annotation space of proteins contains thousands to tens of
thousands of terms, whereas the true annotation of a protein consists of tens
to at most hundreds of terms. Because the task of a prediction algorithm is to
find the most likely annotation, it must devise an efficient procedure to
search through the space of _all_ possible annotations.

Most biomedical ontologies have grown over the years to contain a large number
of terms. Computationally selecting one such “winning” annotation; i.e., a set
of terms, or even providing a short list of most likely annotations, is a
significant challenge [12, 25]. This prediction problem thus belongs to a so-
called extreme classification scenario because the number of possible
(discrete) annotations the algorithm must consider is astronomically large. In
fact, we noticed that it is not even possible to give an exact number of
possible annotations for a protein. Therefore, an answer to such a simple
question (“What is the the number of possible GO annotations a protein can be
assigned?”) requires the development of a practical counting algorithm. The
resulting counts can, in turn, give insight into the nature and the difficulty
of the computational function annotation of biological
macromolecules.333Reasonable approximations can be provided by calculating the
lower and upper bounds, as we have done later in Section 6. Neither of those,
however, provides a full intellectual satisfaction when an exact count can
actually be computed.

It is important to mention that the annotation of biological macromolecules is
one of the most interesting examples of concept annotation, primarily because
of its biomedical significance but also because of the sizes of the available
ontologies. Similar situations, however, arise beyond computational biology,
as in the fields of text mining [8] and computer vision [15].

##  3 Preliminaries

###  3.1 Basic Concepts and Notation

Let $\mathcal{G}=(V,E)$ be a directed graph, where $V$ is a set of vertices
representing concepts and $E\subseteq V\times V$ is a collection of ordered
pairs $(u,v)$ representing directional relationships, $u\to v$, between two
concepts. A sequence of vertices $u_{1},u_{2},\dots,u_{k}$ is called a walk if
$(u_{i},u_{i+1})\in E$ for $i=1,2,\dots,k-1$. A walk of distinct vertices
except for the identical starting and ending vertices is called a cycle. A
directed graph that does not contain cycles is referred to as directed acyclic
graph (DAG).

Given two vertices $u,v\in V$ in a DAG, $u$ is said to be an ancestor of $v$
and $v$ is said to be a descendant of $u$ if there exists a walk from $u$ to
$v$. We denote a set of all ancestors of $v$ as $\mathcal{A}({v})$ and a set
of all descendants of $u$ as $\mathcal{D}({u})$. We next define
$\mathcal{A}^{+}({v})=\\{v\\}\cup\mathcal{A}({v})$ as the set of extended
ancestors of $v$ and $\mathcal{D}^{+}({u})=\\{u\\}\cup\mathcal{D}({u})$ as the
set of extended descendants of $u$. Finally, if $(u,v)\in E$, the vertex $u$
is said to be a parent of $v$, whereas $v$ is said to be a child of $u$. We
denote the set of all parents of $v$ as $\mathcal{P}(v)$ and the set of all
children of $u$ as $\mathcal{C}(u)$.

###  3.2 Transitivity of Relational Ties

When an object is annotated with ontological concepts, it is often considered
that all ancestors of those annotated concepts should be automatically
assigned to the object. For example, annotating the function of a protein with
“enzyme binding” also implicitly annotates it with “protein binding”,
“binding” and, finally, the root term “molecular function”. This type of
reasoning requires all involved relationships between concepts to be
transitive.

Biomedical ontologies, however, usually contain various types of relationships
between concepts, some of which are not transitive. Therefore, we only
consider $\mathsf{is}$-$\mathsf{a}$ and $\mathsf{part}$-$\mathsf{of}$
relationships, both of which maintain transitivity and permit reasoning about
ancestral concepts. It is also worth noting that we define the direction of
edges to be pointing from the general terms to specific so that the depth of a
node aligns with the increasing resolution of the descriptors. We show in
Section 5.2 that the directionality of edges has no impact on the total count.
Throughout this work, we consider an ontology $\mathcal{O}=(V,E)$ to be a DAG,
where edges represent transitive relationships.

###  3.3 Consistent Subgraphs

Let $\mathcal{O}=(V,E)$ be an ontology and $S\subseteq V$ a set of vertices. A
subgraph $(S,E^{S})$ is said to be induced from the original graph
$\mathcal{O}$ by $S$ if $E^{S}$ is the largest subset of pairs $(u,v)$ from
$E$ such that both $u,v\in S$. We denote such vertex-induced subgraph as
${\mathcal{O}}[S]$. We also use $\mathcal{O}[-S]=(V-S,E^{V-S})$ to denote the
subgraph induced by vertices other than $S$.

######  Definition 3.1.

A subgraph ${\mathcal{O}}[S]=(S,E^{S})$ with respect to the original graph
$\mathcal{O}=(V,E)$ is called _consistent_ if $\forall v\in S$, $(u,v)\in
E\implies u\in S$.

##  4 Basic Algorithms

###  4.1 Problem Specification

Given an ontology $\mathcal{O}=(V,E)$, our goal is to develop a practical
algorithm that enumerates all consistent subgraphs of $\mathcal{O}$. We allow
the graph to have more than a single root (a vertex with no incoming edges) as
well as to be disconnected.

An example of the enumeration problem is shown in Figure 2. We generally
observe that the number of consistent subgraphs is bounded from below by
$2^{\ell}$, where $\ell$ is the total number of leaf vertices (those with no
outgoing edges), and from above by $2^{|V|}$. The structure of the graph,
however, determines the exact count and its proximity to either of the bounds.
If the input graph is a chain of $|V|$ vertices ($\ell=1$), the total number
of consistent subgraphs equals $|V|+1$. On the other hand, if the original
graph is a set of $|V|$ disconnected vertices ($\ell=|V|$), there are
$2^{|V|}=2^{\ell}$ consistent subgraphs. This analysis suggests that
enumerating consistent subgraphs has a straightforward intractable solution of
listing all $2^{|V|}$ vertex-induced subgraphs of the ontology and checking
for the consistency of each such subgraph.

Figure 2: Consistent subgraphs of an ontology $\mathcal{O}=(V,E)$ with $|V|=7$
vertices and $|E|=7$ edges. Observe that the reversal of all edges in the
graph would lead to a reversed graph with the same number of consistent
subgraphs (white vertices; Theorem 5.1).

We use $\mathtt{cdag}(\mathcal{O})$ to denote the desired function that takes
a directed acyclic graph $\mathcal{O}$ as input and returns the number of
consistent subgraphs in that graph. We use $\mathtt{ctree}(\mathcal{T})$ and
$\mathtt{cforest}(\mathcal{F})$ for the special cases where the input graph is
a rooted tree $\mathcal{T}$ or a forest $\mathcal{F}$, respectively.

###  4.2 Counting Consistent Subgraphs in Trees

We first discuss a special case where the input graph is a rooted tree; that
is, when each non-root vertex has a single parent. In this case, there exists
a linear algorithm in the number of vertices; see Lemma 1 in [22]. We provide
this solution in Algorithm 1 with a minor modification resulting from the fact
that our algorithm includes an empty tree in the total count. This algorithm
naturally extends to collections of rooted trees. One can enumerate subtrees
for each tree and take the product as the total count. We refer to this
extended algorithm as $\mathtt{cforest}$ (not shown).

Input : A tree $\mathcal{T}_{r}$, rooted at $r$.

Output : The number of consistent subgraphs in $\mathcal{T}_{r}$.

1 Function  _$\mathtt{ctree}(\mathcal{T}_{r})$_

2 if  _$\mathcal{T}_{r}$ is empty_ then

3 return 1;

4

5 else

6 return $1+\prod_{u\in\mathcal{C}(r)}\mathtt{ctree}(\mathcal{T}_{u})$;

7

8 end if

9

10 end

Algorithm 1 Counting the number of consistent subgraphs in rooted trees [22].

Algorithm 1 recursively traverses a tree in a pre-order manner. For any
subtree rooted at vertex $v$, the number of consistent subtrees that contain
$v$ equals the product of all subcounts from its subtrees rooted at each
child. Additionally, we add $1$ for the only consistent subtree that does not
contain $v$; i.e., the empty tree. The recursion terminates at the empty tree
whose count is one.

Input : A directed acyclic graph $\mathcal{O}$.

Output : The number of consistent subgraphs in $\mathcal{O}$.

1 Function  _$\mathtt{cdag}(\mathcal{O})$_

2 if  _$\mathcal{O}$ is a forest_ then

3 return $\mathtt{cforest}(\mathcal{O})$;

4

5 else

6 Pick any vertex $u$ as the pivot;

7 return
$\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({u})])+\mathtt{cdag}({\mathcal{O}}[-\mathcal{A}^{+}({u})])$;

8

9 end if

10

11 end

Algorithm 2 Counting the number of consistent subgraphs in directed acyclic
graphs.

###  4.3 Counting Consistent Subgraphs in Directed Acyclic Graphs

Directed acyclic graphs generalize trees in that they allow for multi-parent
vertices. Such vertices, however, break Algorithm 1 because the recursive
branches are no longer independent. Algorithm 2 circumvents this problem by
recursively decomposing a graph into two strictly smaller subgraphs according
to a selected _pivot_ vertex. We will show in the next section that the number
of consistent subgraphs in the two smaller graphs add up to be the number for
the original graph (Line 2, Algorithm 2). The algorithm continues recursive
enumeration until the graph becomes a forest, in which case it calls
$\mathtt{cforest}$. Figure 3 illustrates the process of graph decomposition
with respect to the pivot vertex $u$. We note that any vertex can serve as
pivot and will discuss the selection of pivots and how they impact the run
time in Sections 5.3 and 6.1.

Figure 3: Illustration of graph decomposition. The enumeration problem of the
original graph from panel (a) is split into two subproblems based on the pivot
vertex $u$; shown in panels (b) and (c). The count in (b) corresponds to the
number of consistent subgraphs in (a) that do not include $u$, while the count
in (c) corresponds to the count of consistent subgraphs in (a) that include
$u$. In panel (a), the set of descendants of $u$ is shaded in orange and the
set of ancestors is shaded in blue.

###  4.4 Correctness and Complexity of the Algorithm

We first observe that the size of the problem in the number of vertices is
guaranteed to decrease during recursive calls, thus ensuring that the
algorithm terminates after a finite number of iterations. Next, we justify the
equation corresponding to the Line 2 in Algorithm 2,

$\displaystyle\mathtt{cdag}(\mathcal{O})=\mathtt{cdag}\left({\mathcal{O}}[-\mathcal{D}^{+}({u})]\right)+\mathtt{cdag}\left({\mathcal{O}}[-\mathcal{A}^{+}({u})]\right).$
(1)

######  Lemma 4.1.

Let $\mathtt{cdag}(\mathcal{O}|\neg u)$ be the number of consistent subgraphs
in $\mathcal{O}$ that do not contain $u$. We have
$\mathtt{cdag}(\mathcal{O}|\neg
u)=\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({u})])$.

###### Proof.

The equal cardinality of the two sets of consistent subgraphs is demonstrated
by showing that both sets are contained in each other. For any $S\subseteq
V-\mathcal{D}^{+}({u})$ that induces a consistent subgraph of
${\mathcal{O}}[-\mathcal{D}^{+}({u})]$, it also induces a unique consistent
subgraph in $\mathcal{O}$. Also, since none of them contains $u$, we have
$\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({u})])\leq\mathtt{cdag}(\mathcal{O}|\neg
u)$. Conversely, for any consistent subgraph induced by $S$ such that $u\notin
S$, we have $\forall v\in\mathcal{D}^{+}({u}),v\notin S$ by the definition of
consistency. Therefore, $S$ also induces a consistent subgraph in
${\mathcal{O}}[-\mathcal{D}^{+}({u})]$. That is,
$\mathtt{cdag}(\mathcal{O}|\neg
u)\leq\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({u})])$.

∎

######  Lemma 4.2.

Let $\mathtt{cdag}(\mathcal{O}|u)$ be the number of consistent subgraphs in
$\mathcal{O}$ that contain $u$. We have
$\mathtt{cdag}(\mathcal{O}|u)=\mathtt{cdag}({\mathcal{O}}[-\mathcal{A}^{+}({u})])$.

###### Proof.

As in Lemma 4.1, for any $S\subseteq V-\mathcal{A}^{+}({u})$ that induces a
consistent subgraph of ${\mathcal{O}}[-\mathcal{A}^{+}({u})]$,
$S\cup\mathcal{A}^{+}({u})$ also induces a unique consistent subgraph (that
contains $u$) in the original graph. That is,
$\mathtt{cdag}({\mathcal{O}}[-\mathcal{A}^{+}({u})])\leq\mathtt{cdag}(\mathcal{O}|u)$.
Also, for any consistent subgraph induced by $S$ and $u\in S$, we have
$\mathcal{A}^{+}({u})\subseteq S$ by the definition of consistency. Note that
the uniqueness of $S$ implies the uniqueness of $S-\mathcal{A}^{+}({u})$. We
can see that the subgraph induced by $S-\mathcal{A}^{+}({u})$ in
${\mathcal{O}}[-\mathcal{A}^{+}({u})]$ is consistent. Given $\forall w\in
S-\mathcal{A}^{+}({u})$, and $(v,w)$ being an edge in
${\mathcal{O}}[-\mathcal{A}^{+}({u})]$, we must have $v\in
S-\mathcal{A}^{+}({u})$ as well, due to the consistency of ${\mathcal{O}}[S]$
with respect to the original graph. That is,
$\mathtt{cdag}(\mathcal{O}|u)\leq\mathtt{cdag}({\mathcal{O}}[-\mathcal{A}^{+}({u})])$.

∎

######  Theorem 4.1.

Given an ontology $\mathcal{O}=(V,E)$ and any $u\in V$, the number of
consistent subgraphs in $\mathcal{O}$ equals the sum of the numbers of
consistent subgraphs in ${\mathcal{O}}[-\mathcal{D}^{+}({u})]$ and
${\mathcal{O}}[-\mathcal{A}^{+}({u})]$.

###### Proof.

Equation 1 holds by combining Lemmas 4.1 and 4.2. ∎

To analyze complexity of the algorithm, let $n$ be the number of vertices in
the graph and $m$ be the number of multi-parent vertices. Assuming a multi-
parent vertex is always selected as pivot, we can express the run time
complexity $T(n)$ via the following recurrence

$T(n)\leq T(n-1)+T(n-3)+f(n),$

where $f(n)$ incorporates the time to select the pivot, split the graph and
add two large integers. Let us further assume that the larger of the two
graphs after decomposition contains $n-n/k$ elements, where $2\leq k\leq n$.
It is now straightforward to show that $T(n)=O(f(n)2^{\min(m,s(k))})$, where
$s(k)=O(n)$ if $k=O(n)$ and $s(k)=O(\log n)$ if $k=O(1)$.

We can now see that the algorithm is exponential in the worst case; however,
it reduces to a polynomial algorithm when $m=O(\log n)$ or when $k=O(1)$.
Assuming linear time to conduct graph decomposition and a constant time for
addition/multiplication, we obtain $T(n)=O(n^{2})$.

##  5 Advanced Algorithms

The run-time of the algorithm heavily depends on the structure of the ontology
and the selection of pivots. Here we discuss several practical considerations
aimed at accelerating Algorithm 2. Once we conclude this discussion, the full
method will be presented in Algorithm 3 (Section 6).

###  5.1 Pruning Branching Components

It is easy to observe that when the ontology consists of multiple connected
components, these components can be independently and, if needed,
simultaneously processed. We take this reasoning a step further to consider a
special scenario of nearly disconnected graphs where (i) the two components
are connected via a single vertex and (ii) all vertices in one component are
descendants of this vertex.

######  Definition 5.1.

Given a graph $\mathcal{O}=(V,E)$ and $u\in V$,
${\mathcal{O}}[\mathcal{D}({u})]$ is called a _branching component_ if
$\forall v\in V-\mathcal{D}^{+}({u})$ and $\forall
w\in\mathcal{D}({u}),(v,w)\notin E$. Vertex $u$ is called a _branching
vertex_.

Figure 4: Illustration of branching components. Panel (a) shows a branching
vertex $u$ that separates the graph into a stem component $\mathcal{O}^{st}$
and a branching component $\mathcal{O}^{br}$. The collection of edges from $u$
to $\mathcal{O}^{br}$ is replaced by a zigzag arrow. Panel (b) shows a
component-wise tree structure.

Figure 4a gives an example in which $u$ is a branching vertex, since the
removal of $u$ disconnects $\mathcal{D}({u})$ (i.e., the branching component,
$\mathcal{O}^{br}$) from the rest of the graph. We refer to the remaining part
of the graph as the stem component, $\mathcal{O}^{st}$. More generally, Figure
4b shows a graph with a component-wise tree structure, where branching
vertices serve as hinges of branching component to their corresponding stems.
We will use $(\mathcal{O}^{st},\mathcal{O}^{br},u)$ to denote the desired
structure.

Given $(\mathcal{O}^{st},\mathcal{O}^{br},u)$, we demonstrate that
$\mathtt{cdag}({\mathcal{O}})$ can be decoupled into two sequential
subproblems: (i) $\mathtt{cdag}(\mathcal{O}^{br})$ and (ii)
$\mathtt{cdag}(\mathcal{O}^{st})$. We use $\varphi(u)$ for the subtotal of
consistent subgraphs in the branching component $\mathcal{O}^{br}$. We also
notice that the entire branching component can be pruned once $\varphi(u)$ is
computed, making $u$ a leaf vertex in $\mathcal{O}^{st}$. Therefore, we modify
the algorithm so as to allow a subtotal count $\varphi(u)$ for every vertex as
if a branching component has been pruned from $u$. Notice that $\varphi(u)=1$
for all intermediate vertices and original leaves.

With the introduction of $\varphi(u)$, the recursive equation in Algorithm 1
becomes

$\displaystyle\mathtt{ctree}(\mathcal{T}_{r})=\varphi(r)+\prod_{u\in\mathcal{C}(r)}\mathtt{ctree}(\mathcal{T}_{u}).$
(2)

Similarly, Equation 1; i.e., Line 2 in Algorithm 2, must be modified to

$\displaystyle\mathtt{cdag}(\mathcal{O}^{st})=\mathtt{cdag}\left({\mathcal{O}^{st}}[-\mathcal{D}^{+}({u})]\right)+\varphi(u)\cdot\mathtt{cdag}\left({\mathcal{O}^{st}}[-\mathcal{A}^{+}({u})]\right),$
(3)

where $\varphi(u)$ accounts for the fact that for any consistent subgraph
$S_{i}$ in the pruned $\mathcal{O}^{br}$ and any consistent subgraph $S_{j}$
in ${\mathcal{O}^{st}}[-\mathcal{A}^{+}({u})]$, ${\mathcal{O}}[S_{i}\cup
S_{j}\cup\mathcal{A}^{+}({u})]$ is a distinct consistent subgraph in
$\mathcal{O}$. The approach naturally extends to multiple (hierachical)
branching components such that we compute the subtotal of consistent subgraphs
within each component and agglomerate them in a reversed topological order.

The pruning operation is preferred before each instance of decomposition for
two main reasons: (i) it divides the problem into smaller non-overlapping
subproblems, while a direct decomposition usually results in substantial
overlapping subproblems; (ii) although a full parallelization over components
is restricted since stem components have to be computed only after all of
their branching components are finished, the unordered components can be
computed simultaneously. For example, as in Figure 4b, $\mathcal{O}^{br}_{1}$,
$\mathcal{O}^{br}_{3}$, $\mathcal{O}^{br}_{4}$ and $\mathcal{O}^{br}_{5}$ can
be computed in parallel.

###  5.2 Reverse Graphs

Let $\mathcal{O}^{R}=(V,E^{R})$ be the reverse graph of $\mathcal{O}$, where
$E^{R}=\\{(u,v)|(v,u)\in E\\}$. We show that the number of consistent
subgraphs in $\mathcal{O}$ equals that in $\mathcal{O}^{R}$.

######  Lemma 5.1.

If ${\mathcal{O}}[S]$ is a consistent subgraph of $\mathcal{O}$,
${\mathcal{O}^{R}}[-S]$ is a consistent subgraph of $\mathcal{O}^{R}$.

###### Proof.

We prove this Lemma by contradiction. For $\forall u\in V-S$ and $\forall
v\in\mathcal{A}^{R}({u})$,444We use $\mathcal{A}^{R}({u})$ and
$\mathcal{D}^{R}({u})$ for ancestors and descendants of $u$ in
$\mathcal{O}^{R}$; $\mathcal{A}^{R}({u})=\mathcal{D}({u})$ and
$\mathcal{D}^{R}({u})=\mathcal{A}({u})$. if $v\notin V-S$, then
$u\in\mathcal{A}({v})\subseteq S$ due to the consistency of
${\mathcal{O}}[S]$. This contradicts $u\in V-S$. Therefore, the assumption
$v\notin V-S$ is false and we have $\forall u\in V-S$, $\forall
v\in\mathcal{A}^{R}({u})\subseteq V-S$. That is, ${\mathcal{O}^{R}}[-S]$ is
consistent. ∎

This Lemma demonstrates that all complementary white vertices in Figure 2 form
consistent subgraphs in the reverse graph.

######  Theorem 5.1.

Given an ontology $\mathcal{O}$,
$\mathtt{cdag}(\mathcal{O})=\mathtt{cdag}(\mathcal{O}^{R})$.

###### Proof.

Given Lemma 5.1, we see that the mapping
$f({\mathcal{O}}[S])={\mathcal{O}^{R}}[-S]$ is a _bijection_ between the two
sets of consistent subgraphs. Therefore, the two sets are of equal
cardinality. ∎

Theorem 5.1 permits graph reversal at any point during the algorithm depending
on which of the graphs is more likely to terminate first. For example, we can
always choose the one with fewer multi-parent vertices so as to greedily
reduce the upper bound of recursive calls. It is worth noting that all the
leaves become roots in the reverse graph. Therefore, in the final algorithm
that incorporates both pruning and reversing modules, we generalize the
algorithm to allow for $\varphi>1$ on roots (branching vertices in the reverse
sense) in order to ensure compatibility.

Having $\varphi(r)>1$ on a root indicates that all the ancestors of $r$ have
been pruned out. For trees (after pruning), we have
${\mathcal{O}}[-\mathcal{D}^{+}({r})]={\mathcal{O}}[\mathcal{A}({r})]$. With
Lemma 4.1 and Theorem 5.1, we have

$\displaystyle\mathtt{cdag}(\mathcal{O}|\neg
r)=\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({r})])=\mathtt{cdag}({\mathcal{O}}[\mathcal{A}({r})])=\mathtt{cdag}({\mathcal{O}^{R}}[\mathcal{D}^{R}(r)])=\varphi(r).$

On the other hand, for any consistent subgraph $S$ containing $r$,
$S-\mathcal{A}^{+}({r})$ induces a consistent subgraph in
${\mathcal{O}}[\mathcal{D}({r})]$ and vice versa; thus,

$\displaystyle\mathtt{cdag}(\mathcal{O}|r)=\prod_{u\in\mathcal{C}(r)}\mathtt{ctree}(\mathcal{T}_{u}).$

Hence, these two subtotals sum to be the total count and Equation 2 remains
unchanged. However, if a root $r$ with $\varphi(r)>1$ is selected to be the
pivot, we have the following equation according to Theorem 5.1 and Equation 3,

$\displaystyle\mathtt{cdag}(\mathcal{O})=\mathtt{cdag}(\mathcal{O}^{R})$
$\displaystyle=\mathtt{cdag}({\mathcal{O}^{R}}[-\mathcal{D}^{R+}(r)])+\varphi(r)\cdot\mathtt{cdag}({\mathcal{O}^{R}}[-\mathcal{A}^{R+}(r)])$
$\displaystyle=\mathtt{cdag}({\mathcal{O}}[-\mathcal{A}^{+}({r})])+\varphi(r)\cdot\mathtt{cdag}({\mathcal{O}}[-\mathcal{D}^{+}({r})]),$

whereas Equation 3 remains unchanged for non-root vertices.

###  5.3 Pivot Selection

As alluded to before, the selection of vertices used for partitioning has the
potential to significantly change the computation time. It is therefore
reasonable to devise a strategy for pivot selection. Besides a random
selection of multi-parent vertices (mpv’s), which aims at directly converting
DAGs into trees one step at a time, we also consider three other pivot
heuristics. The first strategy is to pick a vertex with the maximum degree,
with random selection in case of ties, because decomposing the graph according
to such vertices may increase the chance of having either disconnected
components or branching components. The second strategy selects the pivot so
as to minimize $e-n+r$ over the two subproblems, where $e$, $n$, and $r$ are
the number of edges, vertices, and roots in the two components. We refer to
this quantity as “bound” since it is an upper bound of the number of mpv’s in
the graph (see Supplementary Materials for the proof). Note that it is closely
related to the cyclomatic number of the graph. Finally, the third strategy
simulates a unit network flow for all vertices running in the direction from
leaves to the roots and selects the “bottleneck” vertex; i.e., the one that
maximizes the ratio of the flow in the vertex and the number of its
descendants (see Supplementary Materials for this pivot selection algorithm).
These strategies will be empirically compared in Section 6.

###  5.4 Hashing

It can occur during the recursive procedure that certain subgraphs require
repeated enumeration. In Figure 3, for example, the subgraph $h$-$i$-$j$ is
present in both subproblems shown in Figures 3b-c. Computing the count for
this subgraph would emerge in the Figure 3b subproblem if the ensuing
decomposition were based on vertex $d$, although it would not emerge if the
partitioning were based on vertex $j$. Interestingly, the subgraph $k$-$l$
would be counted twice in the Figure 3b subproblem; i.e., when both
$\mathcal{A}^{+}({d})$ and $\mathcal{D}^{+}({d})$ are removed, and it would
then appear one more time in the Figure 3c subproblem.

To avoid repeated enumeration, whenever a solution to a subproblem is
obtained, the count for this subproblem is stored. Then, during the recursive
calls, we first check if the result is already available before further
calculation. To hash a result, we use the sorted IDs of all vertices in the
subgraph as a key. Obviously, this key is unique because it corresponds to a
vertex-induced subgraph of $\mathcal{O}$. For the pruned subgraph, we store
the key of the subgraph along with the branching vertex. Whenever the ID of
the branching vertex is used to generate a key, the stored key of the
corresponding subgraph is appended to the vertex’s ID with parentheses around
it.

Input : A directed acyclic graph $\mathcal{O}$

Output : The number of consistent subgraphs in $\mathcal{O}$.

1 Function  _$\mathtt{cdag}^{\ast}(\mathcal{O})$_

2 $(count,succeed)\leftarrow\texttt{lookup\\_hash\\_table}(\mathcal{O})$;

3 if  _$succeed$_ then

4 return $count$;

5

6 end if

7 if  _$m(\mathcal{O}^{R}) <m(\mathcal{O})$_ then // check the number of
multi-parent vertices

8 $\mathcal{O}\leftarrow$ reverse($\mathcal{O}$);

9

10 end if

11 $count\leftarrow 1$;

12 foreach  _connected component $\mathcal{O}_{i}$_ do

13 if  _$\mathcal{O}_{i}$ is a tree_ then

14 $count_{i}\leftarrow\mathtt{ctree}^{\ast}(\mathcal{O}_{i})$;

15

16 else

17 $\\{u_{j}\\}\leftarrow\texttt{branching\\_vertices}(\mathcal{O}_{i})$;

18 foreach  _$u_{j}$ in a reversed topological order_ do

19
$\varphi(u_{j})\leftarrow\mathtt{cdag}^{\ast}({\mathcal{O}_{i}}[\mathcal{D}({u_{j}})])$;

20 $\texttt{prune}(\mathcal{D}({u_{j}}))$;

21

22 end foreach

23 $u\leftarrow\texttt{pivot}(\mathcal{O}^{st}_{i})$;

24 if  _$u$ is a root_ then

25
$count_{i}\leftarrow\mathtt{cdag}^{\ast}({\mathcal{O}^{st}_{i}}[-\mathcal{A}^{+}({u})])+\varphi(u)\cdot\mathtt{cdag}^{\ast}({\mathcal{O}^{st}_{i}}[-\mathcal{D}^{+}({u})])$;

26

27 else

28
$count_{i}\leftarrow\mathtt{cdag}^{\ast}({\mathcal{O}^{st}_{i}}[-\mathcal{D}^{+}({u})])+\varphi(u)\cdot\mathtt{cdag}^{\ast}({\mathcal{O}^{st}_{i}}[-\mathcal{A}^{+}({u})])$;

29

30 end if

31

32 end if

33 $count\leftarrow count\cdot count_{i}$;

34

35 end foreach

36 $\texttt{insert\\_hash\\_table}(\mathcal{O},count)$;

37 return $count$;

38

39 end

Algorithm 3 The advanced version of Algorithm 2 with optimization modules.

##  6 Experiments and Results

We empirically evaluate the enumeration procedure from Algorithm 3 and various
practical speedups using randomly generated graphs. We then apply this
algorithm to four biomedical ontologies to gain insight into the sizes of
their concept annotation spaces.

###  6.1 Run-time Evaluation

We generated two sets of graphs to investigate the efficacy of our algorithm.
Each set contained 1000 graphs with either 25 or 100 vertices. To construct
each graph the vertices were added sequentially, with the proposed in-degree
in-deg($v$) of the $k$-th vertex $v$ generated according to a Poisson
distribution with parameter $\lambda$. This vertex then became a child of
$\min$(in-deg($v$), $k-1$) previously generated vertices that were themselves
selected uniformly randomly. The parameter $\lambda$ was selected according to
the $\Gamma(2.0,1.0)$ prior for each new graph and kept constant until the
graph was completed.

With these two sets of simulated graphs, we ran our algorithm with different
modules and pivot selection strategies. In particular, we evaluate pivot
selection based on (i) random selection of vertices, (ii) random selection of
multi-parent vertices, (iii) the degree criterion, (iv) the bound criterion,
and (v) the bottleneck criterion. For each pivoting strategy, we subsequently
add the pruning component, then hashing, and finally the graph reversal. The
criterion for graph reversal was the number of multi-parent vertices; i.e., a
graph will be reversed at any point during the recursive process if the
reversed graph contains fewer multi-parent vertices.

We report the average wall-time and average number of recursive calls over the
two sets of 1000 graphs ($|V|=25$ in Table 1; $|V|=100$ in Table 2). For the
smaller graphs, we also ran a brute-force algorithm that was further
convenient to empirically evaluate the correctness of our algorithm. We see
that simpler schemes perform better on small graphs where the number of
recursive calls per graph has not exceeded a few hundreds. On the other hand,
the advanced techniques show tangible benefits on the larger graphs reducing
the number of recursive calls and total computation time by orders of
magnitude. It is possible to envision other variations that could result in
further speedups; e.g., selecting multi-parent pivots with the highest degree.
These refinements, however, were beyond the scope of this paper.

Table 1: Experiments with simulated graphs with $|V|=25$ vertices. Each field in the table summarizes the per-graph wall-time over a set of 1000 graphs as well as the per-graph number of recursive calls, except for the brute-force method. The columns represent pivot selection strategies: (i) random, (ii) random multi-parent vertex (mpv), (iii) minimum bound, (iv) maximum degree, and (v) bottleneck. The rows represent successive additions of practical modules for speedups: ( ) basic approach from Algorithm 2, ( ) pruning, (  ) pruning and hashing, (  ) pruning, hashing, and graph reversal. brute-force | module | random | random mpv | min. bound | max. degree | bottleneck
---|---|---|---|---|---|---
571ms |  | 22.5ms | (313) | 5.3ms | (39) | 25.7ms | (23) | 1.2ms | (28) | 18.2ms | (47)
| 21.1ms | (97) | 14.3ms | (44) | 26.4ms | (25) | 10.2ms | (28) | 25.3ms | (44)
| 19.6ms | (71) | 14.4ms | (39) | 26.3ms | (23) | 10.2ms | (26) | 24.9ms | (34)
| 19.2ms | (67) | 7.5ms | (28) | 25.5ms | (23) | 7.5ms | (23) | 23.9ms | (31)
Table 2: Experiments with simulated graphs with $|V|=100$ vertices, with rows and columns identical to those in Table 1. The entries with an asterisk indicate that a sample of graphs was considered (instead of a full set of 1000) due to the long run-time. The brute-force algorithm was not considered as it was not feasible to compute the count for even a single graph. module | random | random mpv | min. bound | max. degree | bottleneck
---|---|---|---|---|---
| *3,102s | (119,745,876) | 5.21s | (52,954) | 114s | (25,416) | 9.93s | (101,342) | 122s | (526,925)
| *241s | (1,727,306) | 8.98s | (33,271) | 3.93s | (2,802) | 1.10s | (3,066) | 4.28s | (12,597)
| *132s | (387,910) | 7.35s | (14,913) | 3.67s | (1,111) | 0.92s | (1,107) | 3.08s | (2,052)
| 165s | (508,521) | 4.68s | (9,721) | 3.22s | (1,103) | 0.84s | (1,079) | 2.79s | (2,133)

###  6.2 Consistent Subgraphs in Biomedical Ontologies

We use 02/2017 versions of Gene Ontology (GO) and Human Phenotype Ontology
(HPO) as the target ontologies and compute the number of consistent subgraphs
in each of them. The algorithm is applied to each of the three domains of GO
[1]: (i) molecular function ontology (MFO; 10,789 terms) (ii) biological
process ontology (BPO; 29,575 terms) and (iii) cellular component ontology
(CCO; 4,085 terms). Together with HPO (12,167 terms), these four ontologies
are widely used in annotating functional terms of gene products [17, 11]. We
further define the annotation _level_ for each term in the ontology to be the
length of the longest path to the root. Starting from the root term, we add
more specific terms level-by-level so as to understand how the potential
annotation space grows with increased granularity of functional concepts.

In addition to level-wise full ontologies, we also investigate the “used”
ontologies in which each term was retained only if at least one protein in the
UniProt-GOA [9] and HPO [19] databases has been confidently assigned that term
(confident annotations include all experimental evidence codes as well as
“traceable author statement” and “inferred by curators”). Protein function
annotations were extracted from the 02/13/2017 release of the UniProt-GOA
database, which contains 64,362 proteins with confident MFO annotations,
84,413 proteins with BPO annotations and 79,630 proteins with CCO annotations.
HPO annotations were extracted from the 02/24/2017 release of the HPO database
where 6,411 genes with confident annotations were extracted.

Figure 5 shows the completed counts for both full and used level-wise
ontologies. For each ontology, we additionally compute the lower bound
(generally the larger of $2^{\ell}$ and $2^{r}$, where $r$ is the number of
roots) and estimate the upper bound (we convert a graph into a forest by
keeping only one randomly selected incoming edge for each multi-parent vertex
and then call $\mathtt{cforest}$).

Figure 5: Number of consistent subgraphs in level-wise GO and HPO. In each
panel, solid black lines mark the exact counts for “full” subgraphs and grey
dotted lines mark the exact counts for “used” subgraphs. Colored bars indicate
the estimated upper/lower bounds of the actual counts. The exact integer
counts are available upon request.

Although we were not surprised by the astronomical sizes of concept annotation
spaces, it was interesting to quantify them whenever feasible as well as to
observe an increasing difference between lower and upper bounds (in the 1000s
of orders of magnitude) with the level of the ontology. We also find it
interesting that a large number of ontological terms have never been used (see
Supplementary Materials). These outcomes raise questions regarding the
predictability of ontological annotations as most modern algorithms are asked
to provide accurate deep annotations to be deemed useful. However, annotation
spaces become exceedingly large almost instantaneously with the depth of the
ontology, which presents an immense computational and statistical challenge
for any prediction algorithm. We therefore believe that the balance between
ontology size/complexity and term granularity might become an important topic
for future discussions.

###  6.3 Entropy of Concept Annotation Spaces

The ability to enumerate subgraphs in relatively large ontologies presents an
opportunity to contrast the space of actual ontological annotations in
biological databases with the space of possible ontological annotations. To
investigate this, we first computed the entropy of actual annotations at
different levels in the ontology,

$\displaystyle
H(\mathcal{O}_{lvl})=-\sum_{i}P({\mathcal{O}_{lvl}}[S_{i}])\log_{2}P({\mathcal{O}_{lvl}}[S_{i}]),$

where $\mathcal{O}_{lvl}$ is the truncated ontology as in Section 6.2,
${\mathcal{O}_{lvl}}[S_{i}]$ corresponds to a distinct consistent subgraph
annotation observed at that level and $P({\mathcal{O}_{lvl}}[S_{i}])$ is the
probability that a protein is assigned annotation
${\mathcal{O}_{lvl}}[S_{i}]$. We first enumerated all observed subgraphs from
the UniProt-GOA or HPO database truncated to a particular level, calculated
their relative frequencies, and then plugged these relative frequencies into
the entropy formula above. On the other hand, the maximum entropy was computed
as $\log_{2}\mathtt{cdag}(\mathcal{O}_{lvl})$ by assuming equal probability
for every possible consistent subgraphs.

Figure 6: Ratio of entropies in the four ontologies. Circles with solid lines
show the ratio of observed entropy to the maximum entropy. Dotted lines
correspond to the estimated ratios as the average of the two ratios calculated
by lower/upper bound of the counts. The error bars suggest a possible
placement for the actual ratio.

Figure 6 shows the ratio between the two quantities for levels greater than 1,
suggesting that the world of protein functions, despite great diversity, has
low entropy relative to the possible maximum. Although the currently observed
functional annotations are incomplete, noisy and biased [23, 24, 10], this
suggests considerable departure from the uniform distribution and implies that
an extensive number of possible consistent subgraphs have not been used.

##  7 Related Work

There exists a body of literature in enumerative combinatorics related to our
work. One of the most relevant problems is the enumeration of directed acyclic
graphs with $n$ distinct (labeled) nodes [20]. The resulting count reflects
the size of the structure space of Bayesian networks with $n$ random variables
and, surprisingly, also corresponds to the number of matrices in
$\left\\{0,1\right\\}^{n\times n}$ with all eigenvalues real and positive
[13]. The number of labeled directed acyclic graphs with $n$ nodes does not
have a closed-form solution and is instead available as the A003024 sequence
in the On-Line Encyclopedia of Integer Sequences (OEIS);
https://oeis.org/A003024. The construction was originally proposed by Robinson
[20] and was further investigated by others [26, 21, 6].

The results on rooted labeled trees include both the enumeration of possible
number of trees and also the enumeration of subtrees for a given tree. There
are $n^{n-1}$ labeled rooted trees with $n$ nodes [7] that provide the integer
sequence A000169 in OEIS; https://oeis.org/A000169. The expansion to forests
gives $(n+1)^{n-1}$ using Cayley’s formula [3], as a single root can be added
to connect a forest of unrooted labeled trees into a rooted labeled tree. The
recurrence for the number of subtrees of a given tree was proposed by Ruskey
[22]; see Algorithm 1. The generalization to weighted subtrees was given by
Yan and Yeh [30]. Both algorithms are linear in $n$ assuming constant time
addition and multiplication.

Our work also relates to the research in ontology quality assurance. These
efforts typically include the analysis of irregularities and redundancy in
concept descriptors and graph structure [2, 27, 29]. Our work, primarily the
software we developed, contributes to this area by facilitating the analysis
of the annotation space.

##  8 Conclusions

This work presents a practical algorithm for enumerating consistent subgraphs
of directed acyclic graphs. Although this study largely addresses an
intellectual challenge of efficient subgraph enumeration, it might have
practical utility for the studies of annotation spaces in the biomedical
sciences and beyond; e.g., in text mining [8] and computer vision [15]. As
ontologies are easy to grow and hard to manually interrogate, we provide a
practical tool that can give insights into the complexity of concept
annotation tasks [17, 11]. As such, it may serve as a guide to ontology
developers.

## Acknowledgements

This work has been supported by the National Science Foundation grant
DBI-1458477 and the Indiana University Precision Health Initiative.

## References

  * [1] M. Ashburner, C. A. Ball, et al.  Gene ontology: tool for the unification of biology. The Gene Ontology Consortium.  Nat Genet, 25(1):25–29, 2000.
  * [2] O. Bodenreider.  Strength in numbers: exploring redundancy in hierarchical relations across biomedical terminologies.  AMIA Annu Symp Proc, pages 101–105, 2003.
  * [3] A. Cayley.  A theorem on trees.  Quart J Math, 23:376–378, 1889.
  * [4] W. T. Clark and P. Radivojac.  Information-theoretic evaluation of predicted ontological annotations.  Bioinformatics, 29(13):i53–i61, 2013.
  * [5] I. Friedberg and P. Radivojac.  Community-wide evaluation of computational function prediction.  Methods Mol Biol, 1446:133–146, 2017.
  * [6] I. M. Gessel.  Counting acyclic digraphs by sources and sinks.  Discrete Math, 160:253–258, 1996.
  * [7] J. L. Gross and J. Yellen.  Handbook of graph theory.  CRC Press, Boca Raton, Florida, U.S.A., 2004.
  * [8] M. Grosshans, C. Sawade, et al.  Joint prediction of topics in a URL hierarchy.  In Proceedings of the Joint European Conference on Machine Learning and Knowledge Discovery in Databases, ECML/PKDD 2014, pages 514–529, 2014.
  * [9] R. P. Huntley, T. Sawford, et al.  The GOA database: Gene Ontology annotation updates for 2015.  Nucleic Acids Res, 43(Database issue):D1057–D1063, 2015.
  * [10] Y. Jiang, W. T. Clark, et al.  The impact of incomplete knowledge on the evaluation of protein function prediction: a structured-output learning perspective.  Bioinformatics, 30(17):i609–i616, 2014.
  * [11] Y. Jiang, T. R. Oron, et al.  An expanded evaluation of protein function prediction methods shows an improvement in accuracy.  Genome Biol, 17(1):184, 2016.
  * [12] T. Joachims, T. Finley, and C.-N. J. Yu.  Cutting-plane training of structural SVMs.  Mach Learn, 77(1):27–59, 2009.
  * [13] B. D. McKay, F. E. Oggier, et al.  Acyclic digraphs and eigenvalues of (0,1)-matrices.  J Integer Seq, 7:04.3.3, 2004.
  * [14] Y. Moreau and L. C. Tranchevent.  Computational tools for prioritizing candidate genes: boosting disease gene discovery.  Nat Rev Genet, 13(8):523–536, 2012.
  * [15] Y. Movshovitz-Attias, Q. Yu, et al.  Ontological supervision for fine grained classification of street view storefronts.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2015, pages 1693–1702, 2015.
  * [16] S. Poux and P. Gaudet.  Best practices in manual annotation with the Gene Ontology.  Methods Mol Biol, 1446:41–54, 2017.
  * [17] P. Radivojac, W. T. Clark, et al.  A large-scale evaluation of computational protein function prediction.  Nat Methods, 10(3):221–227, 2013.
  * [18] P. N. Robinson and S. Bauer.  Introduction to bio-ontologies.  CRC Press, Boca Raton, Florida, U.S.A., 2011.
  * [19] P. N. Robinson and S. Mundlos.  The human phenotype ontology.  Clin Genet, 77(6):525–534, 2010.
  * [20] R. W. Robinson.  Counting labeled acyclic digraphs.  In Proceedings of the 3rd Ann Arbor Conference on Graph Theory, pages 239–273. Academic Press, 1971.
  * [21] V. I. Rodionov.  On the number of labeled acyclic digraphs.  Discrete Math, 105:319–321, 1992.
  * [22] F. Ruskey.  Listing and counting subtrees of a tree.  SIAM J Comput, 10(1):141–151, 1981.
  * [23] A. M. Schnoes, S. D. Brown, et al.  Annotation error in public databases: misannotation of molecular function in enzyme superfamilies.  PLoS Comput Biol, 5(12):e1000605, 2009.
  * [24] A. M. Schnoes, D. C. Ream, et al.  Biases in the experimental annotations of protein function and their effect on our understanding of protein function space.  PLoS Comput Biol, 9(5):e1003063, 2013.
  * [25] A. Sokolov and A. Ben-Hur.  Hierarchical classification of gene ontology terms using the GOstruct method.  J Bioinform Comput Biol, 8(2):357–376, 2010.
  * [26] R. P. Stanley.  Acyclic orientations of graphs.  Discrete Math, 5:171–178, 1973.
  * [27] K. Verspoor, D. Dvorkin, et al.  Ontology quality assurance through analysis of term transformations.  Bioinformatics, 25(12):i77–i84, 2009.
  * [28] M. Vihinen.  Variation Ontology for annotation of variation effects and mechanisms.  Genome Res, 24(2):356–364, 2014.
  * [29] G. Xing, G. Q. Zhang, and L. Cui.  FEDRR: fast, exhaustive detection of redundant hierarchical relations for quality improvement of large biomedical ontologies.  BioData Min, 9:31, 2016.
  * [30] W. Yan and Y. N. Yeh.  Enumeration of subtrees of trees.  Theor Comput Sci, 369(1-3):256–268, 2006.