# K-Median Clustering, Model-Based Compressive Sensing, and Sparse Recovery
for Earth Mover Distance††thanks: This research has been supported in part by
the David and Lucille Packard Fellowship, MADALGO (the Center for Massive Data
Algorithmics, funded by the Danish National Research Association) and NSF
grant CCF-0728645. E. Price has been supported in part by an NSF Graduate
Research Fellowship.

Piotr Indyk  Eric Price

(24 April 2011)

###### Abstract

We initiate the study of sparse recovery problems under the Earth-Mover
Distance (EMD). Specifically, we design a distribution over $m\times n$
matrices $A$ such that for any $x$, given $Ax$, we can recover a $k$-sparse
approximation to $x$ under the EMD distance. One construction yields
$m=O(k\log(n/k))$ and a $1+\epsilon$ approximation factor, which matches the
best achievable bound for other error measures, such as the $\ell_{1}$ norm.

Our algorithms are obtained by exploiting novel connections to other problems
and areas, such as streaming algorithms for k-median clustering and model-
based compressive sensing. We also provide novel algorithms and results for
the latter problems.

##  1 Introduction

In recent years, a new “linear” approach for obtaining a succinct approximate
representation of $n$-dimensional vectors (or signals) has been discovered.
For any signal $x$, the representation is equal to $Ax$, where $A$ is an
$m\times n$ matrix, or possibly a random variable chosen from some
distribution over such matrices. The vector $Ax$ is often referred to as the
measurement vector or linear sketch of $x$. Although $m$ is typically much
smaller than $n$, the sketch $Ax$ often contains plenty of useful information
about the signal $x$.

A particularly useful and well-studied problem is that of stable sparse
recovery. The problem is typically defined as follows: for some norm
parameters $p$ and $q$ and an approximation factor $C>0$, given $Ax$, recover
an “approximation” vector $x^{*}$ such that

$\left\lVert x-x^{*}\right\rVert_{p}\leq C\min_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{q}$ (1)

where we say that $x^{\prime}$ is $k$-sparse if it has at most $k$ non-zero
coordinates. Sparse recovery has applications to numerous areas such as data
stream computing [Mut05, Ind07] and compressed sensing [CRT06, Don06], notably
for constructing imaging systems that acquire images directly in compressed
form (e.g., [DDT+08, Rom09]). The problem has been a subject of extensive
study over the last few years, with the goal of designing schemes that enjoy
good “compression rate” (i.e., low values of $m$) as well as good algorithmic
properties (i.e., low encoding and recovery times). It is known by now111In
particular, a random Gaussian matrix [CRT06] or a random sparse binary matrix
([BGI+08], building on [CCFC02, CM04, CM06]) has this property with
overwhelming probability. See [GI10] for an overview. that there exist
matrices $A$ and associated recovery algorithms that produce approximations
$x^{*}$ satisfying Equation (1) with $\ell_{p}=\ell_{q}=\ell_{1}$, constant
approximation factor $C$, and sketch length $m=O(k\log(n/k))$; it is also
known that this sketch length is asymptotically optimal [DIPW10, FPRU10].
Results for other combinations of $\ell_{p}$/$\ell_{q}$ norms are known as
well.

However, limiting the error measures to variants of $\ell_{p}$ norms is quite
inconvenient in many applications. First, the distances induced by $\ell_{p}$
norms are typically only quite raw approximations of the perceptual
differences between images. As a result, in the field of computer vision,
several more elaborate notions have been proposed (e.g., in [RTG00, Low04,
Lyu05, GD05]). Second, there are natural classes of images for which the
distances induced by the $\ell_{p}$ norm are virtually meaningless. For
example, consider images of “point clouds”, e.g., obtained via astronomical
imaging. If we are given two such images, where each point in the second image
is obtained via small random translation of a point in the first image, then
the $\ell_{p}$ distance between the images will be close to the largest
possible, even though the images are quite similar to each other.

Motivated by the above considerations, we initiate the study of sparse
recovery under non-$\ell_{p}$ distances. In particular, we focus on the Earth-
Mover Distance (EMD) [RTG00]. Informally, for the case of two-dimensional
$\Delta\times\Delta$ images (say, $x,y:[\Delta]^{2}\to\mathbb{R}_{+}$) which
have the same $\ell_{1}$ norm, the EMD is defined as the cost of the min-cost
flow that transforms $x$ into $y$, where the cost of transporting a unit of
mass from a pixel $p\in[\Delta]^{2}$ of $x$ to a pixel $q\in[\Delta]^{2}$ of
$y$ is equal to the $\ell_{1}$ distance222One can also use the $\ell_{2}$
distance. Note that the two distances differ by at most a factor of $\sqrt{2}$
for two-dimensional images. between $p$ and $q$. The EMD metric can be viewed
as induced by a norm $\left\lVert\cdot\right\rVert_{EMD}$, such that
$\text{EMD}(x,y)=\left\lVert x-y\right\rVert_{EMD}$; see Section 2 for a
formal definition. Earth-Mover Distance and its variants are popular metrics
for comparing similarity between images, feature sets, etc. [RTG00, GD05].

Results. In this paper we introduce three sparse recovery schemes for the
Earth-Mover Distance. Each scheme provides a matrix (or a distribution of
matrices) $A$, with $m$ rows and $n$ columns for $n=\Delta^{2}$, such that for
any vector $x$, given $Ax$, one can reconstruct a vector $x^{*}$ such that

$\left\lVert x-x^{*}\right\rVert_{EMD}\leq C\min_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}.$ (2)

for some approximation factor $C>0$. We call any recovery scheme satisfying
Equation (2) an _EMD/EMD recovery scheme_. If $A$ is a distribution over
matrices (that is, the scheme is randomized), the guarantee holds with some
probability. The other parameters of the constructions are depicted in Figure
1.

Determinism | Sketch length | Decode time | Approx.
---|---|---|---
Deterministic | $k\log n\log(n/k)$ | $n\log^{O(1)}n$ | $\epsilon$
Deterministic | $k\log(n/k)$ | $n^{O(1)}$ | $\sqrt{\log(n/k)}$
Randomized | $k\log(n/k)$ | $k\log(n/k)$ | $\epsilon$
Figure 1: Our results

In particular, two of our constructions yield sketch lengths $m$ bounded by
$O(k\log(n/k))$, which mimics the best possible bound achievable for sparse
recovery in the $\ell_{1}$ distance [DIPW10]. Note, however, that we are not
able to show a matching lower bound for the EMD case.

#### Connections and applications

What does sparse recovery with respect to the Earth-Mover Distance mean?
Intuitively, a sparse approximation under EMD yields a short “signature” of
the image $x$ that approximately preserves its “EMD properties”. For example,
if $x$ consists of a small number of sparse point clouds (e.g., as in
astronomical imaging), sparse approximation of $x$ will approximately identify
the locations and weights of the clouds. Our preliminary experiments with a
heuristic algorithm for such data [GIP10] show that this approach can yield
substantial improvements over the usual sparse recovery.

Another application [RTG00] stems from the original paper, where such short
signatures were constructed333In fact, the algorithm in [RTG00] vaguely
resembles our approach, in that it uses a kd-tree decomposition to partition
the images. for general images, to extract their color or texture information.
The images were then replaced by their signatures during the experiments,
which significantly reduced the computation time.

The above intuitions can be formalized as follows. Let $x^{\prime}$ be the
minimizer of $\left\lVert x-x^{\prime}\right\rVert_{EMD}$ over all $k$-sparse
vectors. Then one can observe that the non-zero entries of $x^{\prime}$
correspond to the cluster centers in the best $k$-median444For completeness,
in our context the $k$-median is defined as follows. First, each pixel
$p\in[\Delta]^{2}$ is interpreted as a point with weight $x_{p}$. The goal is
to find a set $C\subset[n]^{2}$ of $k$ “medians” that minimizes the objective
function $\sum_{p\in[n]^{2}}\min_{c\in C}\left\lVert
p-c\right\rVert_{1}x_{p}$. clustering of $x$. Moreover, for each such center
$c$, the value of $x^{\prime}_{c}$ is equal to the total weight of pixels in
the cluster centered at $c$. Thus, a solution to the $k$-median problem
provides a solution to our sparse recovery problem as well555If the algorithm
reports both the medians and the weights of clusters..

There has been prior work on the $k$-median problem in the streaming model
under insertions and deletions of points [FS05, Ind04]. Such algorithms
utilize linear sketches, and therefore implicitly provide schemes for
approximating the $k$-median of $x$ from a linear sketch of $x$ (although they
do not necessarily provide the cluster weights, which are needed for the
sparse recovery problem). Both algorithms666The paper [Ind04] claims
$m=k\log^{O(1)}n$. Unfortunately, that is an error, caused by ignoring the
dependencies between the queries and their answers provided by the randomized
data structure MediEval. Fixing this problem requires reducing the probability
of failure of the algorithm so that it is inversely exponential in $k$, which
yields another factor of $k$ in the space bound. yield a method for
approximating the $k$-median from $\Omega(k^{2}\log^{O(1)}n)$ measurements,
with the algorithm of [FS05] providing an approximation factor of
$1+\epsilon$. In contrast, our result achieves an approximation factor of
$1+\epsilon$ with a sketch length $m$ that is as low as $O(k\log(n/k))$.

Thanks to this connection, our results also yield short sketches for the
$k$-median problem. Although the solution $x^{*}$ output by our algorithm does
not have to be $k$-sparse (i.e., we might output more than $k$ medians), one
can post-process the output by computing the best $k$-sparse approximation to
$x^{*}$ using any off-the-shelf (weighted) $k$-median algorithm (e.g.,
[HPM04])). This reduces the number of clusters to $k$, while (by the triangle
inequality of EMD) multiplying the approximation factor by a constant that
depends on the approximation constant of the chosen $k$-median algorithm. See
Appendix C for more details.

#### Techniques

On a high level, our approach is to reduce the sparse recovery problem under
EMD to sparse recovery under $\ell_{1}$. This is done by constructing a linear
mapping $P$ that maps $\mathbb{R}^{[\Delta]^{2}}$ into some space
$\mathbb{R}^{t}$, that has the property that a “good” sparse approximation to
$y=Px$ under $\ell_{1}$ yields a “good” sparse approximation to $x$ under EMD.
777We note that the aforementioned k-median algorithms implicitly rely on some
form of sparse recovery (e.g., see Remark 3.10 in [FS05] or remarks before
Theorem 5 in [Ind04]). However, the bounds provided by those algorithms fall
short of what we aim for. The list of formal constraints that such a mapping
needs to satisfy are given in Section 3. For concreteness, we define one such
mapping below; another one is given in Section 7. Specifically, the _pyramid_
mapping $P$ [IT03, GD05] (building on [Cha02, AV99]) is defined as follows.
First we impose $\log\Delta+1$ nested grids $G_{i}$ on $[\Delta]^{2}$, with
$G=\bigcup G_{i}$. For each _level_ $i=0\ldots l$, $l=\log_{2}\Delta$, the
grid $G_{i}$ is a partition of the image into _cells_ of side length $2^{i}$.
The cells in the grids can be thought of as forming a $4$-ary tree, with each
node $c$ at level $i$ having a set $C(c)$ of children at level $i-1$. For each
$i$, we define a mapping $P_{i}$ such that each entry in $P_{i}x$ corresponds
to a cell $c$ in $G_{i}$, and its value is equal to the sum of coordinates of
$x$ falling into $c$. The final mapping $P$ is defined as

$\displaystyle Px=[2^{0}P_{0}x,2^{1}P_{1}x,\ldots,2^{l}P_{l}x]$ (3)

It is easy to see that, for a vector $x$ that is $k$-sparse, the vector $Px$
is $O(K)$ sparse for $K=kl$. We also show that for any $x$, there exists an
$O(K)$-sparse $y$ such that the difference $\|y-Px\|_{1}$ is comparable to
$\min_{k\text{-sparse }x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}$.
We then find a good approximation $x^{*}$ to $x$ (in the EMD norm) by
“inverting” $P$ on $y$. Since we can recover an $O(K)$-sparse approximation to
$y$ (in the $\ell_{1}$ norm) from a sketch of length $O(K\log(n/K))$, we
obtain the first result from Figure 1.

To improve the sketch length we exploit the particular properties of the
mapping $P$ to recover an $O(K)$-sparse approximation from only $O(K)$
measurements. For any non-negative vector $x$, the coordinates of $Px$ have
the following hierarchical structure: (i) the coordinates are organized into
an $r$-ary tree for $r=4$, and (ii) the value of each internal node is non-
negative and equal to the sum of its children times two. Using one or both of
these properties enables us to reduce the number of measurements.

