Compressed Sensing-Based SARS-CoV-2 Pool Testing

Introduction

In most countries, it seems that the population will have to live with COVID-19 for at least a while. We hope that most readers of this article will require very little convincing of the need to continually monitor the level of infection in the population. However, the monitoring has financial implications, in particular the cost of test kits. This results in a strong desire to reduce the number of tests required to identify infected individuals. One approach is pool testing, where test samples of individuals are pooled together and tested as one sample. If this sample turns out negative, then every individual whose sample is in the pool is declared negative. Otherwise, the group of samples can be subdivided and retested. This has the potential to significantly reduce the number of tests required. It has been intensively studied for COVID-19 (see for example SHB$^{+}$20MNB$^{+}$20), and the United States Federal Drug Agency (FDA) granted pooled testing methods an emergency use authorization in July 2021.

The mathematical field of group testing is concerned with this pooling problem and has therefore also gained interest recently. It was invented by Dorfman Dor43. The methods from classical group testing often suffer from several drawbacks including slowness due to adaptivity of the tests and sensitivity to errors DH06. More on group testing in the next section. Compressed sensing, on the other hand, is a mathematical field concerned with recovering a vector with many zero entries from as few nonadaptive linear measurements as possible Don06. Compressed sensing has achieved several theoretical goals including achieving the optimal number of measurements, independence of the measurements from each other, robustness to noise, and sometimes even error-correcting properties Don06CRT06LDB09. A more detailed discussion on compressed sensing follows after the discussion on group testing.

In this manuscript, we revisit the use of compressed sensing for pool testing. We propose to use compressed sensing instead of classical group testing for viral detection since it tackles some weak points in the method of group testing mentioned above. We by no means claim originality in this; one of the first works is GIS08. However, we advocate for the use of a practical pooling strategy out of the forest of strategies proposed in the literature and the use of a new algorithm with comparable performance to the state-of-the-art in the field but having additional desirable properties that give it an advantage over other algorithms for pool testing.

The use of compressed sensing for viral detection of pools is premised on the fact that it is possible to measure the viral loads of patients. Presently, a standard test for the detection of SARS-CoV-2 is the real-time or quantitative polymerase chain reaction (qPCR).⁠¹ The qPCR measurement are fluorescence intensities at every DNA amplification cycle, and potentially every one of those measurements can be used to compute viral loads.

¹

This should not be confused with a reverse transcription PCR which is often being abbreviated RT-PCR SLB$^{+}$06.

✖

Remark 1.

The main reason for pooling, be it using group testing or compressed sensing, is to reduce the number of tests needed to identify infected individuals. This is why we focus on the number of pools (directly or indirectly in evaluating performance of polling designs and algorithms), typically given in complexity (“Big-O”) notation.

This article is organized as follows. The next two sections give a more detailed introduction to group testing and compressed sensing. Then we discuss the modeling of the pooling process and the identification of infections. We end with a discussion and conclusion.

Group testing

Group testing is concerned with discovering a class of interest from a larger set containing non-members of the class of interest without testing each member of the set. This is done by subdividing the set into subsets and testing each subset as one element. The goal is that with the number of subset/groups being far smaller than the total number of elements in the set, one is able to successfully identify the class of interest. More precisely, group testing seeks to identify $k$ members of a class of interest out of a set of $n$ elements by performing $m$ tests of $m$ groups, where $m\ll n$. In the disease identification setting, we would like to identify $k$ infected individuals from a population of $n$ individuals performing $m$ group (pool) tests, again we require that $m\ll n$. The ratio $k/n$ (denoted as $p$) is known as the prevalence. Usually, group testing works well when $p$ is small. It can be shown by a combinatorial argument that we need $\mathcal{O}\left(k\log _2\left(n/k\right)\right)$ tests for pool testing to work Hwa72. Figure 1 shows an example of pool testing to identify an infected individual. Here each individual participates in only one pool.

Figure 1.

Testing a population of size $n=24$ with only 1 infected individual (in white), i.e., prevalence of $\approx 4\%$. If each individual is tested, it would required 24 tests (i.e., 24 test kits). Using group testing with pool sizes of 6, members of pools 2, 3, and 4 are declared negative; while those in pool 1 are positive. Each individual in pool 1 is then retested, making it a total of 10 tests for group testing as opposed to 24 tests for testing each individual.

