Improved Algorithms For Structured Sparse Recovery

TL;DR

This paper introduces improved, efficient algorithms for structured sparse recovery in compressive sensing, specifically for tree sparsity and CEMD models, enhancing approximation accuracy and reducing computational complexity.

Contribution

It presents the first linear-time approximation algorithms for head- and tail-approximations in structured sparsity models, solving open problems and improving recovery performance.

Findings

Linear time $(1- ext{epsilon})$-approximation for tree sparsity head-approximation

Linear time $(1+ ext{epsilon})$-approximation for tree sparsity tail-approximation

First constant factor approximation for CEMD model head-approximation

Abstract

It is known that certain structures of the signal in addition to the standard notion of sparsity (called structured sparsity) can improve the sample complexity in several compressive sensing applications. Recently, Hegde et al. proposed a framework, called approximation-tolerant model-based compressive sensing, for recovering signals with structured sparsity. Their framework requires two oracles, the head- and the tail-approximation projection oracles. The two oracles should return approximate solutions in the model which is closest to the query signal. In this paper, we consider two structured sparsity models and obtain improved projection algorithms. The first one is the tree sparsity model, which captures the support structure in the wavelet decomposition of piecewise-smooth signals. We propose a linear time $(1-\epsilon)$-approximation algorithm for head-approximation projection and a linear time $(1+\epsilon)$-approximation algorithm for tail-approximation projection. The best previous result is an $\tilde{O}(n\log n)$ time bicriterion approximation algorithm (meaning that their algorithm may return a solution of sparsity larger than $k$) by Hegde et al. Our result provides an affirmative answer to the open problem mentioned in the survey of Hegde and Indyk. As a corollary, we can recover a constant approximate $k$-sparse signal. The other is the Constrained Earth Mover Distance (CEMD) model, which is useful to model the situation where the positions of the nonzero coefficients of a signal do not change significantly as a function of spatial (or temporal) locations. We obtain the first single criterion constant factor approximation algorithm for the head-approximation projection. The previous best known algorithm is a bicriterion approximation. Using this result, we can get a faster constant approximation algorithm with fewer measurements for the recovery problem in CEMD model.

Equations

$$\mathcal{M}=\{x\in {\mathbb{R}}^d\mid \mathbf{supp}(x)\subseteq \Omega \text{ for some }\Omega \in \mathbb{M}\}.$$

$$\Vert x-x_{{\Omega }^{\ast }}{\Vert }_p=\underset{\Omega \in \mathbb{M}}{\min }\Vert x-x_{\Omega }{\Vert }_p$$

$$\Vert x_{\Omega }{\Vert }_p\ge c_H\cdot \Vert x_{{\Omega }^{\ast }}{\Vert }_p$$

$$\Vert x-x_{\Omega }{\Vert }_p\le c_T\cdot \Vert x-x_{{\Omega }^{\ast }}{\Vert }_p,$$

$${\mathcal{T}}_k=\{x\in {\mathbb{R}}^d\mid \mathbf{supp}(x)\subseteq \Omega \text{ for some }\Omega \in {\mathbb{T}}_k\}$$

$$\mathbf{EMD}(A,B)=\underset{\pi :A\to B}{\min }\sum _{a\in A}|a-\pi (a)|,$$

$$\mathbf{EMD}[\Omega ]=\sum _{i=1}^{w-1}\mathbf{EMD}({\Omega }_i,{\Omega }_{i+1}).$$

$${\mathcal{M}}_{k,B}=\{x\in {\mathbb{R}}^d\mid \mathbf{supp}(x)\subseteq \Omega \text{ for some }\Omega \in {\mathbb{M}}_{k,B}\}$$

$$\Vert x_{\hat{\Omega }}{\Vert }_p\ge (1-{ϵ}_1)\underset{\Omega \in {\mathbb{T}}_k}{\max }\Vert x_{\Omega }{\Vert }_p.$$

$$\Vert x-x_{\hat{\Omega }}{\Vert }_p\le (1+{ϵ}_2)\underset{\Omega \in {\mathbb{T}}_k}{\min }\Vert x-x_{\Omega }{\Vert }_p,$$

$$\hat{a}[i_{w-1}]\le \hat{a}[i_w]/(1+\alpha )\le \hat{a}[i_w]/(1+\beta )\le (1+\beta {)}^{\lceil {\log }_{1+\beta }\theta \rceil -1}\le (1+\beta {)}^{\lfloor {\log }_{1+\beta }\theta \rfloor }\le \theta$$

$$s[l_0]\le s^{\prime }[l_0]=\hat{s}[l^{\prime }]=s[l^{\prime }]=(1+\beta )\theta <(1+\beta )s[l+1]\le (1+\beta )s[l_0].$$

$$\sum _{N_{ij}\notin \hat{\Omega }}{x}_{ij}\le \sum _{N_{ij}\notin \hat{\Omega }}\le \sum _{N_{ij}\notin {\Omega }^{\ast }}{\hat{x}}_{ij}\le (1+ϵ)\sum _{N_{ij}\notin {\Omega }^{\ast }}{x}_{ij}=s[k].$$

$$\sum _{i\in (\xi ,\eta ]}\frac{n+1}{2^i}\cdot \frac{2^i}{{ϵ}^{\prime }}=O({ϵ}^{-1}n((\log \log \log n{)}^2+{\log }^2(1/ϵ))).$$

$$\sum _{i\in (\eta ,\log n]}\frac{n+1}{2^i}\cdot \frac{\log n}{{ϵ}_{i-1}{ϵ}_i}=O({ϵ}^{-1}n).$$

$$s[l]\le a[l_1]+b[l_2]\le a^{\prime }[l_1]+b^{\prime }[l_2]=\tilde{s}[l]\le s^{\prime }[l].$$

$$s^{\prime }[l]\le (1+{ϵ}^{\prime })\tilde{s}[l]\le (1+{ϵ}^{\prime })(a^{\prime }[l_1]+b^{\prime }[l_2])\le (1+{ϵ}^{\prime })(1+{ϵ}^{\ast })(a[l_1]+b[l_2])=s[l].$$

$$s^{\prime }[2^i-2]=a^{\prime }[2^{i-1}-1]+b^{\prime }[2^{i-1}-1]\le (1+{ϵ}^{\ast })(a[2^{i-1}]+b[2^{i-1}-1])=(1+{ϵ}^{\ast })s[2^i-2].$$

$$(1+{ϵ}^{\prime })\hat{s}[2^i-2]=(1+{ϵ}^{\prime })s^{\prime }[2^i-2]\le (1+{ϵ}^{\prime })(1+{ϵ}^{\ast })s[2^i-2]\le (1+{ϵ}^{\prime })(1+{ϵ}^{\ast })s[2^i-1]\le \hat{s}[2^i-1].$$

$$O\left(\frac{n^{1-\delta }}{ϵ}\right)\sum _{i=1}^{\eta }O\left(\frac{{\log }^{2}n}{2^i}\cdot 2^{2i-c\sqrt{i}}\right)=O({ϵ}^{-1}{n}^{1-\delta }{\log }^{4}n)=o(n).$$

$$\sum _{i=\eta +1}^{\log n}O\left(\frac{n+1}{2^i}\cdot \frac{\log n}{{ϵ}_i^2}\right)=O({ϵ}^{-2}n/\log n)=o(n).$$

$$\Phi [{\Omega }_{{\check{P}}_t}]+\Phi [{\Omega }_{{\hat{P}}_t}]\ge \Phi [{\Omega }_P],\Delta [{\Omega }_{{\check{P}}_t}]+\Delta [{\Omega }_{{\hat{P}}_t}]\le \Delta [{\Omega }_P].$$

$$\Delta [{\Omega }_{{\check{P}}_0}]=0\le \Delta [{\Omega }_P]/d<\Delta [{\Omega }_P]=\Delta [{\Omega }_{{\check{P}}_h}].$$

$$\Phi [{\Omega }_{{\hat{P}}_{t+1}}]\ge \Phi [{\Omega }_P]-\Phi [{\Omega }_{{\check{P}}_{t+1}}]\ge (1-1/d)\Phi [{\Omega }_P],$$

$$\Delta [{\Omega }_{{\hat{P}}_{t+1}}]\le \Delta [{\Omega }_P]-\Delta [{\Omega }_{{\check{P}}_{t+1}}]\le (1-1/d)\Delta [{\Omega }_P].$$

$$\Phi [{\Omega }_{P^{\prime }}]\ge \Phi [{\Omega }_{{\hat{P}}_{t+1}}]/2(d-1)>\Phi [{\Omega }_P]/2d,\Delta [{\Omega }_{P^{\prime }}]\le \Delta [{\Omega }_{{\hat{P}}_{t+1}}]/(d-1)<\Delta [{\Omega }_P]/d.$$

$$\max \{\Phi [{\Omega }_r^{\prime }],\Phi [{\Omega }_l]\}\ge (\frac{1}{4}-\delta )\mathsf{OPT}.$$

$$\Phi [{\Omega }_{\lambda }]-\lambda \cdot \Delta [{\Omega }_{\lambda }]={\max }_{\Omega \in {\mathbb{M}}_{k,B}}\{\Phi [\Omega ]-\lambda \cdot \Delta [\Omega ]\}$$

$$\Phi [{\Omega }_l]-l\cdot \Delta [{\Omega }_l]$$

$$\Phi [{\Omega }_l]$$

$$=\frac{1}{4}\mathsf{OPT}-\frac{x_{\min }\delta B}{w{h}^2}\ge \frac{1}{4}\mathsf{OPT}-\delta x_{\min }\ge (\frac{1}{4}-\delta )\mathsf{OPT}.$$

Full text

1 Introduction

We consider the robust sparse recovery, an important problem in compressive sensing. The goal of robust sparse recovery is to recover a signal from a small number of linear measurements. Specifically, we call a vector $x\in\mathbb{R}^{n}$ $k$-sparse if it has at most $k$ non-zero entries. The support of $x$, denoted by $\mathbf{supp}(x)\subseteq[n]$, contains the indices corresponding to the nonzero entries in $x$. Given the measurement vector $y=Ax+e$, where $A$ is a measurement matrix, $x$ is a $k$-sparse signal and $e$ is a noise vector, the goal is to find a signal estimate $\hat{x}$ such that $\|x-\hat{x}\|<C\|e\|$ for some constant approximation factor $C>0$.

Structured Sparsity Model: It is well known that in general we need the number of measurements to be $\Omega(k\log(n/k))$ for robust sparse recovery, see [9, 12]. In practice, the support of $x$ usually has some structured constraints, such as tree sparsity and block sparsity, which can help reduce the bound of the number of measurements.

Definition 1** (Structured Sparsity Model [2]).**

Let $\mathbb{M}$ be a family of supports, i.e., $\mathbb{M}=\{\Omega_{1},\Omega_{2},$ $\ldots,\Omega_{L}\}$ where each $\Omega_{i}\subseteq[n]$. Then the corresponding structured sparsity model $\mathcal{M}$ is the set of vectors supported on one of the $\Omega_{i}$:

[TABLE]

To recover such a structured signal $x$ is called structured sparse recovery. Baraniuk et al. [2] provided a general framework called model-based compressive sensing. Their framework depends on a model projection oracle which is defined as follows.

Definition 2** (Model Projection [2]).**

Let $\mathbb{M}$ be a structured sparsity model. A model projection oracle for $\mathbb{M}$ is an algorithm $P(x):\mathbb{R}^{n}\rightarrow\mathbb{M}$ such that $\Omega^{*}=P(x)$ and

[TABLE]

where $x_{\Omega}\in\mathbb{R}^{n}$ is the same as $x$ on the support $\Omega$, and is zero otherwise.

Unfortunately, the best known algorithms for many exact model projection oracles are too slow to be used in practice. Some of the exact model projection oracles are even NP-hard. Recently, Hegde et al. [17] provided a principled method AM-IHT for recovering structured sparse signals. Their framework only requires two approximation oracles called the head- and tail-approximation projection oracles, defined as follows. Let $\Omega^{*}\in\mathbb{M}$ be the optimal support of the model projection oracle as defined in Definition 2.

Definition 3** (Head-Approximation Projection).**

Let $\mathbb{M}$ be a structured sparsity model. A head-approximation oracle for $\mathbb{M}$ is an algorithm $H(x):\mathbb{R}^{n}\rightarrow\mathbb{M}$ such that $H(x)=\Omega$ and

[TABLE]

where $c_{H}\in(0,1]$ is a fixed constant.

Definition 4** (Tail-Approximation Projection).**

Let $\mathbb{M}$ be a structured sparsity model. A tail-approximation oracle for $\mathbb{M}$ is an algorithm $T(x):\mathbb{R}^{n}\rightarrow\mathbb{M}$ such that $T(x)=\Omega$ and

[TABLE]

where $c_{T}\in[1,\infty)$ is a fixed constant.

1.1 Tree Sparsity Model and CEMD Model

Tree Sparsity Model: The tree sparsity model can be used for capturing the support structure of the wavelet decomposition of piecewise-smooth signals and images [4, 8, 15]. In this model, the coefficients of the signal $x$ are arranged as the nodes of a complete $b$-ary tree $T$ rooted at node $N$, and any feasible solution is a subtree which includes the root of $T$ and is of size $k$.

Definition 5** (Tree Sparsity Model).**

Let $T$ be a complete $b$-ary tree with $n$ nodes rooted at node $N$. $\mathbb{T}_{k}(T)=\{\Omega_{1},\Omega_{2},\ldots,\Omega_{L}\}$ is the family of supports where each $\Omega_{i}$ is a subtree of $T$ rooted at $N$ with the number of nodes no more than $k$. We use $\mathbb{T}_{k}$ instead of $\mathbb{T}_{k}(T)$ for short. The tree-structured sparsity model $\mathcal{T}_{k}$ is the set of signals supported on some $\Omega\in\mathbb{T}_{k}$:

[TABLE]

For the tree sparsity model, the head- and tail-approximation projection problems reduces to the following simple-to-state combinatorial problems: for the head-approximation, we want to find a subtree of size $k$ rooted at $N_{1}$ such that the total weight of the subtree is maximized; for the tail-approximation, we want the total weight of the complement of the subtree is minimized. For convenience, we abbreviate them as Tree-Sparsity-Head and Tree-Sparsity-Tail respectively. If our solution is a subtree of size at most $k$, we call it a single-criterion solution. Otherwise if our solution is a subtree of size larger than $k$, we call it a bicriterion solution.

Constrained EMD Model: The CEMD model, introduced by Schmidt et al [24], is particularly useful in 2D image compression and denoising [25, 11]. We first introduce the definition of Earth Mover’s Distance (EMD), also known as the Wasserstein metric or Mallows distance [22].

Definition 6** (EMD).**

The EMD of two finite sets $A,B\subset\mathbb{N}$ with $|A|=|B|$ is defined as

[TABLE]

