This paper concerns a fundamental class of convex matrix optimization problems. It presents the first algorithm that uses optimal storage and provably computes a lowrank approximation of a solution. In particular, when all solutions have low rank, the algorithm converges to a solution. This algorithm, SketchyCGM, modifies a standard convex optimization scheme, the conditional gradient method, to store only a small randomized sketch of the matrix variable. After the optimization terminates, the algorithm extracts a low-rank approximation of the solution from the sketch. In contrast to nonconvex heuristics, the guarantees for Sketchy-CGM do not rely on statistical models for the problem data. Numerical work demonstrates the benefits of SketchyCGM over heuristics. This paper discusses a fundamental class of convex matrix optimization problems with low-rank solutions. We argue that the main obstacle that prevents us from solving these problems at scale is not arithmetic, but storage. We exhibit the first provably correct algorithm for these problems with optimal storage. To explain the challenge, we consider the problem of low-rank matrix completion. Let X ∈ R m×n be an unknown matrix, but assume that a bound r on the rank of X is available, where r min{m, n}. Suppose that we record noisy observations of a subset E of entries from the matrix: The variables ξ ij ∈ R model (unknown) noise. The goal is to approximate the full matrix X . Matrix completion arises in machine learning applications, such as recommendation systems [34] . We can frame the matrix completion problem as a rank-constrained optimization: In general, the formulation (1) is intractable. Instead, we retrench to a tractable convex problem [7, 34] : The Schatten 1-norm · S1 returns the sum of the singular values of its argument; it is an effective proxy for the rank [13] . Adjusting the value of the parameter α modulates the rank of a solution X of (2). If we have enough data and choose α well, we expect that each solution X approximates the target matrix X . The convex problem (2) is often a good model for matrix completion when the number of observations |E| =Õ (r(m + n)), whereÕ suppresses log-like factors; see [7, 34] . We can write a rank-r approximation to a solution X using Θ(r(m + n)) parameters. Thus, we can express the problem and an approximate solution withÕ (r(m + n)) storage. Nevertheless, we need fully mn numbers to express the decision variable X for the optimization problem (2). The cost of storing the decision variable prevents us from solving large-scale instances of (2), even without worrying about arithmetic. This discrepancy raises a question: Is there an algorithm that computes an approximate solution to (2) using the optimal storageÕ (r(m + n))?  Here is another instance of the same predicament. Fix a vector x ∈ C n . Suppose that we acquire d noisy quadratic measurements of x with the form The a i ∈ C n are known measurement vectors, and the ξ i ∈ R model measurement noise. Given the data b and the vectors a i , the phase retrieval problem asks us to reconstruct x up to a global phase shift. Phase retrieval problems are prevalent in imaging science because it is easier to measure the intensity of light than its phase. In practice, the vectors a i are structured because they reflect the physics of the imaging system. For more details, see Appendix B.3 and the examples in [3, 8, 10, 22] . Let us outline a convex approach [3, 8, 10, 22 ] to the phase retrieval problem. The data (3) satisfies Thus, we can formulate phase retrieval as Now, pass to the convex problem We can estimate the parameter α ∈ R + from a i and b; see [39, Sec. II]. To approximate the true vector x , we compute a top eigenvector x of a solution to (5). This procedure is often an effective approach for phase retrieval when the number of measurements d = Θ(n); see [36, Sec. 2.8]. Once again, we recognize a discrepancy. The problem data b ∈ R d and the approximate solution x ∈ C n use storage Θ(n), but the matrix variable in (5) requires Θ(n 2 ) storage. We may ask: Is there an algorithm that computes an approximate solution to (5) using the optimal storage Θ(n)? Matrix completion and phase retrieval are examples of convex low-rank matrix optimization (CLRO) problems. Informally, this class contains convex optimization problems whose decision variable is a matrix and whose solutions are (close to) low rank. These problems often arise as convex relaxations of rankconstrained problems; however, the convex formulations are important in their own right. There has been extensive empirical and theoretical work to validate the use of CLROs in a spectrum of applications. For example, see [7, 8, 13, 16, 22] . Over the last 20 years, optimization researchers have developed a diverse collection of algorithms for CLRO problems. Surprisingly, every extant method lacks guarantees on storage or convergence (or both). Convex optimization algorithms dominated the early literature on algorithms for CLRO. The initial efforts, such as [21] , focused on interior-point methods, whose storage and arithmetic costs are forbidding. To resolve this issue, researchers turned to first-order convex algorithms, including bundle methods [20] , (accelerated) proximal gradient methods [2, 32, 35] , and the conditional gradient method (CGM) [11, 15, 19, 23, 25] . Convex algorithms are guaranteed to solve a CLRO. They come with a complete theory, including rigorous stopping criteria and bounds on convergence rates. They enjoy robust performance in practice. On the other hand, convex algorithms from the literature do not scale well enough to solve large CLRO problems because they operate on and store full-size matrices. The CGM iteration is sometimes touted as a lowstorage method for CLRO [23] . Indeed, CGM is guaranteed to increase the rank of an iterate by at most one per iteration. Nevertheless, the algorithm converges slowly, so intermediate iterates can have very high rank. CGM variants, such as [31, 39] , that control the rank of iterates lack storage guarantees or may not converge to a global optimum. Recently, many investigators have sought recourse in nonconvex heuristics for solving CLROs. This line of work depends on the factorization idea of Burer & Monteiro [6] , which rewrites the matrix variable as a product of two low-rank factors. There are many heuristic procedures, e.g., [4] [5] [6] 24] , that use clever initialization and nonlinear programming schemes in an attempt to optimize the factors directly. The resulting algorithms can have optimal storage costs, and they may achieve a fast rate of local convergence. There has been an intensive effort to justify the application of nonconvex heuristics for CLRO. To do so, researchers often frame unverifiable statistical assumptions on the problem data. For example, in the matrix completion problem (2), it is common to assume that the entries of the matrix are revealed according to some ideal probability distribution [7, 24] . When these assumptions fail, nonconvex heuristics can converge to the wrong point, or they may even diverge. Yurtsever, Udell, Tropp, Cevher Contributions. This paper explains how to extend the convex optimization algorithm CGM to obtain an approximate solution to a class of CLRO problems using optimal storage. Our algorithm operates much like CGM, but it never forms the matrix variable explicitly. Instead, we maintain a small randomized sketch of the matrix variable over the course of the iteration by means of a bespoke sketching method [37] . After the optimization method converges, we extract an approximate solution from the sketch. This technique achieves optimal storage, yet it converges under the same conditions and with the same guarantees as CGM. In summary, this paper presents a solution to the problems posed above: the first algorithm for convex lowrank matrix optimization problems that provably uses optimal storage to compute an approximate solution. We write · for the Euclidean norm, · F for the Frobenius norm, and · S1 for the Schatten 1-norm (aka the trace norm or the nuclear norm). Depending on context, ·, · refers to the Euclidean or Frobenius inner product. The symbol * denotes the conjugate transpose of a vector or matrix, as well as the adjoint of a linear map. The dagger † refers to the pseudoinverse. The symbol [M ] r stands for a best rank-r Frobenius-norm approximation of the matrix M . The function dist F (M ; S) returns the minimum Frobeniusnorm distance from M to a set S. The symbol denotes the semidefinite order. We use the computer science interpretation of the order notation O ,Õ , Ω, Θ. Let us begin with a generalization of the convex matrix completion formulation (2). In §5.5, we return to the psd setting of the phase retrieval problem (5). We consider a convex program with a matrix variable: The linear operator A : R m×n → R d and its adjoint Each coefficient matrix We interpret A X as a set of linear measurements of the matrix X. For example, in the matrix completion problem (2), the image A X lists the entries of X indexed by the set E. The function f : R d → R is convex and continuously differentiable. In many situations, it is natural to regard the objective function as a loss: f (A X) = loss(A X; b) for a vector b ∈ R d of measured data. By choosing the parameter α ∈ R + to be sufficiently small, we can often ensure that each minimizer of (6) is low-rank or close to low-rank. Our goal is to develop a practical algorithm that provably computes a low-rank approximation of a solution to the problem (6). To validate (6) as a model for a given application, one must undertake a separate empirical or theoretical study. We do not engage this question in our work. Suppose that we want to produce a low-rank approximation to a solution of a generic instance of the problem (6). What kind of storage can we hope to achieve? It is clear that we need Θ(r(m + n)) numbers to express a rank-r approximate solution to (6). We must also understand how much extra storage is incurred because of the specific problem instance (A , f ). It is natural to instate a black-box model for the linear map A , its adjoint A * , and the objective function f . For arbitrary vectors u ∈ R m and v ∈ R n and z ∈ R d , assume we have routines that can compute A (uv * ) and u * (A * z) and (A * z)v. We also assume routines for evaluating the function f and its gradient ∇f for any argument z ∈ R d . We may neglect the storage used to compute these primitives. Every algorithm based on these primitives must use storage Ω(m + n + d) just to represent their outputs. Thus, under the black-box model, any algorithm that produces a rank-r solution to a generic instance of (6) must use storage Ω(d + r(m + n)). We say that an algorithm is storage optimal if it achieves this bound. The parameter d often reflects the amount of data that we have acquired, and it is usually far smaller than the dimension mn of the matrix variable in (6). The problems that concern us are data-limited; that is, d mn. This is the situation where a strong structural prior (e.g., low rank or small Schatten 1-norm) is essential for fitting the data. This challenge is common in machine learning problems (e.g., matrix completion for recommendation systems), as well as in scientific applications (e.g., phase retrieval). To the best of our knowledge, no extant algorithm for (6) is guaranteed to produce an approximation of an optimal point and also enjoys optimal storage cost. Convex Low-Rank Matrix Optimization with Optimal Storage To develop our algorithm for the model problem (6), we must first describe a standard algorithm called the conditional gradient method (CGM). Classic and contemporary references include [11, 15, 19, 23, 25] . Here is the CGM algorithm for (6). Start with a feasible point, such as At each iteration t = 0, 1, 2, . . . , compute an update direction H t using the formulas (u t , v t ) = MaxSingVec(A * (∇f (A X t ))); MaxSingVec returns a left/right pair of maximum singular vectors. Update the decision variable: where η t = 2/(t + 2). The convex combination (11) remains feasible for (6) because X t and H t are feasible. CGM is a valuable algorithm for (6) because we can efficiently find the rank-one update direction H t by means of the singular vector computation (10). The weak point of CGM is that the rank of X t typically increases with t, and the peak rank of an iterate X t is often much larger than the rank of the solution of (6). See Figures 6 and 7 for an illustration. The CGM algorithm admits a simple stopping criterion. Given a suboptimality parameter ε > 0, we halt the CGM iteration when the duality gap δ t ≤ ε: Let X be an optimal point for (6). It is not hard to show [23, Sec. 2] that Thus, the condition (12) ensures that the objective value f (A X t ) is ε-suboptimal. The CGM iterates sat- The CGM iteration (9)- (11) requires Θ(mn) storage because it maintains the m × n matrix decision variable X t . We develop a remarkable extension of CGM that provably computes a rank-r approximate solution to (6) with working storage Θ(d + r(m + n)). Our approach depends on two efficiencies: • We use the low-dimensional "dual" variable z t = A X t ∈ R d to drive the iteration. • Instead of storing X t , we maintain a small randomized sketch with size Θ(r(m + n)). It is easy to express the CGM iteration in terms of the "dual" variable z t = A X t . We can obviously rewrite the formula (10) for computing the rank-one update direction H t in terms of z t . We obtain an update rule for z t by applying the linear map A to (11). Likewise, the stopping criterion (12) can be evaluated using z t and H t . Under the black-box model (8), the dual formulation of CGM has storage cost Θ(m + n + d). Yet the dual formulation has a flaw: it "solves" the problem (6), but we do not know the solution! Indeed, we must also track the evolution (11) of the primal decision variable X t . In the next subsection, we summarize a randomized sketching method [37] that allows us to compute an accurate rank-r approximation of X t but operates with storage Θ(r(m + n)). Suppose that X ∈ R m×n is a matrix that is presented to us as a stream of linear updates, as in (11). For a parameter r min{m, n}, we wish to maintain a small sketch that allows us to compute a rank-r approximation of the final matrix X. Let us summarize an approach developed in our paper [37] . Draw and fix two independent standard normal matrices Ω and Ψ where Ω ∈ R n×k with k = 2r + 1; Ψ ∈ R ×m with = 4r + 3. The sketch consists of two matrices Y and W that capture the range and co-range of X: We can efficiently update the sketch (Y , W ) to reflect a rank-one linear update to X of the form Both the storage cost for the sketch and the arithmetic cost of an update are Θ(r(m + n)). The following procedure yields a rank-r approximation X of the matrix X stored in the sketch (15). Yurtsever, Udell, Tropp, Cevher The matrix Q has orthonormal columns that span the range of Y . The extra storage costs of the reconstruction are negligible; its arithmetic cost is Θ(r 2 (m + n)). See [37, §4.2] for the intuition behind this method. It achieves the following error bound. Theorem 1 (Reconstruction error [37, Thm. 5.1]). Fix a target rank r. Let X be a matrix, and let (Y , W ) be a sketch of X of the form (14)-(15) . The procedure (17) yields a rank-r matrixX with Similar bounds hold with high probability. Remarks. The sketch size parameters (k, ) appearing in (14) are recommended to balance storage against reconstruction quality. See [37] and our follow-up work for more details and for other sketching methods. We are now prepared to present SketchyCGM, a storage-optimal extension of the CGM algorithm for the convex problem (6). This method delivers a provably accurate low-rank approximation to a solution of (6). See Algorithm 1 for complete pseudocode. Fix the suboptimality ε and the rank r. Draw and fix standard normal matrices Ω ∈ R n×k and Ψ ∈ R ×m as in (14) . Initialize the iterate and the sketches: At each iteration t = 0, 1, 2, . . . , compute an update direction via Lanczos or via randomized SVD [17]: Set the learning rate η t = 2/(t+2). Update the iterate and the two sketches: The iteration continues until it triggers the stopping criterion: At any iteration t, we can form a rank-r approximate solutionX t of the model problem (6) by invoking the procedure (17) with Y = Y t and W = W t . Suppose that the CGM iteration (9)- (11) generates the sequence (X t : t = 0, 1, 2, . . . ) of decision variables and the sequence (H t : t = 0, 1, 2, . . . ) of update directions. It is easy to verify that the SketchyCGM iteration (18)- (20) maintains the loop invariants z t = A X t and h t = A H t ; In view of the inequality (13) and the invariant (22), the stopping rule (21) ensures that X t is an ε-suboptimal solution to (6) when the iteration halts. Furthermore, Theorem 1 ensures that the computed solutionX t is a near-optimal rank-r approximation of X t at each time t. The total storage cost is Θ(d + r(m + n)) for the dual variable z t , the random matrices (Ω, Ψ), and the sketch (Y , W ). Owing to the black-box assumption (8), the algorithm completes the singular vector computations in (19) with Θ(d + m + n) working storage. At no point during the iteration do we instantiate an m × n matrix! Arithmetic costs are on the same order as the standard version of CGM. SketchyCGM is a provably correct method for computing a low-rank approximation of a solution to (6). Theorem 2. Suppose that the iterates X t from the CGM iteration (9)-(11) converge to a matrix X cgm . LetX t be the rank-r reconstruction of X t produced by SketchyCGM. Then In particular, if rank(X cgm ) ≤ r, then SketchyCGM always works in the fundamental case where each solution of (6) has low rank. Theorem 3. Suppose that the solution set S of the problem (6) contains only matrices with rank r or less. Then SketchyCGM attains E dist F (X t , S ) → 0. Suppose that the optimal point of (6) is stable in the sense that the value of the objective function increases as we depart from optimality. Then the SketchyCGM estimates converge at a controlled rate. Theorem 4. Fix κ > 0 and ν > 0. Suppose the (unique) solution X of (6) has rank(X ) ≤ r and Convex Low-Rank Matrix Optimization with Optimal Storage Algorithm 1 SketchyCGM for model problem (6) Input: Data for (6); suboptimality ε; target rank r Output: Rank-r approximate solutionX = U ΣV * of (6) in factored form for all feasible X. Then we have the error bound where C is the curvature constant [23, Eqn. (3)] of the problem (6). See Appendix A for the proofs of these results. Next, we present a generalization of the convex phase retrieval problem (5). Consider the convex template As before, A : C n×n → C d is a linear map, and f : It is easy to adapt SketchyCGM to handle (24) instead of (6). To sketch the complex psd matrix variable, we follow the approach described in [37, Sec. 7.3]. We also make small changes to the SketchyCGM iteration. Replace the computation (19) with (λ t , u t ) = MinEig(A * (∇f (z t ))); A (αu t u * t ), λ t ≤ 0 0, otherwise. MinEig returns the minimum eigenvalue λ t and an associated eigenvector u t of a conjugate symmetric matrix. This variant behaves the same as SketchyCGM. In this section, we demonstrate that SketchyCGM is a practical algorithm for convex low-rank matrix optimization. We focus on phase retrieval problems because they provide a dramatic illustration of the power of storage-optimal convex optimization. Appendix B adduces additional examples, including some matrix completion problems. To begin, we show that our approach to solving the convex phase retrieval problem (5) has better memory scaling than other convex optimization methods. All five methods solve (24) reliably. The proximal methods (PGM and AT) perform a full eigenvalue decomposition of the iterate at each step, but they can be accelerated by adaptively choosing the number of eigenvectors to compute. The methods based on CGM only need the top eigenvector, so they perform less arithmetic per iteration. To compare the storage costs of the five algorithms, let us consider a synthetic phase retrieval problem. We draw a vector x ∈ C n from the complex standard normal distribution. Then we acquire d = 10n phaseless measurements (3), corrupted with independent Gaussian noise so that the SNR is 20 dB. The measurement vectors a i derive from a coded diffraction pattern; see §B.3.2 for details. We solve the convex problem (5) with α equal to the average of the measurements b i ; see [39, Sec. II]. Then we compute a top eigenvector x of the solution.
In order to overcome such issue, REF recently proposed to combine conditional gradient and sketching techniques in order to maintain a low rank representation of the iterates.
Sketchy Decisions: Convex Low-Rank Matrix Optimization with Optimal Storage
AISTATS
Mathematics Computer Science