The second algorithm from Figure 1 is obtained using the property (i) alone.
Specifically, the problem of recovering a sparse approximation whose support
forms a tree has been well-studied in signal processing (the question is
motivated by an empirical observation that large wavelet coefficients tend to
co-occur in this fashion). In particular, the insightful paper [BCDH10] on
model-based compressive sensing (see Section 5 for an overview) gave a
deterministic scheme that recovers such approximation from a sketch of length
$O(K)$. Although the setup given in that paper is somewhat different from what
we need here, we show that one can modify and re-analyze their scheme to
achieve the desired guarantee. This approach, however, leads to an
approximation factor of $O(\sqrt{\log(n/k)})$.

In order to achieve a constant approximation factor, we employ both properties
(i) and (ii), as well as randomization. Specifically, we recover the tree
coefficients top-down, starting from the root of the tree. This is done in a
greedy manner: we only recurse on the children of nodes that are estimated to
be “heavy”. This first pass identifies a superset $S$ of the locations where
$Px$ is large, but estimates some of the values $(Px)_{S}$ quite poorly. The
set of locations $S$ has $\left|S\right|=O(K)$, so we can recover $(Px)_{S}$
accurately with $O(K)$ measurements using the set query sketches of [Pri11].

Finally, we show that we can achieve the first and second result in Figure 1
by replacing the pyramid mapping by a variant of an even more basic transform,
namely the (two-dimensional) Haar wavelet mapping. Our variant is obtained by
rescaling the original Haar wavelet vectors using exponential weights, to
mimic the pyramid scheme behavior. This result relates the two well-studied
notions (EMD and wavelets) in a somewhat unexpected way. As a bonus, it also
simplifies the algorithms, since inverting the wavelet mapping can now be done
explicitly and losslessly.

##  2 Preliminaries

#### Notation

We use $[n]$ to denote the set $\\{1\ldots n\\}$. For any set $S\subset[n]$,
we use $\overline{S}$ to denote the complement of $S$, i.e., the set
$[n]\setminus S$. For any $x\in\mathbb{R}^{n}$, $x_{i}$ denotes the $i$th
coordinate of $x$, and $x_{S}$ denotes the vector
$x^{\prime}\in\mathbb{R}^{n}$ given by $x^{\prime}_{i}=x_{i}$ if $i\in S$, and
$x^{\prime}_{i}=0$ otherwise. We use $\operatorname{supp}(x)$ to denote the
support of $x$. We use $\mathbb{R}^{[\Delta]^{2}}$ to denote the set of
functions from $[\Delta]\times[\Delta]$ to $\mathbb{R}$; note that
$\mathbb{R}^{[\Delta]^{2}}$ can be identified with $\mathbb{R}^{n}$ since
$n=\Delta^{2}$. We also use $\mathbb{R}_{+}$ to denote $\\{x\in\mathbb{R}\mid
x\geq 0\\}$.

#### EMD

Consider any two non-negative vectors $x,y\in\mathbb{R}_{+}^{[\Delta]^{2}}$
such that $\|x\|_{1}=\|y\|_{1}$. Let $\Gamma(x,y)$ be the set of functions
$\gamma:[\Delta]^{2}\times[\Delta]^{2}\to\mathbb{R}_{+}$, such that for any
$i,j\in[\Delta]^{2}$ we have $\sum_{l}\gamma(i,l)=x_{i}$ and
$\sum_{l}\gamma(l,j)=y_{j}$; that is, $\Gamma$ is the set of possible “flows”
from $x$ to $y$. Then we define

$\text{EMD}^{*}(x,y)=\min_{\gamma\in\Gamma}\sum_{i,j\in[\Delta]^{2}}\gamma(i,j)\|i-j\|_{1}$

to be the min cost flow from $x$ to $y$, where the cost of an edge is its
$\ell_{1}$ distance. This induces a norm $\left\lVert\cdot\right\rVert_{EMD}$
such that $\left\lVert x-y\right\rVert_{EMD}=\text{EMD}^{*}(x,y)$. For general
vectors $w$,

$\left\lVert w\right\rVert_{EMD}=\min_{\begin{subarray}{c}x-y+z=w\\\
\left\lVert x\right\rVert_{1}=\left\lVert y\right\rVert_{1}\\\ x,y\geq
0\end{subarray}}\text{EMD}^{*}(x,y)+D\left\lVert z\right\rVert_{1}$

where $D=2\Delta$ is the diameter of the set $[\Delta]^{2}$. That is,
$\left\lVert w\right\rVert_{EMD}$ is the min cost flow from the positive
coordinates of $w$ to the negative coordinates, with some penalty for
unmatched mass.

#### Signal models

The basic idea of the signal models framework of [BCDH10] is to restrict the
sparsity patterns of the approximations. For some sparsity parameter888We use
$K$ to denote the sparsity in the context of model-based recovery (as opposed
to $k$, which is used in the context of “standard” recovery). $K$ let ${\cal
S}_{K}$ be a family of subsets of $[n]$ such that for each $S\in{\cal S}_{K}$
we have $|S|\leq K$. The family ${\cal S}_{K}$ induces a signal model ${\cal
M}_{K}\subset\mathbb{R}^{n}$ where

${\cal M}_{K}=\\{x\in\mathbb{R}^{n}\mid\operatorname{supp}(x)\subseteq
S\text{\ for some\ }S\in{\cal S}_{K}\\}.$

Note that ${\cal M}_{K}$ is a union of $|{\cal S}_{K}|$ subspaces, each of
dimension at most $K$. The signals in ${\cal M}_{K}$ are called ${\cal
M}_{K}$-sparse.

The following two examples of signal models are particularly relevant to our
paper:

  1. 1.

General $k$-sparse signals, where ${\cal S}_{k}$ contains all $k$-subsets of
$[n]$. In this case the induced signal model (denoted by $\Sigma_{k}$)
contains all $k$-sparse signals.

  2. 2.

Tree sparse signals. In this case, we assume that $n=\frac{c^{l}-1}{c-1}$ for
some (constant) integer $c$ and parameter $l$, and associate each $i\in[n]$
with a node of a full $c$-ary tree $T(c,l)$ of depth $l$. The family ${\cal
S}_{K}$ contains all sets $S$ of size up to $K$ that are connected in $T(c,l)$
and contain the root (so each $S$ corresponds to a graph-theoretic subtree of
$T(c,l)$). The induced signal model is denoted by ${\cal T}^{c}_{K}$, or
${\cal T}_{K}$ for short.999We note that technically this model was originally
defined with respect to the wavelet basis (as opposed to the standard basis
here) and for $c=2$. We adapt that definition to the needs in our paper.

In order to facilitate signal recovery, one often needs to consider the
differences $x-y$ of two signals $x\in{\cal M}$, $y\in{\cal M}^{\prime}$. For
this purpose we define the Minkowski sum of ${\cal M}_{K}$ and ${\cal
M}^{\prime}_{K}$ as ${\cal M}_{K}\oplus{\cal M}^{\prime}_{K}=\\{x+y:x\in{\cal
M}_{K},y\in{\cal M}^{\prime}_{K}\\}$. To simplify the notation, we define
${\cal M}^{(t)}$ to the $t$-wise Minkowski sum of ${\cal M}_{K}$. For all
signal models considered in this paper, we have ${\cal
M}_{K}^{(t)}\subset{\cal M}_{Kt}$.

Restricting sparsity patterns enables to recover sparse approximations from
shorter sketches. We defer a more thorough overview of the results to Section
5.

#### Assumptions

We assume that the sparsity parameters $k$ (and $K$, where applicable) are
smaller than $n/2$. Note that if this assumption does not hold, the problem
becomes trivial, since one can define the measurement matrix $A$ to be equal
to the identity matrix.

##  3 Framework for EMD-sparse recovery

In this section we describe our approach to reducing sparse recovery under EMD
into sparse recovery under $\ell_{1}$. We need the following three components:
(i) a $t\times n$ matrix $B$ (that will be used to map the EMD space into the
$\ell_{1}$ space); (ii) a signal model $\mathcal{M}\subset\mathbb{R}^{t}$; and
(iii) an $\ell_{1}/\ell_{1}$ recovery scheme for ${\cal M}$. The latter
involves an $m\times t$ matrix $A^{\prime}$ (or a distribution over such
matrices) such that, for any $x\in\mathbb{R}^{t}$, given $A^{\prime}x$, one
can recover $x^{*}$ such that

$\left\lVert x-x^{*}\right\rVert_{1}\leq C^{\prime}\min_{x^{\prime}\in{\cal
M}}\left\lVert x-x^{\prime}\right\rVert_{1}$ (4)

for an approximation factor $C^{\prime}$. If $A^{\prime}$ is a distribution
over matrices, we require that the guarantee holds with some constant
probability, e.g., 2/3.

The mapping $B$ must satisfy the following three properties:

  1. A.

(EMD-to-$\ell_{1}$ expansion.) For all $x\in\mathbb{R}^{n}$,

$\left\lVert x\right\rVert_{EMD}\leq\left\lVert Bx\right\rVert_{1}.$

  2. B.

(Model-alignment of EMD with ${\cal M}$.) For all $x\in\mathbb{R}_{+}^{n}$,
there exists a $y\in\mathcal{M}$ with

$\left\lVert y-Bx\right\rVert_{1}\leq\epsilon\min_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}.$

  3. C.

(Invertibility.) There is an efficient algorithm
$\mathcal{B}^{-1}\colon\mathbb{R}^{t}\to\mathbb{R}^{n}$ such that, for some
constant $D$ and all $y\in\mathbb{R}^{t}$,

$\left\lVert y-B\mathcal{B}^{-1}(y)\right\rVert_{1}\leq
D\min_{x\in\mathbb{R}^{n}}\left\lVert y-Bx\right\rVert_{1}.$

######  Lemma 3.1.

Consider $B,A^{\prime},{\cal M}$ satisfying the above properties. Then the
matrix $A=A^{\prime}B$ supports $k$-sparse recovery for EMD (as defined in
Equation (2)) with approximation factor $C=(1+D)C^{\prime}\epsilon$.

###### Proof.

Consider the recovery of any vector $x\in\mathbb{R}_{+}^{n}$. Let

$E=\min_{k\text{-sparse }x^{\prime}}\left\lVert
x-x^{\prime}\right\rVert_{EMD}.$

By Property B, for any $x\in\mathbb{R}^{n}$, there exists a $y\in\mathcal{M}$
with

$\left\lVert y-Bx\right\rVert_{1}\leq\epsilon E.$

Hence our $\ell_{1}/\ell_{1}$ model-based recovery scheme for ${\cal M}$, when
run on $Ax=A^{\prime}Bx$, returns a $y^{*}$ with

$\left\lVert y^{*}-Bx\right\rVert_{1}\leq C^{\prime}\epsilon E.$

Let $x^{*}=\mathcal{B}^{-1}(y^{*})$. We have by Property C that

$\left\lVert y^{*}-Bx^{*}\right\rVert_{1}\leq
D\min_{x^{\prime}\in\mathbb{R}^{n}}\left\lVert
y^{*}-Bx^{\prime}\right\rVert_{1}\leq D\left\lVert
y^{*}-Bx\right\rVert_{1}\leq DC^{\prime}\epsilon E.$

Hence by Property A

$\displaystyle\left\lVert x^{*}-x\right\rVert_{EMD}$
$\displaystyle\leq\left\lVert B(x^{*}-x)\right\rVert_{1}\leq\left\lVert
y^{*}-Bx\right\rVert_{1}+\left\lVert y^{*}-Bx^{*}\right\rVert_{1}$
$\displaystyle\leq(1+D)C^{\prime}\epsilon E$

as desired. ∎

##  4 Pyramid transform

In this section we will show that the pyramid transform $P$ defined in
Equation (3) of Section 1 satisfies properties B and C of Section 3, with
appropriate parameters.

The property A has been shown to hold for $P$ in many other papers (e.g.,
[Cha02, IT03]). The intuition is that the weight of a cell is at least the
Earth-Mover Distance to move all mass in the cell from the center to any
corner of the cell, including the corner that is at the center of the parent
of the cell.

###  4.1 Model-alignment with tree sparsity

In this section we show Property B, where the signal model ${\cal M}$ is equal
to the $K$-tree-sparse model ${\cal T}_{K}$, for $K=O(k\log(n/k))$. In fact,
we show a stronger statement: the trees have their _width_ (the maximum number
of nodes per level) bounded by some parameter $s$. We will exploit the latter
property later in the paper.

######  Lemma 4.1.

For any $x\in\mathbb{R}_{+}^{n}$ there exists a tree $S\subset[t]$ of size $K$
and width $s$ with

$\left\lVert(Px)_{\overline{S}}\right\rVert_{1}\leq\epsilon\min_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}$

for $s=O\left(\frac{1}{\epsilon^{2}}k\right)$ and
$K=O(\frac{1}{\epsilon^{2}}k\log(n/k))$.

###### Proof.