where $\pi$ ranges over all one-to-one mappings from $A$ to $B$.

In the CEMD model, the signal $x\in\mathbb{R}^{n}$ can be interpreted as a matrix $X\in\mathbb{R}^{h\times w}$ with $n=hw$. By this interpretation, the support of $x\in\mathbb{R}^{h\times w}$, denoted by $\mathbf{supp}(x)\subseteq[h]\times[w]$, contains the indices $(i,j)$ ($i\in[h],j\in[w]$) corresponding to the nonzero entries in $x$.

Definition 7** (Support-EMD).**

Consider an $h\times w$ matrix $X$. Let $\Omega\subseteq[h]\times[w]$ be the support of a matrix $X$. Denote $\Omega_{i}$ to be the support of the column $i$ of $X$. Suppose $|\Omega_{i}|=s$ for $i\in[w]$. Then the EMD of the support $\Omega$ (or the support-EMD of $X$) is defined as

[TABLE]

Naturally, we have the following structured sparsity model which contains two constraints: 1) each column is $s(=k/w)$-sparse, 2) the support-EMD is at most $B$.

Definition 8** (Constrained EMD Model [24]).**

Let $\mathbb{M}_{k,B}$ be the family of supports $\{\Omega\subseteq[h]\times[w]\mid\mathbf{EMD}(\Omega)\leq B,\text{ and }|\Omega_{i}|=k/w,\text{ for }i\in[w]\}$. The Constrained EMD (CEMD) model $\mathcal{M}_{k,B}$ is the set of signals supported on some $\Omega\in\mathbb{M}_{k,B}$:

[TABLE]

Consider the CEMD model projection problem. If our solution belongs to $\mathcal{M}_{k,B}$, we call it a single-criterion solution. Otherwise if our solution does not belong to $\mathcal{M}_{k,B}$, i.e., there exists a column of sparsity larger than $s(=k/w)$ or the support-EMD is larger than $B$, we call it a bicriterion solution.

1.2 Our Contributions and Techniques

For both the tree sparsity and CEMD models, we consider the corresponding model-projection problems. We obtain improved approximation algorithms, which have faster running time, and return single-criterion solutions (rather than bicriterion solutions). Consequently, combining with the AM-IHT framework [17], our results implies better structured sparse recovery algorithms, in terms of the number of measurements, the sparsity of the solution, and the running time. We summarize our contributions and main techniques in the following.

Tree Sparsity Model:

Cartis et al. [6] gave an exact tree-sparsity projection algorithm with running time $O(nk)$. For the approximation version, Hegde et al. [15, 16] proposed bicriterion approximation schemes for both head- and tail-approximation tree-sparsity projection problems with running time $\tilde{O}(n\log n)$. Both algorithms achieve constant approximation ratio and output a tree of size at most $2k$. In this paper, we provide the first linear time algorithms for both head- and tail-approximation tree-sparsity projection problems and remove the bicriterion relaxation. This provides an affirmative answer to the open problem in Hegde and Indyk [18], which asks whether there is a nearly-linear time single-criterion approximation algorithm for tree sparsity.

Main Techniques for Tree Sparsity Model: The bottleneck of previous algorithms is computing exact $(\min,+)$-convolutions. Our main technique is to improve the running time of $(\min,+)$-convolutions. In Section 2, we introduce an approach of computing an approximate $(\min,+)$-convolution, called $(\alpha,\beta)$-RS $(\min,+)$-convolution. Instead of maintaining the whole $(\min,+)$-convolution array, we only compute a sparse sequence to approximately represent the whole array. Taking Tree-Sparsity-Tail as an example, we only need to maintain $\tilde{O}(\log n)$ elements in a single node, instead of $k$ elements for the exact $(\min,+)$-convolution. For the computation time, we show that the running time of computing each convolution element can be reduced to $\tilde{O}(1)$, instead of $O(k)$ for the exact $(\min,+)$-convolution. Thus, we only cost $\tilde{O}(\log n)$ to compute our approximate $(\min,+)$-convolution. For Tree-Sparsity-Head, we apply a similar approximate $(\max,+)$-convolution technique, called $(\alpha,\beta)$-RS $(\max,+)$-convolution. Our approximate convolution technique may have independent interest.

In Section 3, we combine the approximate $(\min,+)$-convolution technique and other approaches such as weight discretization, pruning and the lookup table method. Our results can be summarized by the following theorem.

Theorem 9** (Linear time head- and tail-approximation tree-sparsity projection).**

There are linear time algorithms for both head- and tail-approximation tree-sparsity projection problems. Specifically, for any constant $\epsilon_{1}\in(0,1)$, there is an $O(\epsilon_{1}^{-1}n)$ time approximation algorithm that returns a support $\hat{\Omega}\in\mathbb{T}_{k}$ satisfying

[TABLE]

For any constant $\epsilon_{2}\in(0,\infty)$, there is an $O(n+\epsilon_{2}^{-2}n/\log n)$ time approximation algorithm that returns a support $\hat{\Omega}\in\mathbb{T}_{k}$ satisfying

[TABLE]

if $k\leq n^{1-\delta}$ ($\delta\in(0,1)$ is any fixed constant), and there is an $O(\epsilon_{2}^{-1}n(\log\log\log n)^{2})$ time algorithm for general $k$.

Then combining with prior results [2, 16, 17], we provide a more efficient robust sparse recovery algorithm in tree sparsity model as follows. The best prior result can recover an approximate signal $\hat{x}\in\mathcal{T}_{ck}$ for some constant $c>1$ [16] (i.e., the sparsity of their solution is $ck$). In this paper, we improve the constant $c$ to $1$.

Corollary 10**.**

Assume that $k\leq n^{1-\delta}$ ($\delta\in(0,1)$ is any fixed constant). Let $A\in\mathbb{R}^{m\times n}$ be a measurement matrix. Let $x\in\mathcal{T}_{k}$ be an arbitrary signal in the tree sparsity model with dimension $n$, and let $y=Ax+e\in\mathbb{R}^{m}$ be a noisy measurement vector. Here $e\in\mathbb{R}^{m}$ is a noise vector. Then there exists an algorithm to recover a signal approximation $\hat{x}\in\mathcal{T}_{k}$ satisfying $\|x-\hat{x}\|\leq C\|e\|_{2}$ for some constant $C$ from $m=O(k)$ measurements. Moreover, the algorithm runs in $O((n\log n+k^{2}\log n\log^{2}(k\log n))\log\frac{\|x\|_{2}}{\|e\|_{2}})$ time.

CEMD Model: In Section 4, we consider the CEMD model $\mathcal{M}_{k,B}$ and propose the first single-criterion constant factor approximation algorithm for the head-approximation oracle.

Theorem 11**.**

Consider the CEMD model $\mathcal{M}_{k,B}$ with $s=k/w$ sparse for each column and support-EMD $B$. Let $\delta\in(0,1/4)$, $x_{\min}=\min_{|X_{i,j}|>0}|X_{i,j}|^{p}$, and $x_{\max}=\max|X_{i,j}|^{p}$. Let $c=1/4-\delta$. There exists an algorithm running in $O(shn\log\frac{n}{\delta}+\log\frac{x_{\max}}{x_{\min}})$ time, which returns a single-criterion $c^{1/p}$ approximation for the head-approximation projection problem.

Combining with AM-IHT framework [17], we obtain the following corollary which improves the prior result [17] in two aspects:

1. We decrease the total number of measurements from $m=O(k\log(\frac{B}{k}\log\frac{k}{w}))$ to $m=O(k\log(B/k))$. 2) We decrease the running time of the robust sparse recovery from $O(n\log\frac{\|x\|_{2}}{\|e\|_{2}}(k\log n+\frac{kh}{w}(B+\log n+\log\frac{x_{\max}}{x_{\min}})))$ to $O(n\log\frac{\|x\|_{2}}{\|e\|_{2}}(k\log n+\frac{kh}{w}(\log n+\log\frac{x_{\max}}{x_{\min}})))$.

Corollary 12**.**

Let $A\in\mathbb{R}^{m\times n}$ be a measurement matrix. Let $x\in\mathcal{M}_{k,B}$ be an arbitrary signal in the CEMD model with dimension $n=wh$, and let $y=Ax+e\in\mathbb{R}^{m}$ be a noisy measurement vector. Here $e\in\mathbb{R}^{m}$ is a noise vector. Then there exists an algorithm to recover a signal approximation $\hat{x}\in\mathcal{M}_{k,2B}$ satisfying $\|x-\hat{x}\|\leq C\|e\|_{2}$ for some constant $C$ from $m=O(k\log(B/k))$ measurements. Moreover, the algorithm runs in $O(n\log\frac{\|x\|_{2}}{\|e\|_{2}}(k\log n+\frac{kh}{w}(\log n+\log\frac{x_{\max}}{x_{\min}})))$ time, where $x_{\max}=\max|x_{i}|$ and $x_{\min}=\min_{|x_{i}|>0}|x_{i}|$.

1.3 Related Work

In the tree sparsity model, there is an algorithm with running time $O(nk\log n)$ for the exact model projection by dynamic programming. By using a more careful analysis, Cartis et al. [6] improved exact model projection to $O(nk)$. Actually, their dynamic program was based on computing $(\min,+)$-convolutions. The naive algorithm for computing the $(\min,+)$-convolution of arbitrary two length-$n$ arrays requires $O(n^{2})$ time. Williams proposed an improvement algorithm for computing $(\min,+)$-convolutions, and reduced the running time to $O(n^{2}/2^{\Omega(\sqrt{\log n})})$ [26].

For the approximation projection problem, several heuristic algorithms had been proposed, such as CSSA [3], CPRSS [10], optimal-pruning[4]. Hegde et al. [15, 16] improved the running time of tree sparse recovery to $\tilde{O}(n\log n)$.

Schmidt et al. [24] introduced the Constrained Earth Mover’s Distance (CEMD) model. Hegde, Indyk and Schmidt [17] proposed bicriterion approximation algorithms for both head- and tail-approximation projections. How to find a single-criterion approximation algorithm is an open problem mentioned in the survey [18].

Other structured sparsity models also have been studied by researchers. Huang et al. [20] first considered the graph sparsity model, and provided a head-approximation algorithm with a time complexity of $O(n^{c})$, where $c>1$ is a trade-off constant between time and sample complexity. For the tail-approximation projection problem, Hedge et al. [19] proposed a nearly-linear time bicriterion algorithm with the tail-approximation guarantee by modifying the GW scheme [13]. Hedge et al. [14] also studied the $\triangle$-separated model and provided an exact model projection algorithm.

Very recently during SODA17 conference, we knew that, in parallel to our work, Backurs et al. [1] also provided single criteria algorithms for the tree sparsity problem. Their algorithms can handle for more general trees and run in time $n(\log n)^{O(1)}$. Our algorithms only work for b-ary trees, but our running times are much better.

2 Approximate $(\min,+)$-Convolution

In this section, we introduce an approach of computing an approximate $(\min,+)$-convolution, which is useful for the tree sparsity model.

2.1 $(\alpha,\beta)$-RS $(\min,+)$-Convolution

We first introduce a concept called ($\min,+$)-convolution.

Definition 13** ($(\min,+)$-convolution (see e.g. [5, 7])).**

Given two arrays $A=(a[0],a[1],a[2],$ \ldots,$$a[m_{1}]) and $B=(b[0],b[1],b[2],\ldots,b[m_{2}])$, their $(\min,+)$-convolution is the array $S=(s[0],s[1],s[2],\ldots,$ $s[m_{1}+m_{2}])$ where $s[t]=\min_{i=0}^{t}\{a[i]+b[t-i]\},t\in[0,m_{1}+m_{2}]$.

We sketch how to use $(\min,+)$-convolutions in the tree sparsity model. Recall that $T_{ij}$ is the subtree rooted at $N_{ij}$. We maintain an array $S_{ij}=(s[0],s[1],\ldots,s[|T_{ij}|])$ for each node $N_{ij}$. The element $s[l]$ represents the optimal tail value for Tree-Sparsity-Tail on $T_{ij}$, i.e., $s[l]=\min_{\Omega\in\mathbb{T}_{|T_{ij}|-l}(T_{ij})}$ $\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}\setminus\Omega}x_{i^{\prime}j^{\prime}}$ where $\mathbb{T}_{k}(T_{ij})$ is the tree sparsity model defined at the tree $T_{ij}$ (see Definition 5). In fact, the array $S_{ij}$ can be achieved through computing the $(\min,+)$-convolution from the arrays of its two children. 111Note that the value of $s[|T_{ij}|]$ can not be directly obtained by the $(\min,+)$-convolution of the arrays of its two children. In fact, $s[|T_{ij}|]=\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$. Finally, we output the element $s[n-k]$ in the array of the root node $N_{\log(n+1),1}$, which is the optimal tail value of Tree-Sparsity-Tail on $T$.

However, the running time for computing the exact $(\min,+)$-convolution is too long. Instead, we compute an approximate $(\min,+)$-convolution for each node. We first introduce some concepts.

Definition 14** ($\alpha$-RS).**

Given a sequence $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}],\ldots,\hat{a}[i_{m}])$ and a fixed constant $\alpha\in[0,\infty)$. Each element $\hat{a}[i_{v}]\in\hat{A}$ is a real number with an associated index $i_{v}$. If for any $v\in[1,m-1]$, $i_{v+1}>i_{v},\hat{a}[i_{v+1}]\geq(1+\alpha)\hat{a}[i_{v}]\geq 0$, we call the sequence $\hat{A}$ an $\alpha$-representative sequence ($\alpha$-RS).

Definition 15** (Completion of $\alpha$-RS).**

Consider an $\alpha$-RS $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}]\ldots,\hat{a}[i_{m}])$. Define its completion by an array $A^{\prime}=(a^{\prime}[0],a^{\prime}[1],\ldots,a^{\prime}[i_{m}])$ satisfying that: 1) If $0\leq t\leq i_{1},a^{\prime}[t]=a^{\prime}[i_{1}]$; 2) If $i_{v}+1\leq t\leq i_{v+1}\ (1\leq v\leq m-1),a^{\prime}[t]=\hat{a}[i_{v+1}]$.

By the following definition, we show how to use an $\alpha$-RS to approximately represent an array.

Definition 16** (Sequence Approximation).**

Given two non-decreasing arrays $A^{\prime}=(a^{\prime}[0],a^{\prime}[1],\ldots,$ $a^{\prime}[n])$ and $A=(a[0],a[1],\ldots,a[n])$, we say $A^{\prime}$ is an $\alpha$-approximation of $A$ if for any $i$, $a[i]\leq a^{\prime}[i]\leq(1+\alpha)a[i]$. We say an $\alpha$-RS $\hat{A}$ approximates an array $A$ if its completion $A^{\prime}$ is an $\alpha$-approximation of $A$.

