On the optimality of some group testing algorithms

TL;DR

This paper analyzes the optimality of group testing algorithms, showing that the DD algorithm is suboptimal for certain parameters and that SCOMP and LP-based algorithms are at least as effective for large defect rates.

Contribution

The paper provides new bounds on the rate of the DD algorithm and establishes the optimality of SCOMP and LP-based algorithms for specific defect proportions.

Findings

DD is suboptimal for .41 < heta < 0.5

SCOMP and LP algorithms achieve at least the same rate as DD for heta .5

New upper bounds on the rate of the DD algorithm are derived.

Abstract

We consider Bernoulli nonadaptive group testing with $k = \Theta(n^\theta)$ defectives, for $\theta \in (0,1)$. The practical definite defectives (DD) detection algorithm is known to be optimal for $\theta \geq 1/2$. We give a new upper bound on the rate of DD, showing that DD is strictly suboptimal for $\theta < 0.41$. We also show that the SCOMP algorithm and algorithms based on linear programming achieve a rate at least as high as DD, so in particular are also optimal for $\theta \geq 1/2$.

Equations

$$y_t=\{\begin{array}{ll}1& \text{if there exists i∈K with xti=1,}\\ 0& \text{if for all i∈K we have xti=0.}\end{array}$$

$$T\ge {\log }_2\left(\genfrac{}{}{0ex}{}{n}{k}\right)\sim k{\log }_2\frac{n}{k}\sim (1-\theta )k{\log }_{2}n.$$

$$\frac{1}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)}\sum _{\begin{array}{c}\mathcal{K}\subset \{1,\dots ,n\}\\ |\mathcal{K}|=k\end{array}}\mathbb{P}(\hat{\mathcal{K}}\ne \mathcal{K}),$$

$$\frac{\log \left(\genfrac{}{}{0ex}{}{n}{k}\right)}{T}>R$$

$$C=\underset{\nu >0}{\max }\min \left\{h({\mathrm{e}}^{-\nu }),\frac{\nu }{{\mathrm{e}}^{\nu }\ln 2}\frac{1-\theta }{\theta }\right\}.$$