Let $x^{\prime}=\operatorname*{arg\,min}_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}$ be the $k$-medians
approximation of $x$. Consider the cells that contain each point in the
support of $x^{\prime}$. For each such cell at any level $i$, add the
$O(\frac{1}{\epsilon^{2}})$ other cells of the same level within an $\ell_{1}$
distance of $\frac{2}{\epsilon}2^{i}$. The resulting $S$ has
$s=O\left(\frac{1}{\epsilon^{2}}k\right)$ cells per level, and all the
ancestors of any cell in the result also lie in $S$. So $S$ is a tree of width
$s$. It has $O(s)$ elements from the top $\log_{4}s$ levels, and $O(s)$
elements on each of the $\log_{4}t-\log_{4}s$ remaining levels, for a size
$K=O(s\log t/s)$. We will show that
$\left\lVert(Px)_{\overline{S}}\right\rVert_{1}$ is small.

Define $e_{i}$ for $i\in[\Delta]^{2}$ to be the elementary vector with a $1$
at position $i$, so $x_{i}=x\cdot e_{i}$. Suppose that the distance between
$i$ and the nearest center in $x^{\prime}$ is $v_{i}$. Then we have

$\displaystyle\left\lVert(Px)_{\overline{S}}\right\rVert_{1}$
$\displaystyle=\sum_{i\in[\Delta]^{2}}\left\lVert(Px_{i}e_{i})_{\overline{S}}\right\rVert_{1}=\sum_{i\in[\Delta]^{2}}\left\lVert(Pe_{i})_{\overline{S}}\right\rVert_{1}x_{i}$
$\displaystyle\left\lVert x-x^{\prime}\right\rVert_{EMD}$
$\displaystyle=\sum_{i\in[\Delta]^{2}}v_{i}x_{i}.$

so it is sufficient to show
$\left\lVert(Pe_{i})_{\overline{S}}\right\rVert_{1}\leq\epsilon v_{i}$ for any
$i$.

Let $h$ be the highest level such that $e_{i}$ is not contained in a cell at
level $h$ in $S$. If no such $h$ exists,
$\left\lVert(Pe_{i})_{\overline{S}}\right\rVert_{1}=0$. Otherwise,
$v_{i}\geq\frac{2}{\epsilon}2^{h}$, or else $S$ would contain $e_{i}$’s cell
in level $h$. But then

$\displaystyle\left\lVert(Pe_{i})_{\overline{S}}\right\rVert_{1}=\sum_{j=0}^{h}2^{j}$
$\displaystyle<2^{h+1}\leq\epsilon v_{i}$

as desired. ∎

######  Corollary 4.2.

For any $x\in\mathbb{R}_{+}^{n}$, there exists a $y\in\mathcal{T}_{K}$ with

$\left\lVert y-Px\right\rVert_{1}\leq\epsilon\min_{k\text{-sparse
}x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}.$

###  4.2 Invertibility

Given an approximation $b$ to $Px$, we would like to find a vector $y$ with
$\left\lVert b-Py\right\rVert_{1}$ small. Note that this task can be
formulated as a linear program, and therefore solved in time that is
polynomial in $n$. In Appendix A we show a much faster approximate algorithm
for this problem, needed for our fast recovery algorithm:

######  Lemma 4.3.

Given any approximation $b$ to $Px$, we can recover a $y$ in
$O(\left|\operatorname{supp}(b)\right|)$ time with

$\left\lVert Py-Px\right\rVert_{1}\leq 8\left\lVert b-Px\right\rVert_{1}.$

Recall that $P$ has $t=\left\lfloor 4n/3\right\rfloor$ rows. This means
standard $\ell_{1}/\ell_{1}$ $K$-sparse recovery for $Px$ is possible with
$m=O(K\log t/K)=O(\frac{1}{\epsilon^{2}}k\log^{2}(n/k))$. Hence by Lemma 3.1,
using $B=P$ and standard sparse recovery techniques on the model ${\cal
M}=\Sigma_{K}$ gives the first result in Figure 1:

######  Theorem 4.4.

There exists a deterministic EMD/EMD recovery scheme with
$m=O(\frac{1}{\epsilon^{2}}k\log^{2}(n/k))$ and $C=\epsilon$. Recovery takes
$O(n\log^{c}n)$ time for some constant $c$.

##  5 Tree-sparse recovery

To decrease the number of measurements required by our algorithm, we can use
the stronger signal model $\mathcal{T}_{K}$ instead of $\Sigma_{K}$. The paper
[BCDH10] gives an algorithm for model-based sparse recovery of
$\mathcal{T}_{K}$, but their theorem does not give an $\ell_{1}/\ell_{1}$
guarantee. In Appendix B we review the prior work and convert their theorem
into the following:

######  Theorem 5.1.

There exists a matrix $A$ with $O(K)$ rows and a recovery algorithm that,
given $Ax$, returns $x^{*}$ with

$\left\lVert x-x^{*}\right\rVert_{1}\leq
C\sqrt{\log(n/K)}\min_{x^{\prime}\in{\cal T}_{K}}\left\lVert
x-x^{\prime}\right\rVert_{1}$

for some absolute constant $C>1$. As long as the coefficients of $x$ are
integers bounded by $n^{O(1)}$, the algorithm runs in time
$O(K^{2}n\log^{c}n)$ for some constant $c$.

By Lemma 3.1, using this on $B=P$ and ${\cal M}=\mathcal{T}_{K}$ gives the
second result in Figure 1:

######  Theorem 5.2.

There exists a deterministic EMD/EMD recovery scheme with
$m=O(\frac{1}{\epsilon^{2}}k\log(n/k))$ and distortion
$C=O(\epsilon\sqrt{\log(n/k)})$. Recovery takes $O(k^{2}n\log^{c}n)$ time for
some constant $c$.

##  6 Beyond tree sparsity

The previous section achieved $O(\sqrt{\log n})$ distortion deterministically
with $O(k\log(n/k))$ rows. In this section, we improve the distortion to an
arbitrarily small constant $\epsilon$ at the cost of making the algorithm
randomized. To do this, we show that EMD under the pyramid transform is
aligned with a stronger model than just tree sparsity—the model can restrict
the values of the coefficients as well as the sparsity pattern. We then give a
randomized algorithm for $\ell_{1}/\ell_{1}$ recovery in this model with
constant distortion.

######  Definition 6.1.

Define $T_{K}^{s}$ to be the family of sets $S\subseteq[t]$ such that (i) $S$
corresponds to a connected subset of $G$ containing the root and (ii)
$\left|S\cap G_{i}\right|\leq s$ for all $i$. We say that such an $S$ is _$K$
-tree-sparse with width $s$_.

######  Definition 6.2.

Define $\mathcal{M}\subset\mathcal{T}_{K}$ as

$\mathcal{M}=\left\\{y\in\mathbb{R}^{t}\middle|\begin{array}[]{l}\operatorname{supp}(y)\subseteq
S\text{ for some $S\in T_{K}^{s}$, and}\\\ y_{i}\geq 2\left\lVert
y_{C(i)}\right\rVert_{1}\forall i\in[t]\end{array}\right\\}.$

where $s=O(\frac{1}{\epsilon^{2}}k)$ comes from Lemma 4.1.

Note that every $y\in\mathcal{M}$ is non-negative, and
$(Px)_{S}\in\mathcal{M}$ for all $x\in\mathbb{R}_{+}^{n}$. With Lemma 4.1,
this implies:

######  Lemma 6.3.

There is model-alignment of $P$ with ${\cal M}$, i.e., they satisfy Property
B.

We will give a good algorithm for $\ell_{1}/\ell_{1}$ recovery over
$\mathcal{M}$.

###  6.1 Randomized $\ell_{1}/\ell_{1}$ recovery of $\mathcal{M}$

######  Theorem 6.4.

There is a randomized distribution over $m\times t$ matrices $A$ with
$m=O(\frac{1}{\epsilon^{2}}k\log(n/k))$ and an algorithm that recovers $y^{*}$
from $Ay$ in $O(\frac{1}{\epsilon^{2}}k\log(n/k))$ time with

$\left\lVert y^{*}-y\right\rVert_{1}\leq
C\min_{y^{\prime}\in\mathcal{M}}\left\lVert y-y^{\prime}\right\rVert_{1}$

with probability $1-k^{-\Omega(1)}$, for some constant $C$. We assume
$k=\Omega(\log\log n)$.

We will give an algorithm to estimate the support of $y$. Given a sketch of
$y$, it recovers a support $S\in T_{K}^{2s}$ with

$\left\lVert y_{\overline{S}}\right\rVert_{1}\leq
10\min_{y^{\prime}\in\mathcal{M}}\left\lVert y-y^{\prime}\right\rVert_{1}.$

We can then use the set query algorithm [Pri11] to recover a $y^{*}$ from a
sketch of size $O(\left|S\right|)$ with

$\left\lVert y^{*}-y_{S}\right\rVert_{1}\leq\left\lVert
y_{\overline{S}}\right\rVert_{1}.$

Then

$\left\lVert y^{*}-y\right\rVert_{1}\leq\left\lVert
y^{*}-y_{S}\right\rVert_{1}+\left\lVert y-y_{S}\right\rVert_{1}\leq
2\left\lVert y_{\overline{S}}\right\rVert_{1}\leq
20\min_{y^{\prime}\in\mathcal{M}}\left\lVert y-y^{\prime}\right\rVert_{1}.$

as desired. Hence estimating the support of $y$ is sufficient.

###  6.2 Finding a good sparse support $S$ to $y$

Vectors $y^{\prime}\in\mathcal{M}$ have two properties that allow us to find
good supports $S\in T_{K}^{s}$ with constant distortion using only
$O(\left|S\right|)$ rows. First, $\operatorname{supp}(y^{\prime})$ forms a
tree, so the support can be estimated from the top down, level by level.
Second, each coefficient has value at least twice the sum of the values of its
children. This means that the cost of making a mistake in estimating the
support (and hence losing the entire subtree below the missing coefficient) is
bounded by twice the weight of the missing coefficient. As a result, we can
bound the global error in terms of the local errors made at each level.

Of course, $y$ may not be in $\mathcal{M}$. But $y$ is “close” to some
$y^{\prime}\in\mathcal{M}$, so if our algorithm is “robust”, it can recover a
good support for $y$ as well. Our algorithm is described in Algorithm 1.

Algorithm 1 Finding sparse support under $\mathcal{M}$

Definition of sketch matrix $A$. The algorithm is parameterized by a width
$s$. Let $h_{i}$ be a random hash function from $G_{i}$ to $O(s)$ for
$i\in[\log(n/s)]$. Then define $A^{\prime}(i)$ to be the
$O(s)\times\left|G_{i}\right|$ matrix representing $h_{i}$, so
$A^{\prime}(i)_{ab}=1$ if $h_{i}(b)=a$ and $0$ otherwise. Choose $A$ to be the
vertical concatenation of the $A^{\prime}(i)$’s.

Recovery procedure.

$\triangleright$ Find approximate support $S$ to $y$ from $b=Ay$

procedure FindSupport($A$, $b$)

$T_{\log(n/s)}\leftarrow G_{\log(n/s)}$ $\triangleright$
$\left|T_{\log(n/s)}\right|\leq 2s$

for $i=\log(n/s)-1\dotsc 0$ do

$\triangleright$ Estimate $y$ over $C(T_{i+1})$.

$y^{*}_{j}\leftarrow b_{h_{i}(j)}$ for $j\in C(T_{i+1})$.

$\triangleright$ Select the $2s$ largest elements of our estimate.

$T_{i}\leftarrow{\displaystyle\operatorname*{arg\,max}_{\begin{subarray}{c}T^{\prime}\subseteq
C(T_{i+1})\\\ \left|T^{\prime}\right|\leq 2s\end{subarray}}}\left\lVert
y^{*}_{T^{\prime}}\right\rVert_{1}$

end for

$S\leftarrow\displaystyle\bigcup_{i=0}^{\log(n/s)}T_{i}\cup\bigcup_{i\geq\log(n/s)}G_{i}$

end procedure

######  Lemma 6.5.

Algorithm 1 uses a binary sketching matrix of $O(s\log(n/s))$ rows and takes
$O(s\log(n/s))$ time to recover $S$ from the sketch.

###### Proof.

The algorithm looks at $O(\log(n/s))$ levels. At each level it finds the top
$2s$ of $4\times 2s$ values, which can be done in linear time. The algorithm
requires a sketch with $O(\log(n/s))$ levels of $O(s)$ cells each. ∎

The algorithm estimates the value of $y_{C(T_{i+1})}$ by hashing all of
$y_{G_{i}}$ into an $O(s)$ size hash table, then estimating $y_{j}$ as the
value in the corresponding hash table cell. Since $y$ is non-negative, this is
an overestimate. We would like to claim that the $2s$ largest values in our
estimate approximately contain the $s$ largest values in $y_{C(T_{i+1})}$. In
particular, we show that any $y_{j}$ we miss is either (i) not much larger
than $s$ of the coordinates we do output or (ii) very small relative to the
coordinates we already missed at a previous level.

