Sparse Recovery

Model

The focus of this article is recovering a sparse vector $x \in \mathbb{R}^n$ in the streaming setting. Initially, $x = 0$. Then we see a sequence of updates $(i, \Delta)$ each causing the change $x_i \leftarrow x_i + \Delta$. There are three popular studied models

1. insertion only: $\Delta = 1$ for each update.
2. strict turnstile: $\Delta$ can be negative, but it is guaranteed that $x \ge 0$ all the time.
3. general turnstile: $\Delta$ can be negative, and entries of $x$ can be negative.

Let $\text{supp}(x) \subset [n]$ be the set of nonzero indices of $x$. The vector $x$ is called $k$ sparse if $\text{supp}(x) \le k$.

Goal. Given an integer $k$, design a data structure such that, if $x$ is $k$ sparse, the structure returns $x$ exactly; otherwise, the return can be arbitrary.

As we can't actually store real numbers, we use the following assumption:

The updates are integers and the entries of $x$ are bounded by $\pm M$.

With this restriction, it suffices to solve this problem in a finite field $\mathbb{F}_p$ where $p \ge \max \{ n, M \}$. As opposed to the naive solution that keeps all entries of $x$ and uses $O(n \log M)$ bits of memory, we show that the memory can be reduced to $O( k \log (nM) )$ bits.

Algorithm

The problem is easy to solve in the insertion only model. We can use a hash set that keeps the first $k$ occurrence of the indices. If the we detect the $k + 1$ new index, we can just report error.

The scenario is a little more complicated for the other two models. E.g., let $n = 4$ and $k = 1$. Consider the following stream

The vector $x$ is finally $1$-sparse and the only non-zero index is revealed at the end of the stream.

Solution 1

Encode

In the first approach, we use a matrix $H \in \mathbb{F}_p^{m \times n}$ ($m$ to be determined) and encode $x$ as $\mathbb{Enc}(x) = H x$. If the following property holds, we are guaranteed to recover $x$ if it is $k$ sparse.

For every $k$ sparse $x, x' \in \mathbb{F}_p^n$, if $x \neq x'$, then $Hx \neq H x'$.

If $H x \neq Hx'$, then $H (x - x') \neq 0$. Observe that $x - x'$ is a $2k$ sparse vector. This is equivalent to that $H z \neq 0$ for every nonzero $2k$ sparse vector $z \in \mathbb{F}_p^n$.

Such $H$ does exists, with $m = 2k$. Let $a_1, a_2, ..., a_n \in \mathbb{F}_p$ be different integers. Then $H$ is given as

Let $S = \text{supp}(z) \subset [n]$ be a set of size $2k$. Denote $z_S$ be the vector by keeping only the nonzero dimensions of $z$ and $H_S$ be an $2k \times 2k$ matrix by keeping only columns indexed by $S$. Then $H z = H_S z_S$. We show that the rows of $H_S$ are linearly independently, hence its columns. It follows that if $z_S \neq 0 \Longrightarrow H_S z_S \neq 0$.

It is left to show that $\forall c \in \mathbb{F}_p^{2k}$, $c \neq 0 \rightarrow H_S^T c \neq 0$. Define the polynomial in $\mathbb{F}_p$

Then $H_S^T c = ( f(a_i ) )_{i \in [S] }$. As $f(x)$ has at most $2k - 1$ distinct roots, and the $(a_i)_{i \in [S] }$ are different, there exists some $i \in [S]$ such that $f(a_i) \neq 0$.

Decode

Given $y = \mathbb{Enc}(x) = H x$, how can we recover $x$? We can loop over every subset $S \subset [n ]$ of size $2k$ and compute $x' = H_S^{-1} y$. If $x'$ is $k$ sparse, by the property of $H$, it is the unique $k$-sparse vector such that $H_S x' = y$.