A special case is that we say $A^{\prime}$ is a [math]-approximation of $A$ if $A^{\prime}=A$. Figure 1 illustrates these concepts. Now, we are ready to introduce the formal definition of the approximate $(\min,+)$-convolution, called $(\alpha,\beta)$-RS $(\min,+)$-convolution.

Definition 17** ($(\alpha,\beta)$-RS $(\min,+)$-convolution).**

Given two $\alpha$-RSs $\hat{A}$ and $\hat{B}$, suppose $A^{\prime}$ and $B^{\prime}$ are their completions respectively. Suppose the array $S$ is the $(\min,+)$-convolution of $A^{\prime}$ and $B^{\prime}$. We call a sequence $\hat{S}$ an $(\alpha,\beta)$-RS $(\min,+)$-convolution of $\hat{A}$ and $\hat{B}$ if $\hat{S}$ is a $\beta$-RS which approximates the array $S$.

By preserving an $(\alpha,\beta)$-RS $(\min,+)$-convolution instead of an exact $(\min,+)$-convolution, we can reduce the storage space and the computation time.

2.2 A fast algorithm for $(\alpha,\beta)$-RS $(\min,+)$-convolution

Next, we give a simple algorithm RSMinPlus to compute an $(\alpha,\beta)$-RS $(\min,+)$-convolution of two given $\alpha$-RSs $\hat{A}$ and $\hat{B}$.

RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$: We first compute the sum of every pair $(\hat{a}[i_{v}],\hat{b}[j_{t}])$ where $\hat{a}[i_{v}]\in\hat{A},\hat{b}[j_{t}]\in\hat{B}$. Define $\tilde{s}[l_{r}]=\min_{(i_{v},j_{t})\in\Phi_{l_{r}}}\hat{a}[i_{v}]+\hat{b}[j_{t}],$ where $\Phi_{l_{r}}=\{(i_{v},j_{t})\mid i_{v}+j_{t}=l_{r},\hat{a}[i_{v}]\in\hat{A},\hat{b}[j_{t}]\in\hat{B}\}$. Suppose that there are $m$ different elements $\tilde{s}[l_{r}]$. Then we sort $\tilde{s}[l_{r}]$ in the increasing order of the index number $l_{r}$. After sorting, we obtain a monotone increasing array $\tilde{S}=(\tilde{s}[l_{1}],\tilde{s}[l_{2}],\ldots,\tilde{s}[l_{m}]),l_{r}<l_{r+1}\text{ for }r\in[m-1]$. Finally, we construct a $\beta$-RS $\hat{S}$ from $\tilde{S}$ as our solution. Our construction is as follows.

Initially append $\hat{s}[l_{m}]=\tilde{s}[l_{m}]$ to $\hat{S}$. Let $\theta=\hat{s}[l_{m}]/(1+\beta)$. 2. 2.

Sequentially consider all elements in $\tilde{S}$ in decreasing order of the index number. If $\tilde{s}[l_{r}]\leq\theta$, append $\hat{s}[l_{r}]=\tilde{s}[l_{r}]$ to $\hat{S}$. Let $\theta=\hat{s}[l_{r}]/(1+\beta)$. Otherwise, ignore $\tilde{s}[l_{r}]$ and consider the next element $\tilde{s}[l_{r-1}]\in\tilde{S}$. 3. 3.

Return the final sequence $\hat{S}$.

Lemma 18**.**

Suppose $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}],\ldots,\hat{a}[i_{m_{1}}])$ and $\hat{B}=(\hat{b}[j_{1}],\hat{b}[j_{2}],\ldots,\hat{b}[j_{m_{2}}])$. Let $m=\max\{m_{1},m_{2}\}$. RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$ computes an $(\alpha,\beta)$-RS $(\min,+)$-convolution $\hat{S}$ of $\hat{A}$ and $\hat{B}$ in $O(m^{2}\log m)$ time.

Proof.

The correctness is not hard. Let $A^{\prime}$ and $B^{\prime}$ be the completions of $\hat{A}$ and $\hat{B}$ respectively. Let $S^{\prime}$ be the exact $(\min,+)$-convolution of two arrays $A^{\prime}$ and $B^{\prime}$. By the construction of $\tilde{S}$, we have that $S^{\prime}$ is the completion of $\tilde{S}$. Moreover, the $\beta$-RS $\hat{S}$ approximates the array $S^{\prime}$, which proves the correctness. It remains to prove the running time.

By the algorithm RSMinPlus, it takes $m_{1}\cdot m_{2}=O(m^{2})$ time to compute all $\tilde{s}[l_{r}]$. Thus, there are at most $m^{2}$ different elements $\tilde{s}[l_{r}]$ in $\tilde{S}$. Then we need $O(m^{2}\log m)$ time to sort all $\tilde{s}[l_{r}]$ and obtain the array $\tilde{S}$. Finally, scanning all elements in $\tilde{S}$ to construct the $\beta$-RS $\hat{S}$ needs $O(m^{2})$ time. Overall, the runtime of the algorithm is $O(m^{2}\log m)$. ∎

FastRSMinPlus**. A more careful $(\alpha,\beta)$-RS $(\min,+)$-convolution algorithm:** In the tree sparsity model, we have some additional conditions which can improve the running time of the algorithm RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$. We consider the case that $\hat{A}$ are nonnegative sequences with each element $\hat{a}\in\hat{A}$ satisfying that either $\hat{a}=0$ or $\hat{a}\geq 1$. We have the same assumption on $\hat{B}$. Moreover, $\alpha$ and $\beta$ are two constants such that $0<\beta\leq\alpha\leq 1$. The intuition is as follows. Suppose we have just appended some element $\hat{s}[l]$ to the array $\hat{S}$. By the property of $(\alpha,\beta)$-RS $(\min,+)$-convolution, we can safely ignore all $\tilde{s}[l_{r}]$ with $\hat{s}[l]/(1+\beta)<\tilde{s}[l_{r}]\leq\hat{s}[l]$, and focus on finding the largest $\tilde{s}[l_{r}]\in\tilde{S}$ such that $\tilde{s}[l_{r}]\leq\hat{s}[l]/(1+\beta)$ (see the definition of $\tilde{S}$ in RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$).

We first build a hash table $\mathsf{HashB}$. The construction is as follows. Each element in $\mathsf{HashB}$ is a pair $(key,value)$ where the $key$ term is an integer satisfying that $-1\leq key\leq\left\lceil\log_{1+\beta}\hat{b}[j_{m_{2}}]\right\rceil$, and the $value$ term is the largest index $w$ satisfying that $\hat{b}[w]\in\hat{B}$ and $\hat{b}[w]\leq(1+\beta)^{key}$. 222If $\hat{b}[j_{1}]=0$, we define $\mathsf{HashB}(-1)=\hat{b}[j_{1}]$. Otherwise if $\hat{b}[j_{1}]\geq 1$, for a $key$ term, if we have $\hat{b}[j_{1}]>(1+\beta)^{key}$, we ignore this $key$ term when constructing $\mathsf{HashB}$. Symmetrically, we also build a hash table $\mathsf{HashA}$ for the array $\hat{A}$. Let $U=\max\{\hat{a}[i_{m_{1}}],\hat{b}[j_{m_{2}}]\}$. Since both $\hat{A}$ and $\hat{B}$ are increasing sequences, we can construct hash tables $\mathsf{HashA}$ and $\mathsf{HashB}$ in $O(\log_{1+\beta}U)$ time by considering all $key$ terms in increasing order.

Given an element $\hat{s}[l]\in\hat{S}$, we show how to find the largest $\tilde{s}[l_{r}]\in\tilde{S}$ such that $\tilde{s}[l_{r}]\leq\hat{s}[l]/(1+\beta)$ by hash tables. We first reduce this problem to finding the element $\hat{b}[j_{t}]\in\hat{B}$ of the largest index $j_{t}$ for each $\hat{a}[i_{v}]\in\hat{A}$, such that $\hat{a}[i_{v}]+\hat{b}[j_{t}]\leq\hat{s}[l]/(1+\beta)$. A simple scheme is to enumerate all $\hat{a}[i_{v}]\in\hat{A}$, query the hash table $\mathsf{HashB}$, and find the largest $\hat{b}[j_{t}]\in\hat{B}$ such that $\hat{b}[j_{t}]\leq\hat{s}[l]/(1+\beta)-\hat{a}[i_{v}]$. Then among all such $(i_{v},j_{t})$ index pairs, we choose the pair $(i_{v^{*}},j_{t^{*}})$ with the largest sum $i_{v^{*}}+j_{t^{*}}$. We append $\hat{s}[i_{v^{*}}+j_{t^{*}}]=\tilde{s}[i_{v^{*}}+j_{t^{*}}]=\hat{a}[i_{v^{*}}]+\hat{b}[j_{t^{*}}]$ to $\hat{S}$. However, enumerating all elements is not necessary, since we have the following lemma.

Lemma 19**.**

Let $\tau=\left\lceil 1/\alpha\right\rceil$. For any $\hat{a}[i_{v-\tau}],\hat{a}[i_{v}]\in\hat{A}$, we have $\hat{a}[i_{v-\tau}]\leq\hat{a}[i_{v}]/2$. Similarly, for any $\hat{b}[j_{t-\tau}],\hat{b}[j_{t}]\in\hat{B}$, $\hat{b}[j_{t-\tau}]\leq\hat{b}[j_{t}]/2$.

Proof.

W.l.o.g., we only consider the array $\hat{A}$. By Definition 17, we know $(1+\alpha)\hat{a}[i_{v-1}]\leq\hat{a}[i_{v}]$. If $\alpha\geq 1$, $\tau=1$, the lemma is trivially true. Otherwise if $\alpha<1$, we have $\hat{a}[i_{v-\tau}]\leq\hat{a}[i_{v}]/(1+\alpha)^{1/\alpha}\leq\hat{a}[i_{v}]/2$. ∎

Let $\theta=\hat{s}[l]/(1+\beta)$. Assume that $\hat{a}[i_{v}]\in\hat{A}$ (resp. $\hat{b}[j_{t}]\in\hat{B}$) is the largest element such that $\hat{a}[i_{v}]\leq\theta$ (resp. $\hat{b}[j_{t}]\leq\theta$). Hence, both $\hat{a}[i_{v+1}]$ and $\hat{b}[j_{t+1}]$ are at least $\theta$. Therefore, the pair $(i_{v^{*}},j_{t^{*}})$ must satisfy that either $v-\tau\leq v^{*}\leq v$ or $t-\tau\leq t^{*}\leq t$, since $\hat{a}[i_{v-\tau}]+\hat{b}[j_{t-\tau}]\leq\theta$ by Lemma 19. Thus, we only need to consider at most $\tau+1$ elements in $\hat{A}$ or $\hat{B}$. Note that we can directly find such index $i_{v}$ (resp. $j_{t}$) by the following lemma.

Lemma 20**.**