######  Lemma 6.6.

In Algorithm 1, for every level $i$ let $w_{i}=\max_{q\in C(T_{i+1})\setminus
T_{i}}y_{q}$ denote the maximum value that is skipped by the algorithm and let
$f_{i}=\left\lVert y_{G_{i+1}\setminus T_{i+1}}\right\rVert_{1}$ denote the
error from coordinates not included in $T_{i+1}$. Let $c_{i}$ denote the
$s$-th largest value in $y_{T_{i}}$. Then with probability at least
$1-e^{-\Omega(s)}$, $w_{i}\leq\max\\{\frac{f_{i}}{4s},2c_{i}\\}$ for all
levels $i$.

###### Proof.

Define $s^{\prime}=8s\geq\left|C(T_{i+1})\right|$. We make the hash table size
at each level equal to $u=32s^{\prime}$. We will show that, with high
probability, there are at most $s$ coordinates $p$ where $y^{*}_{p}$ is more
than $f_{i}/s^{\prime}$ larger than $y_{p}$. Once this is true, the result
comes as follows: $y^{*}$ is an overestimate, so the top $2s$ elements of
$y^{*}$ contain at least $s$ values that have been overestimated by at most
$f_{i}/s^{\prime}$. Because the algorithm passes over an element of value
$w_{i}$, each of these $s$ values must actually have value at least
$w_{i}-f_{i}/s^{\prime}$. Hence either
$w_{i}<2f_{i}/s^{\prime}=\frac{f_{i}}{4s}$ or all $s$ values are at least
$w_{i}/2$.

To bound the number of badly overestimated coordinates, we split the noise in
two components: the part from $G_{i}\setminus C(T_{i+1})$ and the part from
$C(T_{i+1})$. We will show that, with probability $1-e^{-\Omega(s)}$, the
former is at most $f_{i}/s^{\prime}$ in all but $s/4$ locations and the latter
is zero in all but $3s/4$ locations.

WLOG we assume that the function $h_{i}$ is first fixed for $G_{i}\setminus
C(T_{i+1})$, then randomly chosen for $C(T_{i+1})$. Let $O_{i}\subset[u]$ be
the set of “overflow buckets” $l$ such that the sum $s_{l}=\sum_{p\notin
C(T_{i+1}),h_{i}(p)=l}y_{p}$ is at least $f_{i}/s^{\prime}$. By the definition
of $f_{i}$, $\sum_{l}s_{l}=f_{i}/2$, so

$|O_{i}|/u\leq\frac{f_{i}/2}{f_{i}/s^{\prime}}/u=1/2\frac{s^{\prime}}{32s^{\prime}}=1/64.$

Thus, the probability that a fixed child $q\in C(T_{i+1})$ is mapped to
$O_{i}$ is at most $1/64$. This is independent over $C(T_{i+1})$, so the
Chernoff bound applies. Hence with probability at least $1-e^{-\Omega(s)}$,
the number of $q\in C(T_{i+1})$ mapping to $O_{i}$ is at most twice its
expectation, or $\left|C(T_{i+1})\right|/32=s/4$.

We now bound the collisions within $C(T_{i+1})$. Note that our process falls
into the “balls into bins” framework, but for completeness we will analyze it
from first principles.

Let $Z$ be the number of cells in $C(T_{i+1})$ that collide. $Z$ is a function
of the independent random variables $h_{i}(p)$ for $p\in C(T_{i+1})$, and $Z$
changes by at most $2$ if a single $h_{i}(p)$ changes (because $p$ can cause
at most one otherwise non-colliding element to collide). Hence by McDiarmid’s
inequality,

$\Pr[Z\geq\operatorname{E}[Z]+t]\leq e^{-t^{2}/(2s^{\prime})}$

But we know that the chance that a specific $p$ collides with any of the
others is at most $s^{\prime}/u=1/32$. Hence $\operatorname{E}[Z]\leq
s^{\prime}/32$, and

$\Pr[Z\geq(\frac{1}{32}+\epsilon)s^{\prime}]\leq
e^{-\epsilon^{2}s^{\prime}/2}.$

By setting $\epsilon=2/32$ we obtain that, with probability $1-e^{-\Omega(s)}$
we have that $Z\leq\frac{3s^{\prime}}{32}=3s/4$.

Hence with probability $1-e^{-\Omega(s)}$, only $3s/4$ locations have non-zero
corruption from $C(T_{i+1})$, and we previously showed that with the same
probability only $s/4$ locations are corrupted by $f^{\prime}/s^{\prime}$ from
outside $C(T_{i+1})$. By the union bound, this is true for all levels with
probability at least $1-(\log n)e^{-\Omega(s)}=1-e^{-\Omega(s)}.$ ∎

######  Lemma 6.7.

Let $S$ be the result of running Algorithm 1 on $y\in\mathbb{R}^{t}$. Then

$\left\lVert y_{\overline{S}}\right\rVert_{1}\leq
10\min_{y^{\prime}\in\mathcal{M}}\left\lVert y-y^{\prime}\right\rVert_{1}$

with probability at least $1-e^{\Omega(s)}$.

###### Proof.

From the algorithm definition, $T_{i}=S\cap G_{i}$ for each level $i$. Let
$y^{\prime}\in\mathcal{M}$ minimize $\left\lVert
y-y^{\prime}\right\rVert_{1}$, and let $U=\operatorname{supp}(y^{\prime})$. By
the definition of $\mathcal{M}$, $U\in T_{K}^{s}$.

For each $i$, define $V_{i}=U\cap C(T_{i+1})\setminus T_{i}$ to be the set of
nodes in $U$ that could have been chosen by the algorithm at level $i$ but
were not. For $q\in U\setminus S$, define $R(q)$ to be the highest ancestor of
$q$ that does not lie in $S$; hence $R(q)$ lies in $V_{i}$ for some level $i$.
Then

$\displaystyle\left\lVert
y^{\prime}_{\overline{S}}\right\rVert_{1}=\left\lVert y^{\prime}_{U\setminus
S}\right\rVert_{1}$ $\displaystyle=\sum_{q\in U\setminus S}y^{\prime}_{q}$
$\displaystyle=\sum_{i}\sum_{p\in V_{i}}\sum_{R(q)=p}y^{\prime}_{q}$
$\displaystyle\leq\sum_{i}\sum_{p\in V_{i}}2y^{\prime}_{p}$
$\displaystyle=2\sum_{i}\left\lVert y^{\prime}_{V_{i}}\right\rVert_{1},$ (5)

where the inequality holds because each element of $y^{\prime}$ is at least
twice the sum of its children. Hence the sum of $y^{\prime}$ over a subtree is
at most twice the value of the root of the subtree.

Define the error term $f_{i}=\left\lVert y_{G_{i+1}\setminus
T_{i+1}}\right\rVert_{1}$, and suppose that the statement in Lemma 6.6
applies, as happens with probability $1-e^{\Omega(s)}$. Then for any level $i$
and $p\in V_{i}$, if $c_{i}$ is the $s$th largest value in $y_{T_{i}}$, then
$y_{p}\leq\max\\{f_{i}/4s,2c_{i}\\}$ or $y_{p}\leq\frac{f_{i}}{4s}+2c_{i}$.
Since $y_{T_{i}}$ contains at least $s$ values larger than $c_{i}$, and at
most $\left|U\cap T_{i}\right|=\left|U\cap
C(T_{i+1})\right|-\left|V_{i}\right|\leq s-\left|V_{i}\right|$ of them lie in
$U$, $y_{T_{i}\setminus U}$ must contain at least $\left|V_{i}\right|$ values
larger than $c_{i}$. This, combined with $\left|V_{i}\right|\leq s$, gives

$\displaystyle\left\lVert y_{V_{i}}\right\rVert_{1}\leq f_{i}/4+2\left\lVert
y_{T_{i}\setminus U}\right\rVert_{1}.$ (6)

Combining Equations (5) and (6), we get

$\displaystyle\left\lVert y^{\prime}_{\overline{S}}\right\rVert_{1}$
$\displaystyle\leq
2[\sum_{i}\left\lVert(y-y^{\prime})_{V_{i}}\right\rVert_{1}+\left\lVert
y_{V_{i}}\right\rVert_{1}]$ $\displaystyle\leq
2\left\lVert(y-y^{\prime})_{U}\right\rVert_{1}+\sum_{i}\left(4\left\lVert
y_{T_{i}\setminus U}\right\rVert_{1}+f_{i}/2\right)$ $\displaystyle\leq
2\left\lVert(y-y^{\prime})_{U}\right\rVert_{1}+4\left\lVert y_{S\setminus
U}\right\rVert_{1}+\left\lVert y_{\overline{S}}\right\rVert_{1}/2$
$\displaystyle=2\left\lVert(y-y^{\prime})_{U}\right\rVert_{1}+4\left\lVert(y-y^{\prime})_{S\setminus
U}\right\rVert_{1}+\left\lVert y_{\overline{S}}\right\rVert_{1}/2$
$\displaystyle\leq 4\left\lVert y-y^{\prime}\right\rVert_{1}+\left\lVert
y_{\overline{S}}\right\rVert_{1}/2.$

Therefore

$\displaystyle\left\lVert y_{\overline{S}}\right\rVert_{1}$
$\displaystyle\leq\left\lVert y-y^{\prime}\right\rVert_{1}+\left\lVert
y^{\prime}_{\overline{S}}\right\rVert_{1}$ $\displaystyle\leq 5\left\lVert
y-y^{\prime}\right\rVert_{1}+\left\lVert y_{\overline{S}}\right\rVert_{1}/2$
$\displaystyle\left\lVert y_{\overline{S}}\right\rVert_{1}$ $\displaystyle\leq
10\left\lVert y-y^{\prime}\right\rVert_{1}$

as desired. ∎

###  6.3 Application to EMD recovery

By Lemma 3.1 our $\ell_{1}/\ell_{1}$ recovery algorithm for ${\cal M}$ gives
an $\text{EMD}/\text{EMD}$ recovery algorithm.

######  Theorem 6.8.

Suppose $k=\Omega(\log\log n)$. There is a randomized EMD/EMD recovery scheme
with $m=O(\frac{1}{\epsilon^{2}}k\log(n/k))$, $C=\epsilon$, and success
probability $1-k^{-\Omega(1)}$. Recovery takes
$O(\frac{1}{\epsilon^{2}}k\log(n/k))$ time.

##  7 Wavelet-based method

We can also instantiate the framework of Section 3 using a reweighted Haar
wavelet basis instead of $P$ for the embedding $B$. We will have $\mathcal{M}$
be the tree-sparse model $\mathcal{T}_{O(\frac{1}{\epsilon^{2}}k\log{n/k})}$,
and use the $\ell_{1}/\ell_{1}$ recovery scheme of Section 5.

The details are deferred to Appendix D. We obtain an embedding $W$ defined by
a Haar transform $H$ (after rescaling the rows), and the following theorem:

######  Theorem 7.1.

There exists a matrix $A$ with $O(k\log(n/k))$ rows such that we can recover
$x^{*}$ from $Ax$ with

$\left\lVert x^{*}-x\right\rVert_{EMD}\leq
C\min_{y\in\mathcal{T}_{K}}\left\lVert Wx-y\right\rVert_{1}\leq
C\min_{k\text{-sparse }x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}$

for some distortion $C=O(\sqrt{\log(n/k)})$.

Note that if we ignore the middle term, this gives the same EMD/EMD result as
in Section 5. However the middle term may be small for natural images even if
the right term is not. In particular, it is well known that images tend to be
tree-sparse under $H$.

#### Acknowledgements

The authors would like to thank Yaron Rachlin from Draper Lab for numerous
conversations and the anonymous reviewers for helping clarify the
presentation.