Figure 2.

Testing a population of size $n=25$ with only 1 infected individual (white dot), i.e., prevalence of $4\%$. Individual testing would require 25 tests (i.e., 25 test kits). Using group testing with pool sizes of 5 and each individual participating in two pools, the infected individual is identified by the row and column pools which contain their sample. This pooling strategy requires only 10 tests.

We can have pool testing setups where individuals participate in more than one pool. Figure 2 illustrates this with each individual participating in 2 pools. This can be represented in a two-dimensional (2D) grid; while the example in Figure 1 can be considered one-dimensional (1D). Going from 1D to 2D may reduce the need to retest. Moreover, this can be extended to a $d$-dimensional (dD) setting, where each individual will participate in $d$ pools. This is what the successfully applied hypercube method is about, see MNB$^{+}$20.⁠² The 1D, 2D, …, dD pool testing are instances of the so-called array testing, see PS94.

²

This is the approach being used in Rwanda.

✖

Mathematically, we can represent the population of size $n$ by a binary vector, $\mathbf{x}\in \{0,1\}^n$ where there are $k$ ones for the infected individuals and $n-k$ zeros for the non-infected individuals. Pooling samples of a group together is equivalent to evaluating

$$\begin{equation} \mathbf{a}\circ \mathbf{x}\coloneq \max _{j=1,\dots ,n}\{a_jx_j\}, \tag{1} \end{equation}$$

for some appropriate binary vector $\mathbf{a}\in \{0,1\}^n$. For a matrix $\mathbf{A}\in \{0,1\}^{m\times n}$ whose rows $\mathbf{a}_i$ represents the $m$ groups, we set

$$\begin{equation} \mathbf{y}\coloneq \mathbf{A}\circ \mathbf{x}\coloneq \left[\mathbf{a}_i\circ \mathbf{x}\right]_{i=1,\dots ,m}. \tag{2} \end{equation}$$

The vector $\mathbf{y}$ is referred to as the measurement/observation vector.

The recovery of the infected individuals, equivalently $\mathbf{x}$, from observing $\mathbf{y}$ and knowing $\mathbf{A}$ would lead to solving the following binary system of equations, since both $\mathbf{A}$ and $\mathbf{x}$ are binary.

$$\begin{equation} \mathbf{A}\circ \mathbf{x} = \mathbf{y}. {\tag{3}} \end{equation}$$

The setting described above where the groups are fixed is known as the nonadaptive case of group testing. Alternatively, one may form a group after knowing the outcome of the preceding test(s). In that case we perform adaptive group testing. Nonadaptive approaches are faster but they require more measurements than adaptive approaches. Adaptive methods achieve the optimal $\mathcal{O}\left(k\log _2\left(n/k\right)\right)$ measurements Hwa72; while non-adaptive methods need $\mathcal{O}\left(k^2\log _k\left(n\right)\right)$ measurements AA13.

Many recovery algorithms for group testing have been proposed. One popular such algorithm for nonadaptive group testing is Combinatorial Orthogonal Matching Pursuit (COMP). In the noiseless setting the COMP algorithm assigns a sample as positive (i.e., infected) if and only if all the $k$ tests containing the sample satisfy $y_m>0$. COMP will work perfectly, if the pooling strategy is such that we have “$k$-disjunct” sets of measurements (see PBJ20, Definition 4.1). Another algorithm for the adaptive case, is the binary splitting algorithm by Hwa72.

Compressed sensing

Mathematically, compressed sensing is concerned with the construction of a linear operator and the solution of an underdetermined system of linear equations resulting from the application of the linear operator. More precisely, let $\mathbf{A} \in \mathbb{R}^{m\times n}$ be our linear operator, let $\mathbf{x}\in \mathbb{R}^n$ be the unknown variable vector and let the application of $\mathbf{A}$ result in $\mathbf{y}\in \mathbb{R}^m$. Therefore, with $m\ll n$, we have an underdetermined linear system (note both $\mathbf{A}$ and $\mathbf{x}$ are real):