Either $i_{v}=\mathsf{HashA}(\left\lceil\log_{1+\beta}\theta\right\rceil)$ or $i_{v}=\mathsf{HashA}(\left\lfloor\log_{1+\beta}\theta\right\rfloor)$. Similarly, either $j_{t}=\mathsf{HashB}(\left\lceil\log_{1+\beta}\theta\right\rceil$ or $j_{t}=\mathsf{HashB}(\left\lfloor\log_{1+\beta}\theta\right\rfloor$

Proof.

W.l.o.g., we take $i_{v}$ as example. Let $i_{w}=\mathsf{HashA}(\left\lceil\log_{1+\beta}\theta\right\rceil)$. If $i_{v}\neq i_{w}$, then by the definition of $i_{v}$, we have that $\theta<\hat{a}[i_{w}]\leq(1+\beta)^{\left\lceil\log_{1+\beta}\theta\right\rceil}$. Since $\hat{A}$ is an $\alpha$-RS and $\alpha\geq\beta$, we have

[TABLE]

Thus, we have $i_{v}=i_{w-1}$ by the definition of $i_{v}$. Note that $i_{w-1}=\mathsf{HashA}(\left\lfloor\log_{1+\beta}\theta\right\rfloor)$. We finish the proof. ∎

We summarize our approaches in Algorithm 1. The analysis of the algorithm is as follows.

Lemma 21**.**

Let $U=\max\{\hat{a}[i_{m_{1}}],\hat{b}[j_{m_{2}}]\}$. Let $0\leq\beta\leq\alpha\leq 1$ be two constants. The algorithm $\mathsf{FastRSMinPlus}(\alpha,\beta,\hat{A},\hat{B})$ computes an $(\alpha,\beta)$-RS $(\min,+)$-convolution $\hat{S}$ of $\hat{A}$ and $\hat{B}$ in $O\left(\frac{\log_{1+\beta}U}{\alpha}\right)$ time.

Proof.

We first prove the correctness. Let $A$ and $B$ be the completions of $\hat{A}$ and $\hat{B}$ respectively. Let $S$ be the $(\min,+)$-convolution of $A$ and $B$. Assume that we just append an element $\hat{s}[l]$ to $\hat{S}$ and let $\theta=\hat{s}[l^{\prime}]/(1+\beta)$. Consider the next recursion from Line 5 to Line 20, we append a new element $\hat{s}[l]$ to $\hat{S}$. Initially in Line 6, we find the largest index $j_{t}$ satisfying that $\hat{b}[j_{t}]+\hat{a}[i_{1}]\leq\theta$ by Lemma 20. We first analyse the first loop from Line 8 to Line 12. In Line 10, we find the largest element $\hat{a}[w]$ satisfying that $\hat{a}[w]\leq\theta-\hat{b}[j_{t-\delta}]$ by Lemma 20. Thus, we conclude that $\hat{a}[w]+\hat{b}[j_{t-\Delta}]\leq\theta$ for any $1\leq\Delta\leq\min\{t,\tau\}$. Similarly, we can prove that $\hat{b}[w]+\hat{a}[i_{v-\Delta}]\leq\theta$ for any $1\leq\Delta\leq\min\{v,\tau\}$ in Line 17. Thus, the element $\hat{s}[l]$ always satisfies that $\hat{s}[l]\leq\theta=\hat{s}[l^{\prime}]/(1+\beta)$ during the recursion. Thus, the output $\hat{S}$ is a $\beta$-RS.

On the other hand, let $S^{\prime}$ be the completion of $\hat{S}$. Assume that $s[l]=\hat{a}[i_{v^{*}}]+\hat{b}[j_{t^{*}}]$ is the largest element in $S$ satisfying that $s[l]\leq\theta$. W.l.o.g. we assume that $\theta/2\leq\hat{b}[j_{t^{*}}]\leq\theta-\hat{a}[i_{1}]$ (otherwise $\theta/2\leq\hat{a}[i_{v^{*}}]\leq\theta-\hat{b}[j_{1}]$). We conclude that $t-\tau\leq t^{*}\leq t$ by Lemma 19. Then we must consider the element $\hat{b}[j_{t^{*}}]$ in the loop from Line 8 to Line 12. Note that in Line 10, we find an index $w=i_{v^{*}}$ by Lemma 20. By the updating rules in Line 11-12, we update $\hat{s}[l]=\hat{a}[i_{v^{*}}]+\hat{b}[j_{t^{*}}]$ in Line 12 and append $\hat{s}[l]$ to $\hat{S}$ in Line 20. Considering any element $s^{\prime}[l_{0}]$ with $l+1\leq l_{0}\leq l^{\prime}$, we have that $s^{\prime}[l_{0}]=\hat{s}[l^{\prime}]$ by Definition 15. Moreover, we have the following inequality by the chosen of $l_{r}$,

[TABLE]

Overall, we prove that $S^{\prime}$ is a $\beta$-approximation of $S$ by Definition 16.

Then we analyze the running time. By the definition of $U$, we always have $\theta\leq 2U$. After each iteration, the value $\theta$ decreases by a factor at least $1+\beta$ by the fact that $\hat{S}$ is a $\beta$-RS. Thus, there are at most $\left\lceil\log_{1+\beta}2U\right\rceil$ iterations. For each iteration, we first find $j_{t}$ in $O(1)$ time by Lemma 20. Then we consider at most $\tau+1=\left\lceil 1/\alpha\right\rceil+1$ possible index pairs $(w,j_{t-\Delta})$. We only cost $O(1)$ time for each index pair. For the loop from Line 13 to Line 19, we have the same analysis. Thus, the running time of each iteration is $O(\left\lceil 1/\alpha\right\rceil)$. Overall, the total running time is at most $O(\frac{\log_{1+\beta}U}{\alpha})$. ∎

By Lemma 21, the running time of Algorithm 1 is determined by the term $U=\max\{\hat{a}[i_{m_{1}}],$ $\hat{b}[j_{m_{2}}]\}$. In fact, if $\log_{1+\beta}U$ is larger than the largest index number $M=\max\{i_{m_{1}},j_{m_{2}}\}$ of arrays, we can improve the running time of Algorithm 1 further. The main difference is that we do not use hash tables since it takes $\log_{1+\beta}U$ time for construction. The details are as follows.

Compute the completion $A^{\prime}=(a^{\prime}[0],a^{\prime}[1],\ldots,a^{\prime}[i_{m_{1}}])$ and $B^{\prime}=(b^{\prime}[0],b^{\prime}[1],\ldots,b^{\prime}[j_{m_{2}}])$ of $\hat{A}$ and $\hat{B}$ respectively. 2. 2.

Compute the $(\min,+)$-convolution $S$ of $A^{\prime}$ and $B^{\prime}$ as follows. Let $\tau=\left\lceil 1/\alpha\right\rceil$. Sequentially consider each $L\in[0,i_{m_{1}}+j_{m_{2}}]$ in the increasing order. For a term $L$, find $\hat{a}[i_{v}]\in\hat{A}$ with the largest index satisfying that $i_{v}\leq L$. Similarly, find $\hat{b}[j_{t}]\in\hat{B}$ with the largest index satisfying that $j_{t}\leq L$. 3. 3.

Compute $s[L]=\min\{\min_{0\leq\Delta\leq\min\{v,\tau\}}(a^{\prime}[i_{v-\Delta}]+b^{\prime}[L-i_{v-\Delta}]),\min_{0\leq\Delta\in\leq\min\{t,\tau\}}(a^{\prime}[L-j_{t-\Delta}]+b^{\prime}[j_{t-\Delta}])\}$. 4. 4.

Scan the array $S$ in decreasing order. Construct a $\beta$-RS as in Algorithm RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$

Note that our approach is similar to the iteration of Algorithm RSMinPlus$(\alpha,\beta,\hat{A},\hat{B})$. While we use the properties of $\alpha$-RS during computing $S$, similar to Algorithm 1. The running time of Step 1 is $O(M)$. For each $L$, since we consider $L$ sequentially, it costs $O(1)$ time to find indexes $i_{v}$ and $j_{t}$. Moreover, we cost $O(\tau)$ time to compute $s[L]$. Finally, the running time of Step 4 is $O(M)$. Thus, the total running time is $O(M/\alpha)$. Combining with Lemma 21, we have the following lemma.

Lemma 22**.**

Consider two $\alpha$-RSs $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}],\ldots,\hat{a}[i_{m_{1}}])$ where each $\hat{a}[i_{w}]$ ($1\leq w\leq m_{1}$) satisfies that either $\hat{a}[i_{w}]=0$ or $\hat{a}[i_{w}]\geq 1$, and $\hat{B}=(\hat{b}[j_{1}],\hat{b}[j_{2}],\ldots,\hat{b}[j_{m_{2}}])$ where each $\hat{b}[j_{w}]$ ($1\leq w\leq m_{2}$) satisfies that either $\hat{b}[j_{w}]=0$ or $\hat{b}[j_{w}]\geq 1$. Let $U=\max\{\hat{a}[i_{m_{1}}],\hat{b}[j_{m_{2}}]\}$ and $M=\max\{i_{m_{1}},j_{m_{2}}\}$. Let $0\leq\beta\leq\alpha\leq 1$ be two constants. There exists an algorithm $\mathsf{FastRSMinPlus}(\alpha,\beta,\hat{A},\hat{B})$ computing an $(\alpha,\beta)$-RS $(\min,+)$-convolution $\hat{S}$ of $\hat{A}$ and $\hat{B}$ in $O\left(\min\left\{\frac{\log_{1+\beta}U}{\alpha},\frac{M}{\alpha}\right\}\right)$ time.

3 Tree Sparsity Model

In this section, we discuss the tree sparsity model. We first introduce some essential definitions and techniques such as weight discretization and RS $(\min,+)$-convolution. Using these new techniques, we will give an $O(\epsilon^{-1}n\log n)$ time algorithm. Then we speed up the algorithm to $O(\epsilon^{-1}n(\log\log\log n)^{2})$ time through a faster algorithm for RS $(\min,+)$-convolution and the lookup table method. Our improved algorithm is appropriate for both Tree-Sparsity-Head and Tree-Sparsity-Tail. Moreover, we show that we can obtain a linear time algorithm for Tree-Sparsity-Head by a more careful weight discretization technique. For Tree-Sparsity-Tail, the new weight discretization technique is not suitable. Instead, we give a linear time algorithm for Tree-Sparsity-Tail under the assumption that $k\leq n^{1-\delta}$ ($\delta\in(0,1]$ is a fixed constant) by a pruning technique.

For convenience, we only consider the perfect binary tree in this section. Our algorithm can be naturally extended to the general complete $b$-ary tree sparsity model. We defer the details in Appendix A. In this section, we only consider the $l_{1}$-norm for both Tree-Sparsity-Head and Tree-Sparsity-Tail. Hence, we only consider the case that each node weight $x_{i}\geq 0$. We will see our algorithm can be easily generalized to general $l_{p}$-norm. Again, We defer the details in Appendix A.

We denote the given perfect binary tree by $T$. Consider a node in the tree $T$. Suppose the number of edges on the path between the node and the root is $t$. Define the level of the node by $(\log(n+1)-t)$. For example, all leaves are at level $1$ and the root node is at level $\log(n+1)$. 333Note that $\log(n+1)$ is an integer since $T$ is a perfect binary tree. For each level of $T$, we sort all nodes in the same level by a BFS. We denote by $N_{ij}$ the $j$th node at level $i$. We call the subtree with root $N_{ij}$ and containing all nodes rooted at $N_{ij}$ the largest subtree of $N_{ij}$ and denote it by $T_{ij}$. Note that the left child of $N_{ij}$ is $N_{i-1,2j-1}$ and the right child is $N_{i-1,2j}$.

Assume that each node $N_{ij}$ has a weight $x_{ij}$. Recall that in the tree sparsity model, each support $\Omega\in\mathbb{T}_{k}$ is a subtree of $T$ rooted at the root node $N_{\log(n+1),1}$ with $k$ nodes. For a node $N_{ij}$ and a subtree $\Omega\in\mathbb{T}_{k}$, we use $N_{ij}\notin\Omega$ to denote $N_{ij}\in T\setminus\Omega$. In this section, we first consider the Tree-Sparsity-Tail version. The Tree-Sparsity-Head version is similar to Tree-Sparsity-Tail, and we will show the differences later. We denote the optimal solution of the Tree-Sparsity-Tail problem by $\Omega^{*}$ together with an optimal tail value $\mathsf{OPT}=\sum_{N_{ij}\notin\Omega^{*}}x_{ij}$. We also denote the solution of our algorithm by $\hat{\Omega}$ together with a tail value $\mathsf{SOL}=\sum_{N_{ij}\notin\hat{\Omega}}x_{ij}$. W.l.o.g., we assume that $k\geq\log n$. Otherwise we can safely ignore those nodes $N_{i}$ of depth larger than $k$. We also consider the error parameter $\epsilon>0$ as a constant.

3.1 A Nearly Linear Time Algorithm for Tree-Sparsity-Tail

We first propose a scheme for the tail-approximation projection problem for the general case. We first assume that each node weight $x_{ij}$ is an integer among $[0,\frac{n\log n}{\epsilon}+n]$. Thus there are at most $O(n\log n/\epsilon)$ different weight values. We can remove this assumption by a weight discretization technique, see Appendix B for details. We then introduce a look-up table method, which is inspired by the well known Four Russians Method[21]. Combining FastRSMinPlus and the look-up table method, we give a nearly linear time algorithm for Tree-Sparsity-Tail.

Encoding low levels by the look-up table method: In fact, we can further discretize the weight such that there are at most $O(\log n/\epsilon)$ different discretized weight. Define $\hat{x}_{ij}=(1+\epsilon)^{\left\lceil\log_{1+\epsilon}x_{ij}\right\rceil}$ as the discretized weight of node $N_{ij}$. Therefore, $x_{ij}\leq\hat{x}_{ij}<(1+\epsilon)x_{ij}$. Suppose that $s[k]=\sum_{N_{ij}\notin\Omega^{*}}x_{ij}$ is the optimal tail value for Tree-Sparsity-Tail, where $\Omega^{*}$ is the optimal support using node weights $\{x_{ij}\}$. Suppose that $\hat{\Omega}$ is the optimal support for Tree-Sparsity-Tail using discretized weights $\{\hat{x}_{ij}\}$. We have the following inequality

[TABLE]

Thus, we use the discretized weights in the following. Now we have at most $O(\log_{1+\epsilon}(n\log n/\epsilon))=O(\log n/\epsilon)$ different weights. Consider any node $N_{\xi j}$ at level $\xi=\left\lceil\log\log n-\log(1/\epsilon)-\log\log\log n\right\rceil$. The largest subtree $T_{\xi j}$ rooted at $N_{\xi j}$ has at most $m=\left\lceil\log n/(\epsilon\log\log n)\right\rceil$ nodes. We can compute its exact tail array with running time at most $\sum_{i\in[1,\xi]}(m+1)\cdot 2^{-i}\cdot 2^{2i}=O(\log^{2}n/\epsilon^{2})$ by computing exact $(\min,+)$-convolution level by level. Since we have at most $O(\log n/\epsilon)$ different node weights after discretization, there are at most $O(\log n/\epsilon)^{m}=O(n^{O(\epsilon)})$ possible constructions for $T_{\xi j}$. 444Here, each construction is a weight assignment of all nodes in $T_{\xi j}$. By this observation, we can enumerate all possible constructions and compute the corresponding exact tail array using $O(n^{O(\epsilon)}\log^{2}n/\epsilon^{2})=o(n)$ time and $o(n)$ space. Thus, we encode all possible constructions of subtrees at level $\xi$ into a look-up table. When we need compute the exact tail array of any node at level $\xi$, we search the look-up table and return the array in $O(m)$ time.

Now we are ready to give our algorithm for Tree-Sparsity-Tail. For each node $N_{ij}$, define an array $S_{ij}=(s[0],s[1],s[2],\ldots,s[2^{i}-1])$ to be the exact tail array, where each element $s[l]$ represents the optimal tail value for Tree-Sparsity-Tail on $T_{ij}$, i.e., $s[l]=\min_{\Omega\in\mathbb{T}_{2^{i}-1-l}(T_{ij})}\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}\setminus\Omega}x_{i^{\prime}j^{\prime}}$. In the exact algorithm, we in fact compute the exact tail array $S_{ij}$ for each node $N_{ij}$ through the $(\min,+)$-convolution. Our main technique is to maintain an $\alpha$-RS  $\hat{S}_{ij}$ for each $S_{ij}$. The value of $\alpha$ depends on the level $i$, which will be decided later.

FastTailTree**:** In our algorithm, we use $\mathsf{MinPlus}$ to represent the $O(m^{2}/2^{c\sqrt{\log m}})$ algorithm for exact $(\min,+)$-convolutions mentioned in [26] ($c>0$ is some fixed constant). We divide the whole tree $T$ into three parts as follows.

Step 1:

Let $\xi=\left\lceil(\log\log n-\log(1/\epsilon)-\log\log\log n)\right\rceil,\eta=\left\lceil(\log\log n+\log(1/\epsilon))\right\rceil$. For any node $N_{\xi j}$ at level $\xi$, we use the look-up table method to obtain the exact tail array $S_{\xi j}$. Compute an $\epsilon^{\prime}$-RS $\hat{S}_{\xi j}$ which approximates $S_{\xi j}$, where $\epsilon^{\prime}=\epsilon/(\eta-\xi+1)$.

Step 2:

From level $i=\xi+1$ to level $i=\eta$, we use $\mathsf{FastRSMinPlus}(\epsilon^{\prime},\epsilon^{\prime},\hat{S}_{i-1,2j-1},\hat{S}_{i-1,2j})$ to compute a sequence $\hat{S}_{ij}$ for any node $N_{ij}$, and append $\hat{s}[2^{i}-1]=(1+\epsilon^{\prime})^{i-\xi+1}\cdot\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$ to $\hat{S}_{ij}$. 44footnotetext: We set $\hat{s}[2^{i}-1]$ to be this value, because we want to guarantee that $\hat{S}_{ij}$ is still an $\epsilon^{\prime}$-RS after appending $\hat{s}[2^{i}-1]$. This is convenient for analyzing the algorithm in Theorem 24. Compute an $\epsilon$-RS $\hat{S}_{\eta j}$ which approximates $S_{\eta j}$ for any node $N_{\eta j}$.

Step 3:

From level $i=\eta+1$ to level $i=\log(n+1)$, we use $\mathsf{FastRSMinPlus}(\epsilon_{i-1},\epsilon_{i},\hat{S}_{i-1,2j-1},\hat{S}_{i-1,2j})$ to compute a sequence $\hat{S}_{ij}$ for any node $N_{ij}$, where $\epsilon_{i}=\frac{\epsilon}{3^{(i-\eta)/4}}$ ($\eta\leq i\leq\log(n+1)$). We then append $\hat{s}[2^{i}-1]=(1+\epsilon^{\prime})^{\eta-\xi+1}\cdot\prod_{l=\eta}^{i}(1+\epsilon_{l})\cdot\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$ to $\hat{S}_{ij}$.