## References

  * [AV99] P.K. Agarwal and K. Varadarajan.  Approximation algorithms for bipartite and non-bipartite matching in the plane.  SODA, 1999.
  * [BCDH10] R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde.  Model-based compressive sensing.  IEEE Transactions on Information Theory, 56, No. 4:1982–2001, 2010\.
  * [BGI+08] R. Berinde, A. Gilbert, P. Indyk, H. Karloff, and M. Strauss.  Combining geometry and combinatorics: a unified approach to sparse signal recovery.  Allerton, 2008.
  * [CCFC02] M. Charikar, K. Chen, and M. Farach-Colton.  Finding frequent items in data streams.  ICALP, 2002.
  * [CDDD01] A. Cohen, W. Dahmen, I. Daubechies, and R. DeVore.  Tree approximation and optimal encoding.  Applied and Computational Harmonic Analysis, 2001.
  * [Cha02] M. Charikar.  Similarity estimation techniques from rounding.  In STOC, pages 380–388, 2002.
  * [CIHB09] V. Cevher, P. Indyk, C. Hegde, and RG Baraniuk.  Recovery of clustered sparse signals from compressive measurements.  SAMPTA, 2009.
  * [CM04] G. Cormode and S. Muthukrishnan.  Improved data stream summaries: The count-min sketch and its applications.  Latin, 2004.
  * [CM06] G. Cormode and S. Muthukrishnan.  Combinatorial algorithms for compressed sensing.  Sirocco, 2006.
  * [CRT06] E. J. Candès, J. Romberg, and T. Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Comm. Pure Appl. Math., 59(8):1208–1223, 2006.
  * [DDT+08] M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk.  Single-pixel imaging via compressive sampling.  IEEE Signal Processing Magazine, 2008.
  * [DIPW10] K. Do Ba, P. Indyk, E. Price, and D. Woodruff.  Lower bounds for sparse recovery.  SODA, 2010.
  * [Don06] D. L. Donoho.  Compressed Sensing.  IEEE Trans. Info. Theory, 52(4):1289–1306, Apr. 2006.
  * [FPRU10] S. Foucart, A. Pajor, H. Rauhut, and T. Ullrich.  The gelfand widths of lp-balls for $0<p\leq 1$.  preprint, 2010.
  * [FS05] G. Frahling and C. Sohler.  Coresets in dynamic geometric data streams.  STOC, 2005.
  * [GD05] K. Grauman and T. Darrell.  The pyramid match kernel: Discriminative classification with sets of image features.  ICCV, 2005.
  * [GI10] A. Gilbert and P. Indyk.  Sparse recovery using sparse matrices.  Proceedings of IEEE, 2010.
  * [GIP10] R. Gupta, P. Indyk, and E. Price.  Sparse recovery for earth mover distance.  Allerton, 2010.
  * [HPM04] S. Har-Peled and S. Mazumdar.  Coresets for k-means and k-medians and their applications.  STOC, 2004.
  * [Ind04] P. Indyk.  Algorithms for dynamic geometric problems over data streams.  STOC, 2004.
  * [Ind07] P. Indyk.  Sketching, streaming and sublinear-space algorithms.  Graduate course notes, available at `http://stellar.mit.edu/S/course/6/fa07/6.895/`, 2007.
  * [IT03] P. Indyk and N. Thaper.  Fast color image retrieval via embeddings.  Workshop on Statistical and Computational Theories of Vision (at ICCV), 2003.
  * [Low04] D. Lowe.  Distinctive image features from scale-invariant keypoints.  International Journal of Computer Vision, 60(2):91–110, 2004.
  * [Lyu05] S. Lyu.  Mercel kernels for object recognition with local features.  CVPR, 2005.
  * [Mut05] S. Muthukrishnan.  Data streams: Algorithms and applications).  Foundations and Trends in Theoretical Computer Science, 2005.
  * [NT08] D. Needell and J. A. Tropp.  Cosamp: Iterative signal recovery from incomplete and inaccurate samples.  Arxiv math.NA 0803.2392, 2008.
  * [Pri11] E. Price.  Efficient sketches for the set query problem.  SODA, 2011.
  * [Rom09] J. Romberg.  Compressive sampling by random convolution.  SIAM Journal on Imaging Science, 2009.
  * [RTG00] Y. Rubner, C. Tomasi, and L. Guibas.  The earth mover’s distance as a metric for image retrieval.  International Journal of Computer Vision, 40(2):99–121, 2000.
  * [SDS95] E.J. Stollnitz, A.D. DeRose, and D.H. Salesin.  Wavelets for computer graphics: a primer.  Computer Graphics and Applications, 1995.

##  Appendix A Invertibility of Pyramid Transform

If $b$ were $(Px)_{S}$ for some $S$, then the problem would be fairly easy,
since $b$ tells us the mass $p_{q}$ in cells $q$ (in particular, if $q$ is at
level $i$, $p_{q}=\frac{b_{q}}{2^{i}}$). Define the _surplus_
$s_{q}=p_{q}-\sum_{r\in C(q)}p_{r}$ to be the mass estimated in the cell that
is not found in the cell’s children.

We start from the case when all surpluses are non-negative (as is the case for
$(Px)_{S}$). In this case, we can minimize $\left\lVert b-Py\right\rVert_{1}$
by creating $s_{q}$ mass anywhere in cell $q$.

Algorithm 1 Recovering $y$ from $b$ to minimize $\left\lVert
b-Py\right\rVert_{1}$ when all surpluses are non-negative.

For every cell $q\in G$, let $e_{q}\in\mathbb{R}^{n}$ denote an elementary
unit vector with the $1$ located somewhere in $q$ (for example, at the center
of $q$). Then return

$y=\sum_{q\in G}s_{q}e_{q}.$

######  Lemma A.1.

Suppose $b$ is such that $s_{q}\geq 0$ for all $q\in G$. Let $y$ be the result
of running Algorithm 1 on $b$. Then $y$ minimizes $\left\lVert
b-Py\right\rVert_{1}$.

###### Proof.

The vector $y$ has the property that $(Py)_{q}\geq b_{q}$ for all $q\in G$,
and for the root node $r$ we have $(Py)_{r}=b_{r}$. Because the weights are
exponential in the level value, any $y^{\prime}$ minimizing $\left\lVert
b-Py^{\prime}\right\rVert_{1}$ must have $(Py^{\prime})_{r}\geq b_{r}$, or
else increasing any coordinate of $y^{\prime}$ would decrease $\left\lVert
b-Py^{\prime}\right\rVert_{1}$. But then

$\displaystyle\left\lVert b-Py^{\prime}\right\rVert_{1}$
$\displaystyle=\sum_{i=0}^{\log\Delta}\sum_{q\in
G_{i}}\left|(Py^{\prime})_{q}-b_{q}\right|$
$\displaystyle\geq\sum_{i=0}^{\log\Delta}\sum_{q\in
G_{i}}(Py^{\prime})_{q}-b_{q}$
$\displaystyle=\sum_{i=0}^{\log\Delta}\left(2^{i-\log\Delta}(Py^{\prime})_{r}-\sum_{q\in
G_{i}}b_{q}\right)$
$\displaystyle=(2-2^{-\log\Delta})(Py^{\prime})_{r}-\left\lVert
b\right\rVert_{1}$ $\displaystyle\geq(2-2^{-\log\Delta})b_{r}-\left\lVert
b\right\rVert_{1}.$

Equality holds if and only if $(Py^{\prime})_{q}\geq b_{q}$ for all $q\in G$
and $(Py^{\prime})_{r}=b_{r}$. Since $y$ has these properties, $y$ minimizes
$\left\lVert b-Py\right\rVert_{1}$. ∎

Unfortunately, finding the exact solution is harder when some surpluses
$s_{q}$ may be negative. Then in order to minimize $\left\lVert
b-Py\right\rVert_{1}$ one must do a careful matching up of positive and
negative surpluses. In order to avoid this complexity, we instead find a
greedy 8-approximation. We modify $b$ from the top down, decreasing values of
children until all the surpluses are non-negative.

Algorithm 2 Modifying $b$ to form all non-negative surpluses

Perform a preorder traversal of $G$. At each node $q$ at level $i$, compute
the surplus $s_{q}$. If $s_{q}$ is negative, arbitrarily decrease $b$ among
the children of $q$ by a total of $2^{i-1}\left|s_{q}\right|$, so that $b$
remains non-negative.

######  Lemma A.2.

Suppose we run algorithm 2 on a vector $b$ to get $b^{\prime}$. Then

$\left\lVert b-b^{\prime}\right\rVert_{1}\leq 3\min_{y}\left\lVert
Py-b\right\rVert_{1}.$

###### Proof.

Let $y$ minimize $\left\lVert Py-b\right\rVert_{1}$. As with $Py^{\prime}$ for
any $y^{\prime}$, $Py$ has zero surplus at every node.

At the point when we visit a node $q$, we have updated our estimate of $b$ at
$q$ but not at its children. Therefore if $q$ is at level $i$ we compute
$s_{q}=\frac{1}{2^{i}}b^{\prime}_{q}-\frac{1}{2^{i-1}}\sum_{s\in C(q)}b_{s}$.
Then, because $Py$ has zero surplus,

$\displaystyle\left|s_{q}\right|$
$\displaystyle=\left|\frac{1}{2^{i}}b^{\prime}_{q}-\frac{1}{2^{i}}(Py)_{q}-\frac{1}{2^{i-1}}\sum_{s\in
C(q)}(b_{s}-(Py)_{s})\right|$
$\displaystyle\leq\frac{1}{2^{i}}\left|b^{\prime}_{q}-b_{q}\right|+\frac{1}{2^{i}}\left|b_{q}-(Py)_{q}\right|+\frac{1}{2^{i-1}}\sum_{s\in
C(q)}\left|b_{s}-(Py)_{s}\right|.$

Define $f_{i}=\sum_{q\in G_{i}}\left|b_{q}-(Py)_{q}\right|$ to be the original
$\ell_{1}$ error on level $i$, and $g_{i}=\sum_{q\in
G_{i}}\left|b^{\prime}_{q}-b_{q}\right|$ to be a bound on the amount of error
we add when running the algorithm. Because we only modify values enough to
rectify the surplus of their parent, we have

$\displaystyle g_{i-1}$ $\displaystyle\leq 2^{i-1}\sum_{q\in
G_{i}}\left|s_{q}\right|$ $\displaystyle\leq\sum_{q\in
G_{i}}\frac{1}{2}\left|b^{\prime}_{q}-b_{q}\right|+\frac{1}{2}\left|b_{q}-(Py)_{q}\right|+\sum_{s\in
C(q)}\left|b_{s}-(Py)_{s}\right|$
$\displaystyle\leq\frac{1}{2}g_{i}+\frac{1}{2}f_{i}+f_{i-1}.$

Unrolling the recursion, we get

$\displaystyle g_{i}$ $\displaystyle\leq
f_{i}+\sum_{j=1}^{\log\Delta-i}\frac{1}{2^{j-1}}f_{i+j}$
$\displaystyle\left\lVert
b^{\prime}-b\right\rVert_{1}=\sum_{i=0}^{\log\Delta}g_{i}$
$\displaystyle\leq\sum_{i=0}^{\log\Delta}3f_{i}=3\left\lVert
Py-b\right\rVert_{1}$

as desired. ∎

This lets us prove Lemma 4.3.

  * Lemma 4.3.

Given any approximation $b$ to $Px$, running the previous two algorithms gives
a $y$ with

$\left\lVert Py-Px\right\rVert_{1}\leq 8\left\lVert b-Px\right\rVert_{1}$

in $O(\left|\operatorname{supp}(b)\right|)$ time.

###### Proof.

By running Algorithm 2 on $b$, we get $b^{\prime}$ with $\left\lVert
b-b^{\prime}\right\rVert_{1}\leq 3\left\lVert Px-b\right\rVert_{1}$. Then we
run Algorithm 1 on $b^{\prime}$ to get $y$ that minimizes $\left\lVert
Py-b^{\prime}\right\rVert_{1}$. Then

$\displaystyle\left\lVert Py-Px\right\rVert_{1}$ $\displaystyle\leq\left\lVert
Py-b^{\prime}\right\rVert_{1}+\left\lVert Px-b^{\prime}\right\rVert_{1}$
$\displaystyle\leq 2\left\lVert Px-b^{\prime}\right\rVert_{1}$
$\displaystyle\leq 2(\left\lVert Px-b\right\rVert_{1}+\left\lVert
b^{\prime}-b\right\rVert_{1})$ $\displaystyle\leq 8\left\lVert
Px-b\right\rVert_{1}.$

To bound the recovery time, note that after Algorithm 2 visits a node with
value $0$, it sets the value of every descendant of that node to $0$. So it
can prune its descent when it first leaves $\operatorname{supp}(b)$, and run
in $O(\left|\operatorname{supp}(b)\right|)$ time. Furthermore, this means
$\left|\operatorname{supp}(b^{\prime})\right|\leq\left|\operatorname{supp}(b)\right|$
and $\operatorname{supp}(b^{\prime})$ is a top-down tree. Hence Algorithm 1
can iterate through the support of $b^{\prime}$ in linear time. ∎

##  Appendix B Model-based compressive sensing

In this section we first provide a quick review of model-based sparse
recovery, including the relevant definitions, algorithms and their guarantees.
We then show how to augment the algorithm so that it provides the guarantees
that are needed for our EMD algorithms.

###  B.1 Background

#### Model-based RIP

Given a signal model ${\cal M}_{K}$, we can formulate the ${\cal
M}_{K}$-restricted isometry property (${\cal M}_{K}$-RIP) of an $m\times n$
matrix $A$, which suffices for performing sparse recovery.

######  Definition B.1.

A matrix $A$ satisfies the ${\cal M}_{K}$-RIP with constant $\delta$ if for
any $x\in{\cal M}_{K}$, we have

$(1-\delta)\left\lVert x\right\rVert_{2}\leq\left\lVert
Ax\right\rVert_{2}\leq(1+\delta)\left\lVert x\right\rVert_{2}$

