Low-rank optimization with trace norm penalty

1 Introduction

• The method targets very low-rank solutions (p < 102) in very high-dimensional problems (n > 106), with per-iteration complexity and data storage linear in n.The motivating settings include machine learning and statistical large-scale regression.Source: §1 · p. 2
• The approach uses a low-rank factorization X = UBVT with a positive-definite, potentially non-diagonal scaling matrix B.Allowing B to be non-diagonal makes the factorization non-unique.Source: §1 · p. 2
• The algorithm alternates fixed-rank optimization with rank-one updates, using a Riemannian trust-region method whose iterates converge super-linearly to a local minimum.Rank is incremented to escape local minima; each rank-one update is selected to decrease cost.Source: §1 · p. 2
• For multiple regularization parameters, a predictor-corrector strategy computes the full path efficiently by transforming most fixed-parameter, fixed-rank local minima into global minima for other parameter values.The procedure modifies λ along the way rather than solving each parameter value independently.Source: §1 · p. 2
• For matrix completion, the method uses dense linear algebra and rank-one updates based on dominant singular vectors and values, rather than potentially computing all singular values at each iteration.The authors report that the predictor-corrector path strategy convincingly outperforms warm restarts.Source: §1 · pp. 2–3
• The paper illustrates the algorithm on low-rank matrix completion and multivariate linear regression, with numerical complexity linear in the number of observations.The authors report favorable convergence and precision properties in both applications.Source: §1 · p. 3

2 Relationship between convex program and non-convex formulation

• The rank-p factorization X = UBVT uses orthonormal U and V and a positive-definite scaling B; allowing B to be non-diagonal makes the representation redundant.In this parameterization, the trace norm is Trace(B), so it is differentiable, and the convex problem is recast at fixed rank.Source: §2 · p. 3
• The fixed-rank search space is a product of Stiefel manifolds and the positive-definite cone, with a metric accounting for rotational and scaling invariance.The row and column spaces are represented on the Stiefel factors, while scaling is carried by the positive-definite factor.Source: §2 · p. 3

2.1 First-order optimality conditions

• The paper relates the convex and fixed-rank formulations through their first-order optimality conditions, derived using either KKT conditions or the structured-space gradient.The dual variable is S = GradXf(UBVT).Source: §2.1 · pp. 3–4
• A local minimum of the fixed-rank problem is globally optimal for the convex problem exactly when the dual variable satisfies ∥S∥op = λ; otherwise ∥S∥op ≥ λ.The paper proposes checking global optimality using ∥S∥op − λ ≤ ǫ for a user-defined threshold.Source: §2.1 · p. 4

2.2 Duality gap computation

• The duality gap gives a non-negative measure of distance from the obtained solution to the optimum, beyond the global-optimality criterion alone.The paper presents the Lagrangian dual formulation for computing this gap.Source: §2.2 · pp. 4–5
• A practical dual candidate is the gradient S scaled to meet the operator-norm constraint, M = min{1, λ/∥S∥op}S.For f(X) = ψ(A(X)), the paper also gives a gap expression using the Fenchel conjugate of ψ when ∥A*(M)∥op ≤ λ.Source: §2.2 · p. 5

3 Manifold-based optimization to solve the non-convex problem (2)

• The fixed-rank formulation is optimized with a second-order trust-region method, which the paper says has quadratic convergence under its local assumptions.The trust-region method builds a local quadratic model and accepts or rejects a candidate based on objective decrease.Source: §3 · p. 5
• The factorization has rotational symmetry, so equivalent representations are identified on a quotient manifold where the minima are isolated.The search space is formed from a product of well-studied manifolds.Source: §3 · pp. 5–6

Tangent space of Mp

• Tangent vectors are split into vertical directions along an equivalence class and horizontal directions orthogonal to them under the chosen metric.Projecting onto the horizontal space requires solving a Lyapunov equation of complexity O(p^3).Source: Tangent space of Mp · p. 7

The Riemannian submersion (Mp, g)

• The quotient metric makes the gradient and Hessian computable through horizontal lifts in the total space.The Riemannian Hessian is represented using the directional derivative of the lifted gradient.Source: The Riemannian submersion (Mp, g) · pp. 7–8

Trust-region subproblem and retraction on Mp

• The trust-region subproblem is solved in the horizontal space with truncated conjugate gradients, and a retraction maps the resulting tangent direction back to the manifold.The retraction combines orthogonal-factor extraction and the matrix exponential for the positive-definite cone.Source: Trust-region subproblem and retraction… · pp. 8–9

