An overview of low-rank matrix recovery from incomplete observations

0.78 MB·

by  Mark A. Davenport

An overview of low-rank matrix recovery from incomplete
observations
MarkA.DavenportandJustinRomberg∗
January24,2016
6
1
0
2
Abstract
n
a Low-rankmatricesplayafundamentalroleinmodelingandcomputationalmethodsforsignalprocessing
J
and machine learning. In many applications where low-rank matrices arise, these matrices cannot be fully

4 sampledordirectlyobserved,andoneencounterstheproblemofrecoveringthematrixgivenonlyincomplete
2
and indirect observations. This paper provides an overview of modern techniques for exploiting low-rank

structuretoperformmatrixrecoveryinthesesettings,providingasurveyofrecentadvancesinthisrapidly-
]
T developing field. Specific attention is paid to the algorithms most commonly used in practice, the existing
I theoreticalguaranteesforthesealgorithms,andrepresentativepracticalapplicationsofthesetechniques.
.
s
c
[ 1 Introduction

1
v Low-rankmatricesariseinanincrediblywiderangeofsettingsthroughoutscienceandappliedmathematics.
2 Tonamejustafewexamples,wecommonlyencounterlow-rankmatricesincontextsasvariedas:
2
4 ensemblesofsignals:theoutputofasensorarrayornetwork,acollectionofvideoframes,orasequence
6 • ofsegmentsofalongersignalcanoftenbehighlycorrelatedandrepresentedusingalow-rankmatrix[4,
0
46];
.
1
0 system identification: low-rank (Hankel) matrices correspond to low-order linear, time-invariant sys-
•
6 tems[100];
1
: adjacency matrices: the connectivity structure of many graphs, such as those that arise in manifold
v •
i learningandsocialnetworks,isoftenlowrank[99];
X
r distancematrices: inmanydataembeddingproblems—suchasthosethatariseinthecontextofmulti-
a •
dimensionalscaling[16],sensorlocalization[14,124],nuclearmagneticresonancespectroscopy[122,123],
andothers—thematrixofpairwisedistanceswilltypicallyhavearankdependentonthe(low)dimen-
sionofthespaceinwhichthedatalies;
itemresponsedata:low-rankmodelsarefrequentlyusedinanalyzingdatasetscontainingtheresponses
•
of various individuals to a range of items, such as survey data [135], educational data [90], the data
generatedbyrecommendationsystems[58,118,126],andothers;
machinelearning: low-rankmodelsareubiquitousinmachinelearning,layingthefoundationforboth
•
classicaltechniquessuchasprincipalcomponentanalysis[70,114]aswellasmodernapproachestomulti-
tasklearning[7,107]andnaturallanguageprocessing[15,49];
∗M.D.andJ.R.arewiththeSchoolofElectricalandComputerEngineering,GeorgiaInstituteofTechnology,Atlanta,GA30332USA
(e-mail:{mdav,jrom}@gatech.edu).M.D.issupportedbyNRLgrantN00173-14-2-C001,AFOSRgrantFA9550-14-1-0342,andNSFgrants
CCF-1409406,CCF-1350616,andCMMI-1537261.J.R.issupportedbyONRgrantN00014-11-1-0459,NSFgrantCCF-1422540,andagrant
fromthePackardFoundation.
1
quantumstatetomography: apurequantumstateofN ionscanbedescribedbya2N 2N matrixwith
• ×
rankone[61].
In all of these settings, the matrices we are ultimately interested in can be extremely large. Moreover, as
databecomesincreasinglycheaptoacquire,thepotentialsizewillcontinuetogrow. Thisraisesanumberof
challenges,butoftenakeyobstacleisthatfullyobservingthematrixofinterestcanprovetobeanimpossible
task: it can be prohibitively expensive to fully sample the entire output of a sensor array; we might only be
able to measure the strength of a few connections in a graph; and any particular user of a recommendation
system will provide only a few ratings. In such settings we are left with a highly incomplete set of obser-
vations, and unfortunately, many of the most popular approaches to processing the data in the applications
wherelow-rankmatrices ariseassumethatwe haveafully-sampleddata setandaregenerally notrobustto
missing/incompletedata.Inthesesituationsweareconfrontedwiththeinverseproblemofrecoveringthefull
matrixfromourincompleteobservations.
Whilesuchrecoveryisnotalwayspossibleingeneral,whenthematrixislowrank,itispossibletoexploit
this structure and to perform this kind of recovery in a surprisingly efficient manner. In fact, in recent years
therehasbeentremendousprogressinourunderstandingofhowtosolvesuchproblems.Whilemanyofthese
applications have a relatively long history in which various existing approaches to dealing with incomplete
observationshavebeenindependentlydeveloped,recentadvancesinthecloselyrelatedfieldofcompressive
sensing[23,43,48]haveenabledaburstofprogressinthelastfewyears. Wenowhaveaunifiedframework
whichprovidesastrongbaseoftheoreticalresultsconcerningwhenitispossibletorecoveralow-rankmatrix
fromincompleteobservationsusingefficient,practicalalgorithms[29,31,60,78,115,116].
In this paper we provide a survey of this developing field. We begin with a more formal mathematical
statement of the problem of low-rank matrix recovery in Section 2, followed in Section 3 by an overview of
some of the algorithms most commonly used in practice to solve these kinds of problems. We provide a
brief overview of the existing theoretical guarantees for these algorithms in Sections 4, 5, and 6 for several
concrete observation models with an emphasis on how many observations are required to reliably recover a
low-rank matrix and what additional assumptions are potentially required. Finally, in Section 7 we describe
animportantapplicationofthesetechniquestoanimportantclassofproblemswherewecansolvequadratic
andbilinearsystemsofequationsbyre-castingthemasasimpleproblemoflow-rankmatrixrecovery.
2 The Matrix Recovery Problem
We begin by carefully stating what we mean by low-rank matrix recovery. We observe a matrix X , which we
0
will assume to have size M N and which we can express either exactly or approximately as having rank
×
R M,N. Thismeansthatwecanwrite
(cid:28)
R
X σ u vT
0 ≈ k k k
k=1
(cid:88)
for some scalars σ ,σ ,...,σ 0 and orthonormal vectors u ,u ,...,u RM and v ,v ,...,v RN.
1 2 R 1 2 R 1 2 R
≥ ∈ ∈
The σ can be interpreted as the R largest singular values of X , and the u , v as the corresponding
k 0 k k
singu{lar}vectors. The collection of all such matrices form a union of subspac{es in}R{M}N; each set of vectors
×
u , v defineanR-dimensionalsubspace,andthe σ correspondtoanexpansioninthatsubspace. As
k k k
{ } { } { }
the u , v can vary continuously, this union contains uncountably many such subspaces. However, by
k k
{ } { }
counting the degrees of freedom in the singular value decomposition, we see that the set of rank R matrices
canbeparameterizedbyO(R(M+N))numbers,whichwhenRissmallismuchlessthantheMN requiredto
specifyageneralmatrix. Thissuggeststhatitmightbepossibletoaccuratelyrecoveralow-rankmatrixfrom
relativelyfewmeasurements,andasisnotedabove,recentresultshaveshownthatthisisindeedpossible[29,
31,60,78,115,116].
WewillassumethatratherthanobservingX directlyweinsteadobservey= (X )+zwherezrepresents
0 0
noiseand : RM N RL isalinearmeasurementoperatorthatactsonamatrixAX bytakingstandardinner
× 0
A →
2
productsagainstLpre-definedM N matricesA ,...,A :
1 L
×
M N
y = X ,A +z =trace(ATX )+z = X [m,n]A [m,n]+z . (1)
i (cid:104) 0 i(cid:105) i i 0 i 0 i i
m=1n=1
(cid:88) (cid:88)
Ourdiscussionbelowwillattimesinvolvetheadjointofthisoperator,whichisdefinedas
L
(w)= w A .
∗ i i
A
i=1
(cid:88)
Our survey focuses on three basic variations of this measurement model. The first is taking to be a
A
randomprojection,whereeachoftheA consistingofindependentandidenticallydistributedrandomvariables.
i
Althoughthismodelarisesinonlyalimitednumberofpracticalsituations, thetheoryissostreamlinedthat
itcanbeunderstoodalmostfromfirstprinciples(seeSection4). Foroursecondmodel, returnsasubsetof
A
the entries of the target. Recovering from these samples is known as the matrix completion problem. In this
case, each of the A has exactly one non-zero entry. The analysis for this problem, which we overview in
i
Section5,canalsobeextendedtoobservingasubsetoftheexpansioncoefficientsofX inafixed(andknown)
0
orthobasis. Thethirdmodel,whichwewilldiscussinSection7,isthattheA arerank-1matrices. Theseare
i
encounteredwheneachobservationcanbewrittenasaquadraticorbilinearforminX .
0
Whilethesearethemeasurementmodelsthathavereceivedthemostattentionintheliterature,theyareby
nomeanstheonlyinterestingmodels.Othermodelsinspiredbyapplicationsinimagingandsignalprocessing
havealsoappearedrecentlyintheliterature(seeforexample[4,5,8,86,130]).
3 Algorithms for Matrix Recovery
3.1 Low-rankapproximation
We start by reviewing the classical problem of finding the best low-rank approximation to a given M N
×
matrixX . By“best”,wemeanclosestinthesum-of-squaressense,andweformulatetheproblemas
0
minimize X X 2 subjectto rank(X)=R, (2)
X (cid:107) − 0(cid:107)F
where X X 2 = (X[m,n] X [m,n])2 isthesquareofthestandardFrobeniusnorm,andRisthe
(cid:107) − 0(cid:107)F m,n − 0
desired rank of the approximation. This problem is nonconvex, but is actually easy to solve explicitly using
(cid:80)
thesingularvaluedecomposition(SVD).Inparticular,ifwedecomposeX as
0
K
X =UΣVT = σ u vT,
0 k k k
k=1
(cid:88)
where K = min(M,N), U,V are M K and N K matrices with orthonormal columns u ,...,u and
1 K
× ×
v ,...,v ,andΣisadiagonalK K matrixwithsortedentriesσ σ σ 0,thenthesolutionto
1 K 1 2 K
× ≥ ≥···≥ ≥
(2)isfoundsimplybytruncatingthisexpansion:
R
X = σ u vT.
k k k
k=1
(cid:88)
(cid:99)
This is known as the Eckart-Young theorem; see [69, Chapter 7] for a detailed proof and discussion. For
medium scale problems, computing the SVD to high precision is tractable, with computational complexity
scalingasO(K2max(M,N)).
Computingthebestlow-rankapproximation,then,isakintothresholdingthesingularvalues: wetakethe
matrix,computeitsSVD,keepthelargesingularvalueswhilekillingoffthesmallones,andthenreconstruct.
Avariationoftheprogramabovemakesthisconnectionclearer. TheLagrangianof(2)is
minimize X X 2 +λ rank(X). (3)
X (cid:107) − 0(cid:107)F ·
3
Aswevarytheparameterλabove,thesolutiontotheprogramchanges—infact,thesetofsolutionsproduced
fordifferent0 λ< isexactlythesameasthesetofsolutionsfor(2)producedforall1 R K. Givenλ,
≤ ∞ ≤ ≤
wesolve(3)bycomputingtheSVD,hardthresholdingthesingularvaluesvia
σ , σ γ,
k k
σk(cid:48) = ≥ (4)
(cid:40)0, σk <γ,
withγ =√λ,andthentakingX =UΣ(cid:48)VT.
Acommonvariationofthealgorithmaboveinvolvesreplacingthehardthresholdin(4)withasoftthresh-
old. Inthiscasewestillsetthe(cid:99)singularvaluesthataresmalltozero,butnowthelargevaluesareshrunk:
σ γ, σ γ,
k k
σk(cid:48) = − ≥ (5)
(cid:40)0, σk <γ.
Thisamountstoamoregradualphasingoutofthetermsthatjustcrossthethreshold. Itturnsoutthatthissoft
thresholdingprocesscanalsobeputinvariationalform;whenγ =λ/2theresultoftheprocedureaboveisthe
solutionto
minimize X X 2 +λ X , (6)
X (cid:107) − 0(cid:107)F (cid:107) (cid:107)∗
where X isthenuclearnorm, andisequaltothesumofthesingularvaluesofX. X isalsoknownas
(cid:107) (cid:107)∗ (cid:107) (cid:107)∗
thetracenorm,asitisequaltothetracewhenX issymmetricpositivesemidefinite. Unliketherank, X is
(cid:107) (cid:107)∗
aconvexfunction,andoftenappearsasaconvexproxyforrankinoptimizationproblems[51]. Whileinthe
approximationproblemweareconsideringhere,thesolutionsto(3)and(6)involveverysimilarcomputations,
thiswillnotbethecaseatallwhenweconsiderrecoveryfrompartialobservationsinthenextsection.
In the procedures above, the computational cost is dominated by computing the SVD. For matrices with
M,N ontheorderof100–1000,thereareanumberofexactmethodswithsimilarcomputationalcomplexity—
see[68,Chap.45]foranoverview. WhenX islargebutisverywellapproximatedbyamatrixwithmodest
0
rank,randomizedalgorithmscanbeusedtocomputeanapproximateSVD[64,102,129].
3.2 Low-rankrecoveryandnuclearnormminimization
Thelow-rankapproximationproblemdescribedabovereadilyadmitsastraightforwardsolution. However,in
thissurveyweareconcernedinsteadwiththelow-rankrecoveryproblemwhereweareworkingfrom(possibly
noisy)indirectobservations,y (X ). Inthiscasewewouldideallyliketosolvetheanalogof(2),
0
≈A
minimize y (X) 2 subjectto rank(X)=R.
X (cid:107) −A (cid:107)2
Unfortunately,whereas(2)couldbesolvedviaasimpleSVD,(3.2)isingeneralNP-hard.
In contrast, the nuclear normminimization program remainstractable. In particular, withindirect obser-
vations,wereplace(6)with
minimize (X) y 2+λ X . (7)
X (cid:107)A − (cid:107)2 (cid:107) (cid:107)∗
Thisisanunconstrainedconvexoptimizationprogram,andcanbesolvedinasystematicwayusingaproximal
algorithm[41,112]. Thesolution(s)to(7)willobey,foranyγ >0,thefixedpointcondition[41]
X =prox (X γ ( (X ) y)
(cid:63) γ (cid:63)− A∗ A (cid:63) −
wheretheproximaloperatoris
prox (Z)=argmin X Z 2 +γλ X . (8)
γ X (cid:107) − (cid:107)F (cid:107) (cid:107)∗
Oneclassofmethodsforsolving(7)arebasedontheiteration
Xk+1 =proxγk(Xk−γkA∗(A(Xk)−y), (9)
4
forsomeappropriatelychosensequence γ [19]. Asdiscussedintheprevioussection,thesubproblem(8)is
k
{ }
solvedbysingularvaluesoft-thresholding.ComputingtheSVDoftheX γ ( (X ) y)ateachiteration
k k ∗ k
− A A −
is almost always the dominant cost, as it typically requires significantly more computation than applying
A
and .
∗
A
State-of-the-artmethodsforsolving(7)basedonsingularvaluethresholdingarenottoomuchmorecom-
plicated. Forexample,theFISTAalgorithm[11]modifiesthebasiciterationin(9)throughintelligentchoices
ofthescalingcoefficientγ andbyreplacingX withacarefullychosencombinationofX andX . These
k k k k 1
−
smallchangeshavealmostnoeffectontheamountofcomputationdoneateachiteration,buttheyconverge
insignificantlyfeweriterations.
The simplicity of these proximal-type algorithms makes them very attractive for small to medium sized
problems. However, as the number of rows and columns in the matrix gets to be several thousand, direct
computationoftheSVDbecomesproblematic. Forspeciallystructured ,includingtheimportantcasewhere
A
(X)returnsasubsetoftheentriesofX,fastalgorithmsthattakeadvantageofthisstructuretocomputethe
A
SVDhavebeendevelopedtosolve(7)orcloselyrelatedproblems[19,103].Inmoregeneralsettings,techniques
fromrandomizedlinearalgebrahavebeenappliedtocomputeapproximateSVDs[139].
Storage is also an issue when the target matrix is large. Since we expect the target matrix to have small
rank,wewouldliketosaveonstoragebyrestrictingtheiteratestoalsobelowrank. Wecanreformulatethe
programabovewithX LRT, andoptimizeovertheM RandN RmatricesLandRratherthanthe
≈ × ×
M N matrix X. This reformulation is driven by the fact that the nuclear norm is equal to the minimum
×
Frobeniusnormfactorization[116]:
1
X =min L 2 + R 2 subjectto X =LRT.
(cid:107) (cid:107)∗ L,R 2 (cid:107) (cid:107)F (cid:107) (cid:107)F
(cid:0) (cid:1)
Wecanthenreplace(7)with
λ λ
minimize (LRT) y 2+ L 2 + R 2, (10)
L,R (cid:107)A − (cid:107)2 2(cid:107) (cid:107)F 2(cid:107) (cid:107)F
andifthesolutionto(7)doesindeedhaverankatmostR,itwillalsobethesolutionto(10). Whilethisnew
formulationisnon-convex— (LRT)isacombinationofproductsofunknowns—thereareassurancesthat
A
the local minima in (10) are also global minima if the rank of the true solution is smaller than R [17]. This
technique is sometimes referred to as the Burer-Monteiro heuristic, after the authors of [17,18] who proposed
a version of the above for general semidefinite programming, and it is used in state-of-the-art large scale
implementationsofmatrixrecoveryproblems[117].
Theparameterλin(7)determinesthetrade-offbetweenthefidelityofthesolutiontothemeasurementsy
anditsconformancetothelow-rankmodel. Whenweareveryconfidentinthemeasurements,itmightmake
sensetousethemtodefineasetoflinearequalityconstraints,solving
minimize X subejectto (X)=y. (11)
X (cid:107) (cid:107)∗ A
Theoutputofthisprogramwillmatch(7)asλ 0. SomeoftheanalyticalresultswereviewinSections4and
→
5revealconditionsunderwhich(11)recoversalow-rankmatrixX exactlygivenmeasurementsy = (X )
0 0
A
asconstraints.
3.3 Iterativehardthresholding
Fromthealgorithmicpointofview, iterativehardthresholding(IHT)algorithms[73,131]areverysimilarto
theproximalalgorithmsusedtosolvenuclearnormminimizationintheprevioussection. However,whenthe
targetisverylowrank,theytendtoconvergeextremelyquickly.
Thebasiciterationisasfollows. FromthecurrentestimateX , wefirsttakeastepinthedirectionofthe
k
gradientof (X) y 2,thenprojectontothesetofrankRmatrices:
(cid:107)A − (cid:107)2
Y =X γ ( (X ) y))
k+1 k k ∗ k
− A A −
X =ProjectRank (Y ).
k+1 R k+1
5
The ProjectRank operator computes the top R left and right singular vectors and singular values — when
R
RissmallcomparedtoM andN,thiscanbedoneinsignificantlylesstimethancomputingafullSVD[121],
especially if the operator and its adjoint are structured in a such a way that there is a fast method for
∗
A A
applyingthematricesY toaseriesofvectors.Inthesecases,theintermediatematrixY isnotcomputed
k+1 k+1
explicitly, but each term in the first equation above can be handled efficiently in the SVD computation. The
implementationsofIHTin[73,131]relyonexistingsoftwarepackages[91]todothis.
In contrast to the nuclear norm minimization algorithm above, each of the iterates X IHT produces has
k
a prescribed rank. The storage required for X is roughly R(M +N), as opposed to the MN required for a
k
generalM N matrix. Thisdifferenceiscriticalforlarge-scaleapplications.
×
3.4 Alternatingprojections
Thealternatingprojectionsalgorithmisanotherspaceefficienttechniquewhichstorestheiteratesinfactored
form. Thealgorithmisextraordinarilysimple,andeasytointerpret: lookingforaM N rankRmatrixthat
×
isconsistentwithy:
minimize (X) y 2, subjectto rank(X)=R,
X (cid:107)A − (cid:107)2
isthesameaslookingforaM RmatrixLandaN RmatrixRwhoseproductisconsistentwithy:
× ×
minimize (LRT) y 2. (12)
L,R (cid:107)A − (cid:107)2
Thisoptimizationproblemisstillnon-convex,butwithoneofLorRfixed,itisasimpleleast-squaresproblem.
Thismotivatesthefollowingiteration. GivencurrentestimatesL ,R ,weupdateusing
k k
R =argmin (L RT) y 2,
k+1 R (cid:107)A k − (cid:107)2
(13)
L =argmin (LRT ) y 2.
k+1 L (cid:107)A k+1 − (cid:107)2
EachstepinvolvessolvingalinearsystemofequationswithRM orRN variablesforwhichwecandrawon
well-establishedalgorithmsinnumericallinearalgebra. Itssimplicityandefficiencymakeitoneofthemost
popularmethodsforlarge-scalematrixfactorization[85],andittendstooutperformnuclearnormminimiza-
tion[63,74],especiallyincaseswheretherankRisverysmallcomparedtoN,M.
There are few general convergence guarantees for alternating projections, and the final solution tends to
dependheavilyontheinitializationofLandR.However,guaranteesfortherateofconvergencecanbefound
in[77]andrecentwork[74]hasprovidedsomefirsttheoreticalresultsforconditionsunderwhichtheiterations
aboveconvergetothetruelow-rankmatrix(andamethodforsupplyingareliablestartingpoint). Theseare
discussedfurtherinSections4and5below.
Another advantage alternating projections is that the framework can be extended to handle structure on
one or both of the factors L and R. Typically, this means that the least-squares problems above are either
regularizedorconstrainedinamannerwhichencouragesorenforcesthedesiredstructure. Extensionsofthe
iterations above have been used successfully for problems including non-negative matrix factorization [80],
sparse PCA, where we restrict the number of nonzero terms in L,R [140], and dictionary learning [96,109],
whereLisawell-conditionedmatrixandRissparse. Again, therearefewstrongtheoreticalguaranteesfor
thesealgorithms,withnotableexceptionsintherecentworks[1,95].
Theoptimizationproblemsforalternatingprojections(12)andtheBurer-Monteiroheurtistic(10)aresimilar
(nonlinear)least-squaresproblemsonthematrixfactorsL,R.Thealgorithmsusedtosolvethem,though,have
adistinctdifference: insteadoffixingforonefactorandoptimizingtheotherasin(13)above,solversfor(10)
(e.g.,[117])typicallytakedescentstepsonLandRsimultaneously. Convergenceanalysisforcloselyrelated
localdescentmethodscanbefoundintherecentworks[13,134].
3.5 Otheralgorithmsformatrixrecovery
The methods above are by no means the only algorithms which have been proposed for low-rank matrix
recovery. Weclosethissectionbybrieflymentioningsomeothertechniques.
6
RecentyearshaveseenarenewedinterestinFrank-Wolfe-typealgorithmsforminimizingnormsdefined
by the convex hull of a set of atoms [40,71]. These algorithms are of particular interest for minimizing the
nuclear norm (where the atoms are rank-1 matrices), as they only require computing the leading singular
vector at every iteration, rather than a full SVD [72]. The nuclear norm problems in (7) and (11) can also be
minimizedbysolvingaseriesofweightedleast-squaresproblems,eachofwhichcanbesolvedusingstandard
techniquesforlinearsystemsofequations[54,55,104].
Beyondthenuclearnorm, otherproxiesforrankexist. Themax-norm[127,128]isaconvexfunctionthat
results in a similar optimization program to (7), and is subject to similar heuristics as (10) for storage reduc-
tion[117]. Alternatively,thelogarithmofthedeterminantisanonconvex,butsmooth,proxyforrankwhenX
ispositivesemidefinite(andcanbeappliedtogeneralmatricesbyembeddingtheminaPSDmatrix). Inprac-
tice,locallyminimizingthisfunctionsubjecttoconvexconstraintstendstoproducelow-ranksolutions[52].
Agreedyalgorithmforlow-rankrecoverywaspresentedin[93]. Thisalgorithmalternatesbetweenselect-
ing an estimate for the R-dimensional subspace in which X lives, and projecting onto these subspaces; it is
0
equipped with strong theoretical guarantees. Another class of algorithms evolves the left and right singular
vectorsalongtheGrassmanmanifold[42,78]. Thesealgorithms,whicharespecializedforthematrixcomple-
tionproblem,arecomputationallyefficient,andhaveequallystrongtheoreticalperformanceguarantees.
4 Matrix Recovery from Gaussian Observations
Thetheoryoflow-rankrecoveryiscleanandelegantwhenthemeasurementoperator ()isarandomprojec-
A ·
tion. Applicationswherethisisagoodmodelfortheobservationsarelimited,butlookingatitasanabstract
problemgivesusrealmathematicalinsightaboutwhylow-rankrecoveryworks.
The discussion in this section will center on “Gaussian” (), where the entries of each of the A are
m
A ·
independent and identically distributed normal random variables with zero mean and variance L 1 — this
−
varianceischosensothatE[ (X) 2]= X 2 foranyfixedM N matrixX.
(cid:107)A (cid:107)2 (cid:107) (cid:107)F ×
4.1 Thematrixrestrictedisometryproperty
We first examine the fundamental question of whether we can distinguish different rank-R matrices viewed
through the lens of the operator . One way to formalize this is asking whether preserves the distances
A A
betweenallsuchmatrices;thisiscertainlytrueifthereexistsa0 δ <1suchthat
≤
(1 δ) X X 2 (X ) (X ) 2 (1+δ) X X 2 (14)
− (cid:107) 1− 2(cid:107)F ≤ (cid:107)A 1 −A 2 (cid:107)2 ≤ (cid:107) 1− 2(cid:107)F
forallX ,X ofrankRorsmaller. Thisconditionisknownasthematrixrestrictedisometryproperty(matrix-
1 2
RIP), and is the matrix analog to the restricted isometry property from compressive sensing [43]. The first
immediate consequence is that all matrices of rank R or less have unique images, since if (X ) = (X )
1 2
A A
for X = X , then the lower bound would be violated. Qualitatively, the upper and lower bounds tell us
1 2
(cid:54)
that two rank R matrices are as distinguishable from their measurements as they would be if they were ob-
serveddirectly. With(14)established,wecanuseanynumberoftechniquestorecoveralow-rankmatrixfrom
measurementsthrough ;wediscusssomespecificguaranteesbelow.
A
AGaussian willalmostcertainlyhavethematrixrestrictedisometrypropertywhenthenumberofob-
A
servationsarecommensuratewiththenumberofdegreesoffreedomforanM N matrixwithrankR.Specif-
×
ically,(14)holdswithhighprobabilitywhen1
L (cid:38) R(N +M). (15)
This can be established through relatively simple probabilistic methods. We sketch the argument below; the
detailedproofin[28]smartlycombinesideasfrom[116]and[136].
Tobegin,noticethatitisenoughtoshowthat(1 δ) (X) 2 (1+δ)forallX ofrank2Randunit
− ≤ (cid:107)A (cid:107)2 ≤
Frobeniusnorm. Thentherearethreebasicstepsforestablishing(15).
1WeusethiscompactnotationtomeanthatthereisaconstantsothatL≥Const·R(M+N).
7
1. ForanyarbitraryfixedM N matrixX with X 2 =1,
× (cid:107) (cid:107)F
P (cid:107)A(X)(cid:107)22−1 >t ≤ Ce−cLt2 (16)
fort 1,whereC andcarereasonab(cid:0)l(cid:12)(cid:12)econstantsth(cid:12)(cid:12)atca(cid:1)nbecalculatedexplicitly(standardcalculations
≤ 2
yield C 2 and c 1/8). Since the entries of (X) are independent Gaussian random variables,
≤ ≥ A
(X) 2 isachi-squaredrandomvariablewithLdegreesoffreedom,andtheinequalityabovefollows
(cid:107)A (cid:107)2
fromstandardtailbounds[92].
We want (16) to hold not just for a single matrix, but uniformly over the set of all unit-norm rank 2R
matrices. It is straightforward to get a uniform result over any finite set of such matrices by using a
Q
unionbound:
P sup (X) 2 1 >t Ce cLt2.
(cid:107)A (cid:107)2− ≤ |Q|· −
(cid:18)X∈Q (cid:19)
(cid:12) (cid:12)
Eventhoughthesetofallunit-normr(cid:12)ank2Rmatric(cid:12)esisinfinite,thenextstepshowsusthatitisenough
toconsiderafinitesubset.
2. Let denotethesetofM N matricesX withrank(X) Rand X = 1. Nowlet beanyfinite
F (cid:15)
(cid:15)-apRproximationto —th×ismeansthat andforev≤eryX (cid:107) t(cid:107)hereisanX¯ in Rthatiswithin
(cid:15) (cid:15)
(cid:15): X X¯ (cid:15). TRhenashortargumentRsho⊂wRsthatfor(cid:15)=δ/(4√∈2)R, R
F
(cid:107) − (cid:107) ≤
max (X) 2 1 δ/2
X∈R(cid:15) (cid:107)A (cid:107)2− ≤
(cid:12) (cid:12)
impliesthat (cid:12) (cid:12)
sup (X) 2 1 δ.
(cid:107)A (cid:107)2− ≤
X
∈R
(cid:12) (cid:12)
Thusthesupremumovertheinfinitesetont(cid:12)herightcanb(cid:12)ereplacedwiththemaximumoverthefinite
setontheleft. Usingtheresultfromstep1,wenowhave
P sup (X) 2 1 >δ Ce cLδ2/4.
(cid:107)A (cid:107)2− ≤ |R(cid:15)|· −
(cid:18)X∈R (cid:19)
(cid:12) (cid:12)
(cid:12) (cid:12)
3. Thesizeof canbeestimatedasafunctionof(cid:15). Theboundin[28]reads
(cid:15)
R
R(M+N+1)
9
= eR(M+N+1)log(9/(cid:15)).
(cid:15)
|R |≤ (cid:15)
(cid:18) (cid:19)
Combiningthiswiththechoiceof(cid:15)aboveandpluggingitintotheresultofstep2,tellsusthat(14)will
holdwithprobabilityatleast1 C e c(cid:48)LwhenL(cid:38)δ 2R(M +N).
(cid:48) − −
−
Beyond basic identifiability, the matrix-RIP is also sufficient for concrete algorithms to produce accurate
estimatesoflow-rankmatrices. Forexample,[20,111]refineanargumentin[116]toshowthatnuclearnorm
minimizationwillrecoverrankRX whenδ 0.3in(14).Asimilarresultholdsforiterativehardthresholding
0
≤
techniques [73]and their variants[131]: when δ 1/3, weare guaranteed torecover rank-R X with linear
0
≤
convergence.Thebestknownguaranteesforalternatingminimizationareweaker:[74]statesthatδ Const/R
issufficientforlinearconvergence,whichisachievedforGaussian whenL(cid:38)R3(M +N). ≤
A
The uniformity of the matrix-RIP, that it holds for all pairs of rank R matrices simultaneously, results in
stabilityguaranteesforeachofthesealgorithmswhenthemeasurementsaremadeinthepresenceofnoise,or
thetargetmatrixisonlyapproximatelylowrank.
8
Null( )
S A
X0 0
H
(X )
B TB 0
(a) (b)
Figure 1: (a)Thesolutiontotheoptimizationprogram (11)isX0 whenthenormballB definedin (17)andtheaffine
space in (18) intersect at only one point. (b) The tangent cone TB(X0) from (19) along with a representative element
H ∈TB(X0).
4.2 ConvexgeometryandGaussianwidths
Onestrengthoftheanalysisdescribedabove(basedonthematrix-RIP)isthatitisquitegeneralandcanhelp
in analyzing a variety of different algorithms. However, the resulting bounds are often quite loose. In the
specific case of nuclear norm minimization, however, there are very precise conditions under which it will
succeed(orfail)withhighprobabilitywithaGaussianmeasurementensemble. Theseresultsareofaslightly
differentnaturethanthosebasedonthematrix-RIPintheprevioussection. Theyareweaker,inthattheyonly
speaktotherecoveryofasinglematrix;thetargetX isfixed,theensemble A : (cid:96) = 1 ...,L isgenerated
0 (cid:96)
{ }
independentlyofX ,andtheprobabilitywithwhichX isrecoverediscomputed. Buttheyarealsostronger
0 0
inthesensethattheyaremuchmorepreciseintellinguswhenrecoverywillsucceedorfail. Inthissection,
we will focus on results outlining the conditions for recovering X perfectly in the absence of noise. Given
0
y = (X ),wesolvetheequalityconstrainedproblem(11).
0
A
Recentworks[6,34,110]analyzetheperformanceofthisprogramusingveryintuitivegeometricalprinci-
ples. The target matrix X is a member of two different convex sets; it is in the nuclear norm ball of radius
0
X ,
0
(cid:107) (cid:107)∗
X := X : X X , (17)
0 0
∈B { (cid:107) (cid:107)∗ ≤(cid:107) (cid:107)∗}
anditisintheaffinespaceconsistingofallM N matricesthathavethesamemeasurements,
×
X := X : (X)=y . (18)
0
∈S { A }
By definition, X is the unique solution to (11) if and only if it is the only matrix in the intersection of
0
B
and ; this is illustrated in Figure 1(a). Generating with a Gaussian distribution is the same as choosing
S A
the orientation of uniformly at random. The local geometry of the tip of is determined by the rank of
S B
X , smaller ranks make this point more singular, and hence decrease the probability of an intersection. The
0
dimension of the set is MN L, the same as the null space of ; as L increases, gets smaller, and the
S − A S
probabilityofanintersectiondecreases.
We can make both of these statements precise. The collection of directions that lead into the ball (i.e.,
decreasethenuclearnorm) fromthepointX iscalledthetangentconeof atX :
0 0
B B
(X )= H : X +(cid:15)H X forsome(cid:15)>0 . (19)
0 0 0
TB { (cid:107) (cid:107)∗ ≤(cid:107) (cid:107)∗ }
Askingifthereisabetterfeasiblepointin(11)thanX isexactlythesameasaskingifthesubspace X =
0 0
S−
Null( )(theaffineset shiftedtotheorigin)intersectsthecone (X )atanyplaceotherthantheorigin.
0
A S TB
Whendoesarandomlychosensubspaceintersectaconeonlyattheorigin? Acleananswerisgivenin[34],
which builds on the classic work [59]. This answer depends on the notion of the Gaussian width of a set ,
C
whichisdefinedas
w( )=E sup X,G ,
C (cid:34)X∈C,(cid:107)X(cid:107)F=1(cid:104) (cid:105)(cid:35)
9
whereGisamatrixwhoseentriesareindependentzero-meanGaussianrandomvariableswithunitvariance,
andtheexpectationistakenwithrespecttoG.Thequantityw( )canbeinterpretedastheamountweexpecta
C
set toalignwitharandomlydrawnvector. Itisshownin[34]thatarandomlychosensubspacewillintersect
C
acone onlyattheoriginwithhighprobabilityifthecodimensionofthesubspaceisatleastaslargeasthe
C
squareoftheGaussianwidth:
codim(Null( )) w( (X ))2+1.
0
A ≥ TB
It is clear that for subsets of M N matrices, ω( ) √MN, and if is a D dimensional subspace, a quick
× C ≤ C
calculation shows that w( ) = √D. The codimension of Null( ) is always L, so measuring the Gaussian
C A
width of the tangent cone gives an immediate sufficient condition on the number of measurements needed
foraccuraterecovery. TheboundontheGaussianwidthw( (X ))2 3R(M +N R)givesusthesharp
0
TB ≤ −
sufficientconditionof
L 3R(M +N R)+1. (20)
≥ −
The relation in (20) is very close to being necessary as well. In [6,110], it is shown that if the number of
observationsisnottoofarbeloww( (X ))2,thentheprobabilitythat(11)failstorecoverX isverycloseto
0 0
TB
1. These papers contain an impressive suite of numerical experiments showing that the success or failure of
equality-constrainednuclearnormminimization(11)canbepredictedaccuratelyfromtheparametersM,N,R.
Itshouldbementionedthattheanalysisin[34]and[6]appliestomanydifferenttypesofstructuredrecov-
eryproblemsbasedonconvexoptimization; thelow-rankrecoveryresultsdiscussedaboveareanimportant
specialcase.
5 Matrix Completion
As we have just seen, the theory of low-rank matrix recovery is particularly elegant in the case where the
measurement operator () is a (Gaussian) random projection operator. While this provides some insight
A ·
into the kind of behavior we can hope for in many applications, it is also far from representative of the type
of observations we often encounter in practice. In particular, in many settings of interest the A are highly
i
structured—inthecaseofmatrixcompletion,theA willhaveonlyasinglenonzerovalueof1corresponding
i
totherowandcolumnoftheobservedelement. Anequivalentwaytothinkaboutthistypeofmeasurementis
thatweonlyobservetheentriesofX onasubsetΩofthecompletesetofentries. Thiskindofmeasurement
0
model arises in a variety of practical settings. For example, we might have a large number of questions we
wouldliketopotentiallyaskanumberofusers(asinalarge-scalesurveyorrecommendationsystem)butin
practicewemightonlyexpecttoreceiveresponsestoafewofthesequestionsfromanygivenuser. Similarly,
in many large-scale graphs (such as graphs representing the strength of social connections or the distances
betweensensorsorotheritems)wemightonlybeabletomeasure/observethestrengthofafewconnections
inagraph.
EveniftheelementsofΩarechosenatrandom,thisscenariohassomesignificantdifferencesfromthecase
wheretheA areGaussian. ItisclearthatwecannotexpectthetheorydevelopedusingGaussianwidthstobe
i
ofmuchuse,sinceourobservationsarenotGaussian,butthereisamorefundamentalproblemthatarisesin
thecaseofmatrixcompletion.
5.1 Whichmatricescanbecompleted?
In particular, the most immediate challenge we encounter when developing a theory of matrix completion
is that it is no longer possible to obtain the kind of uniform guarantees that apply to all low-rank matrices
describedinSection4. Toseewhy,consideramatrixwithrankone,butwhereone(orboth)singularvectors
are sparse, meaning that their energy is concentrated on just a few entries. When this occurs, as illustrated
inFigure2, theresultingmatrixwillalsohaveitsenergymostlyconcentratedonjustafewentries, inwhich
casemostentriesgiveusverylittleinformationandtherecoveryproblemishighlyill-posedunlessnearlyall
the entries are observed. (Another way to see this is to realize that if only a few entries of such a matrix are
observed,thematrixisverylikelytoliveinthenullspaceofthemeasurementoperator .) Moregenerally,if
A
10