It is known that random Gaussian matrices with $m=O(k\log(n/k))$ rows satisfy
the $\Sigma_{k}$-RIP (i.e., the “standard” RIP), with very high probability,
and that this bound cannot be improved [DIPW10]. In contrast, it has been
shown that in order to satisfy the ${\cal T}_{K}$-RIP, only $m=O(K)$ rows
suffice [BCDH10]. The intuitive reason behind this is that the number of
rooted trees of size $K$ is $2^{O(K)}$ while the number of sets of size $k$ is
$\binom{n}{k}=2^{\Theta(k\log(n/k))}$.

#### Algorithms

Given a matrix $A$ that satisfies the ${\cal M}_{K}$-RIP, one can show how to
recover an approximation to a signal from its sketch. The specific theorem
(proven in [BCDH10] and re-stated below) considers $\ell_{2}$ recovery of a
“noisy” sketch $Ax+e$, where $e$ is an arbitrary “noise” vector, while
$x\in{\cal M}_{K}$. In the next section we will use this theorem to derive an
$\ell_{1}$ result for a different scenario, where $x$ is an arbitrary vector,
and we are given its exact sketch $Ax$.

######  Theorem B.2.

Suppose that a matrix $A$ satisfies ${\cal M}^{(4)}_{K}$-RIP with constant
$\delta<0.1$. Moreover, assume that we are given a procedure that, given
$y\in\mathbb{R}^{n}$, finds $y^{*}\in{\cal M}_{K}$ that minimizes
$\|y-y^{*}\|_{2}$. Then there is an algorithm that, for any $x\in{\cal
M}_{K}$, given $Ax+e$, $e\neq 0$, finds $x^{*}\in{\cal M}_{K}$ such that

$\|x-x^{*}\|_{2}\leq C\|e\|_{2}$

for some absolute constant $C>1$. The algorithm runs in time
$O((n+T+MM)\log(\|x\|_{2}/\|e\|_{2}))$, where $T$ is the running time of the
minimizer procedure, and $MM$ is the time needed to perform the multiplication
of a vector by the matrix $A$.

Note that the algorithm in the theorem has a somewhat unexpected property: if
the sketch is nearly exact, i.e., $e\approx 0$, then the running time of the
algorithm becomes unbounded. The reason for this phenomenon is that the
algorithm iterates to drive the error down to $\left\lVert e\right\rVert_{2}$,
which takes longer when $e$ is small. However, as long as the entries of the
signals $x,x^{*}$ and the matrix $A$ have bounded precision, e.g., are
integers in the range $1,\dotsc,L$, one can observe that $O(\log L)$
iterations suffice.

The task of minimizing $\|y-y^{*}\|_{2}$ over $y^{*}\in{\cal M}_{K}$ can
typically be accomplished in time polynomial in $K$ and $n$. In particular,
for ${\cal M}_{K}={\cal T}_{K}$, there is a simple dynamic programming
algorithm solving this problem in time $O(k^{2}n)$. See, e.g., [CIHB09] for a
streamlined description of the algorithms for (a somewhat more general)
problem and references. For more mathematical treatment of tree
approximations, see [CDDD01].

The following lemma (from [NT08]) will help us bound the value of $\|e\|_{2}$.

######  Lemma B.3.

Assume that the matrix $A$ satisfies the (standard) $\Sigma_{s}$-RIP with
constant $\delta$. Then for any vector $z$, we have
$\|Az\|_{2}\leq\sqrt{1+\delta}(\|z_{S}\|_{2}+\|z\|_{1}/\sqrt{s})$, where $S$
is the set of the $s$ largest (in magnitude) coefficients of $z$.

For completeness, we also include a proof. It is different, and somewhat
simpler than the original one. Moreover, we will re-use one of the arguments
later.

###### Proof.

We partition the coordinates of $S$ into sets $S_{0},S_{1},$
$S_{2},\dotsc,S_{t}$, such that (i) the coordinates in the set $S_{j}$ are no
larger (in magnitude) than the coordinates in the set $S_{j-1}$, $j\geq 1$,
and (ii) all sets but $S_{t}$ have size $s$. We have

$\displaystyle\|Az\|_{2}$ $\displaystyle\leq$
$\displaystyle\sum_{j=0}^{t}\|Az_{S_{j}}\|_{2}$ $\displaystyle\leq$
$\displaystyle\sqrt{1+\delta}(\|z_{S_{0}}\|_{2}+\sum_{j=1}^{t}\|z_{S_{j}}\|_{2})$
$\displaystyle\leq$
$\displaystyle\sqrt{1+\delta}(\|z_{S_{0}}\|_{2}+\sum_{j=1}^{s}\sqrt{s}(\|z_{S_{j-1}}\|_{1}/s))$
$\displaystyle\leq$
$\displaystyle\sqrt{1+\delta}(\|z\|_{2}+\|z\|_{1}/\sqrt{s})$

∎

###  B.2 New result

We start from the following observation relating general sparsity and tree
sparsity. Consider $k$ and $K$ such that $K=c^{\prime}k\log(n/k)$ for some
constant $c^{\prime}$.

######  Claim B.4.

Assume $n=\frac{c^{l}-1}{c-1}$ for some (constant) integer $c$. Then there
exists a constant $c^{\prime}$ such that $\Sigma_{k}\subset{\cal T}_{K}$.

###### Proof.

It suffices to show that for any $S\subset[n]$ of size $k$ there exists a
rooted connected subset $T$ of $T(c,l)$ of size $K$ such that $S\subset T$.
The set $T$ is equal to $T^{\prime}\cup T^{\prime\prime}$, where (i)
$T^{\prime}$ consist of all nodes in the tree $T(c,l)$ up to level
$\lceil\log_{c}k\rceil$ and (ii) $T^{\prime\prime}$ consists of all paths from
the root to node $i$, for $i\in S$. Note that $|T^{\prime}|=O(k)$, and
$|T^{\prime\prime}\setminus T^{\prime}|=O(k(\log n-\log k))=O(k\log(n/k))$. ∎

This claim is used in the following way. As we will see later, in order to
provide the guarantee for recovery with respect to the model ${\cal T}_{K}$,
we will need to perform the recovery with respect to the model ${\cal
T}_{K}\oplus\Sigma_{k}$. From the claim it follows that we can instead perform
the recovery with respect to the model ${\cal T}^{(2)}_{K}\subset{\cal
T}_{2K}$.

Specifically, we show the following.

######  Theorem B.5.

Suppose that we are given a matrix and minimizer subroutine as in Theorem B.2
for ${\cal T}_{2K}$. Then, for any $x$, given the vector $Ax$, the
approximation $x^{*}$ computed by the algorithm in Theorem B.2 satisfies

$\|x-x^{*}\|_{1}\leq(1+2C\sqrt{(1+\delta)c^{\prime}\log(n/k)})\min_{x^{\prime}\in{\cal
T}_{K}}\|x-x^{\prime}\|_{1}$

###### Proof.

Let $x^{\prime}\in{\cal T}_{K}$ be the minimizer of $\|x-x^{\prime}\|_{1}$.
Let $T$ be a tree of size $K$ such that $x^{\prime}=x_{T}$, and define the
“$\ell_{1}$ approximation error”
$E=\|x-x^{\prime}\|_{1}=\|x_{\overline{T}}\|_{1}$

Let $P\subseteq\overline{T}$ be the set of the $k$ largest (in magnitude)
coordinates of $x_{\overline{T}}$. By Claim B.4 it follows that $P\subseteq
T^{\prime}$, for some $T^{\prime}\in{\cal T}_{K}$. Let $T^{\prime\prime}=T\cup
T^{\prime}$.

We decompose $Ax$ into
$Ax_{T^{\prime\prime}}+Ax_{\overline{T^{\prime\prime}}}=Ax_{T^{\prime\prime}}+e$.
Since $x_{T^{\prime\prime}}\in{\cal T}_{2K}$, by Theorem B.2 we have

$\displaystyle\|x_{T^{\prime\prime}}-x^{*}\|_{2}$ $\displaystyle\leq$
$\displaystyle C\|e\|_{2}$ (7)

Let $T^{*}$ be the support of $x^{*}$. Note that $|T^{*}|\leq 2K$.

Since $A$ satisfies the (standard) RIP of order $k$ with constant
$\delta=0.1$, by Lemma B.3 we have

$\left\lVert e\right\rVert_{2}\leq\sqrt{1+\delta}[\left\lVert
x_{S}\right\rVert_{2}+\left\lVert x_{\overline{T\cup
P}}\right\rVert_{1}/\sqrt{k}]$

where $S\subset\overline{T\cup P}$ is the set of the $k$ largest (in
magnitude) coordinates of $x_{\overline{T\cup P}}$. By the definition of $P$,
every coordinate of $\left|x_{S}\right|$ is not greater than the smallest
coordinate of $|x_{P}|$. By the same argument as in the proof of Lemma B.3 it
follows that $\left\lVert x_{S}\right\rVert_{2}\leq\left\lVert
x_{P}\right\rVert_{1}/\sqrt{k}$, so

$\displaystyle\left\lVert e\right\rVert_{2}$
$\displaystyle\leq\sqrt{(1+\delta)/k}\left\lVert
x_{\overline{T}}\right\rVert_{1}.$ (8)

We have

$\displaystyle\left\lVert x-x^{*}\right\rVert_{1}$
$\displaystyle=\left\lVert(x-x^{*})_{T^{\prime\prime}\cup
T^{*}}\right\rVert_{1}+\left\lVert(x-x^{*})_{\overline{T^{\prime\prime}\cup
T^{*}}}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
x_{T^{\prime\prime}}-x^{*}\right\rVert_{1}+\left\lVert x_{T^{*}\setminus
T^{\prime\prime}}\right\rVert_{1}+\left\lVert(x-x^{*})_{\overline{T^{\prime\prime}\cup
T^{*}}}\right\rVert_{1}$ $\displaystyle=\left\lVert
x_{T^{\prime\prime}}-x^{*}\right\rVert_{1}+\left\lVert
x_{\overline{T^{\prime\prime}}}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
x_{T^{\prime\prime}}-x^{*}\right\rVert_{1}+E$
$\displaystyle\leq\sqrt{4K}\left\lVert
x_{T^{\prime\prime}}-x^{*}\right\rVert_{2}+E$
$\displaystyle\leq\sqrt{4K}C\left\lVert e\right\rVert_{2}+E$
$\displaystyle\leq\sqrt{4K}C\sqrt{(1+\delta)/k}\left\lVert
x_{\overline{T}}\right\rVert_{1}+E$
$\displaystyle=(1+2C\sqrt{(1+\delta)K/k})E$
$\displaystyle=(1+2C\sqrt{(1+\delta)c^{\prime}\log(n/k)})E$

by Equations 7 and 8. ∎

##  Appendix C Strict sparse approximation

In this section we show how to reduce the sparsity of an approximation down to
$k$ for an arbitrary norm $\|\cdot\|$. This reduction seems folklore, but we
could not find an appropriate reference, so we include it for completeness.

Consider a sparse approximation scheme that, given $Ax$, returns (not
necessarily sparse) vector $x^{*}$ such that $\|x^{*}-x\|\leq
C\min_{k\text{-sparse }x^{\prime}}\|x^{\prime}-x\|$; let $x^{\prime}$ be the
the minimizer of the latter expression. Let $\hat{x}$ be the approximately
best k-sparse approximation to $x^{*}$, i.e., such that $\|\hat{x}-x^{*}\|\leq
C^{\prime}\min_{k\text{-sparse }x^{\prime\prime}}\|x^{\prime\prime}-x^{*}\|$;
let $x^{\prime\prime}$ be the minimizer of the latter expression. Note that
since $x^{\prime}$ is $k$-sparse, it follows that
$\|x^{\prime\prime}-x^{*}\|\leq\|x^{\prime}-x^{*}\|$.

######  Claim C.1.

We have

$\|\hat{x}-x\|\leq[(C^{\prime}+1)C+C^{\prime}]\|x^{\prime}-x\|$

###### Proof.

$\displaystyle\|\hat{x}-x\|$ $\displaystyle\leq$
$\displaystyle\|\hat{x}-x^{*}\|+\|x^{*}-x\|$ $\displaystyle\leq$
$\displaystyle C^{\prime}\|x^{\prime\prime}-x^{*}\|+\|x^{*}-x\|$
$\displaystyle\leq$ $\displaystyle C^{\prime}\|x^{\prime}-x^{*}\|+\|x^{*}-x\|$
$\displaystyle\leq$ $\displaystyle
C^{\prime}[\|x^{\prime}-x\|+\|x-x^{*}\|]+\|x^{*}-x\|$ $\displaystyle=$
$\displaystyle(C^{\prime}+1)\|x^{*}-x\|+C^{\prime}\|x^{\prime}-x\|$
$\displaystyle\leq$
$\displaystyle(C^{\prime}+1)C\|x^{\prime}-x\|+C^{\prime}\|x^{\prime}-x\|$
$\displaystyle=$ $\displaystyle[(C^{\prime}+1)C+C^{\prime}]\|x^{\prime}-x\|$