Numerical complexity

• Manifold-related operations have complexity linear in n and m when p is much smaller than both dimensions.The listed operations include metric, projections, retraction, and Lyapunov-equation computations; problem-dependent costs are separate.Source: Numerical complexity · p. 9

4 An optimization scheme to solve convex program (1)

• The algorithm starts at rank 1 and alternates second-order optimization on the fixed-rank manifold with first-order rank-one updates into the next-rank space.The update is a descent step and the new point belongs to Mp+1.Source: §4 · p. 9
• Each rank-one update decreases the objective for sufficiently small positive β when its directions are the dominant left and right singular vectors of the dual variable.The maximum decrease is attained at β = (σ1 − λ)/Lf, where Lf is the gradient Lipschitz constant.Source: §4 · pp. 9–10
• Global convergence follows because both fixed-rank optimization and rank updates are descent iterates, although the algorithm has no theoretical guarantee of stopping at the optimal rank.The authors report empirical convergence at p ≪ min{m, n}, most often at the optimal rank.Source: §4 · pp. 10–11
• The method controls rank at every intermediate iterate more tightly than the cited trace-norm minimization algorithms [CCS10, TY10, MGC11, MHT10].The stopping thresholds for the non-convex and convex problems are chosen separately, so rank increments need not wait for trust-region convergence.Source: §4 · p. 11

5 Regularization path

• Warm restarts ignore the regularization path’s smoothness, whereas the predictor-corrector scheme predicts from two prior solutions on the same rank manifold and corrects by initializing the solver at the prediction.When consecutive solutions belong to different rank manifolds, the scheme uses a warm restart instead.Source: §5 · p. 11
• The predictor-corrector approach outperforms warm restarts in prediction accuracy with minimal extra computations, while potentially making every iterate a global solution for some λ.The full path uses r warm restarts and N − r predictor steps, where r is the number of rank increments and N is the number of λ values.Source: §5 · p. 11, Table 2 · p. 12
• The method assumes a unique solution for every λ; strict convexity of f is sufficient and can be enforced by adding a small multiple of the squared Frobenius norm.The solution path’s smoothness follows from the Implicit Function Theorem; the paper also gives a dual-projection argument.Source: §5 · p. 12

Predictor step on the quotient manifold Mp

• On the quotient manifold, the predictor approximates the geodesic between two same-rank solutions, extends that approximation to estimate the next solution, and uses backtracking line search.The logarithmic mapping has no analytic expression for this case, so the method uses an approximation; the prediction step has complexity O(np^2 + mp^2 + p^3).Source: Predictor step on the quotient… · pp. 12–13
• Figure 2 traces the predictor-corrector solution path, with the blue line denoting the curve of optimal solutions.The predictor is used on a fixed-rank manifold.Source: §5 · p. 11, Figure 2 · p. 13

6.1 Diagonal versus matrix scaling

• In a 200 × 200 rank-5 completion example, constraining the factor B to be diagonal worsens gradient-descent convergence compared with allowing matrix scaling.The experiment fixes λ at 10^-10; the text describes the observation as generic beyond this particular problem.Source: §6.1 · p. 14, Figure 3 · p. 14

6.2 Low-rank matrix completion

• For sparse matrix-completion observations, each trust-region iteration costs O(|W|p + np^2 + mp^2 + p^3), while dominant singular-value and vector computation costs O(|W|); the method assumes the known-entry lower bound is met.Source: §6.2 · p. 15

6.2.2 Regularization path for matrix completion

• The predictor-corrector method computes global minima across λ values and gives better predictions than a pure warm-restart approach.The reported path uses a geometric λ sequence from 1 × 10^3 to 1 × 10^-3 with reduction factor γ = 0.95.Source: §6.2.2 · p. 17, Figure 5 · p. 18

Convergence behavior with varying λ

• Descent-restart + TR outperforms FPCA, SOFT-IMPUTE, and APG across the tested λ values, with the fewest iterations.In a separate convergence test, it sweeps through ranks before trust-region convergence accelerates at the optimal rank; APG exceeds the iteration limit and SOFT-IMPUTE converges to a different solution despite faster initial progress.Source: Convergence behavior with varying λ · p. 19, Figure 6 · p. 20, Convergence test · p. 20

Fixed-rank comparison

