Least Squares and Shrinkage Estimation under Bimonotonicity Constraints

0.42 MB·English

by  Rudolf Beran

Least Squares and Shrinkage Estimation under
Bimonotonicity Constraints
RudolfBeranandLutzDu¨mbgen
UniversityofCaliforniaatDavisandUniversityofBern
September2008,revisedJanuary2009
9
0
0
2 Abstract

n
Inthispaperwedescribeactivesettypealgorithmsforminimizationofasmoothfunction
a
J undergeneralorderconstraints, animportantcasebeingfunctionsonthesetofbimonotone
6 r × s matrices. These algorithms can be used, for instance, to estimate a bimonotone re-
2 gressionfunctionvialeastsquaresor(asmoothapproximationof)leastabsolutedeviations.
Anotherapplicationisshrinkageestimationinimagedenoisingor,moregenerally,regression

]
O problemswithtwoordinalfactorsafterrepresentingthedatainasuitablebasiswhichisin-
dexed by pairs (i,j) ∈ {1,...,r}×{1,...,s}. Various numerical examples illustrate our
C
methods.
.
t
a
t
s
[ Keywords: activesetalgorithm,dynamicprogramming,estimatedrisk,pool-adjacent-violators

3 algorithm,regularization.
v
4
7
9 AMSsubjectclassifications: 62-04,62G05,62G08,90C20,90C25,90C90
0
.
9
0 1 Introduction
8
0
:
v
Monotonicity and other qualitative constraints play an important role in contemporary nonpara-
i
X
metric statistics. One reason for this success is that such constraints are often plausible or even
r
a
justified theoretically, within an appropriate mathematical formulation of the application. More-
over,byimposingshapeconstraintsonecanoftenavoidmoretraditionalsmoothnessassumptions
whichtypicallyleadtoproceduresrequiringthechoiceofsometuningparameter. Agoodstarting
point for statistical inference under qualitative constraints is the monograph by Robertson et al.
[10].
Estimation under order constraints leads often to the following optimization problem: For
1
somedimensionp ≥ 2letQ : Rp → Rbeagivenfunctional. Forinstance,
p
(cid:88)
(1) Q(θ) = w (Z −θ )2
u u u
u=1
withacertainweightvectorw ∈ (0,∞)p andagivendatavectorZ ∈ Rp. Ingeneralweassume
thatQiscontinuouslydifferentiable,strictlyconvexandcoercive,i.e.
Q(θ) → ∞ as(cid:107)θ(cid:107) → ∞,
where (cid:107)·(cid:107) is some norm on Rp. The goal is to minimize Q over the following subset K of Rp:
LetC beagivencollectionofpairs(u,v)ofdifferentindicesu,v ∈ {1,2,...,p},anddefine
K = K(C) = (cid:8)θ ∈ Rp : θ ≤ θ forall(u,v) ∈ C(cid:9).
u v
ThisdefinesaclosedconvexconeinRp containingallconstantvectors.
For instance, if C consists of (1,2), (2,3), ..., (p − 1,p), then K is the cone of all vectors
θ ∈ Rp suchthatθ ≤ θ ≤ ··· ≤ θ . Minimizing(1)overallsuchvectorsisastandardproblem
1 2 p
andcanbesolvedinO(p)stepsviathepool-adjacent-violatorsalgorithm(PAVA).Thelatterwas
introducedinaspecialsettingbyAyeretal.[1]andextendedlaterbynumerousauthors,see[10]
andBestandChakravarti[3].
AssoonasQ(·)isnotoftype(1)orC differsfromtheaforementionedstandardexample,the
minimizationofQ(·)overKbecomesmoreinvolved. HereisanotherexampleforKandC which
isofprimaryinterestinthepresentpaper: Letp = rswithintegersr,s ≥ 2,andidentifyRp with
thesetRr×s ofallmatriceswithr rowsandscolumns. FurtherletK bethesetofallmatrices
r,s
θ ∈ Rr×s suchthat
θ ≤ θ wheneveri < r and θ ≤ θ wheneverj < s.
i,j i+1,j i,j i,j+1
(cid:0) (cid:1)
This corresponds to the set C of all pairs (i,j),(k,(cid:96)) with i,k ∈ {1,...,r} and j,(cid:96) ∈
r,s
{1,...,s} such that either (k,(cid:96)) = (i + 1,j) or (k,(cid:96)) = (i,j + 1). Hence there are #C =
2rs−r−sconstraints.
Minimizingthespecialfunctional(1),i.e.Q(θ) = (cid:80) w (Z −θ )2,overthebimonotone
i,j ij ij ij
coneK isawellrecognizedproblemwithvariousproposedsolutions,see,forinstance,Spouge
r,s
et al. [11], Burdakow et al. [4], and the references cited therein. However, all these algorithms
exploit the special structure of K or (1). For general functionals Q(·), e.g. quadratic functions
r,s
withpositivedefinitebutnon-diagonalhessianmatrix,differentapproachesareneeded.
2
Theremainderofthispaperisorganizedasfollows. InSection2wedescribethebimonotone
regression problem and argue that the special structure (1) is sometimes too restrictive even in
that context. In Section 3 we derive possible algorithms for the general optimization problem
described above. These algorithms involve a discrete optimization step which gives rise to a
dynamic program in case of K = K . For a general introduction to dynamic programming see
r,s
Cormenetal.[6]. Otheringredientsareactivemethodsasdescribedby,forinstance,Fletcher[9],
Best and Chakravarti [3] or Du¨mbgen et al. [8], sometimes combined with the ordinary PAVA in
a particular fashion. It will be shown that all these algorithms find the exact solution in finitely
many steps, at least when Q(·) is an arbitrary quadratic and strictly convex function. Finally, in
Section 4 we adapt our procedure to image denoising via bimonotone shrinkage of generalized
Fourier coefficients. The statistical method in this section was already indicated in Beran and
Du¨mbgen[2]buthasnotbeenimplementedyet,forlackofanefficientcomputationalalgorithm.
2 Least squares estimation of bimonotone regression functions
Supposethatoneobserves(x1,y1,Z1),(x2,y2,Z2),...,(xn,yn,Zn)withrealcomponentsxt,yt
andZt. Thepoints(xt,yt)areregardedasfixedpoints,whichisalwayspossiblebyconditioning,
while
Zt = µ(xt,yt)+εt
foranunknownregressionfunctionµ : R×R → Randindependentrandomerrorsε1,ε2,...,εn
withmeanzero. Insomeapplicationsitisplausibletoassumeµtobebimonotoneincreasing,i.e.
non-decreasinginbotharguments. Thenitwouldbedesirabletoestimateµunderthatconstraint
only. Onepossibilitywouldbetominimize
n
(cid:88)
(Zt−µ(xt,yt))2
t=1
overallbimonotonefunctionsµ. Theresultingminimizerµˆisuniquelydefinedonthefinitesetof
alldesignpoints(xt,yt),1 ≤ t ≤ n.
Foramoredetaileddiscussion,supposethatwewanttoestimateµonafiniterectangulargrid
(cid:8) (cid:9)
(x ,y ) : 1 ≤ i ≤ r,1 ≤ j ≤ s ,
(i) (j)
where x < x < ··· < x and y < y < ··· < y contain at least the different
(1) (2) (r) (1) (2) (s)
elements of {x1,x2,...,xn} and {y1,y2,...,yn}, respectively, but maybe additional points as
3
well. For1 ≤ i ≤ rand1 ≤ j ≤ sletw bethenumberofallt ∈ {1,...,n}suchthat(xt,yt) =
ij
(x ,y ), and let Z be the average of Zt over these indices t. Then (cid:80)n (Zt −µ(xt,yt))2
(i) (j) ij t=1
equals
(cid:88)
Q(θ) = w (Z −θ )2,
ij ij ij
i,j
whereθ = (θ ) standsforthematrix(cid:0)µ(x ,y )(cid:1) ∈ K .
ij i,j (i) (j) i,j r,s
Setting 1: Complete layout. Suppose that w > 0 for all (i,j) ∈ {1,...,r} × {1,...,s}.
ij
Thentheresultingoptimizationproblemispreciselytheonedescribedintheintroduction.
Setting 2a: Incomplete layout and simple interpolation/extrapolation. Suppose that the set
U ofallindexpairs(i,j)withw > 0differsfrom{1,...,r}×{1,...,s}. Then
ij
(cid:88)
Q(θ) = w (Z −θ )2
u u u
u∈U
fails to be coercive. Nevertheless it can be minimized over K with the algorithms described
r,s
later. Letθˇ besuchaminimizer. SinceitisuniquelydefinedonU only, weproposetoreplaceit
withθˆ = 2−1(θ+θ),where
(cid:16) (cid:17)
θ = max (cid:8)θˇ : (i(cid:48),j(cid:48)) ∈ U,i(cid:48) ≤ i,j(cid:48) ≤ j(cid:9)∪{θˇ } ,
ij i(cid:48)j(cid:48) min
(cid:16) (cid:17)
θ = min (cid:8)θˇ : (i(cid:48),j(cid:48)) ∈ U,i ≤ i(cid:48),j ≤ j(cid:48)(cid:9)∪{θˇ } ,
ij i(cid:48)j(cid:48) max
and θˇ and θˇ denote the minimum and maximum, respectively, of {θˇ : u ∈ U}. Note that
min max u
θ and θ belong to K and are extremal in the sense that any matrix θ ∈ K ∩[θˇ ,θˇ ]r×s
r,s r,s min max
withθ = θˇ forallu ∈ U satisfiesnecessarilyθ ≤ θ ≤ θ forall(i,j).
u u ij ij ij
Setting 2b: Incomplete layout and light regularization. Instead of restricting one’s attention
to the index set U, one can estimate the full matrix (cid:0)µ(x ,y )(cid:1) ∈ Rr×s by minimizing a
(i) (j) i,j
suitablypenalizedsumofsquares,
(cid:88)
Q(θ) = w (Z −θ )2+λP(θ),
u u u
u∈U
over K for some small parameter λ > 0. Here P(·) is a convex quadratic function on Rr×s
r,s
suchthatQ(·)isstrictlyconvex. OnepossibilitywouldbeTychonovregularisationwithP(θ) =
4
(cid:80) (θ −θ )2 andacertainreferencevalueθ ,forinstance,θ = (cid:80) w Z (cid:14)(cid:80) w . Inour
i,j ij o o o i,j ij ij i,j ij
particularsettingwepreferthepenalty
(cid:88)
(2) P(θ) = (θ −θ )2,
k(cid:96) ij
((i,j),(k,(cid:96)))∈Cr,s
because it yields smoother interpolations than the recipe for Setting 2a or the Tychonov penalty.
OnecaneasilyshowthattheresultingquadraticfunctionQisstrictlyconvexbutwithnon-diagonal
hessianmatrix. Thusitfulfillsourgeneralrequirementsbutisnotoftype(1).
Note that adding a penalty term such as (2) could be worthwhile even in case of a complete
layout if the underlying function µ is assumed to be smooth. But this leads to the nontrivial task
of choosing λ > 0 appropriately. Here we use the penalty term mainly for smooth interpola-
tion/extrapolationwithλjustlargeenoughtoensureawell-conditionedHessianmatrix. Werefer
tothisas“lightregularization”,andtheexactvalueofλisessentiallyirrelevant.
Example2.1 To illustrate the difference between simple interpolation/extrapolation and light
regularization with penalty (2) we consider just two observations, (x1,y1,Z1) = (2,3,0) and
(x2,y2,Z2) = (6,7,1), and let r = 7, s = 10 with x = i and y = j. Thus w = 0 except
(i) (j) ij
for w = w = 1, while Z = 0 and Z = 1. Any minimizer θˇ of (cid:80) w (Z −θ )2
2,3 6,7 2,3 6,7 u∈U u u u
overK satisfiesθˇ = 0andθˇ = 1,sotherecipeforSetting2ayields
7,10 2,3 6,7

0 ifi ≤ 2,j ≤ 3,


θˆ = 1 ifi ≥ 6,j ≥ 7,
ij

0.5 else.
The left panel of Figure 1 shows the latter fit θˆ, while the right panel shows the regularized fit
based on (2) with λ = 10−4. In these and most subsequent pictures we use a gray scale from
black = 0towhite = 1.
Example2.2 (Binary regression). We generated a random matrix Z ∈ {0,1}r×s with r = 70
rows,s = 100columnsandindependentcomponentsZ ,where
ij
(cid:8) (cid:9)
x +y 1 y ≥ 1/2+cos(πx )/4
(i) (j) (j) (i)
P(Z = 1) = θ = +
ij ij
4 2
with x = (i−0.5)/r and y = (j −0.5)/s. Thereafter we removed randomly all but 700
(i) (j)
of the 7000 components Z . The resulting data are depicted in the upper left panel of Figure 2,
ij
where missing values are depicted grey, while the upper right panel shows the true signal θ. The
5
Figure1: Simpleinterpolation/extrapolationversuslightregularization
lower panels depict the least squares estimator with simple interpolation/extrapolation (left) and
lightregularizationbasedon(2)withλ = 10−4 (right). Notethatbothestimatorsareverysimilar.
Duetothesmallvalueofλ,themaindifferencesoccurinregionswithoutdatapoints.
Thequalityofanestimatorθˆ forθ maybequantifiedbytheaverageabsolutedeviation,
r s
1 (cid:88)(cid:88)
AAD = |θˆ −θ |.
ij ij
rs
i=1 j=1
Fortheestimatorwithsimpleinterpolation/extrapolation,AADturnedouttobe7.5607·10−2,the
estimatorbasedonlightregularizationperformedslightlybetterwithAAD = 7.4039·10−2.
3 The general algorithmic problem
Wereturntothegeneralframeworkintroducedinthebeginningwithacontinuouslydifferentiable,
strictly convex and coercive functional Q : Rp → R and a closed convex cone K = K(C) ∈ Rp
determinedbyacollectionC ofinequalityconstraints.
Beforestartingwithexplicitalgorithms,letuscharacterizethepoint
θˆ = argmin Q(θ).
θ∈K
Itiswell-knownfromconvexanalysisthatapointθ ∈ Kcoincideswithθ if,andonlyif,
(3) ∇Q(θ)(cid:62)θ = 0 ≤ ∇Q(θ)(cid:62)η forallη ∈ K,
where ∇Q(θ) denotes the gradient of Q at θ. This characterization involves infinitely many
inequalities,butitcanbereplacedwithacriterioninvolvingonlyfinitelymanyconstraints.
6
Figure2: Binaryregressionwithincompletelayout
3.1 ExtremaldirectionsofK
NotethatKcontainsallconstantvectorsc1,c ∈ R,where1 = 1 = (1)p . Itcanberepresented
p i=1
asfollows:
Lemma3.1 Define
E = K∩{0,1}p.
Thenanyvectorx ∈ Kmayberepresentedas
(cid:88)
x = min(x)1+ λ e
e
e∈E
(cid:80)
withcoefficientsλ ≥ 0suchthat λ = max(x)−min(x).
e e∈E e
Here min(x) and max(x) denote the minimum and maximum, respectively, of the compo-
nentsofx.
7
Modifiedcharacterizationofθˆ. BymeansofLemma3.1onecaneasilyverifythat(3)isequiv-
alenttothefollowingcondition:
(4) ∇Q(θ)(cid:62)θ = 0 ≤ ∇Q(θ)(cid:62)e foralle ∈ E ∪{−1}.
Thus we have to check only finitely many constraints. Note, however, that the cardinality of E
maybesubstantiallylargerthanthedimensionp,sothatchecking(4)isfarfromtrivial.
ApplicationtoK . ApplyingLemma3.1totheconeK ⊂ Rr×s yieldsthefollowingrepre-
r,s r,s
sentation: With
E = K ∩{0,1}r×s
r,s r,s
anymatrixx ∈ Kmaybewrittenas
(cid:88)
x = a 1 + λ e
o r×s e
e∈Er,s
withcoefficientsa ∈ Randλ ≥ 0,e ∈ E .
o e r,s
There is a one-to-one correspondence between the set E and the set of all vectors e˜ ∈
r,s
{1,2,...,r+s}r withcomponentse˜ < e˜ < ··· < e˜ viathemapping
1 2 r
s
(cid:16) (cid:88) (cid:17)r
e (cid:55)→ i+ e .
ij
i=1
j=1
Sincesuchavectore˜correspondstoasubsetof{1,2,...,r+s}withrelements,weendupwith
(cid:18) (cid:19) (cid:18) (cid:19)
r+s r+s
#E = = .
r,s
r s
HencethecardinalityofE growsexponentiallyinmin(r,s). Nevertheless,minimizingalinear
r,s
functionaloverE ispossibleinO(rs)steps,asexplainedinthenextsection.
r,s
Proof of Lemma 3.1. For x ∈ K let a < a < ··· < a be the different elements of
0 1 m
{x ,x ,...,x },i.e.a = min(x)anda = max(x). Then
1 2 p 0 m
m
(cid:88) (cid:0) (cid:1)p
x = a 1+ (a −a ) 1{x ≥ a } .
0 i i−1 t i t=1
i=1
Obviously,theseweightsa −a arenonnegativeandsumtomax(x)−min(x). Furthermore,
i i−1
onecaneasilydeducefromx ∈ Kthat(cid:0)1{x ≥ a}(cid:1)p belongstoE foranyrealthresholda. (cid:50)
t t=1
8
3.2 AdynamicprogramforE
r,s
Forsomematrixa ∈ Rr×s letL : Rr×s → Rbegivenby
r s
(cid:88)(cid:88)
L(x) = a x .
ij ij
i=1 j=1
TheminimumofL(·)overthefinitesetE maybeobtainedbymeansofthefollowingrecursion:
r,s
For1 ≤ k ≤ r and1 ≤ (cid:96) ≤ sdefine
r s
(cid:110)(cid:88)(cid:88) (cid:111)
H(k,(cid:96)) = min a e : e ∈ E ,e = 1 ,
ij ij r,s k(cid:96)
i=k j=1
r s
(cid:110)(cid:88)(cid:88) (cid:111)
H(k,s+1) = min a e : e ∈ E .
ij ij r,s
i=k j=1
Then
min L(e) = H(1,s+1),
e∈Er,s
and
r s
(cid:88)(cid:88)
H(k,1) = a ,
ij
i=k j=1
s
(cid:16) (cid:88) (cid:17)
H(k,(cid:96)+1) = min H(k,(cid:96)), a +H(k+1,(cid:96)+1)
ij
j=(cid:96)+1
whereweusetheconventionsthatH(k+1,·) = 0and(cid:80)s · = 0. Intherecursionformulafor
j=s+1
H(k,(cid:96)+1),theterm(cid:80)s a +H(k+1,(cid:96)+1)istheminimumofL (e) = (cid:80)r (cid:80)s a e
j=(cid:96)+1 ij k i=k j=1 ij ij
overallmatricese ∈ E withe = 0ande = 1(if(cid:96) < s), whileH(k,(cid:96))istheminimum
r,s k(cid:96) k,(cid:96)+1
ofL (e)overalle ∈ E withe = 1.
k k,s k(cid:96)
Table1providespseudocodeforanalgorithmthatdeterminesaminimizerofL(·)overE .
r,s
3.3 Activesettypealgorithms
Throughout this exposition we assume that minimization of Q over an affine linear subspace of
Rp isfeasible. ThisiscertainlythecaseifQisaquadraticfunctional. IfQistwicecontinuously
differentiablewithpositivedefiniteHessianmatrixeverywhere,thisminimizationproblemcanbe
solvedwitharbitrarilyhighaccuracybyaNewtontypealgorithm.
All algorithms described in this paper alternate between two basic procedures which are
described next. In both procedures θ ∈ K is replaced with a vector θ ∈ K such that
new
Q(θ ) < Q(θ)unlessθ = θ.
new new
9
Algorithme ← DynamicProgram(a)
b ← (cid:0)(cid:80)s a (cid:1)
j=(cid:96) k,j k≤r,(cid:96)≤s+1
(cid:0) (cid:1)
H ← 0
k≤r+1,(cid:96)≤s+1
fork ← r downto1do
H ← H +b
k,1 k+1,1 k,1
for(cid:96) ← 1tosdo
(cid:0) (cid:1)
H ← min H ,b +H
k,(cid:96)+1 k,(cid:96) k,(cid:96)+1 k+1,(cid:96)+1
endfor
endfor
(cid:0) (cid:1)
e ← 0
k≤r,(cid:96)≤s
k ← 1,(cid:96) ← s
whilek ≤ r and(cid:96) ≥ 1do
ifH = H then
k,(cid:96)+1 k,(cid:96)
(e )r ← (1)r
i,(cid:96) i=k i=k
(cid:96) ← (cid:96)−1
else
k ← k+1
endif
endwhile.
Table1: MinimizingalinearfunctionaloverE
r,s
Basicprocedure1: Checkingoptimalityofθ ∈ K
Supposethatθ ∈ Ksatisfiesalreadythethefollowingtwoequations:
(5) ∇Q(θ)(cid:62)θ = 0 = ∇Q(θ)(cid:62)1.
Accordingto(3),thisvectorisalreadythesolutionθˆ if,andonlyif,∇Q(θ)(cid:62)e ≥ 0foralle ∈ E.
Thuswedetermine
∆ ∈ argmin ∇Q(θ)(cid:62)e
e∈E
anddothefollowing: If∇Q(θ)(cid:62)∆ ≥ 0,weknowthatθ = θˆ andstopthealgorithm. Otherwise
wedetermine
t = argmin Q(θ+t∆) > 0
o
t∈R
andreplaceθ with
θ := θ+t ∆.
new o
This vector θ lies in the cone K, too, and satisfies the inequality Q(θ ) < Q(θ). Then we
new new
proceedwithbasicprocedure2.
10