##  Appendix D Wavelet-based method

We start by recalling the definition of the non-standard two-dimensional Haar
wavelet basis (see [SDS95] for an overview). Let $H\in\mathbb{R}^{n\times n}$
be the matrix with rows corresponding to the basis vectors. We will define $H$
in terms of the grids $G_{i}$. The first row of $H$ has all coordinates equal
to $1/n$. The rest of $H$ consists of three rows for each cell $C\in G_{i}$
for $i\geq 1$. For each cell $C$, the corresponding rows contain zeros outside
of the coordinates corresponding to $C$. The entries corresponding to $C$ are
defined as follows: (i) one row has entries equal to $2^{-i}$ for each entry
corresponding to the left half of $C$ and equal to $-2^{-i}$ for each entry
corresponding to the right half of $C$; (ii) the second row has entries equal
to $2^{-i}$ for the top half of $C$ and to $-2^{-i}$ for the bottom half; (ii)
and the third row has entries equal to $2^{-i}$ for the top left and bottom
right quadrants of $C$, and equal to $-2^{-i}$ for the other two quadrants.

We define $W$ to transform into the same basis as $H$, but with rescaled basis
vectors. In particular, the basis vectors from level $i$ are smaller by a
factor of $2^{2i-2}$, so the non-zero entries have magnitude $2^{2-3i}$. This
is equivalent to changing the coefficients of the corresponding rows of $W$ to
be $2^{i-2}$ rather than $2^{-i}$. Similarly, we rescale the all-positive
basis vector to have coefficients equal to $1/n^{3}$. Then $W=DH$ for some
diagonal matrix $D$.

This rescaling is such that the columns of $W^{-1}$, call them $v_{i}$, all
have $\left\lVert v_{i}\right\rVert_{EMD}=1$. This is because the min-cost
matching moves each of $2^{2i}/2$ coefficients by $2^{i}/2$. So we have

$\displaystyle\left\lVert x\right\rVert_{EMD}$
$\displaystyle=\left\lVert\sum(Wx)_{i}v_{i}\right\rVert_{EMD}$
$\displaystyle\leq\sum\left\lVert(Wx)_{i}v_{i}\right\rVert_{EMD}$
$\displaystyle=\sum\left|(Wx)_{i}\right|=\left\lVert Wx\right\rVert_{1},$

which is Property A of the framework.

Property C is easy since $W$ has a known inverse (namely $H^{T}D^{-1}$),
giving $\left\lVert y-WW^{-1}y\right\rVert_{1}=0$ for all $y$. All that
remains to show is Property B.

######  Lemma D.1.

For all $x\in\mathbb{R}_{+}^{n}$, there exists a
$y\in\mathcal{T}_{O(\frac{1}{\epsilon^{2}}k\log(n/k))}$ with

$\left\lVert y-Wx\right\rVert_{1}\leq\epsilon\min_{k\text{-sparse
}x_{k}}\left\lVert x-x_{k}\right\rVert_{EMD}.$

###### Proof.

We will show this using Lemma 4.1 as a black box. We know there exists a
support $S$ of $Px$ corresponding to a tree of grid cells such that

$\left\lVert(Px)_{\overline{S}}\right\rVert_{1}\leq\epsilon\min_{k\text{-sparse
}x_{k}}\left\lVert x-x_{k}\right\rVert_{EMD}.$

Let $S^{\prime}$ be a support of $Wx$ that contains the all-constant basis
vector as well as, for each cell $C\in G_{i}$ in $S$ with $i\geq 1$, the three
coefficients in $Wx$ corresponding to $C$. Then $S^{\prime}$ is also a tree.

For any cell $C\in G_{i}$, let $u$ be the row in $P$ corresponding to $C$ and
$v$ be any of the three rows in $W$ corresponding to $C$. Then

$\left\lVert v\right\rVert_{\infty}=2^{i-2}=\frac{1}{4}\left\lVert
u\right\rVert_{\infty}.$

So the only difference between $v$ and $u$ is that (i) $v$ has one fourth the
magnitude in each coefficient and (ii) some coefficients of $v$ are negative,
while all of $u$ are positive. Hence for positive $x$, $\left|v\cdot
x\right|\leq\frac{1}{4}\left|u\cdot x\right|$. This gives

$\left\lVert(Wx)_{\overline{S}^{\prime}}\right\rVert_{1}\leq\frac{3}{4}\left\lVert(Px)_{\overline{S}}\right\rVert_{2}\leq\frac{3}{4}\epsilon\min_{k\text{-sparse
}x_{k}}\left\lVert x-x_{k}\right\rVert_{EMD}.$

as desired. ∎

######  Theorem D.2.

This gives

$\left\lVert x^{*}-x\right\rVert_{EMD}\leq
C\min_{y\in\mathcal{T}_{K}}\left\lVert Wx-y\right\rVert_{1}\leq
C\min_{k\text{-sparse }x^{\prime}}\left\lVert x-x^{\prime}\right\rVert_{EMD}$

for some distortion $C=O(\sqrt{\log(n/k)})$.

∎


AI

    Artificial Intelligence
HCI

    Human-Computer Interaction
IRR

    Inter-Rater Reliability
LLM

    Large Language Model
NLP

    Natural Language Processing
QDA

    Qualitative Data Analysis

S1pt

# Decoding Complexity: Exploring Human-AI Concordance in Qualitative Coding

Elisabeth Kirsten  0009-0003-0680-8916<EMAIL_ADDRESS>Max Planck
Institute for Security and PrivacyBochumGermany ,  Annalina Buckmann
0000-0002-7959-9743 Ruhr University BochumBochumGermany ,  Abraham Mhaidli
0000-0002-9519-245X Max Planck Institute for Security and PrivacyBochumGermany
and  Steffen Becker  0000-0001-7526-5597 Ruhr University BochumBochumGermany
Max Planck Institute for Security and PrivacyBochumGermany

###### Abstract.

Qualitative data analysis provides insight into the underlying perceptions and
experiences within unstructured data. However, the time-consuming nature of
the coding process, especially for larger datasets, calls for innovative
approaches, such as the integration of Large Language Models (LLMs). This
short paper presents initial findings from a study investigating the
integration of LLMs for coding tasks of varying complexity in a real-world
dataset. Our results highlight the challenges inherent in coding with
extensive codebooks and contexts, both for human coders and LLMs, and suggest
that the integration of LLMs into the coding process requires a task-by-task
evaluation. We examine factors influencing the complexity of coding tasks and
initiate a discussion on the usefulness and limitations of incorporating LLMs
in qualitative research.

Large Languge Models, Qualitative Data Analysis, Generative Pre-trained
Transformer

††ccs: Human-centered computing HCI design and evaluation methods††ccs:
Computing methodologies Natural language processing

##  1\. Introduction

Qualitative Data Analysis (QDA) is an approach used in Human-Computer
Interaction (HCI) research to analyze unstructured data, such as interview
data or responses to open-ended survey questions. A fundamental step in QDA is
coding, the systematic process of identifying, annotating, and categorizing
data segments based on recurring themes and patterns, and assigning according
categories or “codes” (Braun and Clarke, 2006; Flick, 2022; Bergin, 2018;
Kuckartz, 2019; Adams et al., 2008). Coding data can be a time-consuming
process, particularly when analyzing large volumes of data (e. g., tens of
thousands of open-ended survey responses). The emergence of Generative
Artificial Intelligence, and in particular Large Language Models (LLMs) (e.
g., GPT-4 (Achiam et al., 2023), Gemini (Team et al., 2023), or Llama 2
(Touvron et al., 2023)) could significantly speed up the coding process, yet
their value and validity needs to be evaluated. Large Language Models have
demonstrated remarkable performance in annotation tasks in zero-shot or few-
shot learning scenarios where no or little labeled data is given (Törnberg,
2023; Gilardi et al., 2023; Ziems et al., 2023). As these models become
integral to HCI research, it is imperative to critically evaluate their
potential benefits and risks, as well as methodological challenges associated
with the integration of LLMs. In this workshop submission, we present the
preliminary results of a study scrutinizing how LLMs perform at applying tags
for QDA in coding tasks of different complexity. Specifically, we explore LLM-
assisted QDA of a real-world dataset to understand how LLMs perform, how they
compare to human coders, and how their performance is affected by the type of
coding (e. g., semantic vs. latent themes (Braun and Clarke, 2006)). First, we
develop a strategy for prompting LLMs to apply tags for QDA. We then implement
this strategy with GPT-3.5 (Brown et al., 2020) and GPT-4 (Achiam et al.,
2023) on interview data and compare the models’ performance with human coding
in three tasks that require different levels of interpretation.

##  2\. Research Method

To evaluate the performance of LLMs in QDA, we provide the model with segments
of qualitative data from interview transcripts111The data for our experiments
is drawn from a separate interview study (n=47) focused on users’ perceptions
of digital security and privacy. The interview data is in German, and we
conduct all experiments using German data segments and prompts., a human-
generated codebook, and the instruction to assign zero, one, or more tags to
each data segment using codes from the codebook. We then calculate the Inter-
Rater Reliability (IRR) to evaluate the agreement of the resulting coded
datasets between LLMs and human coders. We analyze three distinct themes in
the data, each of which presents unique challenges and considerations for
coding.

### Coding Tasks

The models and two human coders annotate the dataset using the same human-
generated codebook per task222To facilitate the coding process and adhere to
the models’ token length limits, two researchers extracted relevant segments
from the interview transcripts for all three coding tasks.. We hypothesize
that the difficulty of assigning tags to interview segments is influenced by
several factors, such as the length of the segments, the codebook and its
length, and the background knowledge and judgment required to assign
appropriate categories. In this context, Braun and Clarke (2006) classify
possible themes as semantic/latent or a combination of both. Semantic themes
can be identified from the surface meanings of the data by simply looking at
what a participant has said. Latent themes go beyond the surface meaning of
the data and require interpretation of the underlying ideas and assumptions.
To analyze how the similarity of results between human coders and LLMs is
affected by coding for semantic versus latent themes, we perform coding for
three different coding tasks:

  1. (1)

Task A (Internet-connected devices): This coding task involves identifying the
Internet-connected devices that participants use. We expect that the
assignment of semantic codes for Internet-connected devices is straightforward
for both humans and LLMs, and requires only the identification of the entities
as communicated by the interviewees. (average segment length: 118 words;
codebook length: 18 codes)

  2. (2)

Task B (Apps, programs, services, and use cases): This coding task focuses on
apps, programs, and services participants use on their Internet-connected
devices and for what purpose. This introduces a layer of complexity, as
participants may articulate their interactions in different ways, e. g., by
enumerating individual apps, grouping applications, or explicitly describing
their use cases. This variability introduces a hierarchy in the data that
demands both semantic and latent coding. (average segment length: 274 words;
codebook length: 24 codes)

  3. (3)

Task C (Trusted sources): This coding task explores participants’ practices
and sources when seeking guidance on digital security and privacy. This task
goes beyond capturing the semantic content of the data, as it requires
examining underlying ideas and assumptions. (average segment length: 469
words; codebook length: 32 codes)

### Experimental Design

For our experiments, we use OpenAI’s GPT-3.5 (Brown et al., 2020) and GPT-4
(Achiam et al., 2023) models via the API service. To achieve more consistent
and less random completions, we set the temperature parameter to $0$. To
enable performance comparisons between humans and LLMs, we implement both
coding approaches in an annotation task, where the models and human coders
code all data segments using a human-generated codebook as a common benchmark
for evaluation. We design a prompt for each coding task to instruct the models
to perform the coding given an interview segment and the codebook, following
the scheme in Figure 1. We experiment with two different models and prompt
engineering techniques, including few-shot learning (OpenAI, 2023; White et
al., 2023; Brown et al., 2020):

  * •

GPT-3.5 Turbo vs. GPT-4: We compare the performance of OpenAI’s GPT-3.5 and
GPT-4 models. While GPT-4 offers broader general knowledge via a larger
context window and can follow complex instructions, GPT-3.5 is offered at a
lower price, and versions of this model are also accessible via OpenAI’s free
AI system ChatGPT.

  * •

Zero-shot vs. One-shot vs. Few-shot: Providing the model with expected inputs
and outputs can increase performance, and teach the model to adopt a specific
style (Brown et al., 2020). In the zero-shot setting, we do not include
examples in the prompt. In the one-shot setting, we provide the model with one
exemplary data segment and expected output. In the few-shot setting, we
provide three examples333As examples, we use synthetic interview data
generated with GPT-4 and annotated by human coders..

### Evaluation of Agreement

For each coding task, we calculate the Inter-Rater Reliability (IRR) between
the two human coders as well as between the final, unified human annotations
and both models via Cohen’s Kappa (Cohen, 1960; McHugh, 2012). While the IRR
between human coders is based on the 80% of annotations that were not used to
generate the codebook, we compute the IRR between models and humans based on
all annotations.

