Improved Algorithms For Structured Sparse Recovery

0.32 MB·

by  Lingxiao Huang

Improved Algorithms For Structured Sparse Recovery
LingxiaoHuang1, YifeiJin1, Jian Li 1, and HaitaoWang 2
∗ †
1TsinghuaUniversity,Beijing100084,China
7
2Utah StateUniversity,Logan,UT 84322,USA
1
0
2
January 23,2017
n
a
J

0
Abstract
2
Itisknownthatcertainstructuresofthesignalinadditiontothestandardnotionofsparsity

] (calledstructuredsparsity)canimprovethesamplecomplexityinseveralcompressivesensing
T
applications.Recently,Hegdeetal.[17]proposedaframework,calledapproximation-tolerant
I
. model-basedcompressivesensing,forrecoveringsignalswithstructuredsparsity.Theirframe-
s
c workrequirestwooracles,the head-andthetail-approximationprojectionoracles. Thetwo
[ oracles should return approximate solutions in the model which is closest to the query sig-

1 nal. In this paper, we consider two structured sparsity models and obtain improved projec-
v tionalgorithms. Thefirstoneisthetreesparsitymodel,whichcapturesthesupportstructure
4 in the wavelet decompositionof piecewise-smoothsignals and images. We proposea linear
7
time (1 − ǫ)-approximation algorithm for head-approximation projection and a linear time
6
(1+ǫ)-approximationalgorithmfortail-approximationprojection.Thebestpreviousresultis
5
0 an O˜(nlogn) time bicriterion head-approximation(tail-approximation)algorithm (meaning
. thattheiralgorithmmayreturna solutionofsparsitylargerthank)byHegdeetal[16]. Our
1
resultprovidesan affirmativeanswertothe openproblemmentionedin thesurveyofHegde
0
7 andIndyk[18]. Asa corollary,we canrecovera constantapproximatek-sparsesignal. The
1 otheristheConstrainedEarthMoverDistance(CEMD)model,whichisusefultomodelthe
: situationwherethepositionsofthenonzerocoefficientsofasignaldonotchangesignificantly
v
i as a function of spatial (or temporal) locations. We obtain the first single criterion constant
X
factorapproximationalgorithmforthehead-approximationprojection[17]. Thepreviousbest
r knownalgorithmisabicriterionapproximation.Usingthisresult,wecangetafasterconstant
a
approximationalgorithmwithfewermeasurementsfortherecoveryprobleminCEMDmodel.
1 Introduction
Weconsidertherobustsparserecovery,animportantproblemincompressivesensing. Thegoalof
robust sparse recovery istorecoverasignal fromasmallnumberoflinearmeasurements. Specif-
∗[email protected]
†[email protected]
1
ically, we call a vector x ∈ Rn k-sparse if it has at most k non-zero entries. The support of x,
denoted bysupp(x) ⊆ [n],contains theindices corresponding tothe nonzero entries inx. Given
the measurement vector y = Ax + e, where A is a measurement matrix, x is a k-sparse signal
and eisanoise vector, the goal istofindasignal estimate xˆ such that kx−xˆk < Ckekforsome
constant approximation factorC > 0.
StructuredSparsityModel: Itiswellknownthatingeneralweneedthenumberofmeasurements
to be Ω(klog(n/k)) for robust sparse recovery, see [9,12]. In practice, the support of x usually
hassomestructuredconstraints, suchastreesparsityandblocksparsity,whichcanhelpreducethe
boundofthenumberofmeasurements.
Definition 1 (Structured Sparsity Model [2]). Let M be a family of supports, i.e., M = {Ω ,Ω ,
1 2
...,Ω } whereeach Ω ⊆ [n]. Thenthe corresponding structured sparsity model Misthe setof
L i
vectorssupported ononeoftheΩ :
i
M= {x ∈ Rd | supp(x) ⊆ ΩforsomeΩ ∈ M}.
To recover such a structured signal x is called structured sparse recovery. Baraniuk et al. [2]
providedageneralframeworkcalledmodel-basedcompressivesensing. Theirframeworkdepends
onamodelprojection oraclewhichisdefinedasfollows.
Definition 2 (Model Projection [2]). Let M be a structured sparsity model. A model projection
oracleforMisanalgorithm P(x) : Rn → MsuchthatΩ = P(x)and
∗
kx−x k = minkx−x k
Ω∗ p Ω p
Ω M
∈
wherex ∈ Rn isthesameasxonthesupport Ω,andiszerootherwise.
Ω
Unfortunately, thebestknownalgorithmsformanyexactmodelprojectionoraclesaretooslow
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.
Theirframeworkonly requires twoapproximation oracles called thehead-andtail-approximation
projection oracles, definedasfollows. LetΩ ∈ Mbetheoptimalsupport ofthemodelprojection
∗
oracleasdefinedinDefinition2.
Definition 3 (Head-Approximation Projection). Let M be a structured sparsity model. A head-
approximation oracleforMisanalgorithm H(x) : Rn → MsuchthatH(x)= Ωand
kx k ≥ c ·kx k
Ω p H Ω∗ p
wherec ∈ (0,1]isafixedconstant.
H
Definition 4 (Tail-Approximation Projection). Let M be a structured sparsity model. A tail-
approximation oracleforMisanalgorithm T(x) : Rn → MsuchthatT(x) = Ωand
kx−x k ≤ c ·kx−x k ,
Ω p T Ω∗ p
wherec ∈ [1,∞)isafixedconstant.
T
2
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 ofthesignal xarearranged asthenodes ofacomplete b-ary treeT rooted atnodeN,
andanyfeasiblesolution isasubtree whichincludes therootofT andisofsizek.
Definition 5 (Tree Sparsity Model). Let T be a complete b-ary tree with n nodes rooted at node
N. T (T) = {Ω ,Ω ,...,Ω } is the family of supports where each Ω is a subtree of T rooted
k 1 2 L i
at N with the number of nodes no more than k. We use T instead of T (T) for short. The
k k
tree-structured sparsity modelT isthesetofsignals supported onsomeΩ ∈ T :
k k
T = {x ∈ Rd |supp(x)⊆ ΩforsomeΩ ∈ T }
k k
For the tree sparsity model, the head- and tail-approximation projection problems reduces to
thefollowingsimple-to-statecombinatorialproblems: forthehead-approximation, wewanttofind
a subtree of size k rooted at N such that the total weight of the subtree is maximized; for the
1
tail-approximation, we want the total weight of the complement of the subtree is minimized. For
convenience, weabbreviate them asTree-Sparsity-Headand Tree-Sparsity-Tailrespectively. If
our solution is a subtree of size at most k, we call it a single-criterion solution. Otherwise if our
solution isasubtreeofsizelargerthank,wecallitabicriterion solution.
Constrained EMD Model: The CEMD model, introduced by Schmidt et al [24], is particularly
useful in 2D image compression and denoising [11,25]. Wefirst introduce the definition of Earth
Mover’sDistance(EMD),alsoknownastheWasserstein metricorMallowsdistance [22].
Definition6(EMD). TheEMDoftwofinitesetsA,B ⊂ Nwith|A|= |B|isdefinedas
EMD(A,B) = min |a−π(a)|,
π:A B
→ aXA
∈
whereπ rangesoverallone-to-one mappingsfromAtoB.
IntheCEMDmodel,thesignalx ∈RncanbeinterpretedasamatrixX ∈ Rh wwithn= hw.
×
By this interpretation, the support of x ∈ Rh w, denoted by supp(x) ⊆ [h]×[w], contains the
×
indices(i,j)(i ∈ [h],j ∈ [w])corresponding tothenonzeroentries inx.
Definition7(Support-EMD). Consider anh×w matrix X. LetΩ ⊆ [h]×[w] bethe support of
amatrixX. DenoteΩ tobethesupportofthecolumniofX. Suppose|Ω | = sfori ∈[w]. Then
i i
theEMDofthesupportΩ(orthesupport-EMD ofX)isdefined as
w 1
−
EMD[Ω] = EMD(Ω ,Ω ).
i i+1
Xi=1
3
Naturally, wehave the following structured sparsity model which contains twoconstraints: 1)
eachcolumniss(= k/w)-sparse, 2)thesupport-EMDisatmostB.
Definition8(ConstrainedEMDModel[24]). LetM bethefamilyofsupports{Ω ⊆ [h]×[w] |
k,B
EMD(Ω) ≤ B, and|Ω | = k/w, fori ∈ [w]}. TheConstrained EMD(CEMD)model M is
i k,B
thesetofsignalssupported onsomeΩ ∈ M :
k,B
M = {x ∈Rd |supp(x) ⊆ ΩforsomeΩ ∈ M }
k,B k,B
Consider the CEMD model projection problem. If our solution belongs to M , we call it
k,B
a single-criterion solution. Otherwise if our solution does not belong to M , i.e., there exists
k,B
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 OurContributions andTechniques
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
summarizeourcontributions andmaintechniques inthefollowing.
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
withrunning timeO˜(nlogn). Bothalgorithms achieve constant approximation ratio andoutput a
treeofsizeatmost2k. Inthispaper, weprovidethefirstlineartimealgorithmsforbothhead-and
tail-approximation tree-sparsity projection problems and remove the bicriterion relaxation. This
provides anaffirmativeanswer tothe open problem inHegdeand Indyk [18], whichasks whether
thereisanearly-linear timesingle-criterion approximation algorithm fortreesparsity.
MainTechniquesforTreeSparsityModel: Thebottleneck ofprevious algorithms iscomputing
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 (α,β)-RS (min,+)-convolution. Instead of maintaining the whole (min,+)-
convolution array,weonlycomputeasparsesequence toapproximately represent thewholearray.
Taking Tree-Sparsity-Tail as an example, we only need to maintain O˜(logn) elements in a sin-
gle 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 O˜(1),
instead of O(k) for the exact (min,+)-convolution. Thus, we only cost O˜(logn) to compute
our approximate (min,+)-convolution. For Tree-Sparsity-Head, we apply a similar approxi-
mate (max,+)-convolution technique, called (α,β)-RS (max,+)-convolution. Our approximate
convolution technique mayhaveindependent interest.
4
InSection3,wecombinetheapproximate(min,+)-convolutiontechniqueandotherapproaches
suchasweightdiscretization,pruningandthelookuptablemethod. Ourresultscanbesummarized
bythefollowingtheorem.
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. Specifi-
cally, foranyconstant ǫ ∈ (0,1), thereisanO(ǫ 1n)timeapproximation algorithm thatreturns
1 −1
asupport Ωˆ ∈T satisfying
k
kx k ≥ (1−ǫ )maxkx k .
Ωˆ p 1 Ω∈Tk Ω p
Forany constant ǫ ∈ (0,∞), there isan O(n+ǫ 2n/logn)timeapproximation algorithm that
2 −2
returnsasupport Ωˆ ∈T satisfying
k
kx−x k ≤ (1+ǫ ) min kx−x k ,
Ωˆ p 2 Ω∈Tk Ω p
if k ≤ n1 δ (δ ∈ (0,1) is any fixed constant), and there is an O(ǫ 1n(logloglogn)2) time
− −2
algorithm forgeneralk.
Thencombiningwithpriorresults[2,16,17],weprovideamoreefficientrobustsparserecovery
algorithmintreesparsitymodelasfollows. Thebestpriorresultcanrecoveranapproximatesignal
xˆ ∈ T for some constant c > 1 [16] (i.e., the sparsity of their solution is ck). In this paper, we
ck
improvetheconstant cto1.
Corollary 10. Assume that k ≤ n1 δ (δ ∈ (0,1) is any fixed constant). Let A ∈ Rm n be a
− ×
measurement matrix. Letx ∈ T be an arbitrary signal in the tree sparsity model withdimension
k
n, and let y = Ax + e ∈ Rm be a noisy measurement vector. Here e ∈ Rm is a noise vector.
Then there exists an algorithm to recover a signal approximation xˆ ∈ T satisfying kx − xˆk ≤
k
Ckek for some constant C from m = O(k) measurements. Moreover, the algorithm runs in
2
O((nlogn+k2lognlog2(klogn))log kxk2)time.
e 2
k k
CEMD Model: In Section 4, we consider the CEMD model M and propose the first single-
k,B
criterion constantfactorapproximation algorithm forthehead-approximation oracle.
Theorem 11. Consider the CEMD model M with s = k/w sparse for each column and
k,B
support-EMD B. Let δ ∈ (0,1/4), x = min |X |p, and x = max|X |p. Let
c = 1/4−δ. Thereexistsanalgorithmmriunnning in|OX(i,sj|h>n0logi,nj +log xmmaxa)xtime,whichir,jeturns a
δ xmin
single-criterion c1/p approximation forthehead-approximation projection problem.
Combining withAM-IHT framework [17], weobtain the following corollary which improves
the prior result [17] in two aspects: 1) We decrease the total number of measurements from
m = O(klog(B log k)) to m = O(klog(B/k)). 2) We decrease the running time of the robust
k w
sparserecoveryfromO(nlog kxk2(klogn+kh(B+logn+log xmax)))toO(nlog kxk2(klogn+
kh(logn+log xmax))). kek2 w xmin kek2
w xmin
5
Corollary 12. LetA ∈ Rm n beameasurement matrix. Letx ∈ M beanarbitrary signal in
× k,B
the CEMD model with dimension n = wh, and let y = Ax+e ∈ Rm be a noisy measurement
vector. Heree ∈ Rm is anoise vector. Then there exists analgorithm to recover asignal approx-
imation xˆ ∈ M satisfying kx−xˆk ≤ Ckek forsomeconstant C from m = O(klog(B/k))
k,2B 2
measurements. Moreover, the algorithm runs in O(nlog kxk2(klogn + kh(logn + log xmax)))
time,wherex = max|x |andx = min |x |. kek2 w xmin
max i min |xi|>0 i
1.3 Related Work
In the tree sparsity model, there is an algorithm with running time O(nklogn) 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 com-
puting (min,+)-convolutions. The naive algorithm for computing the (min,+)-convolution of
arbitrary two length-n arrays requires O(n2) time. Williams proposed an improvement algorithm
forcomputing(min,+)-convolutions, andreducedtherunning timetoO(n2/2Ω(√logn))[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
timeoftreesparserecoverytoO˜(nlogn).
Schmidtetal.[24]introducedtheConstrainedEarthMover’sDistance(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 mentionedinthesurvey[18].
Other structured sparsity models also have been studied by researchers. Huang etal. [20]first
considered the graph sparsity model, and provided a head-approximation algorithm with a time
complexity of O(nc), 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 △-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(logn)O(1). Our algorithms only work for b-ary
trees,butourrunning timesaremuchbetter.
2 Approximate (min,+)-Convolution
In this section, we introduce an approach of computing an approximate (min,+)-convolution,
whichisusefulforthetreesparsitymodel.
6
2.1 (α,β)-RS(min,+)-Convolution
Wefirstintroduce aconceptcalled(min,+)-convolution.
Definition 13 ((min,+)-convolution (see e.g. [5,7])). Given two arrays A = (a[0],a[1],a[2],
...,a[m ]) and B = (b[0],b[1],b[2],...,b[m ]), their (min,+)-convolution is the array S =
1 2
(s[0],s[1],s[2],...,s[m +m ])wheres[t]= mint {a[i]+b[t−i]},t ∈ [0,m +m ].
1 2 i=0 1 2
We sketch how to use (min,+)-convolutions in the tree sparsity model. Recall that T is
ij
the subtree rooted at N . We maintain an array S = (s[0],s[1],...,s[|T |]) for each node
ij ij ij
N . The element s[l] represents the optimal tail value for Tree-Sparsity-Tail on T , i.e., s[l] =
ij ij
minΩ∈T|Tij|−l(Tij) Ni′j′∈Tij\Ωxi′j′ where Tk(Tij) is the tree sparsity model defined at the tree
T (see DefinitionP5). In fact, the array S can be achieved through computing the (min,+)-
ij ij
convolution from the arrays of its two children. 1 Finally, we output the element s[n−k] in the
arrayoftherootnodeN ,whichistheoptimaltailvalueofTree-Sparsity-TailonT.
log(n+1),1
However,therunning timeforcomputing theexact(min,+)-convolution istoolong. Instead,
we compute an approximate (min,+)-convolution for each node. We first introduce some con-
cepts.
Definition 14 (α-RS). Given a sequence Aˆ = (aˆ[i ],aˆ[i ],...,aˆ[i ]) and a fixed constant α ∈
1 2 m
[0,∞). Each element aˆ[i ] ∈ Aˆ is a real number with an associated index i . If for any v ∈
v v
[1,m −1], i > i ,aˆ[i ] ≥ (1+α)aˆ[i ] ≥ 0, we call the sequence Aˆ an α-representative
v+1 v v+1 v
sequence (α-RS).
Definition 15 (Completion of α-RS). Consider an α-RS Aˆ = (aˆ[i ],aˆ[i ]...,aˆ[i ]). Define its
1 2 m
completion by an array A = (a[0],a [1],...,a[i ]) satisfying that: 1) If 0 ≤ t ≤ i ,a[t] =
′ ′ ′ ′ m 1 ′
a[i ];2)Ifi +1≤ t ≤ i (1 ≤ v ≤ m−1),a [t]= aˆ[i ].
′ 1 v v+1 ′ v+1
Bythefollowingdefinition, weshowhowtouseanα-RStoapproximately represent anarray.
Definition16(SequenceApproximation). Giventwonon-decreasingarraysA = (a[0],a [1],...,
′ ′ ′
a[n]) and A = (a[0],a[1],...,a[n]), we say A is an α-approximation of A if for any i, a[i] ≤
′ ′
a[i] ≤ (1 + α)a[i]. We say an α-RS Aˆ approximates an array A if its completion A is an α-
′ ′
approximation ofA.
A special case is that we say A is a 0-approximation of A if A = A. Figure 1 illustrates
′ ′
theseconcepts. Now,wearereadytointroducetheformaldefinitionoftheapproximate(min,+)-
convolution, called(α,β)-RS(min,+)-convolution.
Definition 17 ((α,β)-RS (min,+)-convolution). Given two α-RSs Aˆand Bˆ, suppose A and B
′ ′
are their completions respectively. Suppose the array S is the (min,+)-convolution of A and
′
B . We call a sequence Sˆ an (α,β)-RS (min,+)-convolution of Aˆ and Bˆ if Sˆ is a β-RS which
′
approximates thearrayS.
1Notethat thevalue of s[|T |] can not be directlyobtained by the (min,+)-convolution of the arrays of itstwo
ij
children.Infact,s[|Tij|]=PNi′j′∈Tijxi′j′.
7
aˆ(8)
value
a(4) a(5) a(6) a(7) a(8)
′ ′ ′ ′ ′
a(8)
a(6) a(7)
a(5)
a(4)
aˆ(3)
a(2) a(3)
′ ′
aˆ(1)
aˆ(0) a(1) a(3)
′
a(0) a(2)
′
a(1)
a(0) n = 8 k
Figure 1: The figure illustrates the concepts α-RS and its completion. Here, Aˆ =
(aˆ[0],aˆ[1],aˆ[3],aˆ[8]) is an α-RS. The array (a[0] = aˆ[0],a[1],...,a[8]) is the completion of
′ ′ ′
Aˆ. Bythisfigure,wecanseethattheα-RSAˆapproximates thearrayA = (a[0],a[1],...,a[8]).
By preserving an (α,β)-RS (min,+)-convolution instead of an exact (min,+)-convolution,
wecanreducethestorage spaceandthecomputation time.
2.2 A fastalgorithmfor(α,β)-RS(min,+)-convolution
Next, we give a simple algorithm RSMinPlusto compute an (α,β)-RS (min,+)-convolution of
twogivenα-RSsAˆandBˆ.
RSMinPlus(α,β,Aˆ,Bˆ): We first compute the sum of every pair (aˆ[i ],ˆb[j ]) where aˆ[i ] ∈
v t v
Aˆ,ˆb[j ] ∈ Bˆ. Define s˜[l ] = min aˆ[i ] + ˆb[j ], where Φ = {(i ,j ) | i + j =
t r (iv,jt)∈Φlr v t lr v t v t
l ,aˆ[i ] ∈ Aˆ,ˆb[j ] ∈ Bˆ}. Suppose that there are m different elements s˜[l ]. Then we sort s˜[l ] in
r v t r r
the increasing order ofthe index number l . Aftersorting, weobtain amonotone increasing array
r
S˜ = (s˜[l ],s˜[l ],...,s˜[l ]),l < l forr ∈ [m−1]. Finally, weconstruct aβ-RSSˆfrom S˜as
1 2 m r r+1
oursolution. Ourconstruction isasfollows.
1. Initially appendsˆ[l ]= s˜[l ]toSˆ. Letθ = sˆ[l ]/(1+β).
m m m
2. SequentiallyconsiderallelementsinS˜indecreasingorderoftheindexnumber. Ifs˜[l ]≤ θ,
r
append sˆ[l ] = s˜[l ] to Sˆ. Let θ = sˆ[l ]/(1+β). Otherwise, ignore s˜[l ] and consider the
r r r r
nextelements˜[l ]∈ S˜.
r 1
−
3. Returnthefinalsequence Sˆ.
Lemma 18. Suppose Aˆ = (aˆ[i ],aˆ[i ],...,aˆ[i ]) and Bˆ = (ˆb[j ],ˆb[j ],...,ˆb[j ]). Let m =
1 2 m1 1 2 m2
max{m ,m }. RSMinPlus(α,β,Aˆ,Bˆ)computesan(α,β)-RS(min,+)-convolution SˆofAˆand
1 2
Bˆ inO(m2logm)time.
8
Proof. Thecorrectness isnothard. LetA andB bethecompletions ofAˆandBˆ respectively. Let
′ ′
S bethe exact (min,+)-convolution oftwoarrays A and B . Bytheconstruction ofS˜, wehave
′ ′ ′
thatS isthecompletion ofS˜. Moreover, theβ-RSSˆapproximates thearrayS ,whichprovesthe
′ ′
correctness. Itremainstoprovetherunningtime.
BythealgorithmRSMinPlus,ittakesm ·m =O(m2)timetocomputealls˜[l ]. Thus,there
1 2 r
areatmost m2 different elements s˜[l ]inS˜. Thenweneed O(m2logm)timetosort alls˜[l ]and
r r
obtainthearrayS˜. Finally,scanning allelementsinS˜toconstruct theβ-RSSˆneedsO(m2)time.
Overall,theruntimeofthealgorithm isO(m2logm).
FastRSMinPlus. Amorecareful(α,β)-RS (min,+)-convolution algorithm: Inthetreespar-
sity model, we have some additional conditions which can improve the running time of the algo-
rithm RSMinPlus(α,β,Aˆ,Bˆ). Weconsider the case that Aˆare nonnegative sequences witheach
element aˆ ∈ Aˆsatisfying that either aˆ = 0oraˆ ≥ 1. Wehave thesame assumption on Bˆ. More-
over, α and β are two constants such that 0 < β ≤ α ≤ 1. The intuition is as follows. Suppose
we have just appended some element sˆ[l] to the array Sˆ. By the property of (α,β)-RS (min,+)-
convolution, wecansafelyignorealls˜[l ]withsˆ[l]/(1+β) < s˜[l ] ≤ sˆ[l],andfocusonfindingthe
r r
largests˜[l ]∈ S˜suchthats˜[l ] ≤ sˆ[l]/(1+β)(seethedefinitionofS˜inRSMinPlus(α,β,Aˆ,Bˆ)).
r r
Wefirstbuildahashtable HashB. Theconstruction isasfollows. EachelementinHashBisa
pair (key,value) where the key term is an integer satisfying that −1 ≤ key ≤ log ˆb[j ] ,
1+β m2
l m
and the value term is the largest index w satisfying that ˆb[w] ∈ Bˆ and ˆb[w] ≤ (1 + β)key. 2
Symmetrically, we also build a hash table HashA for the array Aˆ. Let U = max{aˆ[i ],ˆb[j ]}.
m1 m2
Since both Aˆand Bˆ are increasing sequences, we can construct hash tables HashA and HashB in
O(log U)timebyconsidering allkey termsinincreasing order.
1+β
Given an element sˆ[l] ∈ Sˆ, we show how to find the largest s˜[l ] ∈ S˜ such that s˜[l ] ≤
r r
sˆ[l]/(1+β) by hash tables. We first reduce this problem to finding the element ˆb[j ] ∈ Bˆ of the
t
largest index j for each aˆ[i ] ∈ Aˆ, such that aˆ[i ]+ˆb[j ] ≤ sˆ[l]/(1 + β). A simple scheme is
t v v t
to enumerate all aˆ[i ] ∈ Aˆ, query the hash table HashB, and find the largest ˆb[j ] ∈ Bˆ such that
v t
ˆb[j ] ≤ sˆ[l]/(1+β)−aˆ[i ]. Thenamongallsuch(i ,j )indexpairs,wechoosethepair(i ,j )
t v v t v∗ t∗
withthelargestsumi +j . Weappendsˆ[i +j ]= s˜[i +j ] = aˆ[i ]+ˆb[j ]toSˆ. However,
v∗ t∗ v∗ t∗ v∗ t∗ v∗ t∗
enumerating allelementsisnotnecessary, sincewehavethefollowinglemma.
Lemma19. Letτ = ⌈1/α⌉. Foranyaˆ[i ],aˆ[i ] ∈ Aˆ,wehaveaˆ[i ] ≤ aˆ[i ]/2. Similarly, for
v τ v v τ v
anyˆb[j ],ˆb[j ] ∈ Bˆ,ˆb[j ]≤ˆb[j ]/2. − −
t τ t t τ t
− −
Proof. W.l.o.g.,weonlyconsider thearrayAˆ. ByDefinition17,weknow(1+α)aˆ[i ] ≤aˆ[i ].
v 1 v
−
If α ≥ 1, τ = 1, the lemma is trivially true. Otherwise if α < 1, we have aˆ[i ] ≤ aˆ[i ]/(1+
v τ v
α)1/α ≤ aˆ[i ]/2. −
v
Let θ = sˆ[l]/(1 + β). Assume that aˆ[i ] ∈ Aˆ (resp. ˆb[j ] ∈ Bˆ) is the largest element such
v t
that aˆ[i ] ≤ θ (resp. ˆb[j ] ≤ θ). Hence, both aˆ[i ]andˆb[j ]are atleast θ. Therefore, the pair
v t v+1 t+1
2Ifˆb[j ]=0,wedefineHashB(−1)=ˆb[j ].Otherwiseifˆb[j ]≥1,forakeyterm,ifwehaveˆb[j ]>(1+β)key,
1 1 1 1
weignorethiskeytermwhenconstructingHashB.
9
(i ,j ) must satisfy that either v −τ ≤ v ≤ v or t−τ ≤ t ≤ t, since aˆ[i ]+ˆb[j ] ≤ θ
v∗ t∗ ∗ ∗ v τ t τ
byLemma19. Thus,weonlyneedtoconsideratmostτ +1elementsinAˆorBˆ−. Notethat−wecan
directlyfindsuchindexi (resp. j )bythefollowinglemma.
v t
Lemma 20. Either i = HashA( log θ ) or i = HashA( log θ ). Similarly, either j =
v 1+β v 1+β t
HashB( log θ orj = HashB(cid:6)( log (cid:7)θ (cid:4) (cid:5)
1+β t 1+β
(cid:6) (cid:7) (cid:4) (cid:5)
Proof. W.l.o.g., we take i as example. Let i = HashA( log θ ). If i 6= i , then by the
v w 1+β v w
definition of iv, wehave that θ < aˆ[iw] ≤ (1+β)⌈log1+βθ⌉(cid:6). Since Aˆ(cid:7)is an α-RS and α ≥ β, we
have
aˆ[iw 1] ≤ aˆ[iw]/(1+α) ≤ aˆ[iw]/(1+β) ≤ (1+β)⌈log1+βθ⌉−1 ≤ (1+β)⌊log1+βθ⌋ ≤ θ
−
Thus,wehavei = i bythedefinition ofi . Notethati = HashA( log θ ). Wefinish
v w 1 v w 1 1+β
− −
theproof. (cid:4) (cid:5)
Wesummarizeourapproaches inAlgorithm 1. Theanalysis ofthealgorithm isasfollows.
Lemma 21. Let U = max{aˆ[i ],ˆb[j ]}. Let 0 ≤ β ≤ α ≤ 1 be two constants. The algo-
m1 m2
rithm FastRSMinPlus(α,β,Aˆ,Bˆ) computes an (α,β)-RS (min,+)-convolution Sˆ of Aˆand Bˆ in
log U
O 1+β time.
α
(cid:16) (cid:17)
Proof. Wefirstprovethecorrectness. LetAandBbethecompletionsofAˆandBˆ respectively. Let
Sbethe(min,+)-convolutionofAandB. Assumethatwejustappendanelementsˆ[l]toSˆandlet
θ = sˆ[l]/(1+β). ConsiderthenextrecursionfromLine5toLine20,weappendanewelementsˆ[l]
′
toSˆ. InitiallyinLine6,wefindthelargestindexj satisfying thatˆb[j ]+aˆ[i ] ≤ θbyLemma20.
t t 1
Wefirstanalyse thefirstloopfromLine8toLine12. InLine10,wefindthelargest elementaˆ[w]
satisfyingthataˆ[w] ≤ θ−ˆb[j ]byLemma20. Thus,weconcludethataˆ[w]+ˆb[j ] ≤ θforany
t δ t ∆
1 ≤ ∆ ≤ min{t,τ}. Similarl−y,wecanprovethatˆb[w]+aˆ[i ] ≤ θforany1 ≤−∆ ≤ min{v,τ}
v ∆
−
inLine17. Thus,theelementsˆ[l]alwayssatisfiesthatsˆ[l] ≤ θ = sˆ[l]/(1+β)duringtherecursion.
′
Thus,theoutputSˆisaβ-RS.
On the other hand, let S′ be the completion of Sˆ. Assume that s[l] = aˆ[iv∗] +ˆb[jt∗] is the
largest element in S satisfying that s[l] ≤ θ. W.l.o.g. we assume that θ/2 ≤ ˆb[j ] ≤ θ −aˆ[i ]
t∗ 1
(otherwise θ/2 ≤ aˆ[i ] ≤ θ−ˆb[j ]). Weconclude that t−τ ≤ t ≤ t by Lemma 19. Then we
v∗ 1 ∗
mustconsidertheelementˆb[j ]intheloopfromLine8toLine12. NotethatinLine10,wefindan
t∗
indexw = i byLemma20. BytheupdatingrulesinLine11-12,weupdatesˆ[l] = aˆ[i ]+ˆb[j ]
v∗ v∗ t∗
inLine12andappendsˆ[l]toSˆinLine20. Consideringanyelements[l ]withl+1 ≤ l ≤ l ,we
′ 0 0 ′
havethats[l ] = sˆ[l ]byDefinition15. Moreover,wehavethefollowinginequality bythechosen
′ 0 ′
ofl ,
r
s[l ]≤ s[l ]= sˆ[l] = s[l]= (1+β)θ < (1+β)s[l+1] ≤(1+β)s[l ].
0 ′ 0 ′ ′ 0
Overall,weprovethatS isaβ-approximation ofS byDefinition16.
′
10