Remark: The proposed method is exponential respect to $k$. However, there is an algorithm (called syndrome decoding) to perform the recovery in $O(k^2 poly(\log p))$ time.

Space and Time

Space. The algorithm uses $O(k \log (nM))$ bits.

This is because we don't store $H$. For the finite field $\mathbb{F}_p$, there exists at least one primitive element whose powers give all nonzero elements. We can pick such a primitive element $a$, and set $a_1 = a$, $a_2 = a^2$, ..., $a_n = a^n$. Storing $y = Hx$ requires only $O(k \log (nM))$ bits.

Time. Processing each update $(i, \Delta)$ in the stream takes $O(k)$ time.

Solution 2

We first show how to recover an $1$ sparse vector when negative updates are allowed, then use it as a black box to recover $k$ sparse vector.

1 sparse recovery

We show an algorithm that successes with probability $n / p$.

First, we pick a non-zero integer $z$ uniformly at random from $\mathbb{F}_p$. Let

The algorithm consists of 3 parts.

Initialize:
$\qquad$ $w \leftarrow \{ 0 \}^3$
$\qquad$ Sample $z$ uniformly from $1, ..., p - 1$.

Updates: $(i, \Delta)$
$\qquad$ $w \leftarrow w + \Delta \cdot H_i$, where $H_i$ is the $i$th column of $H$.

Output:
$\qquad$ If $w_1 = 0$ or $w_3 \neq w_1 \cdot z^{w_2 / w_1}$:
$\qquad\qquad$ Return $x$ is not $1$ sparse;
$\qquad$ Else:
$\qquad\qquad$ Return $j \doteq w_2 / w_1$

Failure Analysis

Theorem. The algorithm fails with probability $n / p$.

Observe that

If $w_1 = 0$, either $x = 0$ or there is cancellation. Conditioned on $w_1 \neq 0$, if $x$ is $1$ sparse, then

Hence, $j = w_2 / w_1$ is exactly the non-zero index, and $w_3 = w_1 \cdot z^{j}$.

If $x$ is not $1$ sparse, then

is a polynomial of degree at most $n$. It has at most $n$ distinct roots with one of them being $0$. The probability that $z$ being a root is at most $(n - 1) / (p - 1) \le n / p$.

Space and Time

Space. The algorithm uses $O(\log (nM))$ bits.

Time. Processing each update $(i, \Delta)$ takes $O(1)$ time.

$k$ sparse recovery

To apply the $1$ sparse recovery data structure, we partition $x$ into $2k$ subsets randomly. In particular, we use a pairwise hash function that maps each $i \in [n]$ to $[2k]$ uniformly at random. For $\ell \in [2k]$, let $S_\ell$ be the set of indices mapped to $\ell$:

Then we run $1$ sparse recovery on each partition. If some partition contains only one nonzero index, then we recover it successfully.

Theorem. If $x$ is $k$ sparse, an nonzero dimension $i$ is recovered with probability $1 / 2$.

Proof. It suffices to prove that $i$ is the only nonzero index in $S_{h (i) }$. For $j \in \text{supp}(x)$, define the indicator variable $X_j$ which equals $1$ if $h(j) = h(i)$ and equals 0 otherwise. By Markov's inequality and linearity of expectation,

Remark. The Markov's inequality is rather tight here. Even if we have a fully independent hash function $h$, the probability that none of the nonzero indices collides with $i$ is given by

$\square$

Corollary. Running $\log (k / \delta)$ independent copies of the above structure reduce the failure probability to $\delta$.

Space and Time

Space. The algorithm uses $O(k (\log (k / \delta)) \log (nM))$ bits.

Time. Processing each update $(i, \Delta)$ takes $O(\log (k / \delta))$ time.

Reference

[1]. Jelani Nelson, Chapter 4, Sketching Algorithms, 2020.

[2]. Moses Charikar, $\ell_{0}$ Sampling and Sketching Distances, CS369G: Algorithmic Techniques for Big Data, 2016