• On 50000 × 50000 rank-5 completion, the trust-region scheme is competitive with LMaFit but performs more poorly than LRGeom, despite similar asymptotic convergence rates.LMaFit has better per-iteration time complexity, but its convergence suffers on this large-scale problem.Source: Figure 9 · p. 22, Fixed-rank comparison · p. 22

Comparison with APGL

• On large-scale completion, the framework competes favorably with APGL for very low ranks and higher accuracy, while APGL performs better as rank increases and at moderate accuracy.Traversing ranks one by one is costly; the paper identifies computing the full regularization path as the setting where this approach is most effective.Source: Figure 10 · p. 23, Comparison with APGL · p. 23

Comments on matrix completion algorithms

• Across random-example benchmarks, Descent-restart + TR performs favorably and is especially suited to low-rank solutions and regularization-path computation.The summary identifies strict movement in low-rank space as enabling a predictor-corrector path beyond standard warm restarts.Source: Comments on matrix completion algorithms · p. 24

6.3 Multivariate linear regression

• For multivariate regression, trace-norm regularization models low-rank coefficients when responses are related; the per-iteration cost includes O(nqp + nkp + kp^2 + nk).The framework is described as applicable to other smooth loss functions, beyond the quadratic loss considered here.Source: §6.3 · p. 24
• Low-rank factorization reduces per-iteration numerical complexity from O(q^2k) to O(q^2p + qkp), which is quadratic and linear in n.The comparison is with using the full matrix W, whose per-iteration numerical complexity is O(q^2k).Source: §6.3 · p. 25

Fenchel dual and duality gap for multivariate linear regression

• For W=UBV^T, the duality gap is f(W)+λ∥W∥∗+ψ∗(M), and computing ψ∗(M) costs O(nk).The passage says the cost of computing ψ∗(M) is of the same order as computing the factored objective.Source: Fenchel dual and duality gap for… · p. 25

6.3.1 Regularization path for multivariate linear regression

• The regularization-path experiment uses synthetic Gaussian data, four SNR values, and results averaged over five random 70/30 train-test splits.The design matrix has size 5000 × 120, and the ground-truth coefficient matrix has size 120 × 100 and rank 5; ranks are restricted to below 30 along the path.Source: §6.3.1 · p. 25, Figure 11 · p. 26

• The framework uses increasing fixed-rank problems, a numerically tractable duality gap, and rank increments that decrease cost.The duality gap monitors convergence to the original problem.Source: §7 · p. 26
• A second-order trust-region algorithm on a quotient search space guarantees quadratic convergence with proper parameter tuning.Source: §7 · p. 26
• The predictor-corrector approach has superior performance to the warm-restart approach.Source: §7 · p. 27
• Applications to low-rank matrix completion and multivariate linear regression produced encouraging numerical results.Source: §7 · p. 27

• The gradient component for U is SVB.Source: §A.1 · p. 27
• The gradient component for B is U^T SV + λI.Source: §A.1 · p. 27
• The gradient component for V is S^T UB.Source: §A.1 · p. 27

• The trace-norm subdifferential consists of UV^T plus an orthogonal component whose operator norm is at most one.Source: §A.2 · p. 27
• At a stationary point, Sym(U^T SV) + λI = 0 implies λ = ∥Sym(U^T SV)∥op ≤ ∥S∥op.Source: §A.2 · p. 27
• Equality holds when U and V span dominant row and column subspaces of S, yielding the global optimality condition.Source: §A.2 · p. 27

• The proof assumes f is smooth and its gradient is Lipschitz continuous with constant Lf.Source: §A.3 · p. 27
• The rank-one update X+ = X − βuv^T is bounded using the trace-norm increase and the leading singular value of S.The objective bound includes the term β(σ1 − λ − Lf ...).Source: §A.3 · p. 28
• With Lf > 0, the maximizing step size βmax is zero only at optimality.Source: §A.3 · p. 28

A.4 Proof of Proposition (2.3)

• Introducing Z with the constraint Z = X separates the trace-norm penalty in the reformulated optimization problem.The Lagrangian uses M as the dual variable for this constraint.Source: §A.4 · p. 28
• The derivation uses the operator norm as the dual norm of the trace norm and uses the Fenchel conjugate of f.These yield the stated dual function under ∥M∥op ≤ λ.Source: §A.4 · p. 28
• The dual formulation maximizes −f*(M) subject to ∥M∥op ≤ λ.The dual function reduces to −f*(M) when the operator norm constraint holds.Source: §A.4 · p. 28