Step 4:

Let $\hat{s}[L]\in\hat{S}_{\log(n+1),1}$ be the smallest element of index $L$ satisfying that $L\geq n-k$. Output $\hat{\Omega}\leftarrow$ FindTree($L,T$). Here, $\mathsf{FindTree}$ is a backtracking process with running time $O(n)$ which obtains a feasible solution $\hat{\Omega}$. We defer the details in Algorithm 3 in Appendix A.

Before analyzing FastTailTree, we give some intuitions about why we compute $(\alpha,\beta)$-RS $(\min,+)$-convolutions. Note that the weight $x_{i}$ of each node is an integer at most $O(n\log n/\epsilon)$. Thus, the maximum value in each $S_{ij}$ is at most $O(n^{2}\log n/\epsilon)$. In our algorithm, we use a sequence $\hat{S}_{ij}$ to approximate $S_{ij}$. By Definition 16 and 17, the number of elements in $\hat{S}_{ij}$ is at most $\log_{(1+\epsilon_{i})}(n^{2}\log n/\epsilon)=O(\log n/\epsilon_{i})$, which means that the size of $\hat{S}_{ij}$ is sublinear on $n$. Thus, if the level $i$ is high enough, the array $\hat{S}_{ij}$ maintains much fewer elements than $S_{ij}$, and can be constructed faster. By Lemma 22, we have the following corollary.

Corollary 23**.**

At Step 2, $\mathsf{FastRSMinPlus}(\epsilon^{\prime},\epsilon^{\prime},\hat{S}_{i-1,2j-1},\hat{S}_{i-1,2j})$ can be computed in time $O(\frac{2^{i-1}}{\epsilon^{\prime}})$. At Step 3, $\mathsf{FastRSMinPlus}(\epsilon_{i-1},\epsilon_{i},\hat{S}_{i-1,2j-1},\hat{S}_{i-1,2j})$ can be computed in time $O(\frac{\log n}{\epsilon_{i-1}\epsilon_{i}})$.

Proof.

For each node $N_{ij}$, the largest index of $\hat{S}_{i-1,2j-1}$ or $\hat{S}_{i-1,2j}$ is at most $M\leq 2^{i-1}$. On the other hand, the maximum value $U$ in $\hat{S}_{i-1,2j-1}$ or $\hat{S}_{i-1,2j}$ is at most $O(n^{2}\log n/\epsilon)$. By Lemma 22, we prove the corollary. ∎

Now we are ready to give the following main theorem.

Theorem 24**.**

Algorithm FastTailTree is a $(1+\epsilon)$-approximation algorithm with running time $O(\epsilon^{-1}n(\log\log\log n)^{2})$ for Tree-Sparsity-Tail.

Proof.

We first prove the running time. There are $(n+1)/2^{\xi}$ nodes at level $\xi$. We need $O(2^{\xi})$ time to compute each exact tail array through searching the look-up table and need $O(2^{\xi})$ time to compute an $\epsilon^{\prime}$-RS. Thus, the runtime of Step 1 is $O(2^{-\xi}(n+1)\cdot 2^{\xi})=O(n)$. Considering Step 2, the running time for each node $N_{ij}$ is at most $O(2^{i}/\epsilon^{\prime})$. Thus, the total time of Step 2 is

[TABLE]

For Step 3, the running time for each node $N_{ij}$ is at most $O(\log n/(\epsilon_{i}\epsilon_{i-1}))$ by Corollary 23. Thus, the total time of Step 3 is

[TABLE]

The running time of Step 4 is $O(n)$. Overall, the total running time is $O(n+\epsilon^{-1}n(\log\log\log n)^{2}+\epsilon^{-1}n)=O(\epsilon^{-1}n(\log\log\log n)^{2})$.

Then we prove the correctness by showing that our solution $\hat{\Omega}$ is a $(1+\epsilon)$-approximation for the optimal solution $\Omega^{*}$. We first prove by induction that for each node $N_{ij}$ at level $\xi\leq i\leq\eta$, the array $\hat{S}_{ij}$ is a \bigl{(}(1+\epsilon^{\prime})^{i-\xi+1}-1\bigr{)}-RS which approximates the exact tail array $S_{ij}$. The base case at level $\xi$ is true since each sequence $\hat{S}_{\xi j}$ is an $\epsilon^{\prime}$-RS of the exact tail array $S_{\xi j}$ at Step 1. Then we suppose that for level $i-1$ ($\xi+1\leq i\leq\eta$), any sequence $\hat{S}_{i-1,j}$ is an \epsilon^{*}=\bigl{(}(1+\epsilon^{\prime})^{i-\xi}-1\bigr{)}-RS which approximates the array $S_{i-1,j}$. We consider an arbitrary node $N$ and its sequence $\hat{S}$ at level $i$. Let $S^{\prime}_{1}$ be the completion of the sequence $\hat{S}_{1}$ maintained by $N$’s left child. Let $S^{\prime}_{2}$ be the completion of the sequence $\hat{S}_{2}$ maintained by $N$’s right child. Let $S_{1}$ and $S_{2}$ be the exact tail arrays of $N$’s left and right children respectively. Let $S$ be the exact $(\min,+)$-convolution of $S_{1}$ and $S_{2}$, i.e., $S$ is the exact tail array of $N$ without the last term $s[2^{i}-1]$. By induction, we know that the two arrays $S^{\prime}_{1}$ and $S^{\prime}_{2}$ are $\epsilon^{*}$-approximations of $S_{1}$ and $S_{2}$ respectively. Let $\tilde{S}$ be the exact $(\min,+)$-convolution of $S^{\prime}_{1}$ and $S^{\prime}_{2}$. Let $S^{\prime}$ be the completion of $\hat{S}$ (without the element $\hat{s}[2^{i}-1]$). By Definition 16 and 17, we have that $S^{\prime}$ is an $\epsilon^{\prime}$-approximation of $\tilde{S}$.

Consider any element $\tilde{s}[l]\in\tilde{S}$ such that $\tilde{s}[l]=a^{\prime}[l_{1}]+b^{\prime}[l_{2}]$ for $l_{1}+l_{2}=l,a^{\prime}[l_{1}]\in S^{\prime}_{1},b^{\prime}[l_{2}]\in S^{\prime}_{2}$. By induction, we have that $a[l_{1}]\leq a^{\prime}[l_{1}]$ for $a[l_{1}]\in S_{1}$ and $b[l_{2}]\leq b^{\prime}[l_{2}]$ for $b[l_{2}]\in S_{2}$. Therefore, we have that

[TABLE]

The last inequality follows from the fact that $S^{\prime}$ is an $\epsilon^{\prime}$-approximation of $\tilde{S}$. On the other hand, consider any element $s[l]\in S$ such that $s[l]=a[l_{1}]+b[l_{2}]$ for $l_{1}+l_{2}=l,a[l_{1}]\in S_{1},b[l_{2}]\in S_{2}$. By induction, we have that $a^{\prime}[l_{1}]\leq(1+\epsilon^{\prime})a[l_{1}]$ for $a^{\prime}[l_{1}]\in S^{\prime}_{1}$ and $b^{\prime}[l_{2}]\leq(1+\epsilon^{*})b[l_{2}]$ for $b^{\prime}[l_{2}]\in S^{\prime}_{2}$. Thus, we conclude that

[TABLE]

By the above argument, $S^{\prime}$ is a \bigl{(}(1+\epsilon^{\prime})(1+\epsilon^{*})-1\bigr{)}-approximation of $S$. More specifically, we have the following inequality

[TABLE]

Now we consider the element $\hat{s}[2^{i}-1]$ appended to $\hat{S}$ at Step 2. On one hand, since the exact tail value $s[2^{i}-1]=\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$, we have that $s[2^{i}-1]\leq\hat{s}[2^{i}-1]\leq(1+\epsilon^{\prime})(1+\epsilon^{*})s[2^{i}-1]$. On the other hand, we have

[TABLE]

The last inequality follows from the fact that $\hat{s}[2^{i}-1]=(1+\epsilon^{\prime})(1+\epsilon^{*})s[2^{i}-1]$. Thus, we conclude that $\hat{S_{ij}}$ is still an $\epsilon^{\prime}$-RS approximating the exact tail array $S_{ij}$, which proves the induction.

By a similar reduction, we can prove that $\hat{S}_{\log(n+1),1}$ is a \bigl{(}(1+\epsilon^{\prime})^{\eta-\xi+1}\cdot\prod_{l=\eta}^{\log(n+1)}(1+\epsilon_{l})-1\bigr{)}-RS which approximates $S_{\log(n+1),1}$. Overall, the approximation ratio for the root array $\hat{S}_{\log(n+1),1}$ is $1+O(\epsilon)$. Therefore, let $S^{\prime}_{\log(n+1),1}$ be the completion of the sequence $\hat{S}_{\log(n+1),1}$ maintained in the root node. Let $s^{\prime}[n-k]\in S^{\prime}_{\log(n+1),1}$, we have that $\hat{s}[L]\leq(1+\epsilon_{\log(n+1)})s^{\prime}[n-k]\leq(1+O(\epsilon))s[n-k]$ for $s[n-k]\in S_{\log(n+1),1}$. By using a small enough value $\theta(\epsilon)$ to replace $\epsilon$, we can guarantee that the value $\hat{s}[L]$ is a $(1+\epsilon)$-approximation tail value for Tree-Sparsity-Tail.

∎

3.2 A Linear Time Algorithm for Tree-Sparsity-Tail if $k\leq n^{1-\delta}$

For a special case that $k\leq n^{1-\delta}$ for some fixed constant $\delta\in(0,1]$, we can further improve the running time to linear for Tree-Sparsity-Tail. Note that in practice, this is a reasonable assumption which generalizes the assumption $k\leq n^{1/2-\delta}$ in the previous work [16]. Our main approach is to show that we can safely ignore many nodes at low levels.

We divide the tree into two parts. Let $\eta=\left\lceil 2\log\log n\right\rceil$. The first part is from level $1$ to $\eta$ and the second part is from level $(\eta+1)$ to $\log(n+1)$. For the second part, we still use FastRSMinPlus algorithm to maintain an approximate tail array. The difference is that for the first part, we show that we only need to consider at most $O(n^{1-\delta}/\epsilon)$ nodes. Recall that $T_{ij}$ is the perfect binary subtree rooted at $N_{ij}$, and $u_{ij}=\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$ is the total subtree weight of $N_{ij}$. Note that there are at most \bigl{(}(n+1)/\log^{2}n\bigr{)} nodes at level $\eta$. Let $u$ be the $\left\lceil(1+\epsilon)n^{1-\delta}/\epsilon\right\rceil$-largest total subtree weight among these nodes $\{N_{\eta j}\}_{j}$. We argue that we can safely ignore all subtrees $T_{\eta j}$ if its corresponding total subtree weight $u_{\eta j}<u$. The details can be found in Algorithm 2.

Theorem 25**.**

Algorithm 2 is a $(1+\epsilon)$-approximation algorithm with running time $O(n+\epsilon^{-2}n/\log n)$ for Tree-Sparsity-Tail if $k\leq n^{1-\delta}$.

Proof.

We first prove the correctness. Let $\mathcal{C}=\{T_{\eta j}:\ u_{\eta j}\geq u\}$ be the collection of those subtrees with total subtree weight at least $u$. Let $\bar{\mathcal{C}}=\{T_{\eta j}\}_{j}\setminus\mathcal{C}$ be the complement of $\mathcal{C}$. We argue that the influence caused by deleting the subtrees in $\bar{\mathcal{C}}$ in Step 2 is negligible. Let $\Omega^{*}$ be the optimal support with the optimal tail value $\mathsf{OPT}=\sum_{N_{i}\notin\Omega^{*}}x_{i}$. Let $\check{\Omega}$ be the optimal subtree of the case, in which we delete all subtrees in $\bar{\mathcal{C}}$. Let $\check{\mathsf{OPT}}=\sum_{N_{i}\notin\check{\Omega}}x_{i}$. Note that our algorithm obtains a $(1+\epsilon)$-approximation $\hat{\Omega}$ of $\check{\Omega}$ following from the analysis in Theorem 24 and Corollary 23.

Thus, we only need to prove that $\check{\mathsf{OPT}}\leq(1+\epsilon)\mathsf{OPT}$. By the assumption that $k\leq n^{1-\delta}$, $\Omega^{*}$ contains at most $n^{1-\delta}$ nodes in $\bar{\mathcal{C}}$, which have a total weight at most $n^{1-\delta}u$. It means that $\check{\mathsf{OPT}}-\mathsf{OPT}\leq n^{1-\delta}u$. On the other hand, there are at least $n^{1-\delta}/\epsilon$ subtrees in $\mathcal{C}$ that do not intersect $\Omega^{*}$, since $\Omega^{*}$ can contain at most $n^{1-\delta}$ nodes in $\mathcal{C}$. Thus, we have that $\mathsf{OPT}\geq n^{1-\delta}u/\epsilon$. Hence, $\check{\mathsf{OPT}}-\mathsf{OPT}\leq\epsilon\mathsf{OPT}$ which proves the correctness.

Then we analyze the running time. It costs $O(n)$ time to compute all $u_{\eta j}$ and $u$ in Step 2. For each node $N_{ij}$ at level $1\leq i\leq\eta$, it costs $O(2^{2i-c\sqrt{i}})$ time to compute $\hat{S_{ij}}$ using the procedure MinPlus. Among each subtree in $\mathcal{C}$, the number of nodes at level $1\leq i\leq\eta$ is $2^{\eta-i}=O(\log^{2}n/2^{i})$. On the other hand, there are at most $O(n^{1-\delta}/\epsilon)$ trees in $\mathcal{C}$. Thus, the total running time from Step 3 to Step 7 is

[TABLE]

Considering Step 8, it costs $O(2^{\eta})$ time for each node $N_{\eta j}$. Thus, the total running time for Step 8 is $O(\frac{n^{1-\delta}}{\epsilon}\cdot 2^{\eta})=o(n)$. By Corollary 23, the construction time of all $\hat{S}_{ij}$ at level $\eta+1\leq i\leq\log(n+1)$ from Step 9-12 is

[TABLE]

Finally, the backtracking process FindTree$(L,T)$ in Step 14 costs $O(n)$ time . Therefore, the total running time of Algorithm 2 is $O(n+\epsilon^{-2}n/\log n)$. ∎

3.3 A Linear Time Algorithm for Tree-Sparsity-Head