$$C=\{\begin{array}{ll}1& \text{for }\theta \le \frac{1}{3},\\ 0.531\frac{1-\theta }{\theta }& \text{for }\theta >0.359,\end{array}$$

$$R_{\mathrm{DD}}^{\ast }\ge \frac{1}{\mathrm{e}\ln 2}\min \left\{1,\frac{1-\theta }{\theta }\right\},$$

$$R_{\mathrm{DD}}^{\ast }=\frac{1}{\mathrm{e}\ln 2}\frac{1-\theta }{\theta }=0.531\frac{1-\theta }{\theta }\text{for }\theta \ge \frac{1}{2},$$

$$\begin{split}R^{*}_{\mathrm{DD}}\leq\max_{\nu>0}\min\bigg{\{}\mathrm{e}^{-\nu}\log\mathrm{e}^{\nu}+\nu\mathrm{e}^{-\nu}\log\frac{1}{1-\mathrm{e}^{-\nu}},\\ \frac{\nu}{\mathrm{e}^{\nu}\ln 2}\frac{1-\theta}{\theta}\bigg{\}}.\end{split}$$

$${\theta }^{\ast }=\frac{1}{2-\ln (1-{\mathrm{e}}^{-1})}=\mathrm{0.407.}$$

$$R$$

$$=q_0\log \frac{1}{q_0}+q_{1+}\log \frac{1}{q_{1+}}$$

$$=q_0\log \frac{1}{q_0}+(q_1+q_{2+})\log \frac{1}{q_{1+}}.$$

$$R_{\mathrm{DD}}^{\ast }\le q_0\log \frac{1}{q_0}+q_1\log \frac{1}{q_{1+}}.$$

$$q_0=(1-p{)}^k\sim {\left(1-\frac{\nu }{k}\right)}^k\to {\mathrm{e}}^{-\nu },$$

$$q_1=kp(1-p{)}^{k-1}\sim k\frac{\nu }{k}{\left(1-\frac{\nu }{k}\right)}^{k-1}\to \nu {\mathrm{e}}^{-\nu },$$

$$q_{1+}=1-q_0\sim 1-{\mathrm{e}}^{-\nu }.$$

$$R_{\mathrm{DD}}^{\ast }\le {\mathrm{e}}^{-\nu }\log {\mathrm{e}}^{\nu }+\nu {\mathrm{e}}^{-\nu }\log \frac{1}{1-{\mathrm{e}}^{-\nu }}.$$

$$R^{*}_{\mathrm{DD}}\leq\max_{\nu>0}\min\bigg{\{}\mathrm{e}^{-\nu}\log\mathrm{e}^{\nu}+\nu\mathrm{e}^{-\nu}\log\frac{1}{1-\mathrm{e}^{-\nu}},\\ h(\mathrm{e}^{-\nu})\,,\ \frac{\nu}{\mathrm{e}^{\nu}\ln 2}\frac{1-\theta}{\theta}\bigg{\}}.$$

$$h({\mathrm{e}}^{-\nu })=q_0\log \frac{1}{q_0}+q_{1+}\log \frac{1}{q_{1+}}$$

$${\mathrm{e}}^{-1}\log \mathrm{e}+{\mathrm{e}}^{-1}\log \frac{1}{1-{\mathrm{e}}^{-1}}<\frac{1}{\mathrm{e}\ln 2}\frac{1-\theta }{\theta }.$$

$$R_{\mathrm{COMP}}^{\ast }\le q_0\log \frac{1}{q_0}.$$

$$R_{\mathrm{COMP}}^{\ast }\le {\mathrm{e}}^{-1}\log \mathrm{e}=\frac{1}{\mathrm{e}\ln 2}.$$

$$R_{\mathrm{COMP}}^{\ast }=\frac{1}{\mathrm{e}\ln 2}(1-\theta ),$$

$$R_{\text{A}}^{\ast }\ge R_{\mathrm{DD}}^{\ast }\ge \min \left\{\frac{1}{\mathrm{e}\ln 2},\frac{1}{\mathrm{e}\ln 2}\frac{1-\theta }{\theta }\right\},$$

$$R_{\mathrm{SCOMP}}^{\ast }\ge R_{\mathrm{DD}}^{\ast }\ge \min \left\{\frac{1}{\mathrm{e}\ln 2},\frac{1}{\mathrm{e}\ln 2}\frac{1-\theta }{\theta }\right\}.$$

$$\sum _{i=1}^n{z}_i$$

$$\sum _{i=1}^n{x}_{ti}{z}_i\ge 1\text{when }{y}_t=1,$$

$$\sum _{i=1}^n{x}_{ti}{z}_i=0\text{when }{y}_t=0,$$

$$z_i\in \{0,1\}.$$

$$R_{\mathrm{LiPo}}^{\ast }\ge \frac{1}{\frac{8}{3}{\text{e}}^2\ln 2}\frac{1-\theta }{1+\theta }\approx 0.0732\frac{1-\theta }{1+\theta }.$$

$$R_{\mathrm{LP}}^{\ast }\ge R_{\mathrm{DD}}^{\ast }\ge \min \left\{\frac{1}{\mathrm{e}\ln 2},\frac{1}{\mathrm{e}\ln 2}\frac{1-\theta }{\theta }\right\}.$$

Full text

I Introduction

Group testing is a sparse inference problem that involves detecting the defective items in a population through the means of pooled tests. By testing a ‘pool’ – that is, a subset – of items, we discover whether that pool is free of defective items or whether it contains at least one defective item. After a sequence of such tests, we must infer (or ‘detect’) the true set of defective items. The goal is to find the defective set with low probability of error using as few tests as possible.

This naturally splits into two problems: first, designing the testing strategy, and second, detecting the defective items, given the test design and the results of the test.

For the design, here we consider Bernoulli nonadaptive testing. In nonadaptive testing, the test pools are all defined in advance, so the tests can be carried out in parallel. Bernoulli nonadaptive testing is the simplest form of nonadaptive testing, where each item is placed in each test independently with some fixed probability $p$. Bernoulli nonadaptive testing is simple, order optimal for all $\theta\in[0,1)$ (see [1, 2, 3, 4]), even has optimal constants for $\theta\in[0,1/3)$ (see [4, 5, 6]), and has been studied by many authors (see [1, 7, 2, 3, 4, 5, 6, 8, 9, 10, 11] for just a few examples).

For the detection problem, we consider three detection algorithms: one called DD (for ‘definite defectives’), one called SCOMP (which can be seen as a Sequential version of the earlier COMP algorithm), and approaches based on linear programming. In contrast to the theoretical algorithms often used to prove achievability results, these three algorithms are practical, in that they are computationally efficient and do not require knowledge of the true number of defectives. It is important to know when these practical algorithms are as good as the theoretically optimal but impractical detection algorithms. Here, we study this question in the context of Bernoulli nonadaptive testing. It is known that DD is optimal for denser cases, where we have more defectives. Our first result here is that DD is strictly suboptimal for sparser cases. Our second result is that SCOMP and linear programming approaches are also optimal in denser cases.

We begin with some notation. There are $n$ items, of which $k$ are defective, and we use $T$ tests. We assume that defectivity is rare, in that $k=o(n)$, and concentrate on the case where $k$ scales like $k=\Theta(n^{\theta})$ for some sparsity parameter $\theta\in(0,1)$. (The $\theta=0$ case behaves much the same, but one has to be a little careful over limiting arguments, as $k$ does not tend to $\infty$ as $n$ grows. We omit this case for brevity, although the results of this paper do carry over.) We assume that the true defective set $\mathcal{K}$ is chosen uniformly at random for all subsets of size $k$.

A useful way to keep track of the design of the test pools is to write $x_{ti}=1$ to denote that item $i\in\{1,2,\dots,n\}$ is in the pool for test $t\in\{1,2,\dots,T\}$, and $x_{ti}=0$ to denote that item $i$ is not in the pool for test $t$.

Definition 1

In Bernoulli nonadaptive group testing we set $x_{ti}=1$ with probability $p$ and $x_{ti}=0$ with probability $1-p$, independently over $t\in\{1,2,\dots,T\}$ and $i\in\{1,2,\dots,n\}$, for some fixed $p\in[0,1]$.

Recall that in standard group testing, a test outcome is positive if at least one item from $\mathcal{K}$ is in the test, and negative if no items from $\mathcal{K}$ are in the test. That is, writing $y_{t}$ for the outcome of test $t$, we have

[TABLE]

In general, since there are $\binom{n}{k}$ possible defective sets and $2^{T}$ possible outcomes from $T$ tests, we see that we have a lower bound on the number of tests required to reconstruct the defective set of

[TABLE]

One way of interpreting this is to say that specifying the defective set requires $\log_{2}\binom{n}{k}$ bits of information, and each test, with a yes–no answer, can impart at most $1$ bit of information. (From now on, $\log$ always means $\log_{2}$.) Given this, for a scheme that uses $T$ tests, we can regard the ‘bits learned per test’ $\log\binom{n}{k}/T$ as the rate of group testing. (This terminology follows that of Baldassini, Johnson, and Aldridge [12].)

Definition 2

The average error probability is defined as follows: if our estimate of the defective set $\mathcal{K}$ is $\hat{\mathcal{K}}$, then the error probability is

[TABLE]

where the probability is over the random choice of test design.

Definition 3

Given a detection algorithm, we say a rate $R$ is achievable if, for any $\epsilon>0$, then for $n$ sufficiently large and $k=k(n)=\Theta(n^{\theta})$, there exists a Bernoulli parameter $p=p(n)$ and a number of tests $T=T(n)$ with rate

[TABLE]

and Bernoulli nonadaptive testing having error probability less than $\epsilon$.

The capacity of Bernoulli nonadaptive group testing (that is, the maximum achievable rate for any detection algorithm) is

[TABLE]

Achievability was shown by Scarlett and Cevher [4] and the converse by Aldridge [6]. It is easy to check that for most $\theta$ this can be simplified to

[TABLE]

so the $\max$-$\min$ is only needed in the small interval $\theta\in(1/3,0.359)$.

However, the achievability result of Scarlett and Cevher [4] does not give a practical decoding algorithm to achieve the rate (1).

Aldridge, Baldassini and Johnson [2] showed that a simple algorithm called DD (for ‘definite defectives’) is optimal for $\theta\geq 1/2$. We define the DD algorithm in Section II. Also, here ‘optimal’ is used throughout to mean ‘achieves the capacity (1).’ One might wonder if DD is also optimal for $\theta<1/2$. Our first main result mostly answers this in the negative: we give a new upper bound for the rate of DD, showing that for $\theta<0.407$ DD is strictly suboptimal (Theorem 4).

Two decoding algorithms seem to perform better than DD in practice – thus far without theoretical justification. These are the SCOMP algorithm of Aldridge, Baldassini, and Johnson [2], and approaches based on linear programming (see, for example, [10, 7]). We define these algorithms in Section III-B and III-C. Figure 1 shows a simulation of all three algorithms, demonstrating the superiority of SCOMP and linear programming. This theoretical performance suggests that the rate of these algorithms should be at least as high as DD, so in particular they should also be optimal for $\theta>1/2$. Our second main result proves this (Theorems 8 and 9).

The rate bounds implied by our main results are illustrated in Figure 2.

II DD is suboptimal for $\theta<0.42$

The DD algorithm of Aldridge, Baldassini, and Johnson [2] works as follows:

Any item in a negative test is definitely nondefective. The remaining items we call ‘possible defectives’. 2. 2.

If a (necessarily positive) test contains exactly one possible defective, then that item is in fact definitely defective (‘DD’). 3. 3.

The definitely defective items are declared to be defective. All other items – the definitely nondefective items and any remaining possible defectives – are declared to be nondefective.

DD is ‘practical’ in the sense that it requires only $O(nT)=O(kn\log n)=O(n^{1+\theta}\log n)$ work and storage – which is what is required just to store and read the test design. In comparison, the decoder of Scarlett and Cevher [4] involves solving the maximum likelihood problem, which is in NP in general. Further, DD does not require the detector to know $k$ exactly or even approximately, while the maximum likelihood decoder requires the exact value of $k$.

It was shown in [2] that, with Bernoulli tests, the DD algorithm can achieve a rate of

[TABLE]

where $1/\mathrm{e}\ln 2\approx 0.531$.

This can be compared with the capacity for Bernoulli testing (1). Since (1) gives an upper bound of $(1/\mathrm{e}\ln 2)(1-\theta)/\theta$ for $\theta>0.359$, and (2) gives the same rate for $\theta\geq 1/2$, we see that for $\theta\geq 1/2$ we have the equality

[TABLE]

and DD is optimal in this region. (This was first noted in [6].)

Here we give an upper bound for DD, which shows that it is strictly suboptimal for some smaller values of $\theta$.

Theorem 4

The rate of DD for Bernoulli nonadaptive group testing is bounded above by

[TABLE]

In particular, DD is strictly suboptimal for $\theta<\theta^{*}$, where

[TABLE]

Although the bound (4) appears complicated, it is easy to check that for $\theta>\theta^{*}$ the second minimand dominates and gives the same bound $0.531(1-\theta)/\theta$ as the standard capacity bound (1). Further, one can check numerically that for $\theta<0.357$ the first minimand of (4) dominates and gives a bound of $0.853$. It is only for $\theta\in(0.357,0.407)$ that the complicated expression (4) is necessary. Whether or not DD is optimal in the range $\theta\in(0.407,1/2)$ remains an open question.

The above bounds are illustrated in Figure 2.

Proof:

The key point is that DD can make use of tests with no defectives (to eliminate ‘possible defectives’), and tests with one defective (if there are no other possible defectives in the test, the defective is ‘definitely defective’), but does not make use of any tests containing two or more defectives (they necessarily contain at least two possible defectives).

Suppose we use independent tests (as we do with a Bernoulli design). Write $q_{d}$ for the probability a test contains exactly $d$ defectives and $q_{d+}$ for the probability a test contains at least $d$ defectives. Then a standard information argument gives a bound in terms of the entropy

[TABLE]

(Note that when $q_{0}=q_{1+}=1/2$, this gives the standard counting bound $R\leq 1$.) While it is helpful to think of this in terms of ‘bits learned per tests’, it can be made rigorous using the concept of typical sets – see, for example, [6, Lemma 1] – which tells us that with high probability we see one of $2^{Th(q_{0})}$ ‘typical’ outcomes.

Since DD does make use of those tests containing $2$ or more defectives, we get a tighter bound

[TABLE]

This is because we learn at most $\log 1/q_{0}$ bits with probability $q_{0}$ from the negative tests, and from positive tests we learn $\log 1/q_{1+}$ bits only with probability $q_{1}$, since DD gets no information from a (positive) test with probability $q_{2+}$.

Now recall we are using a Bernoulli test design. Setting the parameter as $p\sim\nu/k$ (as we must to get a nonzero rate [6]), we have

[TABLE]

as $k\to\infty$. Hence the bound (5) becomes

[TABLE]

Combined with the capacity bound (1), this gives

[TABLE]

Note, however, that the second minimand

[TABLE]

is always dominated by the first minimand, so the second minimand can be deleted.

For the final part of the theorem, we note that second minimand in (4) is maximised at $\nu=1$, and this dominates until

[TABLE]

By rearranging, and noting that $\log\mathrm{e}=1/\ln 2$, we see that this is precisely when $\theta<\theta^{*}$. ∎

We remark that a similar proof technique could be used to bound the rate of the COMP algorithm of Chan et al [1]. This algorithm makes no use of of any positive tests, so the equivalent to the bound (5) would be

[TABLE]

Taking $q_{0}\sim\mathrm{e}^{-\nu}$ as above, and optimising at $\nu=1$ gives the bound

[TABLE]

Compared to the actual maximum achievable rate of COMP [1, 6], which is

[TABLE]

we see that such a bound is tight for $\theta=0$ but loose for $\theta>0$.

III Algorithms that are optimal for $\theta\geq 1/2$

We have seen that the DD algorithm with Bernoulli testing achieves the rate (2). Together with the capacity bound (1), we see that DD is optimal for Bernoulli testing for $\theta\geq 1/2$.

As we saw in Figure 1, the SCOMP algorithm [2] and approaches based on linear programming [10, 7] outperform DD in simulations. In this section, we prove that the SCOMP algorithm and LP approaches achieve at least as high a rate as DD, so in particular satisfy the bound (2) and are also optimal for $\theta\geq 1/2$.

III-A Property $P$

The crucial step to showing an algorithm outperforms DD is the a property we call ‘Property $P$’. First recall the following definition.

Definition 5

A set $\mathcal{L}\subset\{1,2,\dots,n\}$ is satisfying if no negative test contains any item from $\mathcal{L}$ and every positive test contains at least one item from $\mathcal{L}$.

So a satisfying set is one that ‘fits’ the observations. In particular, the true defective set $\mathcal{K}$ is always a satisfying set, although their may be others.

Definition 6

A detection algorithm A has Property $P$ if, whenever the DD algorithm outputs a satisfying set, then A also outputs that same satisfying set.

Lemma 7

Suppose the detection algorithm A satisfies Property $P$. Then A succeeds whenever DD succeeds. In particular, with Bernoulli nonadaptive group testing, we have

[TABLE]

and A is optimal for $\theta\geq 1/2$.

Proof:

Note that DD is not guaranteed to output a satisfying set. If the output of DD is not satisfying, then clearly the algorithm is in error. If the output of DD is satisfying, then it is necessarily the unique smallest satisfying set (since each definite defective appears uniquely in some test, so must be in the true defective set), and the algorithm may be successful.

Since algorithm A has Property $P$, in this latter case, A also outputs this unique smallest satisfying set. Hence, whenever DD succeeds, A will succeed also.

The second part of the lemma then follows immediately from (2). ∎

III-B SCOMP

SCOMP is an algorithm due to Aldridge, Baldassini, and Johnson [2] that builds a satisfying set by starting from the DD set of definite defectives and sequentially adding new items until a satisfying set is reached. (The name comes from ‘Sequential COMP’, as it can be viewed as a sequential version of the COMP algorithm of Chan et al [1].)

Although recalling the precise details of how SCOMP extends from DD to a satisfying set is not required for the proof of following result, we define it here for completeness.

Any item in a negative test is definitely nondefective. The remaining items we call ‘possible defectives’. 2. 2.

If a (necessarily positive) test contains exactly one possible defective, then that item is in fact definitely defective (‘DD’). 3. 3.

The definitely defective items are declared to be defective, we call these $\hat{\mathcal{K}}$, and definitely nondefective items are declared to be nondefective. The other possible defectives are not yet declared either way. 4. 4.

Any positive test is called ‘unexplained’ if it does not contain any items from $\hat{\mathcal{K}}$. Add to $\hat{\mathcal{K}}$ the possible defective not in $\hat{\mathcal{K}}$ that is in the most unexplained tests, marking the corresponding tests as no longer unexplained. (Ties may be broken arbitrarily.) 5. 5.

Repeat step 4 until no tests remain unexplained. The estimate of the defective set is $\hat{\mathcal{K}}$.

Note that a satisfying set leaves no unexplained tests, and that any set containing no definite nondefectives and leaving no unexplained tests is satisfying. Note also that the set of all possible defectives is satisfying, so the SCOMP algorithm does indeed terminate.

Theorem 8

SCOMP with Bernoulli tests can achieve the rate

[TABLE]

In particular, SCOMP is optimal for $\theta\geq 1/2$.

Proof:

Following Lemma 7, we have wo show that SCOMP satisfies Property $P$.

As explained above, the SCOMP algorithm begins with the DD algorithm, then takes further steps to ensure a satisfying set is reached. However, if DD already provides a satisfying set, no further steps are taken, and SCOMP halts and outputs the DD set. Hence Property $P$ is satisfied. ∎

III-C Linear programming approaches

LP algorithms are based on the fact choosing the smallest satisfying set as an estimate is known to be optimal [6], but is likely to be impractical, as it requires solving an instance of the NP-hard set cover problem. Specifically, the smallest satisfying set is related to the solution to the [math]–$1$ linear program

[TABLE]

(Recall that $x_{ti}$ indicates if item $i$ is in test $t$, and $y_{t}$ is the outcome of test $t$.) The set $\hat{\mathcal{K}}=\{i:z_{i}=1\}$ is the smallest satisfying set.

LP approaches attempt to estimate the defective set via the relaxed version of the [math]–$1$ problem, where the $z_{i}$s can be any non-negative real numbers. Linear programs like this can be solved quickly using, for example, the simplex algorithm.

There are various heuristics for how to turn the optimal solution $(z_{i})$ to the relaxed program into an estimate of the defective set. For the purposes of Theorem 9, it suffices that in the event that all $z_{i}$s are [math] or $1$ then the heuristic chooses $\hat{\mathcal{K}}=\{i:z_{i}=1\}$ (as any sensible heuristic must).

For example, one could consider the following very crude method: if there is any $i$ with $z_{i}\neq 0,1$, declare a global error; else estimate $\hat{\mathcal{K}}=\{i:z_{i}=1\}$ to be the defective set. Malioutov and Malyutov [10] suggest using the estimate $\hat{\mathcal{K}}=\{i:z_{i}>0\}$, and show strong performance on simulated problems. (This was the ‘Malioutov rule’ used for the simulations in Figure 1.) Another suggestion could be to estimate $\hat{\mathcal{K}}=\{i:z_{i}\geq 1/2\}$, or to place each item $i$ in $\hat{\mathcal{K}}$ independently with probability $z_{i}$.

Along these lines, the ‘LiPo’ algorithm of Chan et al [7] assumes the detector knows $k$ exactly, so the LP relaxation can be instead phrased as a feasibility problem. They show that LiPo achieves a nonzero rate for all $\alpha\in(0,1)$, and specifically

[TABLE]

However, by using Property $P$, we can show a higher rate, which is optimal for $\theta>1/2$.

Theorem 9

Any LP approach as described above with Bernoulli tests can achieve the rate

[TABLE]

In particular, LP approaches are optimal for $\theta\geq 1/2$.

Proof:

Following Lemma 7, we must show that Property $P$ is satisfied.

Note that any item $i$ in a negative test (a definite nondefective) will have $z_{i}=0$ to satisfy the second constraint of the linear program. Further, if a test $t$ contains only one possible (and thus definite) defective $i$, the LP solution must have $z_{i}\geq 1$ to satisfy the corresponding constraint, and will choose $z_{i}=1$ to minimise $\sum_{i}z_{i}$. Finally, if these definite defectives form a satisfying set, all constraints are satisfied, and the LP will set all other $z_{i}$s to [math], to minimise $\sum_{i}z_{i}$. Hence Property $P$ holds, and we are done. ∎

Acknowledgments

The author thanks Jonathan Scarlett for useful discussions.