$$\begin{equation} \mathbf{A}\mathbf{x} = \mathbf{y}. {\tag{4}} \end{equation}$$

In compressed sensing parlance, $\mathbf{A}$ is called the sensing/measurement matrix, $\mathbf{y}$ is the measurement vector, and $\mathbf{x}$ is the signal of interest. Linear algebra tells us that there are infinitely many solutions to 4. However, it is possible to obtain a unique solution if we make some assumptions on $\mathbf{A}$ and $\mathbf{x}$. The main assumption on $\mathbf{x}$ is that it has some simplicity/redundancy. The simplicity of $\mathbf{x}$ is either sparsity or compressibility. By sparsity of $\mathbf{x}$ we mean it has few nonzero elements. Let $\mathcal{S}$ be the support of $\mathbf{x}$, i.e., $\mathcal{S} = \{\, j\mathrel{:}x_j\neq 0, \ \forall j\in 1, \ldots , n \,\}$ and let $|\mathcal{S}|\leq k$, then $\mathbf{x}$ is said to be $k$-sparse. Conventionally, the $\|\mathbf{x}\|_0$ (the $\ell _0$-norm of $\mathbf{x}$)⁠³ is used to denote the sparsity of $\mathbf{x}$. Precisely,

³

This is not a norm because $\|\lambda \mathbf{x}\|_0 = \|\mathbf{x}\|_0$ for all scalars $\lambda \neq 0$ and all vectors $\mathbf{x}$.

✖

$$\begin{equation} \|\mathbf{x}\|_0 \coloneq \sum _{j=1}^{n}\mathbb{I}_{\mathcal{S}} \left(j\right), {\tag{5}} \end{equation}$$

where $\mathbb{I}$ is the indicator/characteristic function. On the other hand, $\mathbf{x}$ is said to be $k$-compressible when $\mathbf{x}$ can be approximated quite well by a $k$-sparse vector. Furthermore, $\mathbf{x} = \Psi \mathbf{z}$ is $k$-sparse ($k$-compressible) in a basis $\Psi$ when $\mathbf{z}$ is $k$-sparse ($k$-compressible). We denote the restriction of $\mathbf{x}$ on its support as $\mathbf{x}_{\mathcal{S}} \in \mathbb{R}^k$ and the restriction of $\mathbf{x}$ on the complement of the support as $\mathbf{x}_{\mathcal{S}^c} \in \mathbb{R}^{n-k}$.

The assumption on $\mathbf{A}$ is that it is an information-preserving projection or a bi-Lipschitz linear metric space embedding of all $k$-sparse vectors into $\mathbb{R}^m$ (henceforth referred to as a stable linear embedding). The conditions for $\mathbf{A}$ to fulfill the information preservation requirement include the following.

(i)

Restricted Isometry Property (RIP). A more general definition of the RIP, i.e., $\ell _p$-norm-restricted isometry property ($\mathrm{RIP}_{p}$) for $p\geq 1$ is the following.

Definition 2.

Matrix $\mathbf{A}$ has $\mathrm{RIP}_{p}$ of order $k$, with constants $\delta _k<1$, if for all $k$-sparse $\mathbf{x}$, it satisfies

$$\begin{equation} (1-\delta _k)\|\mathbf{x}\|_p^p \leq \|\mathbf{A}\mathbf{x}\|_p^p \leq (1 + \delta _k)\|\mathbf{x}\|_p^p. {\tag{6}} \end{equation}$$

For small $\delta _k$, we have $\mathbf{A}$ being a near isometry and this implies that it is information preserving. Note that $p=2$ is typically written as RIP (without the subscript $2$).

(ii)

Nullspace Property (NSP), here defined in the noiseless settings.

Definition 3.

Matrix $\mathbf{A}$ has the nullspace property of order $k$, if for any $\mathbf{v} \in \operatorname {Null}(\mathbf{A})$ and any set $\mathcal{S}\subset \{1,\dots ,n\}$ with $|\mathcal{S}|\leq k$, we have

$$\begin{equation} \|\mathbf{v}_{\mathcal{S}}\|_1 \leq \|\mathbf{v}_{\mathcal{S}^c}\|_1. {\tag{7}} \end{equation}$$

