Table Of ContentRandomized Extended Kaczmarz for Solving Least Squares
Anastasios Zouzias∗ Nikolaos M. Freris†
January 8,2013
3
1
0
2 Abstract
n
Wepresentarandomizediterativealgorithmthatexponentiallyconvergesinexpectationtothemini-
a
mumEuclideannormleastsquaressolutionofagivenlinearsystemofequations. Theexpectednumber
J
of arithmetic operations required to obtain an estimate of given accuracy is proportional to the square
5
conditionnumberofthesystemmultipliedbythenumberofnon-zerosentriesoftheinputmatrix. The
proposedalgorithmisanextensionoftherandomizedKaczmarzmethodthatwasanalyzedbyStrohmer
]
A andVershynin.
N
. 1 Introduction
h
t
a
m TheKaczmarzmethodisaniterativeprojectionalgorithmforsolvinglinearsystemsofequations[Kac37].
Due to its simplicity, the Kaczmarz method has found numerous applications including image recon-
[
struction, distributed computation and signal processing to name a few [FCM+92, Her80, Nat01, FZ12],
3 see [Cen81] for more applications. The Kaczmarz method has also been rediscovered in the field of im-
v
agereconstructionandcalledART(AlgebraicReconstructionTechnique)[GBH70],seealso[CZ97,Her80]
0
for additional references. It has been also applied to more general settings, see [Cen81, Table 1] and
7
7 [Tom55,McC75]fornon-linearversionsoftheKaczmarzmethod.
5 LetA ∈ Rm×n and b ∈ Rm. Throughout the paperallvectorsareassumed to becolumn vectors. The
. Kaczmarzmethodoperatesasfollows: Initially,itstartswithanarbitraryvectorx(0) ∈Rn.Ineachiteration,
5
0 the Kaczmarzmethod goes through the rowsof A in a cyclic manner1 and for eachselected row, sayi-th
2 rowA(i),itorthogonallyprojectsthecurrentestimatevectorontotheaffinehyperplanedefinedbythei-th
1 constraintofAx = b, i.e.,{x| A(i), x = b }whereh·, ·iistheEuclideaninnerproduct. Moreprecisely,
i
:
v assumingthatthei -throwhasbeenselectedatk-thiteration,thenthe(k+1)-thestimatevectorx(k+1) is
k (cid:10) (cid:11)
i inductivelydefinedby
X
r x(k+1) :=x(k)+λ bik − A(ik), x(k) A(ik)
a k (cid:10)A(ik) 2 (cid:11)
2
where λ ∈ R are the so-called relaxationparametersan(cid:13)d k·k(cid:13)denotes the Euclideannorm. The original
k (cid:13) 2(cid:13)
Kaczmarzmethod correspondstoλ = 1forallk ≥ 0andallother settingofλ ’sareusuallyreferredas
k k
therelaxedKaczmarzmethodintheliterature[Cen81,Gal03].
Kaczmarz proved that this process converges to the unique solution for square non-singular matri-
ces[Kac37],butwithoutanyattempttoboundtherateofconvergence.Boundsontherateofconvergenceof
theKaczmarzmethodaregivenin[McC75],[Ans84]and[Gal03,Theorem4.4,p.120].Inaddition,anerror
analysisoftheKaczmarzmethodunderthefiniteprecisionmodelofcomputationisgivenin[Kni93,Kni96].
∗Anastasios Zouzias is with the Department of Computer Science at the University of Toronto, Canada. E-mail:
[email protected]. Part ofthiswork wasdonewhile theauthorwasvisitingtheDepartmentofComputerScience at
PrincetonUniversity.
†NikolaosM.FrerisiswithIBMResearch-Zu¨rich,Sa¨umerstrasse4,8803Ru¨schlikon,Switzerland.E-mail:[email protected]
1Thatis,selectingtheindicesoftherowsfromthesequence1,2,...,m,1,2,....
1
Nevertheless, the Kaczmarz method converges even if the linear system Ax = b is overdetermined
(m> n)andhasnosolution. InthiscaseandprovidedthatAhasfullcolumnrank,theKaczmarzmethod
converges to the least squares estimate. This was first observed by Whitney and Meany [WM67] who
proved that the relaxed Kaczmarz method converges provided that the relaxation parameters are within
[0,2]andλ →0,seealso[CEG83,Theorem1],[Tan71]and[HN90]foradditionalreferences.
k
Intheliteraturetherewasempiricalevidencethatselectingtherowsnon-uniformlyatrandommaybe
moreeffectivethanselectingtherowsviaKaczmarz’scyclicmanner [HM93,FCM+92]. Towardsexplain-
ing such an empirical evidence, Strohmer and Vershynin proposed a simple randomized variant of the
Kaczmarzmethodthathasexponentialconvergenceinexpectation[SV09]assuming thatthelinearsystem
issolvable;seealso[LL10]forextensionstolinearconstraints. Arandomizediterativealgorithmthatcom-
putesasequenceofrandomvectorsx(0),x(1),...issaidtoconvergeinexpectationtoavectorx∗ ifandonly
ifE x(k)−x∗ 2 → 0ask → ∞,wheretheexpectationistakenovertherandomchoicesofthealgorithm.
2
Soonafter[SV09],NeedellanalyzedthebehavioroftherandomizedKaczmarzmethodforthecaseoffull
(cid:13) (cid:13)
colu(cid:13)mnrankl(cid:13)inearsystemsthatdonothaveanysolution [Nee10]. Namely, Needellprovedthatthe ran-
domizedKaczmarzestimatevectoris(inthelimit)within afixeddistancefromtheleastsquaressolution
andalsothatthisdistanceisproportionaltothedistanceofbfromthecolumnspaceofA. Inotherwords,
NeedellprovedthattherandomizedKaczmarzmethodiseffectiveforleastsquaresproblemswhoseleast
squareserrorisnegligible.
In this paper we present a randomized iterative least squares solver (Algorithm 3) that converges in
expectationtotheminimumEuclideannormsolutionof
min kAx−bk . (1)
x∈Rn 2
Theproposedalgorithmisbasedon[SV09,Nee10]andinspiredby[Pop99]. Morepreciselytheproposed
algorithmcanbethoughtofasarandomizedvariantofPopa’sextendedKaczmarzmethod[Pop99],there-
forewenameditasrandomizedextendedKaczmarz.
Organization of the paper In Section 2, we briefly discuss related work on the design of deterministic
and randomized algorithms for solving least squares problems. In Section 3, we present a randomized
iterative algorithm for projecting a vector onto a subspace (represented as the column space of a given
matrix) which maybe of independent interest. In addition, we discuss the convergence propertiesof the
randomizedKaczmarzalgorithmforsolvablesystems(Section3.2)andrecallitsanalysisfornon-solvable
systems(Section3.3). InSection4,wepresentandanalyzetherandomizedextendedKaczmarzalgorithm.
Finally,inSection5weprovideanumericalevaluationoftheproposedalgorithm.
2 Least squares solvers
Inthissectionwegiveabriefdiscussiononleastsquaressolversincludingdeterministicdirectanditerative
algorithmstogetherwithrecentlyproposedrandomizedalgorithms. Foradetaileddiscussionondetermin-
isticmethods, thereaderisreferredto[Bj96]. Inaddition,weplaceourcontribution incontextwithprior
work.
Deterministicalgorithms Intheliterature,severalmethodshavebeenproposedforsolvingleastsquares
problems of the form (1). Here we briefly describe a representative sample of such methods including
theuseofQRfactorizationwithpivoting,theuseofthesingularvaluedecomposition(SVD)anditerative
methodssuchasKrylovsubspacemethodsappliedonthenormalequations[Saa03].LAPACKprovidesro-
bustimplementationsofthefirsttwomethods;DGELSYusesQRfactorizationwithpivotingandDGELSD
usesthesingularvaluedecomposition[ABD+90]. Fortheiterativemethods,LSQRisequivalenttoapply-
ingtheconjugategradientmethodonthenormalequations[PS82]anditisarobustandnumericallystable
method.
2
Randomized algorithms To the best of our knowledge, most randomized algorithms proposed in the
theoreticalcomputerscienceliteratureforapproximatelysolvingleastsquaresaremainlybasedonthefol-
lowinggenerictwostepprocedure:firstrandomly(andefficiently)projectthelinearsystemintosufficiently
manydimensions,andsecondreturnthesolutionofthedown-sampledlinearsystemasanapproximation
tothe originaloptimalsolution [DMM06,Sar06,CW09,NDT09,MZ11,DMMS11], seealso[CW12]. Con-
centrationof measureargumentsimply thatthe optimal solution of the down-sampled system isclose to
theoptimalsolutionoftheoriginalsystem. Theaccuracyoftheapproximatesolutionusingthisapproach
depends on the sample size and to achieve relative accuracy ε, the sample size should depend inverse
polynomially on ε. This makes these approachesunsuitable for the high-precision regime of error that is
consideredhere.
A different approach is the so called randomized preconditioning method, see [RT08, AMT10]. The
authorsof[AMT10]implementedBlendenpik,ahigh-precisionleastsquaressolver. Blendenpikconsistsof
twosteps. Inthefirststep,theinputmatrixisrandomlyprojectedandaneffectivepreconditioningmatrix
isextractedfromtheprojectedmatrix.Inthesecondstep,aniterativeleastsquaressolversuchastheLSQR
algorithm of Paige and Saunders[PS82] is appliedon the preconditioned system. Blendenpikis effective
foroverdeterminedandunderdeterminedproblems.
AparalleliterativeleastsquaressolverbasedonnormalrandomprojectionscalledLSRNwasrecently
implemented by Meng, Saundersand Mahoney [MSM11]. LSRN consists of two phases. In the first pre-
conditioning phase, the original system is projected using random normal projection from which a pre-
conditioner is extracted. In the second step, an iterative method such as LSQR or the Chebyshev semi-
iterativemethod[GV61]isappliedonthepreconditionedsystem. Thisapproachisalsoeffectiveforover-
determined and under-determined least squares problems assuming the existence of a parallelcomputa-
tionalenvironment.
2.1 Relation with ourcontribution
InSection5,wecomparetherandomizedextendedKaczmarzalgorithmagainstDGELSY,DGELSD,Blenden-
pik. LSRN[MSM11]didnotperformwellunderasetupinwhichnoparallelizationisallowed, sowedo
notincludeLSRN’sperformance. ThenumericalevaluationofSection5indicatesthattherandomizedex-
tendedKaczmarziseffectiveonthecaseofsparse,well-conditionedandstronglyrectangular(bothoverde-
terminedandunderdetermined)leastsquaresproblems,seeFigure1. Moreover,therandomizedextended
Kaczmarzalgorithm has also comparable performance with LAPACK’s routine for the dense random in-
putmatrices, see Figure 2 (notice that the proposedalgorithm almost matchesBlendenpik’s performance
fortheunderdeterminedcase,seeFigure2(b)).Ontheotherhand,apreconditionedversionoftheproposed
algorithmdoesnotperformwellunderthecaseofill-conditionedmatrices,seeFigure3.
3 Background
PreliminariesandNotation Foranintegerm,let[m]:={1,...,m}. Throughoutthepaperallvectorsare
assumedtobecolumnvectors. WedenotetherowsandcolumnsofAbyA(1),...,A(m) andA ,...,A ,
(1) (n)
respectively(bothviewedascolumnvectors). R(A)denotesthecolumnspaceofA,i.e.,R(A) :={Ax|x∈
Rn}andR(A)⊥ denotestheorthogonalcomplementofR(A). Givenanyb ∈ Rm,wecanuniquelywriteit
asbR(A) +bR(A)⊥,wherebR(A) istheprojectionofbontoR(A). kAkF := mi=1 nj=1|aij|2 andkAk2 :=
maxx=6 0kAxk2/kxk2 denotes the Frobenius norm and spectral norm, respqePctivelyP. Let σ1 ≥ σ2 ≥ ... ≥
σrank(A) be the non-zero singular values of A. We will usually refer to σ1 and σrank(A) as σmax and σmin,
respectively.TheMoore-Pensorepseudo-inverseofAisdenotedbyA†[GL96]. Recallthat A† =1/σ .
2 min
Foranynon-zerorealmatrixA,wedefine
(cid:13) (cid:13)
(cid:13) (cid:13)
κ2(A):=kAk2 A† 2. (2)
F F 2
(cid:13) (cid:13)
(cid:13) (cid:13)
3
RelatedtothisisthescaledsquareconditionnumberintroducedbyDemmelin[Dem88],seealso[SV09].
It is easy to check that the above parameter κ2(A) is related with the condition number of A, κ2(A) :=
F
σ2 /σ2 ,viatheinequalities: κ2(A) ≤ κ2(A) ≤ rank(A)·κ2(A). Wedenotebynnz(·)thenumberofnon-
max min F
zeroentries of its argumentmatrix. We define the average row sparsity and averagecolumnsparsity of A by
R andC ,respectively,asfollows:
avg avg
m n
Ravg := qinnz A(i) and Cavg := pjnnz A(j)
Xi=1 (cid:16) (cid:17) Xj=1 (cid:0) (cid:1)
wherep := A 2/kAk2 foreveryi ∈ [n]andq := A(i) 2/kAk2 foreveryi ∈ [m]. Thefollowing fact
j (j) 2 F i 2 F
willbeusedextensivelyinthepaper.
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
Fact1. LetAbeanynon-zerorealm×nmatrixandb∈Rm. DenotebyxLS :=A†b. ThenxLS =A†bR(A).
Wefrequentlyusetheinequality1−t≤exp(−t)foreveryt≤1.Weconcludethissectionbycollectinga
fewbasicfactsfromprobabilitytheorythatwillbefrequentlyused. ForanyrandomvariableX,wedenote
its expectation by E[X] or EX. If X is a non-negative random variable, Markov’s inequality states that
P(X >t) ≤ t−1E[X]. LetX andY betworandomvariables,thenE[X +Y] = E[X]+E[Y]. Wewillrefer
tothisfactaslinearityofexpectation.LetE ,E ,...,E beasetofeventsdefinedoversomeprobabilityspace
1 2 l
holding with probabilities p ,p ,...p respectively, then P(E ∪E ∪...∪E ) ≤ l p . We refer to this
1 2 l 1 2 l i=1 i
factasunionbound.
P
3.1 RandomizedApproximateOrthogonal Projection
Algorithm1RandomizedOrthogonalProjection
1: procedure(A,b,T) ⊲A∈Rm×n,b∈Rm,T ∈N
2: Initializez(0) =b
3: fork =0,1,2,...,T −1 do
4: Pickjk ∈[n]withprobabilitypj := A(j) 22/kAk2F, j ∈[n]
A A⊤
5: Setz(k+1) = Im− k(Ajk) (kjk2) z(k) (cid:13)(cid:13) (cid:13)(cid:13)
(cid:18) (jk) 2 (cid:19)
6: endfor
7: Outputz(T)
8: endprocedure
Inthissectionwepresentarandomizediterativealgorithm(Algorithm1)that,givenanyvectorb∈Rm
andalinearsubspaceofRmrepresentedasthecolumnspaceofagivenmatrixA,approximatelycomputes
theorthogonalprojectionofbontothecolumnspaceofA(denotedbybR(A),bR(A) = AA†b),see[CRT11]
foradifferentapproach.
Algorithm1isiterative. Initially, itstartswithz(0) = b. Atthe k-thiteration, the algorithmrandomly
selects a column A of A for some j, and updates z(k) by projecting it onto the orthogonal complement
(j)
of the space of A . The claimis that randomly selecting the columns of A with probabilityproportional
(j)
to their square norms implies that the algorithm converges to bR(A)⊥ in expectation. After T iterations,
the algorithmoutputsz(T) andbyorthogonality b−z(T) servesasanapproximationforbR(A). Thenext
theoremboundstheexpectedrateofconvergenceforAlgorithm1.
Theorem 2. Let A ∈ Rm×n, b ∈ Rm and T > 1 be the input to Algorithm 1. Fix any integer k > 0. In exact
arithmetic,afterkiterationsofAlgorithm1itholdsthat
E z(k)−bR(A)⊥ 22 ≤ 1− κ21(A) k bR(A) 22.
(cid:13) (cid:13) (cid:18) F (cid:19) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13)
4
Moreover, each iteration of Algorithm 1 requires in expectation (over the random choices of the algorithm) at most
5C arithmeticoperations.
avg
kA⊤z(k)k
Remark1. AsuggestionforastoppingcriterionforAlgorithm1istoregularlycheck: 2 ≤εforsomegiven
kAk kz(k)k
F 2
accuracyε > 0. Itiseasytoseethatwheneverthiscriterionissatisfied,itholdsthat bR(A)⊥ −z(k) 2/ z(k) 2 ≤
εκF(A),i.e.,b−z(k) ≈bR(A). (cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
A A⊤
WedevotetherestofthissubsectiontoproveTheorem2. DefineP(j):=I − (j) (j) foreveryj ∈[n].
m kA k2
ObservethatP(j)P(j) = P(j),i.e.,P(j)isaprojectormatrix. LetX bearandomvaria(jb)le2 over{1,2,...,n}
thatpicksindexj withprobability A 2/kAk2. ItisclearthatE[P(X)]=I −AA⊤/kAk2. Laterwewill
(j) 2 F m F
makeuseofthefollowingfact.
(cid:13) (cid:13)
(cid:13) (cid:13)
Fact3. ForeveryvectoruinthecolumnspaceofA,itholds I − AA⊤ u ≤ 1− σm2in kuk .
m kAk2F 2 kAk2F 2
(cid:13)(cid:16) (cid:17) (cid:13) (cid:16) (cid:17)
Definee(k) :=z(k)−bR(A)⊥ foreveryk ≥0. Adirect(cid:13)(cid:13)calculationimp(cid:13)(cid:13)liesthat
e(k) =P(j )e(k−1).
k
Indeed,e(k) =z(k)−bR(A)⊥ =P(jk)z(k−1)−bR(A)⊥ =P(jk)(e(k−1)+bR(A)⊥)−bR(A)⊥ =P(jk)e(k−1)using
thedefinitionsofe(k),z(k),e(k−1) andthefactthatP(jk)bR(A)⊥ = bR(A)⊥ foranyjk ∈ [n]. Moreover,itis
easytoseethatforeveryk ≥ 0e(k) isinthecolumnspaceofA,sincee(0) = b−bR(A)⊥ = bR(A) ∈ R(A),
e(k) =P(j )e(k−1) andinadditionP(j )isaprojectormatrixforeveryj ∈[n].
k k k
LetX ,X ,...beasequenceofindependentandidenticallydistributedrandomvariablesdistributedas
1 2
X. Foreaseofnotation, wedenotebyE [·] = E [·|X ,X ,...,X ],i.e.,theconditionalexpectation
k−1 Xk 1 2 k−1
conditionedonthefirst(k−1)iterationofthealgorithm. Itfollowsthat
2 2
E e(k) =E P(X )e(k−1) = E P(X )e(k−1), P(X )e(k−1)
k−1 k−1 k k−1 k k
2 2
(cid:13) (cid:13) (cid:13) (cid:13) D E
(cid:13) (cid:13) =E (cid:13)e(k−1), P(X(cid:13))P(X )e(k−1) = e(k−1), E [P(X )]e(k−1)
(cid:13) (cid:13) k−1(cid:13) k(cid:13) k k−1 k
D AA⊤ E D σ2 2 E
≤ e(k−1) I − e(k−1) ≤ 1− min e(k−1)
2(cid:13) m kAk2! (cid:13) kAk2! 2
(cid:13) (cid:13) (cid:13) F (cid:13)2 F (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
whereweusedlinearityofe(cid:13)xpecta(cid:13)tio(cid:13)(cid:13)n,thefactthatP(·)isa(cid:13)(cid:13)projectormatrix,C(cid:13)auchy-(cid:13)Schwarzinequality
andFact3. Repeatingthesameargumentk−1timeswegetthat
2 1 k 2
E e(k) ≤ 1− e(0) .
2 κ2(A) 2
(cid:13) (cid:13) (cid:18) F (cid:19) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
Notethate(0) =b−bR(A)⊥ =bR(A) t(cid:13)ocon(cid:13)clude. (cid:13) (cid:13)
Step5 canbe rewrittenasz(k+1) = z(k) − A , z(k) / A 2 A . Ateveryiteration, the inner
(jk) (jk) 2 (jk)
pinroedxpuectctaantidontheeauchpditaetreatfiroonmrezq(ku)irtoesza(tk+m1)osrteq(cid:16)u(cid:10)njir=e1a5tpmjnonszt(cid:11)5An(nj(cid:13)(cid:13))z(cid:0)=A(j5(cid:13)(cid:13)kC)(cid:1)a(cid:17)vogpoepreartiaotniosnfso.rsomejk ∈ [n];hence
P (cid:0) (cid:1)
3.2 RandomizedKaczmarz
StrohmerandVershyninproposedthefollowingrandomizedvariantofKaczmarzalgorithm(Algorithm2),
see [SV09] for more details. The following theorem is a restatementof the main result of [SV09] without
imposingthefullcolumnrankassumption.
5
Algorithm2RandomizedKaczmarz[SV09]
1: procedure(A,b,T) ⊲A∈Rm×n,b∈Rm
2: Setx(0) tobeanyvectorintherowspaceofA
3: fork =0,1,2,...,T −1 do
4: Pickik ∈[m]withprobabilityqi := A(i) 22/kAk2F,i∈[m]
5: Setx(k+1) =x(k)+ bik−hx(k),A(ik)iA(cid:13)(ik) (cid:13)
kA(ik)k2 (cid:13) (cid:13)
2
6: endfor
7: Outputx(T)
8: endprocedure
Theorem4. LetA ∈ Rm×n,b ∈ Rm andT > 1betheinputtoAlgorithm2. AssumethatAx = bhasasolution
anddenotex :=A†b. Inexactarithmetic,Algorithm2convergestox inexpectation:
LS LS
k
2 1 2
E x(k)−x ≤ 1− x(0)−x ∀k >0. (3)
LS 2 κ2(A) LS 2
(cid:13) (cid:13) (cid:18) F (cid:19) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
Remark2. Theabovetheorem(cid:13) hasbeenp(cid:13)rovedin[SV09]forth(cid:13)ecaseoffu(cid:13)llcolumnrank. Also,therateofexpected
convergencein[SV09] is1−1/κ2(A) whereκ2(A) := kAk2/σ (A⊤A). Noticethatifrank(A) < n, then
F min(m,n)
κ2(A)isinfinitewhereasκ2(A)isbounded.
F
e e
We devote the restof this subsection to prove Theorem4 following [SV09]. The proof is based on the
e
following two elementary lemmas which both appeared in [SV09]. However, in our setting, the second
lemmaisnotidenticaltothatin[SV09]. WedeferredtheirproofstotheAppendix.
Lemma5(Orthogonality). AssumethatAx=bhasasolutionandusethenotationofAlgorithm2,thenx(k+1)−
x isperpendiculartox(k+1)−x(k)foranyk ≥0. Inparticular,inexactarithmeticitholdsthat x(k+1)−x 2 =
LS LS 2
x(k)−xLS 22− x(k+1)−x(k) 22. (cid:13)(cid:13) (cid:13)(cid:13)
(cid:13) The abo(cid:13)ve le(cid:13)mma provide(cid:13)sa formula for the error at each iteration. Ideally, we seek to minimize the
(cid:13) (cid:13) (cid:13) (cid:13)
errorateachiterationwhichisequivalenttomaximizing x(k+1)−x(k) overthechoiceoftherowprojec-
2
tionsofthealgorithm. ThenextlemmasuggeststhatbyrandomlypickingtherowsofAreducestheerror
(cid:13) (cid:13)
inexpectation. (cid:13) (cid:13)
Lemma 6 (Expected Error Reduction). Assume that Ax = b has a solution. Let Z be a random variable over
kA(i)k2
[m] with distribution P(Z =i) = 2 and assume that x(k) is a vector in the row space of A. If x(k+1) :=
kAk2
F
x(k)+ bZ−hx(k),A(Z)iA(Z) (inexactarithmetic),then
kA(Z)k2
2
2 1 2
E x(k+1)−x ≤ 1− x(k)−x . (4)
Z LS 2 κ2(A) LS 2
(cid:13) (cid:13) (cid:18) F (cid:19)(cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
Theorem4followsbyiteratin(cid:13)gLemma6,(cid:13)wegetthat (cid:13) (cid:13)
2 1 k 2
E x(k+1)−x ≤ 1− x(0)−x .
LS 2 κ2(A) LS 2
(cid:13) (cid:13) (cid:18) F (cid:19) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
3.3 RandomizedKaczmarzApplied toNoisy LinearSystems
The analysis of Strohmer and Vershynin is basedon the restrictive assumption that the linear system has
asolution. Needellmadeastepfurtherandanalyzedthemoregeneralsettinginwhichthelinearsystem
6
does not have any solution and A has full column rank [Nee10]. In this setting, it turns out that the ran-
domizedKaczmarzalgorithmcomputesanestimatevectorthatiswithinafixeddistancefromthesolution;
thedistanceisproportionaltothenormofthe“noisevector”multipliedbyκ2(A)[Nee10]. Thefollowing
F
theorem is a restatement of the main result in [Nee10] with two modifications: the full column rank as-
sumptionontheinputmatrixisdroppedandtheadditivetermγofTheorem2.1in[Nee10]isimprovedto
kwk2/kAk2. Theonlytechnicaldifferenceherefrom[Nee10]isthatthefullcolumnrankassumptionisnot
2 F
necessary,sowedefertheprooftotheAppendixforcompleteness.
Theorem 7. Assume that the system Ax = y has a solution for some y ∈ Rm. Denote by x∗ := A†y. Let xˆ(k)
denotethek-thiterateoftherandomizedKaczmarzalgorithmappliedtothelinearsystemAx =bwithb:=y+w
foranyfixedw ∈Rm,i.e.,runAlgorithm2withinput(A,b). Inexactarithmetic,itfollowsthat
2 1 2 kwk2
E xˆ(k)−x∗ ≤ 1− E xˆ(k−1)−x∗ + 2. (5)
2 κ2(A) 2 kAk2
(cid:13) (cid:13) (cid:18) F (cid:19) (cid:13) (cid:13) F
(cid:13) (cid:13) (cid:13) (cid:13)
Inparticular, (cid:13) (cid:13) (cid:13) (cid:13)
2 1 k 2 kwk2
E xˆ(k)−x∗ ≤ 1− x(0)−x∗ + 2.
(cid:13) (cid:13)2 (cid:18) κ2F(A)(cid:19) (cid:13) (cid:13)2 σm2in
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
4 Randomized Extended Kaczmarz
Given any least squares problem, Theorem 7 with w = bR(A)⊥ tells us that the randomized Kaczmarz
algorithmworkswellforleastsquareproblemswhoseleastsquareserrorisveryclosetozero,i.e.,kwk ≈0.
2
Roughlyspeaking,inthiscasetherandomizedKaczmarzalgorithmapproachestheminimumℓ -normleast
2
squaressolution uptoanadditiveerrorthatdependsonthedistancebetweenbandthecolumnspaceof
A.
In the present paper, the main observation is that it is possible to efficiently reduce the norm of the
“noisy” part of b, bR(A)⊥ (using Algorithm 1) and then apply the randomized Kaczmarz algorithm on
a new linear system whose right hand side vector is now arbitrarily close to the column space of A, i.e.,
Ax ≈ bR(A). This idea together with the observation that the least squares solution of the latter linear
system is equal (in the limit) to the least squares solution of the original system (see Fact 1) implies a
randomizedalgorithmforsolvingleastsquares.
NextwepresenttherandomizedextendedKaczmarzalgorithmwhichisaspecificcombinationofthe
randomizedorthogonalprojectionalgorithmtogetherwiththerandomizedKaczmarzalgorithm.
4.1 The algorithm
Wedescribearandomizedalgorithmthatconvergesinexpectationtotheminimumℓ -normsolutionvector
2
x (Algorithm3). Theproposedalgorithmconsistsoftwocomponents. Thefirstcomponentconsistingof
LS
Steps 5 and 6 is responsible to implicitly maintain an approximation to bR(A) formed by b − z(k). The
second component, consisting of Steps 4 and 7, applies the randomized Kaczmarz algorithm with input
A and the currentapproximationb−z(k) of bR(A), i.e., appliesthe randomized Kaczmarzon the system
Ax= b−z(k). Sinceb−z(k) convergestobR(A),x(k) willeventuallyconvergetotheminimumEuclidean
normsolutionofAx=bR(A) whichequalstoxLS =A†b(seeFact1).
Thestopping criterionof Step8 wasdecidedbasedonthe following analysis. Assumethatthe termi-
nation criteria are met for some k > 0. Let z(k) = bR(A)⊥ +w for some w ∈ R(A) (which holds by the
definitionofz(k)). Then,
A⊤z(k) 2 = A⊤(bR(A)⊥ +w) 2 = A⊤w 2 ≥σmin(A) z(k)−bR(A)⊥ 2.
(cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13)
7
Algorithm3RandomizedExtendedKaczmarz(REK)
1: procedure(A,b,ε) ⊲A∈Rm×n,b∈Rm,ε>0
2: Initializex(0) =0andz(0) =b
3: fork =0,1,2,... do
4: Pickik ∈[m]withprobabilityqi := A(i) 22/kAk2F,i∈[m]
5: Pickjk ∈[n]withprobabilitypj := (cid:13)(cid:13)A(j)(cid:13)(cid:13)22/kAk2F, j ∈[n]
6: Setz(k+1) =z(k)− hAk(Ajk),zk(k2)iA(jk) (cid:13)(cid:13) (cid:13)(cid:13)
(jk) 2
7: Setx(k+1) =x(k)+ bik−zi(kk)−hx(k),A(ik)iA(ik)
kA(ik)k2
2
8: Checkevery8min(m,n)iterationsandterminateifitholds:
Ax(k)−(b−z(k)) A⊤z(k)
2 ≤ε and 2 ≤ε.
(cid:13)(cid:13) kAkF x(k) 2 (cid:13)(cid:13) kA(cid:13)(cid:13)k2F x(k(cid:13)(cid:13)) 2
(cid:13) (cid:13) (cid:13) (cid:13)
9: endfor (cid:13) (cid:13) (cid:13) (cid:13)
10: Outputx(k)
11: endprocedure
Byre-arrangingtermsandusingthesecondpartoftheterminationcriterion,itfollowsthat z(k)−bR(A)⊥ 2 ≤
εkAk2F x(k) . Now, (cid:13) (cid:13)
σmin 2 (cid:13) (cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) A(x(k) −x ) ≤ Ax(k)−(b−z(k)) + b−z(k)−Ax
LS LS
2 2 2
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) ≤(cid:13)(cid:13)εkAkF x(k) + b(cid:13)(cid:13)R(A)⊥(cid:13)(cid:13)−z(k) (cid:13)(cid:13)
2 2
(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13)kAk2 (cid:13)(cid:13)
≤εkAk (cid:13)x(k)(cid:13) +(cid:13)ε F x(k) ,(cid:13)
F 2 σmin 2
(cid:13) (cid:13) (cid:13) (cid:13)
whereweusedthetriangleinequality,thefirstp(cid:13)(cid:13)artof(cid:13)(cid:13)thetermina(cid:13)(cid:13)tion(cid:13)(cid:13)ruletogetherwithbR(A) =AxLS and
theabovediscussion. Now,sincex(k),x ∈R(A⊤),itfollowsthat
LS
x(k)−x
LS 2 ≤εκ (A)(1+κ (A)). (6)
x(k) F F
(cid:13) 2(cid:13)
(cid:13) (cid:13)
Equation(6)demonstratesthattheforw(cid:13)arde(cid:13)rrorofREKafterterminationisbounded.
(cid:13) (cid:13)
4.2 Rate of convergence
ThefollowingtheoremboundstheexpectedrateofconvergenceofAlgorithm3.
Theorem8. AfterT > 1iterations, inexactarithmetic,Algorithm3withinputA(possiblyrank-deficient)andb
computesavectorx(T) suchthat
⌊T/2⌋
2 1
E x(T)−x ≤ 1− 1+2κ2(A) kx k2.
LS 2 κ2(A) LS 2
(cid:13) (cid:13) (cid:18) F (cid:19)
Proof. Forthesakeofnota(cid:13)(cid:13)tion, setα (cid:13)(cid:13)= 1−1/κ2F(A)anddeno(cid:0)tebyEk[·] :(cid:1)= E[·|i0,j0,i1,j1,...,ik,jk], i.e.,
theconditionalexpectationwithrespecttothefirstkiterationsofAlgorithm3. ObservethatSteps5and6
areindependentfromSteps4and7ofAlgorithm3,soTheorem 2impliesthatforeveryl ≥0
2
E z(l)−bR(A)⊥ 2 ≤αl bR(A) 22 ≤ bR(A) 22. (7)
(cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
8
Fixaparameterk∗ := ⌊T/2⌋. Afterthek∗-thiterationofAlgorithm3,itfollowsfromTheorem7(Inequal-
ity(5))that
E(k∗−1) x(k∗)−xLS 22 ≤α x(k∗−1)−xLS 22+ (cid:13)bR(A)⊥k−Akz2(k∗−1)(cid:13)22.
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) F (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
Indeed,therandomizedKa(cid:13)czmarzalgo(cid:13)rithm(cid:13)isexecutedwi(cid:13)thinput(A,b−z(k∗−1))andcurrentestimate
vectorx(k∗−1). Sety=bR(A)andw=bR(A)⊥−z(k∗−1)inTheorem7andrecallthatxLS =A†b=A†bR(A) =
A†y.
Now, averagingtheaboveinequalityoverthe randomvariablesi1,j1,i2,j2,...,ik∗−1,jk∗−1 andusing
linearityofexpectation,itholdsthat
E x(k∗)−x 2 ≤αE x(k∗−1)−x 2+ E bR(A)⊥ −z(k∗−1) 22 (8)
LS 2 LS 2 (cid:13) kAk2 (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) F (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) 2
(cid:13) (cid:13) ≤αE(cid:13)x(k∗−1)−x (cid:13)2+ bR(A) 2 byIneq.(7)
LS 2 (cid:13) kAk2(cid:13)
(cid:13) (cid:13) (cid:13) F(cid:13)
≤...≤(cid:13)(cid:13) αk∗ x(0)−x(cid:13)(cid:13) 2+k∗−2αl bR(A) 22, (repeattheabovek∗−1times)
LS 2 (cid:13) kAk2(cid:13)
(cid:13) (cid:13) Xl=0 (cid:13) F(cid:13)
≤kx k2+ ∞(cid:13)(cid:13) αl bR(A(cid:13)(cid:13)) 22, sinceα<1andx(0) =0.
LS 2 (cid:13) kAk2(cid:13)
Xl=0 (cid:13) F(cid:13)
Simplifyingtherighthandsideusingthefactthat ∞ αl = 1 =κ2(A),itfollows
l=0 1−α F
E x(k∗)−xLS 22 ≤PkxLSk22+ bR(A) 22/σm2in. (9)
(cid:13) (cid:13)
Moreover,observethatforeveryl(cid:13)(cid:13)≥0 (cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13)
E bR(A)⊥ −z(l+k∗) 22 ≤αl+k∗ bR(A) 22 ≤αk∗ bR(A) 22. (10)
Nowforanyk>0,similar(cid:13)(cid:13)(cid:13)considerationsas(cid:13)(cid:13)(cid:13)Ineq.(8)im(cid:13)(cid:13)pliest(cid:13)(cid:13)hat (cid:13)(cid:13) (cid:13)(cid:13)
E x(k+k∗)−x 2 ≤αE x(k+k∗−1)−x 2+ E bR(A)⊥ −z(k−1+k∗) 22
LS 2 LS 2 (cid:13) kAk2 (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) F (cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) ≤...≤(cid:13)(cid:13) αkE x(k∗)−x(cid:13)(cid:13) 2+k−1α(k−1)−lE bR(A)⊥ −z(l+k∗) 22 (byinduction)
LS 2 (cid:13) kAk2 (cid:13)
(cid:13) (cid:13) Xl=0 (cid:13) F (cid:13)
≤αkE x(k∗)(cid:13)(cid:13)−x 2+ α(cid:13)(cid:13)k∗ bR(A) 22 k−1αl (byIneq.(10))
LS 2 (cid:13)kAk2 (cid:13)
(cid:13) (cid:13) (cid:13) F (cid:13) Xl=0
≤αk k(cid:13)(cid:13)xLSk22+ bR(cid:13)(cid:13)(A) 22/σm2in +αk∗ bR(A) 22/σm2in (byIneq.(9))
=αk(cid:16)kxLSk22+(α(cid:13)(cid:13)k+αk(cid:13)(cid:13)∗) bR(A(cid:17)) 22/σm2(cid:13)(cid:13)in (cid:13)(cid:13)
≤αkkxLSk22+(αk+αk∗)κ(cid:13)(cid:13)2(A)k(cid:13)(cid:13)xLSk22 since bR(A) 2 ≤σmaxkxLSk2
≤αk∗(1+2κ2(A))kx k2. (cid:13) (cid:13)
LS 2 (cid:13) (cid:13)
Toderivethelastinequality,considertwocases. IfT iseven,setk = k∗,otherwisesetk = k∗+1. Inboth
cases,(αk+αk∗)≤2αk∗.
9
4.3 Theoretical boundson time complexity
In this section, we discuss the running time complexity of the randomized extended Kaczmarz (Algo-
rithm 3). Recall that REK is a Las-Vegas randomized algorithm, i.e., the algorithm always outputs an
“approximatelycorrect”leastsquaresestimate(satisfying(6))butitsrunnning timeisarandomvariable.
Givenanyfixedaccuracyparameterε>0andanyfixedfailureprobability0<δ <1weboundthenumber
ofiterationsrequiredbythealgorithmtoterminatewithprobabilityatleast1−δ.
Lemma9. Fixanaccuracyparameter0<ε<2andfailureprobability0<δ <1. Inexactarithmetic,Algorithm3
terminatesafteratmost
32(1+2κ2(A))
T∗ :=2κ2(A)ln
F δε2
(cid:18) (cid:19)
iterationswithprobabilityatleast1−δ.
Proof. Denoteα:=1−1/κ2(A)fornotationalconvenience. Itsufficestoprovethatwithprobabilityatleast
F
1−δtheconditionsofStep8ofAlgorithm3aremet. Insteadofprovingthis,wewillshowthat:
1. Withprobabilityatleast1−δ/2: (b−z(T∗))−bR(A) 2 ≤ε bR(A) 2/4.
2. Withprobabilityatleast1−δ/2: (cid:13)(cid:13)x(T∗)−xLS 2 ≤εkx(cid:13)(cid:13)LSk2/4(cid:13)(cid:13). (cid:13)(cid:13)
Laterwe prove thatItems (1) and (2)i(cid:13)mply the Le(cid:13)mma. First we prove Item (1). By the definition of the
(cid:13) (cid:13)
algorithm,
P (b−z(T∗))−bR(A) 2 ≥ε bR(A) 2/4 =P bR(A)⊥ −z(T∗) 22 ≥ε2 bR(A) 22/16
(cid:16)(cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (cid:17)≤ 16(cid:18)E(cid:13)(cid:13)(cid:13) z(T∗)−bR(A)⊥(cid:13)(cid:13)(cid:13) 22 (cid:13)(cid:13) (cid:13)(cid:13) (cid:19)
2
(cid:13)(cid:13)ε2 bR(A) 2 (cid:13)(cid:13)
≤16αT∗/ε(cid:13)2 ≤δ/2(cid:13)
(cid:13) (cid:13)
thefirstequalityfollowssinceb−bR(A) = bR(A)⊥,thesecondinequalityisMarkov’sinequality,thethird
inequalityfollowsbyTheorem2,andthelastinequalitysinceT∗ ≥κ2(A)ln(32).
F δε2
Now,weproveItem(2):
16E x(T∗)−x 2
P x(T∗)−x ≤εkx k /4 ≤ LS 2
LS 2 LS 2 (cid:13)ε2kx k2 (cid:13)
(cid:16)(cid:13) (cid:13) (cid:17) (cid:13) LS 2 (cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) ≤16α⌊T∗/2⌋(1+2κ2(A))/ε2 ≤δ/2.
thefirstinequalityisMarkov’sinequality,thesecondinequalityfollowsbyTheorem8,andthelastinequal-
ityfollowsprovidedthatT∗ ≥2κ2(A)ln 32(1+2κ2(A))
F δε2
Aunionboundonthecomplemento(cid:16)ftheabovet(cid:17)woevents(Item(1)and(2))impliesthatbothevents
happenwithprobabilityatleast1−δ. NowweshowthatconditioningonItems(1)and(2),itfollowsthat
REKterminatesafterT∗iterations,i.e.,
A⊤z(k)
Ax(T∗)−(b−z(T∗)) ≤εkAk x(T∗) and 2 ≤ε.
2 F 2 kA(cid:13)k2 x(k(cid:13))
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) F (cid:13) 2
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
Westartwiththefirstcondition. First,usingtriangleinequalityandItem2,(cid:13)itfoll(cid:13)owsthat
x(T∗) ≥kx k − x −x(T∗) ≥(1−ε/4)kx k . (11)
LS 2 LS LS 2
2 2
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
10