##  3\. Results: Agreement between Humans and LLMs

Table 1 presents the IRR between human coders’ and models’ annotations. For
tasks A and B, human coders achieved almost perfect agreement and for Task C
they reached substantial agreement. GPT-4 consistently outperforms its
predecessor in all three tasks. It achieves almost perfect agreement with
human annotations on Task A in all settings – comparable to inter-human
agreement. In contrast, with zero or one example, GPT-3.5 only achieves
substantial agreement in Task A, which can be raised to almost perfect
agreement when three examples are provided. In Task B, both models show
moderate to substantial agreement. However, in Task C, the agreement for
GPT-3.5 is only fair for the one-shot and few-shot settings, while GPT-4
achieves moderate agreement across all settings. As the complexity of the
coding tasks increases from Task A to Task C, there is both a notable decrease
in inter-human agreement and a widening gap between inter-human and model-
human scores.

Table 1. IRR in the form of Cohen’s $\kappa$ between the two human coders
(top) and between both models and human coders (bottom). Scores for Cohen’s
$\kappa$ range from $-1$ to $1$, and can be characterized as indicating no
agreement ($\kappa<0.0$), slight agreement ($0.0\leq\kappa\leq 0.2$), fair
agreement ($0.2<\kappa\leq 0.4$), moderate agreement ($0.4<\kappa\leq 0.6$),
substantial agreement ($0.6<\kappa\leq 0.8$), or almost perfect agreement
($0.8<\kappa\leq 1$) (Landis and Koch, 1977; McHugh, 2012).

|  Task A (Internet-connected devices) |  Task B (Apps, programs, services) |  Task C (Trusted sources)
---|---|---|---
inter-human | 0.97 | 0.83 | 0.63
| Zero-Shot | One-Shot | Few-Shot | Zero-Shot | One-Shot | Few-Shot | Zero-Shot | One-Shot | Few-Shot
gpt-3.5-turbo | 0.74 | 0.78 | 0.85 | 0.60 | 0.58 | 0.56 | 0.43 | 0.36 | 0.34
gpt-4 | 0.97 | 0.94 | 0.95 | 0.75 | 0.74 | 0.75 | 0.56 | 0.55 | 0.58

Providing the model with examples did not significantly improve IRR. However,
it played a critical role in mitigating certain challenges, particularly for
GPT-3.5, given its inherent limitations, such as restricted input and output
sizes. In particular, few-shot learning helped reduce the number of formatting
errors that required manual correction. It also played a key role in reducing
the number of codes introduced by the model that did not exist in the given
codebook, commonly referred to as hallucinations (Ji et al., 2023). Notably,
GPT-3.5 introduced as many as $47$ new, incorrect codes for Task B – an amount
that could be reduced by more than half with the provision of examples. In
contrast, GPT-4 generated only two incorrect codes in the same setting.

##  4\. Discussion

### Inter-human vs. model-human performance

Although inter-human agreement is in one case equivalent to, and generally
higher than model-human agreement, our research highlights the shared
challenges faced by human coders and LLMs when confronted with increasingly
complex QDA tasks. Our quantitative findings underscore the inherent
difficulty in coding with longer codebooks and contexts, which necessitates
capturing nuanced signals and latent themes. Thus, we advocate for the
evaluation of LLMs on a task-specific basis, recognizing that not all tasks
are uniformly compatible with LLMs. Various factors, including data collection
methods (e. g., survey data, interview data) and segmentation techniques,
contribute to the complexity of QDA. Our ongoing work delves into these
dimensions to gain a comprehensive understanding of the challenges and nuances
associated with using LLMs for qualitative coding.

### Model choice and few-shot learning

GPT-4 outperformed GPT 3.5 on all tasks, indicating an improved ability to
understand and encode qualitative data. While GPT-4 shows higher agreement
with human coders, GPT-3.5 offers a considerably lower cost, with input tokens
priced 60 times less in our experiments (OpenAI, 2024). This makes GPT-3.5
particularly suitable for rapid prototyping, e. g., when engineering prompts
for QDA. While few-shot learning did not improve the models’ performance, it
could play a crucial role in mitigating hallucinations and formatting errors,
especially for less capable models such as GPT-3.5. However, researchers will
need to consider whether to use synthetic examples or to manually annotate a
subset of data to provide examples. This decision may impact performance,
calling for a careful analysis of the choice of examples for few-shot
learning.

### Methodological challenges

While IRR provides a quantitative measure of agreement, it falls short of
capturing the depth and subtleties of qualitative analysis (McDonald et al.,
2019). In our future work, we aim to conduct a qualitative assessment of our
results, along with a detailed error analysis. It is critical to recognize
that the human benchmark is not infallible, as human coders are susceptible to
subjectivity and potential bias (Seaman, 1999; Ortloff et al., 2023). This
raises the normative question of whether to hold LLMs to the standard of human
coding or to explore alternative scales and benchmarks that might better
capture the strengths and limitations of these models.

### Risks and limitations

While LLMs exhibit remarkable performance in processing and generating
language (Brown et al., 2020; Ziems et al., 2023), there are inherent
limitations in their ability to understand and encode complex contexts,
especially in tasks that require a deeper understanding, potentially leading
to oversights and misinterpretations. Moreover, the role and integration of
LLMs in QDA requires a normative discussion. For some, the researcher’s
subjectivity is integral to QDA (Wa-Mbaleka, 2020), where they essentially
serve as the measurement instrument (Merriam and Tisdell, 2015). In this view,
coding involves becoming familiar with the data, and its meaning is derived
through human interpretation. Using LLMs for coding could disrupt this
process, removing human involvement in text interpretation. Thus, aside from
determining the suitability of LLMs for specific coding tasks, it is crucial
to reconsider the fundamental purpose of QDA. In addition, the current lack of
transparency in LLMs makes it difficult to track the reasoning behind the
results generated. Small changes in prompt phrasing can result in completely
different outputs (Jones and Steinhardt, 2022). Memory issues and limited
context windows raise concerns about the consistency and reproducibility of
model outputs, as the model does not retain information about its coding
decisions between segments. In addition, the use of LLMs involves the
processing of significant amounts of data, raising questions about the
confidentiality and protection of sensitive information of participants,
especially when using ChatGPT as a non-API consumer application. As LLMs
become more integrated into qualitative research, ethical considerations
regarding these challenges and their potential impact on research findings
become paramount. Continued attention to ethical guidelines and responsible
use of LLMs is essential to maintaining the integrity and reliability of
qualitative research.

## References

  * (1)
  * Achiam et al. (2023) Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al. 2023\.  GPT-4 Technical Report.  _arXiv preprint arXiv:2303.08774_ (2023).
  * Adams et al. (2008) Anne Adams, Peter Lunt, and Paul Cairns. 2008.  A qualitative approach to HCI research.  (2008).
  * Bergin (2018) Tiffany Bergin. 2018\.  An Introduction to Data Analysis: Quantitative, Qualitative and Mixed Methods.  _An Introduction to Data Analysis_ (2018), 1–296.
  * Braun and Clarke (2006) Virginia Braun and Victoria Clarke. 2006.  Using Thematic Analysis in Psychology.  _Qualitative research in psychology_ 3, 2 (2006), 77–101.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. 2020\.  Language models are few-shot learners.  _Advances in neural information processing systems_ 33 (2020), 1877–1901.
  * Cohen (1960) Jacob Cohen. 1960\.  A Coefficient of Agreement for Nominal Scales.  _Educational and psychological measurement_ 20, 1 (1960), 37–46.
  * Flick (2022) Uwe Flick. 2022\.  An Introduction to Qualitative Research.  _An introduction to qualitative research_ (2022), 1–100.
  * Gilardi et al. (2023) Fabrizio Gilardi, Meysam Alizadeh, and Maël Kubli. 2023\.  ChatGPT Outperforms Crowd-Workers for Text-Annotation Tasks.  _arXiv preprint arXiv:2303.15056_ (2023).
  * Ji et al. (2023) Ziwei Ji, Nayeon Lee, Rita Frieske, Tiezheng Yu, Dan Su, Yan Xu, Etsuko Ishii, Ye Jin Bang, Andrea Madotto, and Pascale Fung. 2023.  Survey of Hallucination in Natural Language Generation.  _Comput. Surveys_ 55, 12 (2023), 1–38.
  * Jones and Steinhardt (2022) Erik Jones and Jacob Steinhardt. 2022.  Capturing Failures of Large Language Models via Human Cognitive Biases.  _Advances in Neural Information Processing Systems_ 35 (2022), 11785–11799.
  * Kuckartz (2019) Udo Kuckartz. 2019\.  Qualitative Text Analysis: A Systematic Approach.  _Compendium for early career researchers in mathematics education_ (2019), 181–197.
  * Landis and Koch (1977) J Richard Landis and Gary G Koch. 1977.  The Measurement of Observer Agreement for Categorical Data.  _biometrics_ (1977), 159–174.
  * McDonald et al. (2019) Nora McDonald, Sarita Schoenebeck, and Andrea Forte. 2019\.  Reliability and Inter-Rater Reliability in Qualitative Research: Norms and Guidelines for CSCW and HCI Practice.  _Proc. ACM Hum.-Comput. Interact. (CSCW)_ (2019).
  * McHugh (2012) Mary L McHugh. 2012\.  Interrater reliability: the kappa statistic.  _Biochemia medica_ 22, 3 (2012), 276–282.
  * Merriam and Tisdell (2015) Sharan B Merriam and Elizabeth J Tisdell. 2015.  _Qualitative Research: A Guide to Design and Implementation_.  John Wiley & Sons.
  * OpenAI (2023) OpenAI. 2023.  _Prompt Engineering Guide_.  Retrieved 2024-02-02 from https://platform.openai.com/docs/guides/prompt-engineering
  * OpenAI (2024) OpenAI. 2024.  _Pricing_.  Retrieved 2024-02-23 from https://openai.com/pricing
  * Ortloff et al. (2023) Anna-Marie Ortloff, Matthias Fassl, Alexander Ponticello, Florin Martius, Anne Mertens, Katharina Krombholz, and Matthew Smith. 2023.  Different Researchers, Different Results? Analyzing the Influence of Researcher Experience and Data Type During Qualitative Analysis of an Interview and Survey Study on Security Advice. In _Proceedings of the 2023 CHI Conference on Human Factors in Computing Systems_. 1–21.
  * Seaman (1999) C.B. Seaman. 1999\.  Qualitative methods in empirical studies of software engineering.  _IEEE Transactions on Software Engineering_ 25, 4 (1999), 557–572.  https://doi.org/10.1109/32.799955
  * Team et al. (2023) Gemini Team, Rohan Anil, Sebastian Borgeaud, Yonghui Wu, Jean-Baptiste Alayrac, Jiahui Yu, Radu Soricut, Johan Schalkwyk, Andrew M Dai, Anja Hauth, et al. 2023\.  Gemini: A Family of Highly Capable Multimodal Models.  _arXiv preprint arXiv:2312.11805_ (2023).
  * Törnberg (2023) Petter Törnberg. 2023\.  ChatGPT-4 Outperforms Experts and Crowd Workers in Annotating Political Twitter Messages with Zero-Shot Learning.  _arXiv preprint arXiv:2304.06588_ (2023).
  * Touvron et al. (2023) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al. 2023\.  Llama 2: Open Foundation and Fine-Tuned Chat Models.  _arXiv preprint arXiv:2307.09288_ (2023).
  * Wa-Mbaleka (2020) Safary Wa-Mbaleka. 2020\.  The Researcher as an Instrument. In _Computer Supported Qualitative Research_ , António Pedro Costa, Luís Paulo Reis, and António Moreira (Eds.). Springer International Publishing, Cham, 33–41.
  * White et al. (2023) Jules White, Quchen Fu, Sam Hays, Michael Sandborn, Carlos Olea, Henry Gilbert, Ashraf Elnashar, Jesse Spencer-Smith, and Douglas C Schmidt. 2023.  A Prompt Pattern Catalog to Enhance Prompt Engineering with ChatGPT.  _arXiv preprint arXiv:2302.11382_ (2023).
  * Ziems et al. (2023) Caleb Ziems, William Held, Omar Shaikh, Jiaao Chen, Zhehao Zhang, and Diyi Yang. 2023\.  Can Large Language Models Transform Computational Social Science?  _arXiv preprint arXiv:2305.03514_ (2023).

##  Appendix A PROMPT DESIGN

Figure 1. Prompt template for coding Internet-connected devices with optional
examples. We modify the highlighted text sections to change the coding task.

  *[QDA]: Qualitative Data Analysis
  *[HCI]: Human-Computer Interaction
  *[LLMs]: Large Language Model
  *[LLM]: Large Language Model
  *[IRR]: Inter-Rater Reliability