If $\mathbf{A}$ satisfies this property, then there exists a unique $k$-sparse vector solving 4.

(iii)

$M^+$ criterion is defined as follows.

Definition 4.

Matrix $\mathbf{A}$ obeys the $M^+$ criterion with vector $\mathbf{u}$ and constant $\kappa$, if $\mathbf{A}^T\mathbf{u} > 0$ and $\kappa \coloneq \max _{i\in [n]} \left|\left(\mathbf{A}^T\mathbf{u}\right)_i\right|\times \max _{i\in [n]} \left|\left(\left(\mathbf{A}^T\mathbf{u}\right)_i\right)^{-1}\right|$.

Note that $\kappa$ is actually a condition number of the diagonal matrix with diagonal $\mathbf{A}^T\mathbf{u}$.

An important goal of compressed sensing is to make projections by $\mathbf{A}$ to map to spaces with very small dimensions (i.e., small $m$) such that recovery is possible to a certain error tolerance. This translates to finding matrices $\mathbf{A}$ with small $m$, which have a good restricted isometry property or a good nullspace property. It is known that matrices can have $\mathrm{RIP}_{p}$ with $p=1$ or $p=2$ or the nullspace property of order $k$, if the number of measurements, $m = \mathcal{O}\left(k\log \left(n/k\right)\right)$ CDD09. This is referred to as an optimal sampling rate.

The second part of the compressed sensing problem is the reconstruction of $\mathbf{x}$ from the projection $\mathbf{A}\mathbf{x}$. The system in 4 is underdetermined, which means $\mathbf{A}$ is not invertible. The aim of getting the sparsest solution that is faithful to the data leads to what is called the $\ell _0$ problem, i.e.:

$$\begin{equation} \hat{\mathbf{x}} = \operatorname {argmin}_{\mathbf{z}} \|\mathbf{z}\|_0 \quad \text{such that} \quad \mathbf{A}\mathbf{z} = \mathbf{y}. {\tag{8}} \end{equation}$$

This is a combinatorial and nonconvex problem, which has been shown to be NP-hard to solve. However, many discrete (also referred to as “greedy”) algorithms have been proposed for this problem with provable recovery guarantees. These include iterative hard thresholding (IHT), orthogonal matching pursuit (OMP), and compressive sampling matching pursuit (CoSaMP) FR13.

On the other hand, the $\ell _0$ problem can be relaxed to a convex one. A typical case is the Basis Pursuit (BP) algorithm that solves

$$\begin{equation} \hat{\mathbf{x}} = \operatorname {argmin}_{\mathbf{z}} \|\mathbf{z}\|_1 \quad \text{such that} \quad \mathbf{A}\mathbf{z} = \mathbf{y}. {\tag{9}} \end{equation}$$

It has been well established that the solution of BP coincides with the solution of the $\ell _0$ problem under most of the conditions on $\mathbf{A}$ and $\mathbf{x}$ stated above. The recovery of all compressed sensing algorithms are expected to be stable and obey an instance optimality ($\ell _p/\ell _q$-approximation) guarantee according to which any solution $\hat{\mathbf{x}}$ satisfies

$$\begin{equation} \|\hat{\mathbf{x}} - \mathbf{x}\|_p \leq C\frac{\sigma _k(\mathbf{x})_q}{\sqrt {k}}, \quad \text{for }1\leq q \leq p \leq 2, \tag{10} \end{equation}$$

where $C$ is an absolute constant independent of $\mathbf{x}$ and