Now we consider the Tree-Sparsity-Head version. Recall that our goal is to find a subtree $\Omega\in\mathbb{T}_{k}$ such that $\sum_{N_{ij}\in\Omega}x_{ij}\geq(1-\epsilon)\sum_{N_{ij}\in\Omega^{*}}x_{ij}$, where $\Omega^{*}$ is the optimal solution of the Tree-Sparsity-Head problem. In this subsection, we denote $\mathsf{OPT}_{H}=\sum_{N_{i}\in\Omega^{*}}x_{i}$ to be the optimal head value for Tree-Sparsity-Head. Our framework is similar to the framework for Tree-Sparsity-Tail. We again assume that each node weight $x_{ij}$ is an integer among $[0,O(n\log n/\epsilon)]$. Thus there are at most $O(n\log n/\epsilon)$ different weight values. Similar to Tree-Sparsity-Tail, we can remove this assumption by a weight discretization technique, see Appendix B for details. By this assumption, we still construct a dynamic program for Tree-Sparsity-Head. However, our techniques and definitions have some differences. We then show the differences in details in the following.

Approximate $(\max,+)$-Convolution: At first, we introduce another concept called $(\max,+)$-convolution which is similar to $(\min,+)$-convolution (see Definition 13).

Definition 26** ($(\max,+)$-convolution).**

Given two arrays $A=(a[0],a[1],a[2],\ldots,a[m_{1}])$ and $B=(b[0],b[1],b[2],\ldots,b[m_{2}])$, their $(\max,+)$-convolution is the array $S=(s[0],s[1],s[2],\ldots,$ $s[m_{1}+m_{2}])$ where $s[t]=\max_{i=0}^{t}\{a[i]+b[t-i]\},t\in[0,m_{1}+m_{2}]$.

The only difference from $(\min,+)$-convolution is that $s_{t}=\max_{i}\{a_{i}+b_{t-i}\}$. In fact, these two definitions are equivalent. Suppose that $-S=(-s[0],-s[1],\ldots,-s[m_{1}+m_{2}])$ is the $(\min,+)$-convolution of $-A=(-a[0],-a[1],\ldots,$ $-a[m_{1}])$ and $-B=(-b[0],-b[1],\ldots,-b[m_{2}])$. Then $S=(s[0],s[1],\ldots,s[m_{1}+m_{2}])$ is exactly the $(\max,+)$-convolution of two arrays $A$ and $B$.

For each node $N_{ij}$ on the tree $T$, we define $S_{ij}=(s[0],s[1],\ldots,s[2^{i}-1])$ to be the head array of $N_{ij}$, where each element $s[l]$ represents the optimal head value for Tree-Sparsity-Head on $T_{ij}$, i.e., $s[l]=\max_{\Omega\in\mathbb{T}_{l}(T_{ij})}\sum_{N_{i^{\prime}j^{\prime}}\in T_{ij}}x_{i^{\prime}j^{\prime}}$. In fact, the array $S_{ij}$ can be achieved through computing the $(\max,+)$-convolution from the arrays $S_{i-1,2j-1}$ and $S_{i-1,2j}$ of its two children.666Note that for each element $s[l]\in S_{ij}$ $(1\leq l\leq 2^{i}-1)$, we have that $s[l]=x_{ij}+\max_{t}\{a[t]+b[l-1-t]\}$ where $a[t]\in S_{i-1,2j-1}$ and $b[l-1-t]\in S_{i-1,2j}$. Similar to Tree-Sparsity-Tail, our key approach is to maintain a head sequence $\hat{S}_{ij}$ as an approximation of $S_{ij}$ which reduces the running time. We first introduce some concepts to describe $\hat{S}_{ij}$.

Definition 27** (Head-Completion of $\alpha$-RS).**

Consider an $\alpha$-RS $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}]\ldots,\hat{a}[i_{m}])$. Define its head-completion of cardinality $M$ by an array $A^{\prime}=(a^{\prime}[0],a^{\prime}[1],\ldots,a^{\prime}[M])$ satisfying that: 1) If $0\leq t\leq i_{1}-1,a^{\prime}[t]=0$; 2) If $i_{v}\leq t\leq i_{v+1}-1\ (1\leq v\leq m-1),a^{\prime}[t]=\hat{a}[i_{v}]$; 3) If $i_{m}\leq t\leq M,a^{\prime}[t]=\hat{a}[i_{m}]$.

Definition 28** (Head-Sequence Approximation).**

Given two $n$-length non-decreasing arrays $A^{\prime}=(a^{\prime}[0],a^{\prime}[1],\ldots,a^{\prime}[n])$ and $A=(a[0],a[1],\ldots,a[n])$, we say $A^{\prime}$ is an $\alpha$-head-approximation of $A$ if for any $i$, $(1-\alpha)a[i]\leq a^{\prime}[i]\leq a[i]$. We say an $\alpha$-RS $\hat{A}$ head-approximates an array $A$ if its head-completion $A^{\prime}$ of cardinality $n$ is an $\alpha$-head-approximation of $A$.

Figure 2 illustrates these definitions. Note that the above definitions have some differences from in Tree-Sparsity-Tail. By comparing Figure 1 and 2, we can see the differences. Now we are ready to define the concept of $(\alpha,\beta)$-RS $(\max,+)$-convolution.

Definition 29** ($(\alpha,\beta)$-RS $(\max,+)$-convolution).**

Given two $\alpha$-RSs $\hat{A}$ and $\hat{B}$, suppose $A^{\prime}$ and $B^{\prime}$ are their head-completions of cardinality $M_{1}$ and $M_{2}$ respectively. Suppose the array $S$ is the $(\max,+)$-convolution of $A^{\prime}$ and $B^{\prime}$. We call $\hat{S}$ an $(\alpha,\beta)$-RS $(\max,+)$-convolution of $\hat{A}$ and $\hat{B}$ if $\hat{S}$ is a $\beta$-RS which head-approximates the array $S$.

Similar to Lemma 22, we have the following lemma.

Lemma 30**.**

Consider two $\alpha$-RSs $\hat{A}=(\hat{a}[i_{1}],\hat{a}[i_{2}],\ldots,\hat{a}[i_{m_{1}}])$ where each $\hat{a}[i_{w}]$ ($1\leq w\leq m_{1}$) satisfies that either $\hat{a}[i_{w}]=0$ or $\hat{a}[i_{w}]\geq 1$, and $\hat{B}=(\hat{b}[j_{1}],\hat{b}[j_{2}],\ldots,\hat{b}[j_{m_{2}}])$ where each $\hat{b}[j_{w}]$ ($1\leq w\leq m_{2}$) satisfies that either $\hat{b}[j_{w}]=0$ or $\hat{b}[j_{w}]\geq 1$. Let $U=\max\{\hat{a}[i_{m_{1}}],\hat{b}[j_{m_{2}}]\}$ and $M=\max\{i_{m_{1}},j_{m_{2}}\}$. Let $0\leq\beta\leq\alpha\leq 1$ be two constants. There exists an algorithm $\mathsf{FastRSMaxPlus}(\alpha,\beta,\hat{A},\hat{B})$ computing an $(\alpha,\beta)$-RS $(\max,+)$-convolution $\hat{S}$ of $\hat{A}$ and $\hat{B}$ in $O(\min\{\frac{\log_{1+\beta}U}{\alpha},\frac{M}{\alpha}\})$ time.

Using the same scheme as Algorithm FastTailTree in Section 3.1, we can design a $(1-\epsilon)$-approximation algorithm for Tree-Sparsity-Head with running time $O(\epsilon^{-1}n(\log\log\log n)^{2})$. One difference is that we compute $\hat{S}_{ij}$ by an approximate $(\max,+)$-convolution scheme FastRSMaxPlus by Lemma 30. The other difference is that after we compute the sequence $\hat{S}_{\log(n+1),1}$ for the root node, we find the largest element $\hat{s}[L]\in\hat{S}_{\log(n+1),1}$ of index $L$ satisfying that $L\leq k$ and return a solution $\hat{\Omega}$ by a backtracking process.

A linear time algorithm for Tree-Sparsity-Head: In fact, we can improve the running time to linear by some additional properties of Tree-Sparsity-Head. Let $\xi=\lceil\log\log n-\log(1/\epsilon)$ $-\log\log\log n\rceil$, $\eta=\left\lceil\log\log n+\log(1/\epsilon)\right\rceil$. As in Algorithm FastTailTree, the time cost of the second part (i.e., from level $\xi+1$ to level $\eta$) is the bottleneck. Fortunately for Tree-Sparsity-Head, we can speed up the second part by a new weight discretization technique. Recall that $\epsilon$ is a constant number.

Then we show how to compute an array $\hat{S}_{\eta j}$ as an $\epsilon$-head-approximation of $S_{\eta j}$ for all nodes $N_{\eta j}$ in linear time. We first divide the array $\hat{S}_{\eta j}$ into two sub-arrays. One sub-array consists of the first $\left\lceil 2\log\log n\right\rceil$ elements $\hat{s}[l]$ $(0\leq l\leq\left\lceil 2\log\log n\right\rceil-1$. The other sub-array consists of the remaining elements $\hat{s}[l]$ $(\left\lceil 2\log\log n\right\rceil\leq l\leq 2^{\eta}-1)$. In the following, we show how to compute these two sub-arrays respectively.

Case 1, $\hat{s}[l]\in\hat{S}_{\eta j},0\leq l\leq\left\lceil 2\log\log n\right\rceil-1$: We still compute $\hat{S}_{\xi j}$ through the look-up table method as in Step 1 of Algorithm FastTailTree. Then for any node $N_{ij}$ at level $\xi+1\leq i\leq\eta$, we construct a sub-array $\hat{S}_{ij}$ by computing an exact $(\max,+)$-convolution from its two children, while we only compute the first $\left\lceil 2\log\log n\right\rceil$ elements $\hat{s}[l]\in\hat{S}_{ij}$ $(0\leq l\leq\left\lceil 2\log\log n\right\rceil-1)$.

Case 2, $\hat{s}[l]\in\hat{S}_{\eta j},\left\lceil 2\log\log n\right\rceil\leq l\leq 2^{\eta}-1$: We consider a more careful weight discretization. Consider the perfect binary subtree $T_{\eta j}$ rooted at some node $N_{\eta j}$. Let $N_{\max}\in T_{\eta j}$ be the node of the largest weight $x_{\max}$. Consider an element $s[l]\in S_{\eta j}$ ($l\geq\left\lceil 2\log\log n\right\rceil\geq\eta$) representing the optimal head value of sparsity $l$ for Tree-Sparsity-Head on $T_{\eta j}$. Then we have that $s[l]\geq x_{\max}$, since there exists a subtree rooted at $N_{\eta j}$ with $l$ nodes and containing node $N_{\max}$. 777Note that this property is only satisfied in Tree-Sparsity-Head. We define the new discretized weight for each node $N_{i^{\prime}j^{\prime}}\in T_{\eta j}$ to be $\hat{x}_{i^{\prime}j^{\prime}}=\left\lfloor\frac{x_{i^{\prime}j^{\prime}}\log^{2}n}{\epsilon x_{\max}}\right\rfloor$. After weight discretizing, each node weight in $T_{\eta j}$ is an integer among the range $\left[0,\left\lfloor\log^{2}n/\epsilon\right\rfloor\right]$.

Based on these node weights $\{\hat{x}_{i^{\prime}j^{\prime}}\}$, we again use the look-up table method to compute all arrays $\hat{S}_{\xi j}$ at level $\xi$. Similar to Step 2 of Algorithm FastTailTree, we compute an $(\epsilon^{\prime},\epsilon^{\prime})$-RS $(\max,+)$-convolution from level $\xi+1$ to level $\eta$, where $\epsilon^{\prime}=\epsilon/(\eta-\xi)$. Now for a node $N_{\eta j}$, we obtain an approximate sequence and we compute its head-completion $\hat{S}_{\eta j}$ of cardinality $2^{\eta}-1$. Finally for each element $\hat{s}[l]\in\hat{S}_{\eta j},\left\lceil 2\log\log n\right\rceil\leq l\leq 2^{\eta}-1$, we multiply it by a normalization factor $\epsilon x_{\max}/\log^{2}n$.

Overall, we combine the above two sub-arrays, and obtain an approximate array $\hat{S}_{\eta j}$. We will prove that $\hat{S}_{\eta j}$ is an $\epsilon$-head-approximation of $S_{\eta j}$. Then we compute an $\epsilon$-RS which head-approximates $\hat{S}_{\eta j}$ for each node $N_{\eta j}$. Finally we use the similar technique as in Step 3 of Algorithm FastTailTree. For any node $N_{ij}$ at level $\eta+1\leq i\leq\log(n+1)$, we use Algorithm $\mathsf{FastRSMaxPlus}(\epsilon_{i-1},\epsilon_{i},\hat{S}_{i-1,2j-1},\hat{S}_{i-1,2j})$ to compute a sequence $\hat{S}_{ij}$, where $\epsilon_{i}=3^{-(i-\eta)/4}\epsilon$ ($\eta\leq i\leq\log(n+1)$).

Theorem 31**.**

There is a $(1-\epsilon)$-approximation algorithm with running time $O(\epsilon^{-1}n)$ for Tree-Sparsity-Head.

Proof.

We first consider the running time. For Case 1, the running time for computing all $\hat{S}_{\eta j}$ is $o(n)$. For each node $N_{ij}$ at level $\xi+1\leq i\leq\eta$, since we only compute $\left\lceil 2\log\log n\right\rceil$ elements, the running time for constructing the sub-array is $O(\log^{2}\log n)$. Note that there are at most O(n/2^{\xi})=O\bigl{(}n\log\log n/(\epsilon\log n)\bigr{)} nodes. Thus, the total running time for Case 1 is $o(n)$. For Case 2, using the lookup table method costs $o(n)$ time. For each node $N_{ij}$ at level $\xi+1\leq i\leq\eta$, the running time for computing an $(\epsilon^{\prime},\epsilon^{\prime})$-RS $(\min,+)$-convolution is $O(\log\log n/\epsilon^{\prime 2})$. by Lemma 30. Thus, the total running time for this case is $\sum_{i\in(\xi,\eta]}2^{-i}n\cdot\log\log n/\epsilon^{\prime 2}=o(n)$. For those nodes $N_{ij}$ at level $\eta+1\leq i\leq\log(n+1)$, by the same analysis in Theorem 24, the total running time is $O(\epsilon^{-1}n)$. Overall, the running time is $O(\epsilon^{-1}n)$.

Then we prove the approximation ratio. For Case 1, by Definition 26, each element $\hat{s}[l]\in\hat{S}_{\xi j}$ $(0\leq l\leq\left\lceil 2\log\log n\right\rceil-1)$ satisfies that $(1-\epsilon)s[l]\leq\hat{s}[l]\leq s[l]$. For Case 2, we can show that the new weight discretization scheme leads to a $(1-\epsilon)$-approximation following from the same argument as in Lemma 40. Then by the same argument as in Theorem 24, we have that the processes $\mathsf{FastRSMaxPlus}$ from level $\eta$ to $\log(n+1)$ compute an $(1-\epsilon)$-head-approximation array for Tree-Sparsity-Head. Thus, the total approximation ratio is $(1-\epsilon)$. ∎

Combining Theorem 24,  25 and 31, we obtain Theorem 9.

3.4 Compressive Sensing Recovery

By Theorem 9, we can obtain a faster tree sparse recovery algorithm by the framework AM-IHT in [17]. The framework AM-IHT is an iterative scheme. In each iteration, we need to complete two matrix multiplications, a head-approximation, and a tail-approximation projections.

Theorem 10 (Restated).

Assume that $k\leq n^{1-\delta}$ ($\delta\in(0,1)$ is any fixed constant). Let $A\in\mathbb{R}^{m\times n}$ be a measurement matrix. Let $x\in\mathcal{M}_{k}$ be an arbitrary signal in the tree sparsity model with dimension $n$, and let $y=Ax+e\in\mathbb{R}^{m}$ be a noisy measurement vector. Here $e\in\mathbb{R}^{m}$ is a noise vector. Then there exists an algorithm to recover a signal approximation $\hat{x}\in\mathcal{M}_{k}$ satisfying $\|x-\hat{x}\|\leq C\|e\|_{2}$ for some constant $C$ from $m=O(k)$ measurements. Moreover, the algorithm runs in $O((n\log n+k^{2}\log n\log^{2}(k\log n))\log\frac{\|x\|_{2}}{\|e\|_{2}})$ time.

Proof.

Our theorem is very similar to Theorem 3 in [16]. The only difference is that we output a solution $\hat{x}\in\mathbb{M}_{k}$ instead of $\hat{x}\in\mathbb{M}_{ck}$ for some constant $c>1$. That is because the final solution is obtained from a tail oracle, and our tail oracle is a single-criterion oracle. ∎

4 CEMD Model

In this section, we discuss another structured sparsity model known as the Constrained EMD model [24].

4.1 A Single-Criterion Approximation Algorithm for Head-Approximation Projection

We develop a single-criterion constant approximation algorithm for the head approximation projection in the CEMD model, improving the result in [17] which relaxes the support space to $\Omega\in\mathbb{M}_{k,B\log k}$. We first use an EMD flow network [17], and similarly obtain two supports $\Omega_{l}$ and $\Omega_{r}$. Then from these two supports, we construct a single-criterion constant factor approximate solution. Formally speaking, given an arbitrary signal $x$, we want to find a support $\hat{\Omega}\in\mathbb{M}_{k,B}$ such that $\sum_{x_{i,j}\in\hat{\Omega}}|x_{i,j}|^{p}\geq c\cdot\max_{\Omega\in\mathbb{M}_{k,B}}\sum_{x_{i,j}\in\Omega}|x_{i,j}|^{p}$ for some fixed constant $c\in(0,1]$

Step 1: Constructing an EMD flow network: We first recall the EMD flow network construction defined in [17]. See Figure 3 as an example.

Definition 32** (EMD flow network).**

For a given signal $X$, sparsity $k$, and a parameter $\lambda>0$, the flow network $G_{X,k,\lambda}$ is defined as follows:

Each entry $x_{i,j}\in X$ corresponds to a node $v_{i,j}$ for $i\in[h],j\in[w]$. Additionally, add a source node $\mu$ and a sink node $\nu$. 2. 2.

Add an edge from every $v_{i_{1},j}$ to every $v_{i_{2},j+1}$ for $i_{1},i_{2}\in[h],j\in[w-1]$. Moreover, add an edge from the source to every $v_{i,1}$ and from every $v_{i,w}$ to the sink. 3. 3.

The capacity on every edge and node (except source and sink) is 1. 4. 4.

The cost of node $v_{i,j}$ is $-|x_{i,j}|^{p}$. The cost of an edge from $v_{i_{1},j}$ to $v_{i_{2},j+1}$ is $\lambda|i_{1}-i_{2}|$. The cost of the source, the sink, and each edge incident to the source or sink is 0. 5. 5.

Both the supply at the source and the demand at the sink are $s(=\frac{k}{w})$.

Since all edge capacities, the source supply, and the sink demand are integers, by Theorem 9.10 in [23], we know that $G_{X,k,\lambda}$ always has an integer min-cost max-flow. Note that this integer min-cost max-flow must be a set of disjoint paths through the network $G_{X,k,\lambda}$, and it corresponds to a support in $X$. For a flow network $G_{X,k,\lambda}$, we denote the support of this integer min-cost max-flow by $\Omega_{\lambda}=\mathsf{MinCostFlow}(G_{X,k,\lambda})$. Thus, for any $\lambda$, a solution of the min-cost max-flow problem on $G_{X,k,\lambda}$ reveals a subset $S$ of nodes that corresponds to a support $\Omega_{\lambda}$ satisfying the following two properties: 1) in each column, $\Omega_{\lambda}$ has exactly $s$ indices; 2) $\Omega_{\lambda}$ is the support which minimizes $-\sum_{x_{i,j}\in\Omega}|x_{i,j}|^{p}+\lambda\mathbf{EMD}[\Omega]$ (also equivalent to maximize $\sum_{x_{i,j}\in\Omega}|x_{i,j}|^{p}-\lambda\mathbf{EMD}[\Omega]$).

