Low-rank approximation

In mathematics, low-rank approximation is a minimization problem, in which the cost function measures the fit between a given matrix (the data) and an approximating matrix (the optimization variable), subject to a constraint that the approximating matrix has reduced rank. The problem is used for mathematical modeling and data compression. The rank constraint is related to a constraint on the complexity of a model that fits the data. In applications, often there are other constraints on the approximating matrix apart from the rank constraint, e.g., non-negativity and Hankel structure.
Low-rank approximation is closely related to:
Definition.
Given
Basic low-rank approximation problem.
The unstructured problem with fit measured by the Frobenius norm, i.e., 
has analytic solution in terms of the singular value decomposition of the data matrix. The result is referred to as the matrix approximation lemma or Eckart–Young–Mirsky theorem. Let
be the singular value decomposition of formula_8 and partition formula_9, formula_10, and formula_11 as follows:
where formula_13 is formula_14, formula_15 is formula_16, and formula_17 is formula_18. Then the rank-formula_4 matrix, obtained from the truncated singular value decomposition
is such that
The minimizer formula_22 is unique if and only if formula_23.
Proof of Eckart–Young–Mirsky theorem (for spectral norm).
Let formula_24 be a real (possibly rectangular) matrix with formula_25. Suppose that 
is the singular value decomposition of formula_27. Recall that formula_28 and formula_11 are orthogonal matrices, and formula_30 is an formula_31 diagonal matrix with entries formula_32 such that formula_33.
We claim that the best rank formula_34 approximation to formula_27 in the spectral norm, denoted by formula_36, is given by
where formula_38and formula_39 denote the formula_40th column of formula_9 and formula_11, respectively.
First, note that we have
Therefore, we need to show that if formula_44 where formula_45 and formula_46 have formula_34 columns then formula_48.
Since formula_46 has formula_34 columns, then there must be a nontrivial linear combination of the first formula_51 columns of formula_11, i.e.,
such that formula_54. Without loss of generality, we can scale formula_55 so that formula_56 or (equivalently) formula_57. Therefore,
The result follows by taking the square root of both sides of the above inequality.
Proof of Eckart–Young–Mirsky theorem (for Frobenius norm).
Let formula_24 be a real (possibly rectangular) matrix with formula_25. Suppose that
is the singular value decomposition of formula_27.
We claim that the best rank formula_34 approximation to formula_27 in the Frobenius norm, denoted by formula_65, is given by
where formula_38 and formula_39 denote the formula_40th column of formula_9 and formula_11, respectively.
First, note that we have
Therefore, we need to show that if formula_44 where formula_45 and formula_46 have formula_34 columns then
By the triangle inequality with the spectral norm, if formula_78 then formula_79. Suppose formula_80 and formula_81 respectively denote the rank formula_34 approximation to formula_83 and formula_84 by SVD method described above. Then, for any formula_85
Since formula_87, when formula_88 and formula_89 we conclude that for formula_90
Therefore,
as required.
Weighted low-rank approximation problems.
The Frobenius norm weights uniformly all elements of the approximation error formula_93. Prior knowledge about distribution of the errors can be taken into account by considering the weighted low-rank approximation problem
where formula_95 vectorizes the matrix formula_27 column wise and formula_97 is a given positive (semi)definite weight matrix.
The general weighted low-rank approximation problem does not admit an analytic solution in terms of the singular value decomposition and is solved by local optimization methods, which provide no guarantee that a globally optimal solution is found.
In case of uncorrelated weights, weighted low-rank approximation problem also can be formulated in this way: for a non-negative matrix formula_97 and a matrix formula_27 we want to minimize formula_100 over matrices, formula_101, of rank at most formula_4.
Entry-wise formula_103 low-rank approximation problems.
Let formula_104. For formula_105, the fastest algorithm runs in formula_106 time. One of the important ideas been used is called Oblivious Subspace Embedding (OSE), it is first proposed by Sarlos.
For formula_107, it is known that this entry-wise L1 norm is more robust than the Frobenius norm in the presence of outliers and is indicated in models where Gaussian assumptions on the noise may not apply. It is natural to seek to minimize formula_108. For formula_109 and formula_110, there are some algorithms with provable guarantees.
Distance low-rank approximation problem.
Let formula_111 and formula_112 be two point sets in an arbitrary metric space. Let formula_113 represent the formula_114 matrix where formula_115. Such distances matrices are commonly computed in software packages and have applications to learning image manifolds, handwriting recognition, and multi-dimensional unfolding. In an attempt to reduce their description size, one can study low rank approximation of such matrices.
Distributed/Streaming low-rank approximation problem.
The low-rank approximation problems in the distributed and streaming setting has been considered in.
Image and kernel representations of the rank constraints.
Using the equivalences
and
the weighted low-rank approximation problem becomes equivalent to the parameter optimization problems
and
where formula_120 is the identity matrix of size formula_4.
Alternating projections algorithm.
The image representation of the rank constraint suggests a parameter optimization method in which the cost function is minimized alternatively over one of the variables (formula_122 or formula_123) with the other one fixed. Although simultaneous minimization over both formula_122 and formula_123 is a difficult biconvex optimization problem, minimization over one of the variables alone is a linear least squares problem and can be solved globally and efficiently.
The resulting optimization algorithm (called alternating projections) is globally convergent with a linear convergence rate to a locally optimal solution of the weighted low-rank approximation problem. Starting value for the formula_122 (or formula_123) parameter should be given. The iteration is stopped when a user defined convergence condition is satisfied.
Matlab implementation of the alternating projections algorithm for weighted low-rank approximation:
function [dh, f] = wlra_ap(d, w, p, tol, maxiter)
[m, n] = size(d); r = size(p, 2); f = inf;
for i = 2:maxiter
 % minimization over L
 bp = kron(eye(n), p);
 vl = (bp' * w * bp) \ bp' * w * d(:);
 l = reshape(vl, r, n);
 % minimization over P
 bl = kron(l', eye(m));
 vp = (bl' * w * bl) \ bl' * w * d(:);
 p = reshape(vp, m, r);
 % check exit condition
 dh = p * l; dd = d - dh;
 f(i) = dd(:)' * w * dd(:);
 if abs(f(i - 1) - f(i)) < tol, break, end
end 
Variable projections algorithm.
The alternating projections algorithm exploits the fact that the low rank approximation problem, parameterized in the image form, is bilinear in the variables formula_122 or formula_123. The bilinear nature of the problem is effectively used in an alternative approach, called variable projections.
Consider again the weighted low rank approximation problem, parameterized in the image form. Minimization with respect to the formula_123 variable (a linear least squares problem) leads to the closed form expression of the approximation error as a function of formula_122
The original problem is therefore equivalent to the nonlinear least squares problem of minimizing formula_133 with respect to formula_122. For this purpose standard optimization methods, e.g. the Levenberg-Marquardt algorithm can be used.
Matlab implementation of the variable projections algorithm for weighted low-rank approximation:
function [dh, f] = wlra_varpro(d, w, p, tol, maxiter)
prob = optimset(); prob.solver = 'lsqnonlin';
prob.options = optimset('MaxIter', maxiter, 'TolFun', tol); 
prob.x0 = p; prob.objective = @(p) cost_fun(p, d, w);
[p, f ] = lsqnonlin(prob); 
[f, vl] = cost_fun(p, d, w); 
dh = p * reshape(vl, size(p, 2), size(d, 2));
function [f, vl] = cost_fun(p, d, w)
bp = kron(eye(size(d, 2)), p);
vl = (bp' * w * bp) \ bp' * w * d(:);
f = d(:)' * w * (d(:) - bp * vl);
The variable projections approach can be applied also to low rank approximation problems parameterized in the kernel form. The method is effective when the number of eliminated variables is much larger than the number of optimization variables left at the stage of the nonlinear least squares minimization. Such problems occur in system identification, parameterized in the kernel form, where the eliminated variables are the approximating trajectory and the remaining variables are the model parameters. In the context of linear time-invariant systems, the elimination step is equivalent to Kalman smoothing.
A Variant: convex-restricted low rank approximation.
Usually, we want our new solution not only to be of low rank, but also satisfy other convex constraints due to application requirements. Our interested problem would be as follows,
This problem has many real world applications, including to recover a good solution from an inexact (semidefinite programming) relaxation. If additional constraint formula_136 is linear, like we require all elements to be nonnegative, the problem is called structured low rank approximation. The more general form is named convex-restricted low rank approximation.
This problem is helpful in solving many problems. However, it is challenging due to the combination of the convex and nonconvex (low-rank) constraints. Different techniques were developed based on different realizations of formula_136. However, the Alternating Direction Method of Multipliers (ADMM) can be applied to solve the nonconvex problem with convex objective function, rank constraints and other convex constraints, and is thus suitable to solve our above problem. Moreover, unlike the general nonconvex problems, ADMM will guarantee to converge a feasible solution as long as its dual variable converges in the iterations