$$\begin{equation} \sigma _k(\mathbf{x})_q \coloneq \min _{ {k-\text{sparse } \mathbf{x}'}} \|\mathbf{x} - \mathbf{x}'\|_q, \tag{11} \end{equation}$$

that is the best $k$-term approximation of $\mathbf{x}$. There is a more general definition of the instance optimality encompassing the noise setting too FR13.

We conclude this section by mentioning compressed sensing with non-negativity constraints. This problem is typically formulated in the following way.

$$\begin{equation} \hat{\mathbf{x}} = \operatorname {argmin}_{\mathbf{z}\geq 0} \|\mathbf{A}\mathbf{z} - \mathbf{y}\|. {\tag{12}} \end{equation}$$

In the case relevant to the focus of this manuscript, which is pool testing, the norm in 12 is taken to be $\|\cdot \|_1$ (i.e., the $\ell _1$-norm). In this setting 12 is known as the non-negative least absolute deviation (NNLAD) problem. The authors of PBJ21 proposed an efficient and tuning-free algorithm (dubbed NNLAD) for this problem.

Viral Detection in Pooled Tests

Testing design

Suppose we want to find $k$ individuals infected with a virus among $n$ individuals. According to the information theoretic lower bound we require at least $\log \left({\binom{n}{k}}\right)$ tests to find the infected individuals. The binary splitting algorithm proposed by Hwa72 finds the $k$ infected individuals with a number of tests in $\mathcal{O}\left(k\log _2\left({{n}/{k}}\right)\right)$. However, since the tests are adaptive, each subsequent test is designed depending on the outcome of previous tests, and thus each test has to wait for the result of the previous one. If it takes a long time to perform a test, it might be desirable to perform multiple tests at once. In such cases nonadaptive methods are preferable. There are many deterministic nonadaptive methods for group testing. Most of these prove their results using disjunct matrices. For instance, in AA13 a nonadaptive group testing method is presented whose number of tests is in $\mathcal{O}\left(k^2{\log _2(n)}/{\log _2(k)}\right)$.

We propose to use compressed sensing with additional nonnegativity constraints to solve this problem. We collect specimens from $n$ individuals arranged in a vector $\mathbf{x}\in \mathbb{R}^n$, and the amount of virus in the specimen of the $j$th individual is denoted by the non-negative quantity $x_j\in \mathbb{N} \cup \left\{0\right\} \subset \mathbb{R}$. We assume that viruses are evenly distributed in each specimen, meaning that if we take $\alpha \in \left[0,1\right]$ of the volume of the specimen of the $i$th individual, it will contain roughly $\alpha x_j$ viruses. Since, $k$ individuals are infected we have $\left\|\mathbf{x}\right\|_0=k$. Let the sample of the $i$th test contain a fraction $a_{ij}$ of the amount of specimen of the $j$th individual. The sample of the $i$th test thus contains, up to rounding errors, the amount of viruses

$$\begin{equation*} \sum _{j\in \left[n\right]}a_{ij}x_j=\left(\mathbf{A}\mathbf{x}\right)_i \eqcolon y_i, \end{equation*}$$

where $\mathbf{A}$ is an $m\times n$ matrix with entries $a_{ij}\in [0,1]$ with column sums of at most one. The number of tests is thus $m$. Often one makes the assumption that $\mathbf{x}$ is a random vector with i.i.d. elements, for instance a Poisson random variable multiplied by a Bernoulli random variable with parameter $p={k}/{n}$. Instead we take it here as a deterministic unknown. It is assumed that qPCR can be used to generate an estimate $y_i$ of the amount of virus in the $i$th test $\left(\mathbf{A}\mathbf{x}\right)_i$. This procedure is not accurate and errors $e_i\coloneq y_i-\left(\mathbf{A}\mathbf{x}\right)_i$ might occur. We try to recover the amount of virus in the specimen of the individuals according to

$$\begin{align} \mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{e}, {\tag{13}} \end{align}$$

with a possibly small $\left\|\mathbf{x}\right\|_0$ and $m$ as small as possible. As discussed above, this is exactly a compressed sensing problem where, due to the nature of the problem, $\mathbf{x}$ is non-negative and exactly $k$-sparse (compressed sensing guarantees often also work for compressible vectors which are well approximated by sparse vectors). The theory of compressed sensing states that there exist matrices and efficient decoders that allow recovery of $\mathbf{x}$ if $m$ is $\mathcal{O}\left(k\log {\left({n}/{k}\right)}\right)$ Don06CRT06. It remains to design a suitable measurement/sensing matrix and a reconstruction algorithm.

Measurement matrices

Recall from the discussion on compressed sensing above that there exist matrices that allow recovery of $\mathbf{x}$ from 13 if $m$ is in $\mathcal{O}\left(k\log {\left({n}/{k}\right)}\right)$, referred to as the optimal scaling. The testing design described above is a special compressed sensing task where the elements of the matrix $\mathbf{A}$ and the unknown sparse vector $\mathbf{x}$ are non-negative. Non-negative $\mathbf{A}$ does not have $\mathrm{RIP}_2$ in the optimal regime and instead one has to resort to $\mathrm{RIP}_1$ or use tools like the NSP 7. For example, it is known that matrices with independent and uniformly distributed entries on $\left\{0,1\right\}$ achieve optimal scaling KJ18 using the NSP, while adjacency matrices of expander graphs have the optimal scaling using $\mathrm{RIP}_1$. However, deterministic (i.e., non-random) construction of such matrices is an open problem.

As far as we know, the best deterministic construction (in terms of optimal scaling) of binary compressed sensing matrices (expander matrices) is GUV09. This construction has a near optimal scaling but is difficult to implement. Precisely, for any $\alpha >0$ there exists a constant $C_\alpha$ such that the number of rows of the matrix is:

$$\begin{equation} m\leq C_\alpha k^{1+\alpha }\left(\log _2(n)\log _2(k)\right)^{2+\frac{2}{\alpha }}. \tag{14} \end{equation}$$

The downside of this result is that the constant $C_\alpha$ is rather large for small $\alpha$.

For ease of implementation (i.e., practical reasons), we propose using the sub-optimal matrices with explicit constructions in LV20 and KS64, i.e., adjacency matrices of low-girth left-regular bipartite graphs and $k$-disjunct matrices respectively, in which the number of measurements $m$ is in the order of $k\sqrt {n}$. Actually, we established the equivalence of the two constructions in PBJ20, Theorem 4.3. Such $k$-disjunct matrices do not achieve the optimal rate for fixed $k$ as $n\rightarrow \infty$ according to DR82. However, for viral detection we have a fixed prevalence $p={k}/{n}$ in mind. The prevalence might be larger in the beginning of a pandemic and smaller when healthcare professionals are testing regularly and asymptotically, but for each of these applications we consider a fixed $p={k}/{n}$. In this scenario our result achieves the rate

$$\begin{equation} m = pn^\frac{3}{2}+n^\frac{1}{2}, \tag{15} \end{equation}$$

which can outperform results that achieve optimality according to DR82. For instance the construction of AA13, Corollary 1 achieves a number of measurements

$$\begin{equation} m \leq 2k^2\frac{\log _2(n)}{\log _2(k)} = 2p^2n^2\frac{\log _2(n)}{\log _2(n)+\log _2(p)}. {\tag{16}} \end{equation}$$

For a suitably chosen $n$, this upper bound and other constructions can be outperformed by the construction we propose.

Determining infections

The classical decoding procedure for disjunct matrices (known as COMP) iterates over all $j\in \left[n\right]$ and, if for a fixed $j$ and for all $i$ with $a_{ij}\neq 0$, the $i$th test is positive, it declares the $j$th individual as infected. This process is simple and, if $\mathbf{A}$ is scaled to be $k$-disjunct with column sums of $1$, there is no noise and there are no more than $k$ infected individuals, then it is guaranteed to find all infected individuals KS64. However, every false negative test will result in at least one individual that is falsely flagged as not infected. Thus, this decoding procedure is incredibly sensitive to noise. There are extensions to these matrices which tolerate a fixed number of errors under restrictive conditions AA13.

Compressed sensing, on the other hand, not only detects infected individuals but also estimates the viral load, which may have further benefits to the medical practitioners. The noise in 13 is nonzero in general, since the estimate is affected by some errors including, rounding errors and inaccuracy of the qPCR. In GAR$^{+}$20 the noise is modeled as a heavy-tailed random variable depending on the unknown $\mathbf{A}\mathbf{x}$.

Therefore, it is difficult to apply recovery methods from compressed sensing for independent additive noise out of the box. Parameter tuning, using, e.g., cross validation, is often crucial for most of these methods in these noise settings. Interestingly, non-negativity helps with such noise models. Combining the heavy-tailed noise model with the non-negativity of the viral load data (i.e., $\mathbf{x}$), we recommend using the parameter tuning free non-negative least absolute deviation (NNLAD), proposed in PBJ21, for recovery which is any minimizer

$$\begin{align} \hat{\mathbf{x}} = \operatorname *{argmin}_{\mathbf{z}\geq 0}\left\|\mathbf{A}\mathbf{z}-\mathbf{y}\right\|_1. \tag{NNLAD} \end{align}$$

In PBJ21, the authors showed that for certain matrices $\mathbf{A}$, for example what they refer to as random walk matrices of lossless expander graphs (which can be considered to include disjunct matrices), the convex NNLAD approach is indeed sparsity promoting and allows for an estimate of the form

$$\begin{equation} \left\|\mathbf{x}-\hat{\mathbf{x}}\right\|_1\leq C\left\|\mathbf{e}\right\|_1, \tag{17} \end{equation}$$

for some constant $C$ independent of $\mathbf{e},\mathbf{x},\mathbf{y}$, and $\hat{\mathbf{x}}$.

Finally, given a certain threshold $\epsilon$, one declares the individual $n$ to be infected if $\hat{x}_j>\epsilon$, for small $\epsilon >0$. If the noise level $\left\|\mathbf{e}\right\|_1$ is small enough, this method will guarantee that these are exactly the infected individuals. Thus, in this sense compressed sensing gives a nonadaptive testing procedure which provably finds $k$ infected individuals among $n$ with the scaling-optimal number of tests, and small errors would have no effect on the test result. The guarantees in PBJ21 follow from a self-regularization feature in the non-negative case which has been worked out already for other cases, like the non-negative least squares in KJ18. Note that a comprehensive comparative study of the performance of NNLAD vis-a-vis other traditional compressed sensing algorithms was conducted in PBJ21.

The take-home message of this section is that we suggest a quasi-optimal pooling procedure with an efficient noise-robust recovery algorithm. This is a more practical setup that can be implemented in a straight forward way at a medical lab doing COVID testing. In particular, for the measurement matrix/design we trade-off optimal scaling with ease of use; while for the decoder, i.e., NNLAD, we trade-off ease of decoding for noise robustness (when compared to COMP) and for no parameter tuning (when compared to many compressing decoders).

Discussion

Now we try to put the above theoretical results into the right perspective by doing a bit more detailed comparison of our proposed approach to other approaches proposed for group/pool testing. We also discuss some empirical results and then draw conclusions.

Nonadaptive methods

Consider the nonadaptive method from AA13 with number of tests less than $p^2n^2\frac{\log _2\left(n\right)}{\log _2\left(n\right)+\log _2\left(p\right)}$, see 16. Some empirical example comparisons are made in Table 1. In many scenarios, our method requires a similar number of tests as other methods from nonadaptive group testing and sometimes even requires fewer.

Table 1.

A comparison of the method of AA13 to ours looking at different prevalences and population sizes. We denote prevalence by $p$, number of infected individuals by $k$, population size by $n$, pool size by $Q$ and number of tests by $m$. Note that the rates ($m/n$) are based on approximations of the exact expected number of tests.

Work	$p$ ($k/n$)	$n$	$Q$	$m/n$
AA13	0.01	900	unknown	0.5573
Ours	0.01	900	30	0.3333
AA13	0.001	10000	unknown	0.0800
Ours	0.001	10000	100	0.1100

Adaptive methods

Adaptive methods outperform non-adaptive methods in general. Although nonadaptive, we compare our method to other adaptive methods in Table 2 below. Many adaptive group testing methods require half as many tests as our method, but ours has the advantage that it only requires one stage and can thus be performed faster.

Table 2.

A comparison of our method to other methods including Dor43MNB$^{+}$20 for different prevalences and population sizes. We denote prevalence by $p$, number of infected individuals by $k$, population size by $n$, pool size by $Q$ and number of tests by $m$. Note: (i) that the rates ($m/n$) are based on approximations of the exact expected number of tests; (ii) where these approximations can be computed without involving $n$, we put independent (abbreviated indep.) of n under the $n$ column.

Work	$p$ ($k/n$)	$n$	$Q$	$m/n$
Dor43	0.01	indep. of n	10	0.2000
MNB$^{+}$20	0.01	indep. of n	35	0.1168
Ours	0.01	900	30	0.3333
Dor43	0.001	indep. of n	32	0.0633
MNB$^{+}$20	0.001	indep. of n	350	0.0179
Ours	0.001	10000	100	0.1100

Error correcting properties

We demonstrate an advantage of the NNLAD approach over other compressed sensing decoders in a short experiment with synthetic data. One common error model is to use multiplicative noise, i.e., $y_m= g_m\times (\mathbf{A}\mathbf{x})_m$ for some random variable $g_m$. This yields that the magnitude of certain noise components are significantly larger than others and you end up with a peaky noise. This motivates us to model the additive noise vector $\mathbf{e}\coloneq \mathbf{y}-\mathbf{A}\mathbf{x}$ approximately as a sparse random vector. Following PBJ20, Theorem 4.5 for group sizes of $31$ and $k=7$, we are guaranteed to detect up to $7$ infected people among $961$ by using $248$ tests. However, empirically the identification will succeed even if more than $7$ individuals are infected and even if multiple measurements are corrupted. In Figure 3, we vary the prevalence $p=\frac{\left\|\mathbf{x}\right\|_0}{n}$ and the fraction of corrupted measurements $p_e=\frac{\left\|\mathbf{e}\right\|_0}{m}$ and plot the probability that the NNLAD estimator $\hat{\mathbf{x}}$ is sufficiently close to the true signal $\mathbf{x}$ in the $\ell _1$-norm.

Figure 3.

A phase-transition plot of the probability of recovery as a function of the prevalence and the number of corrupted measurements for pool sizes of $31$ and $k=7$.

We see that as guaranteed by PBJ20, Theorem 4.5 for $p=\frac{\left\|\mathbf{x}\right\|_0}{n}\leq \frac{7}{31^2}\approx 0.0073$ and $p_e=0$, the recovery succeeds. But empirically the recovery also succeeds for $p\leq 0.08$ and $p_e=0$, i.e., 10 times higher than guaranteed. This suggests that $k$ and thus $m$ could be reduced and still recover whenever $p\leq 0.02$ for instance, which might be sufficient for a prevalence of $p=0.01$.

Further the NNLAD seems to successfully recover even in the presence of sparse noise, i.e., when $p_e$ is small. This is not surprising as the NNLAD minimizes the $\ell _1$ norm of all possible noises and $\ell _1$-minimization is sparsity promoting CT05. Recovery seems to succeed whenever $p_e\leq 0.06$ and $\frac{4}{3}p_e+p\leq 0.08$. This error correcting property gives the NNLAD advantages over other compressed sensing decoders in the presence of heavy outliers.

However, the error correcting properties cannot be guaranteed uniformly for all $\mathbf{e}$ with $p_e\leq 0.06$. If the noise components are active exactly on the support of a column of $\mathbf{A}$ where the signal $\mathbf{x}$ is nonvanishing, the measurement might correspond to a different signal with the same support. Hence, recovery might fail as soon as $\left\|\mathbf{e}\right\|_0$ is at least as large as the number of nonzero entries in a column of $\mathbf{A}$, in this case $k+1=8$. However, there are $\binom{m}{k+1}=\binom{248}{8}\approx 3.16\cdot 10^{14}$ different possible supports of a $k+1$-sparse noise, but only $n=961$ of these appear as columns of the matrix. Thus, if the noise support is drawn uniformly at random, such an event is highly unlikely. Thus, we see that the NNLAD is empirically correcting more errors than expected.

Conclusion

We have explained how compressed sensing can be used to solve the viral detection problem. It generates a nonadaptive testing procedure. Further, we have presented a construction of design matrices that can be used for classical nonadaptive group testing and compressed sensing based viral detection. The construction requires roughly as many tests as other methods from nonadaptive group testing and can possibly outperform those. Adaptive group testing methods still require fewer tests but can only be performed sequentially which is a critical problem in a pandemic with an exponential growth rate. We have proposed to use the NNLAD as a compressed sensing decoder, since, compared to other compressed sensing decoders, it is robust against heavy-tailed noise and does not requires knowledge of noise level. In particular, it also has certain error correcting properties. Lastly our method also computes the viral load of infected individuals.