For convenience, we define $\Phi[\Omega]$ to be the head value $\sum_{x_{i,j}\in\Omega}|x_{i,j}|^{p}$ of support $\Omega$ and denote the $\mathbf{EMD}[\Omega]$ by $\Delta[\Omega]$. In [17], we can obtain the following theorem by this flow network.

Theorem 33** (Theorem 34 and 36 in [17]).**

Let $\delta>0$, $x_{\min}=\min_{|X_{i,j}|>0}|X_{i,j}|^{p}$, and $x_{\max}=\max_{|X_{i,j}|>0}|X_{i,j}|^{p}$. There exists an algorithm running in $O(snh(\log\frac{n}{\delta})+\log\frac{x_{\max}}{x_{\min}})$ time, which returns two solutions $\Omega_{l}=\mathsf{MinCostFlow}(G_{X,k,l})$, and $\Omega_{r}=\mathsf{MinCostFlow}(G_{X,k,r})$. We have that $l,r\geq 0$, $l-r\leq\frac{\delta x_{\min}}{wh^{2}}$, and $\Delta[\Omega_{l}]\leq B\leq\Delta[\Omega_{r}]$.

Then we show how to construct a single-criterion solution by $\Omega_{l}$ and $\Omega_{r}$.

Step 2: Constructing a single-criterion solution: By Theorem 33, assume that we have two solutions $\Omega_{l}$ and $\Omega_{r}$ now. We want to construct a single-criterion solution which is also a constant approximation. Note that $\Omega_{l}\in\mathbb{M}_{k,B}$ and $\Omega_{r}$ may not be in $\mathbb{M}_{k,B}$. We first construct a single-criterion solution $\Omega_{r}^{\prime}$ based on $\Omega_{r}$ such that $\Phi[\Omega_{r}^{\prime}]\geq\Phi[\Omega_{r}]\cdot(2(\left\lfloor\Delta[\Omega_{r}]/B\right\rfloor+1))^{-1}$. We need the following lemma for preparation.

Lemma 34**.**

Given any path $P$ on the flow network $G_{X,k,\lambda}$ from source to sink, let $\Omega_{P}$ be the support of $P$. Let $d\geq 1$ be some positive integer. There exists an $O(n)$ time algorithm which finds another path $P^{\prime}$ with support $\Omega_{P^{\prime}}$ satisfying that $\Delta[\Omega_{P^{\prime}}]\leq\Delta[\Omega_{P}]/d$, and $\Phi[\Omega_{P^{\prime}}]\geq\Phi[\Omega_{P}]/2d$.

Proof.

W.l.o.g., assume that the lowest node on path $P$ is at row 1. Consider the row $L_{t}$ which separates the lowest $t$ rows and the upper $h-t$ rows. Row $L_{t}$ decomposes the path $P$ into two paths $\check{P}_{t}$ and $\hat{P}_{t}$. Specifically, for any edge $(v_{i_{1},j},v_{i_{2},j+1})$, we add two edges in $\check{P}_{t}$ and $\hat{P}_{t}$ respectively as follows.

• •

If $i_{1}>t$ and $i_{2}>t$, we add the edge $(v_{i_{1},j},v_{i_{2},j+1})$ in $\hat{P}_{t}$ and add the edge $(v_{t,j},v_{t,j+1})$ in $\check{P}_{t}$. Similarly, if $i_{1}\leq t$ and $i_{2}\leq t$, we add the edge $(v_{t+1,j},v_{t+1,j+1})$ in $\hat{P}_{t}$ and add the edge $(v_{i_{1},j},v_{i_{2},j+1})$ in $\check{P}_{t}$.

• •

If $i_{1}>t$ and $i_{2}\leq t$, we add the edge $(v_{i_{1},j},v_{t+1,j+1})$ in $\hat{P}_{t}$ and add the edge $(v_{t,j},v_{i_{2},j+1})$ in $\check{P}_{t}$. Similarly, if $i_{1}\leq t$ and $i_{2}>t$, we add the edge $(v_{t+1,j},v_{i_{2},j+1})$ in $\hat{P}_{t}$ and add the edge $(v_{i_{1},j},v_{t,j+1})$ in $\check{P}_{t}$.

See Figure 4 as an example. Suppose for an edge $(v_{i_{1},j},v_{i_{2},j+1})$ in $P$, we add an edge $(v_{\hat{i}_{1},j},v_{\hat{i}_{2},j+1})$ in path $\hat{P}_{t}$ and an edge $(v_{\check{i}_{1},j},v_{\check{i}_{2},j+1})$. It is not difficult to check that $|i_{2}-i_{1}|\geq|\hat{i}_{1}-\hat{i}_{2}|+|\check{i}_{1}-\check{i}_{2}|$. Moreover, $\Omega_{\hat{P}_{t}}\cup\Omega_{\check{P}_{t}}\supset\Omega_{P}$. Thus,

[TABLE]

Also observe that as $t$ increases, $\Delta[\Omega_{\check{P}_{t}}]$ is non-decreasing.

If $\Delta[\Omega_{P}]$, then the path $P$ itself satisfies the lemma. Thus, we assume that $\Delta[\Omega_{P}]>0$. We then prove the lemma by induction on $d$. If $d=1$, the path $P$ itself satisfies the lemma. Suppose the lemma is true for any positive integer no more than $d-1$. Now we consider the integer $d$.

We first find the highest row $L_{t}$ such that $\Delta[\Omega_{\check{P}_{t}}]\leq\Delta[\Omega_{P}]/d$, and $\Delta[\Omega_{\check{P}_{t+1}}]>\Delta[\Omega_{P}]/d$. Note that such an index $t$ must exist, since

[TABLE]

Note that row $L_{t}$ is a path with $\Delta[\Omega_{l_{t}}]=0$ where $\Omega_{l_{t}}$ is the support of $L_{t}$. We distinguish three cases.

If $\Phi[\Omega_{\check{P}_{t}}]\geq\Phi[\Omega_{P}]/2d$, then the path $P_{t}$ satisfies the lemma. 2. 2.

If $\Phi[\Omega_{\check{P}_{t}}]<\Phi[\Omega_{P}]/2d$ and $\Phi[{\Omega_{L_{t+1}}}]\geq\Phi[\Omega_{P}]/2d$, then the path $L_{t+1}$ satisfies the lemma. 3. 3.

If $\Phi[\Omega_{\check{P}_{t}}]<\Phi[\Omega_{P}]/2d$ and $\Phi[{\Omega_{L_{t+1}}}]<\Phi[\Omega_{P}]/2d$, then we have that $\Phi[{\Omega_{\check{P}_{t+1}}}]<\Phi[\Omega_{P}]/d$ and $\Delta[{\Omega_{\check{P}_{t+1}}}]>\Delta[\Omega_{P}]/d$. Thus, according to Inequalities 1, we have that

[TABLE]

and

[TABLE]

By induction, we can find a path $P^{\prime}$ from path $\hat{P}_{t+1}$, such that

[TABLE]

By the above discussion, we prove the lemma. ∎

Note that $\Omega_{r}$ consists of $s$ disjoint paths. According to Lemma 34, we can construct a single-criterion solution $\Omega_{r}^{\prime}$ as follows.

Corollary 35**.**

Let $d=\left\lfloor\Delta[\Omega_{r}]/B\right\rfloor$. We can construct a support $\Omega_{r}^{\prime}\in\mathcal{M}_{k,B}$ such that $\Phi[\Omega_{r}^{\prime}]\geq\Phi[\Omega_{r}]/2(d+1)$ in $O(ns)$ time.

We next compare $\Phi[\Omega_{r}^{\prime}]$ with $\Phi[\Omega_{l}]$. If $\Phi[\Omega_{r}^{\prime}]>\Phi[\Omega_{l}]$, then we output $\Omega_{r}^{\prime}$ as our solution. Otherwise, we output $\Omega_{l}$ as our solution. By the following lemma, we show that our solution is a constant approximation.

Lemma 36**.**

Suppose $\mathsf{OPT}=\max_{\Omega\in\mathbb{M}_{k,B}}\Phi[\Omega]$. Then we have that

[TABLE]

Proof.

Recall that $\mathsf{MinCostFlow}$ solves the min-cost max-flow $\Omega_{\lambda}$ of the graph $G_{X,k,\lambda}$, i.e.,

[TABLE]

The value of objective is no less than [math] for any $\lambda$ because there exists some support $\Omega$ such that $\Delta[\Omega]=0$ and $\Phi[\Omega]\geq 0$ for any $\Omega$.

We get $\Omega_{l}$ and $\Omega_{r}$ from $\mathsf{MinCostFlow}$ algorithm for $\lambda$ equaling to $l$ and $r$ respectively. Thus, we have $\Phi[\Omega_{r}]-r\cdot\Delta[\Omega_{r}]\geq 0$. Moreover, $\Phi[\Omega_{r}]\geq\mathsf{OPT}$. Suppose $\Phi[\Omega^{*}]=\mathsf{OPT}$. If $\Phi[\Omega_{r}]<\mathsf{OPT}$, changing $\Omega_{r}$ to $\Omega^{*}$ would increase the objective $\Phi[\Omega_{r}]-r\cdot\Delta[\Omega_{r}]$ since $\Delta[\Omega_{r}]\geq B$, which yields a contradiction.

Assume that $\Delta[\Omega_{r}]\in[dB,(d+1)B)$ for some positive integer $d\geq 1$. We distinguish three cases.

If $d=1$, by Corollary 35, $\Omega_{r}^{\prime}$ satisfies $\Phi[\Omega_{r}^{\prime}]\geq\Phi[\Omega_{r}]/2\geq\mathsf{OPT}/2$ and $\Delta[\Omega_{r}^{\prime}]\leq\Delta[\Omega_{r}]/2\leq B$. 2. 2.

If $d\geq 2$ and $\Phi[\Omega_{r}]\geq 3d\mathsf{OPT}/4$, by Corollary 35, $\Omega_{r}^{\prime}$ satisfies $\Phi[\Omega_{r}^{\prime}]\geq\Phi[\Omega_{r}]/2(d+1)\geq 3d\mathsf{OPT}/8(d+1)\geq\mathsf{OPT}/4$ and $\Delta[\Omega_{r}^{\prime}]\leq\Delta[\Omega_{r}]/(d+1)\leq B$. 3. 3.

If $d\geq 2$ and $\Phi[\Omega_{r}]<3d\mathsf{OPT}/4$, we have that $\mathsf{OPT}>4rB/3$ since $\Phi[\Omega_{r}]\geq r\cdot\Delta[\Omega_{r}]\geq rdB$. Then we have the following inequalities.

[TABLE]

Here, the first inequality follows from the fact

[TABLE]

Besides, $l-r\leq\frac{x_{\min}\delta B}{wh^{2}}$ follows from Theorem 33.

Overall, we prove the lemma. ∎

Combining Theorem 33 and Lemma 36, we have the following theorem.

Theorem 11 (Restated).

Consider the CEMD model $\mathcal{M}_{k,B}$ with $s=k/w$ sparse for each column and support-EMD $B$. Let $\delta\in(0,1/4)$, $x_{\min}=\min_{|X_{i,j}|>0}|X_{i,j}|^{p}$, and $x_{\max}=\max_{|X_{i,j}|>0}|X_{i,j}|^{p}$. Let $c=1/4-\delta$. There exists an algorithm running in $O(shn\log\frac{n}{\delta}+\log\frac{x_{\max}}{x_{\min}})$ time, which returns a single-criterion $c^{1/p}$ approximation for head-approximation projection.

Note that the exponent $1/p$ of $c$ comes from $l_{p}$-norm.

4.2 Compressive Sensing Recovery

Similar to tree sparsity model, our head oracle in Theorem 33 can also lead to a model-based compressive sensing recovery algorithm, combining with AM-IHT and the tail oracle in [17]. We summarize our result as follows.

Theorem 12 (Restated).

Let $A\in\mathbb{R}^{m\times n}$ be a measurement matrix. Let $x\in\mathbb{M}_{k,B}$ be an arbitrary signal in the CEMD model with dimension $n=wh$, and let $y=Ax+e\in\mathbb{R}^{m}$ be a noisy measurement vector. Here $e\in\mathbb{R}^{m}$ is a noise vector. Then there exists an algorithm to recover a signal approximation $\hat{x}\in\mathbb{M}_{k,2B}$ satisfying $\|x-\hat{x}\|\leq C\|e\|_{2}$ for some constant $C$ from $m=O(k\log(B/k))$ measurements. Moreover, the algorithm runs in $O(n\log\frac{\|x\|_{2}}{\|e\|_{2}}(k\log n+\frac{kh}{w}(\log n+\log\frac{x_{\max}}{x_{\min}})))$ time, where $x_{\max}=\max|x_{i}|$ and $x_{\min}=\min_{|x_{i}|>0}|x_{i}|$.

Proof.

Our theorem is very similar to Theorem 37 in [17] except two improvements. The first improvement is that we reduce the number of measurements. That is because the head oracle in [17] outputs a solution in $\mathbb{M}_{k,\gamma B}$, and the number $m$ of measurements has the following bound:

[TABLE]

In [17], $\gamma=O(\log(k/w))$. In contrast, our head oracle confirms that $\gamma=1$ by Theorem 33.

The second improvement is the running time. There are $O(\log\frac{\|x\|_{2}}{\|e\|_{2}})$ iterations in the framework AM-IHT. In each iteration, we need to complete two matrix multiplications, a head-approximation, and a tail-approximation. In [17], the time complexity of a head oracle is $O(\frac{nkhB}{w})$, while the time complexity of our head oracle is exactly the same as the tail oracle in [17] by Theorem 33. ∎

Appendix A Some missing details

Extend the Algorithm for $l_{1}$-norm to $l_{p}$-norm: It is not different to extend our results to the $l_{p}$-norm for both Tree-Sparsity-Head and Tree-Sparsity-Tail. The only difference is that we compute the $l_{p}$-norm weight $|x_{i}|^{p}$ for each node $N_{i}\in T$ at the beginning. Then we run our algorithms for both Tree-Sparsity-Head and Tree-Sparsity-Tail using these $l_{p}$-norm weight $|x_{i}|^{p}$. We will obtain a $(1+\epsilon)^{1/p}$-approximation for Tree-Sparsity-Head and a $(1-\epsilon)^{1/p}$-approximation for Tree-Sparsity-Tail respectively. Thus, we only need to set the value of $\epsilon$ to be $O(p\epsilon)$ instead.

Extend the Algorithm for Binary Tree to $b$-ary: On the other hand, we can extend our algorithms to $b$-ary trees. Note that our algorithms are based on $(\min,+)$-convolutions (or $(\max,+)$-convolutions). Consider any node $N$. We want to compute an approximate sequence $\hat{S}$. In a $b$-ary tree, each node $N$ has $b$ children. Denote them by $N_{1},N_{2},\ldots,N_{b}$. We compute $\hat{S}$ by the following iterations.

[TABLE]

The iteration takes time $b$ times as much as before for a binary tree. Since we assume $b$ is a constant integer, it does not affect the time complexity asymptotically.

FindTree Algorithm: We give the $\mathsf{FindTree(L,T)}$ as follows. By the backtracking process $\mathsf{FindTree}(L,T)$, we obtain a support $\hat{\Omega}$ with a tail value at most $\hat{s}[L]$ since each element $\hat{s}[L]$ is at least as large as the exact tail value $s[L]$ by the algorithm. On the other hand, $|\hat{\Omega}|\leq k$ since $L\geq n-k$.

Then, we analyse the running time of the backtracking process $\mathsf{FindTree}(L,T)$. In fact, for each element $\hat{s}\in\hat{S}$ of index $L$, we can save the two indices $L_{1}$ and $L_{2}$ satisfying the condition in Line 3 of Algorithm 3, during constructing $\hat{S}$ in Algorithm FastTailTree. Thus, we only cost $O(1)$ time for each node in $\mathsf{FindTree}(L,T)$. Then the running time of $\mathsf{FindTree}(L,T)$ is $O(n)$.

Appendix B Weight Discretization in the Tree Sparse Model

Weight Discretization for Tree-Sparsity-Tail: We first introduce a linear time $O(\log n)$-approximation algorithm for Tree-Sparsity-Tail, which offers a criterion to discretize the weight.

For each node $N_{ij}$, denote the largest subtree rooted at $N_{ij}$ by $T_{ij}$. Compute the subtree weight $u_{ij}=\sum_{i,j:N_{ij}\in T_{ij}}x_{ij}$ of all nodes in the subtree $T_{ij}$. Let $u$ be the $k$th largest weight among $\{u_{ij}\}_{i,j}$. 2. 2.

Add all nodes with $u_{ij}>u$ into $\Omega$ directly. Then do a BFS (breath-first-search) on tree $T$ and add all nodes with $u_{ij}=u$ into $\Omega$ until $|\Omega|=k$. Denote $\sum_{N_{ij}\notin\Omega}x_{ij}$ by $W$ and return $W$.

We have $W\leq\log n\cdot\mathsf{OPT}$, which means that $W$ is a $\log n$-approximation for Tree-Sparsity-Tail.

Lemma 37**.**

The above algorithm is a $\log n$-approximation algorithm with running time $O(n)$ for Tree-Sparsity-Tail.

Proof.

Observe that for any two nodes $N_{ij}$ and $N_{i^{\prime}j^{\prime}}$, if $N_{ij}$ is the ancestor of $N_{i^{\prime}j^{\prime}}$, we have $u_{ij}\geq u_{i^{\prime}j^{\prime}}$. Combining this fact and the BFS procedure, we have that the support $\Omega$ is a subtree rooted at $N_{\log(n+1),1}$. Then we analyze the time complexity and approximation ratio.

The weight $u_{ij}$ is the summation of the weights of its left subtree, right subtree and itself. We compute $u_{ij}$ from leafs to root. Hence, it takes $O(1)$ time to compute each $u_{ij}$. Constructing $\Omega$ needs $O(n)$ time since we only do a BFS. Thus, the total running time is $O(n)$.

Finally we prove the approximation ratio. Recall that $\Omega^{*}$ is the optimal subtree rooted at $N_{\log{n+1},1}$. We have the following inequality.

[TABLE]

The last inequality follows from the fact that the algorithm selects the $k$ nodes with the largest weight $u_{ij}$. Note that each node appears in at most $\log n$ different subtrees $T_{ij}$ except the root node, we have

[TABLE]

∎

Next, we show how to discretize the weights. For each node $N_{ij}$, if its node weight $x_{ij}\in[0,W]$, we define $\hat{x}_{ij}=\left\lceil\frac{x_{ij}n\log n}{\epsilon W}\right\rceil$ to be the discretized weight. Otherwise if $x_{ij}>W$, we define $\hat{x}_{ij}=\left\lceil\frac{n\log n}{\epsilon}\right\rceil+n$. By this discretization, we have that each $\hat{x}_{ij}$ is an integer among $[0,\frac{n\log n}{\epsilon}+n]$. Assume that $\hat{\Omega}$ is the optimal solution for Tree-Sparsity-Tail based on the discretized weights $\{\hat{x}_{ij}\}_{ij}$, together with a tail value $\widehat{\mathsf{OPT}}=\sum_{N_{ij}\notin\hat{\Omega}}\hat{x}_{ij}$. Denote $\mathsf{OPT}^{\prime}=\sum_{N_{ij}\notin\hat{\Omega}}x_{ij}$ to be the tail value of $\hat{\Omega}$ based on the original weights $\{x_{ij}\}_{i,j}$. We prove the following lemma.

Lemma 38**.**

$\mathsf{OPT}^{\prime}\leq(1+\epsilon)\mathsf{OPT}$.

Proof.

Since $W\geq\mathsf{OPT}$, those nodes of weight larger than $W$ must be in the optimal solution $\Omega^{*}$. Hence, we have that

[TABLE]

By the construction of $\hat{\Omega}$, we have that

[TABLE]

By the above inequality, we conclude that all nodes of weight larger than $W$ are also in the solution $\hat{\Omega}$. Thus, for any node $N_{ij}\notin\Omega^{*}\cup\hat{\Omega}$, we have that $x_{ij}\leq W$ and $\hat{x}_{ij}=\left\lceil\frac{x_{ij}n\log n}{\epsilon W}\right\rceil\leq\left\lceil\frac{n\log n}{\epsilon}\right\rceil$ . By this observation, we have the following inequality.

[TABLE]

Here the fourth inequality follows from the fact that $\frac{W}{\log n}\leq\mathsf{OPT}$ by Lemma 37. ∎

By Lemma 38, we know that the influence caused by the weight discretization is negligible. Note that all nodes of weight larger than $W$ are in the solution $\hat{\Omega}$. W.l.o.g., we assume that each node is of weight $x_{ij}\leq W$ and $\hat{x}_{ij}\leq\left\lceil\frac{n\log n}{\epsilon}\right\rceil$. From now on, we focus on the discretized weight $\{\hat{x}_{ij}\}_{ij}$. For convenience, we use $x_{ij}$ to represent $\hat{x}_{ij}$.

Weight Discretization for Tree-Sparsity-Head: In order to discretize the weight, we still need to introduce a linear time $O(1/\log n)$-approximation algorithm for Tree-Sparsity-Head as follows.

Let $Q$ be the collection of $\left\lfloor k/\log n\right\rfloor$ nodes with the largest node weights (breaking ties arbitrarily). 2. 2.

For each node $N_{ij}\in Q$, append to the solution $\hat{\Omega}$ all nodes on the path from $N_{ij}$ to the root node. Let $W=\sum_{N_{ij}\in\hat{\Omega}}x_{ij}$.

Then we prove the head value $W\geq\mathsf{OPT}_{H}/3\log n$.

Lemma 39**.**

$\hat{\Omega}$* is an $(1/3\log n)$-approximation for Tree-Sparsity-Head with running time $O(n)$.*

Proof.

Note that the number of nodes in $\hat{\Omega}$ is at most $\log n\cdot\left\lfloor k/\log n\right\rfloor\leq k$. Thus $\hat{\Omega}$ is a feasible solution. On the other hand, assume that the minimum weight of nodes in $Q$ is $w$. Then the head value $W$ is at least $W=\sum_{N_{ij}\in\hat{\Omega}}x_{ij}\geq w\cdot\left\lfloor k/\log n\right\rfloor$, while the optimal solution $\mathsf{OPT}_{H}$ is at most $\sum_{N_{ij}\in\Omega^{*}}x_{ij}<\sum_{N_{ij}\in Q}x_{ij}+k\cdot w\leq W+k\cdot w$. So we can conclude that $W\leq\mathsf{OPT}_{H}<(2\log n+1)W$.

Consider the running time. We cost $O(n)$ time to construct the collection $Q$, and cost $O(|\hat{\Omega}|)=O(k)$ time to construct $\hat{\Omega}$. Overall, the running time is $O(n)$. ∎

Next, we show how to discretize the node weights by the criterion $W$. We define $\hat{x}_{i}$ to be $\left\lfloor\frac{kx_{i}}{\epsilon W}\right\rfloor$. Assume the optimal solution for Tree-Sparsity-Head based on the discretized node weights is $\hat{\Omega}$. Denote $\mathsf{OPT}_{H}^{\prime}=\sum_{N_{ij}\in\hat{\Omega}}x_{ij}$. We next analyze the difference between two solutions $\Omega^{*}$ and $\hat{\Omega}$. We have the following lemma.

Lemma 40**.**

$\mathsf{OPT}_{H}^{\prime}\geq(1-\epsilon)\mathsf{OPT}_{H}$.

Proof.

By the definition of $\mathsf{OPT}_{H}^{\prime}$, we have that

[TABLE]

The second inequality follows from the definition of $\hat{\Omega}$, and the last inequality follows from the fact that $W\leq\mathsf{OPT}_{H}$.

∎

By Lemma 40, we know that the loss caused by the weight discretization is negligible. From now on, we focus on the discretized weights $\{\hat{x}_{ij}\}_{ij}$. For convenience, we use $x_{ij}$ to represent $\hat{x}_{ij}$. In the following, we only consider the case that the weight of each node is at most $(3\log n\cdot W)$. Thus, the weight of each node is an integer among the range $[0,\left\lfloor 3k\log n/\epsilon\right\rfloor]$. Note that by Lemma 39, each node with weight at least $(3\log n\cdot W)$ must appear in the optimal solution $\Omega^{*}$. We can directly append such nodes and the nodes on the path from such nodes to root to our solution. Suppose the number of these nodes are $k^{\prime}$. The problem is reduced to find the $(k-k^{\prime})$ nodes with maximum head value among the remaining nodes which can be solved by the same method.