Low-rank approximation

(Redirected from Eckart-Young theorem)

In mathematics, low-rank approximation refers to the process of approximating a given matrix by a matrix of lower rank. More precisely, it 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 numerous other techniques, including principal component analysis, factor analysis, total least squares, latent semantic analysis, orthogonal regression, and dynamic mode decomposition.

Definition

edit

Given

  • structure specification  ,
  • vector of structure parameters  ,
  • norm  , and
  • desired rank  ,
 

Applications

edit

Basic low-rank approximation problem

edit

The unstructured problem with fit measured by the Frobenius norm, i.e.,

 

has an 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. This problem was originally solved by Erhard Schmidt[1] in the infinite dimensional context of integral operators (although his methods easily generalize to arbitrary compact operators on Hilbert spaces) and later rediscovered by C. Eckart and G. Young.[2] L. Mirsky generalized the result to arbitrary unitarily invariant norms.[3] Let

 

be the singular value decomposition of  , where   is the   rectangular diagonal matrix with the singular values  . For a given  , partition  ,  , and   as follows:

 

where   is  ,   is  , and   is  . Then the rank-  matrix, obtained from the truncated singular value decomposition

 

is such that

 

The minimizer   is unique if and only if  .

Proof of Eckart–Young–Mirsky theorem (for spectral norm)

edit

Let   be a real (possibly rectangular) matrix with  . Suppose that

 

is the singular value decomposition of  . Recall that   and   are orthogonal matrices, and   is an   diagonal matrix with entries   such that  .

We claim that the best rank-  approximation to   in the spectral norm, denoted by  , is given by

 

where  and   denote the  th column of   and  , respectively.

First, note that we have

 

Therefore, we need to show that if   where   and   have   columns then  .

Since   has   columns, then there must be a nontrivial linear combination of the first   columns of  , i.e.,

 

such that  . Without loss of generality, we can scale   so that   or (equivalently)  . 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)

edit

Let   be a real (possibly rectangular) matrix with  . Suppose that

 

is the singular value decomposition of  .

We claim that the best rank   approximation to   in the Frobenius norm, denoted by  , is given by

 

where   and   denote the  th column of   and  , respectively.

First, note that we have

 

Therefore, we need to show that if   where   and   have   columns then

 

By the triangle inequality with the spectral norm, if   then  . Suppose   and   respectively denote the rank   approximation to   and   by SVD method described above. Then, for any  

 

Since  , when   and   we conclude that for  

 

Therefore,

 

as required.

Weighted low-rank approximation problems

edit

The Frobenius norm weights uniformly all elements of the approximation error  . Prior knowledge about distribution of the errors can be taken into account by considering the weighted low-rank approximation problem

 

where   vectorizes the matrix   column wise and   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:[4][5] for a non-negative matrix   and a matrix   we want to minimize   over matrices,  , of rank at most  .

Entry-wise Lp low-rank approximation problems

edit

Let  . For  , the fastest algorithm runs in   time.[6][7] One of the important ideas been used is called Oblivious Subspace Embedding (OSE), it is first proposed by Sarlos.[8]

For  , 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  .[9] For   and  , there are some algorithms with provable guarantees.[10][11]

Distance low-rank approximation problem

edit

Let   and   be two point sets in an arbitrary metric space. Let   represent the   matrix where  . 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,[12][13] one can study low rank approximation of such matrices.

Distributed/Streaming low-rank approximation problem

edit

The low-rank approximation problems in the distributed and streaming setting has been considered in.[14]

Image and kernel representations of the rank constraints

edit

Using the equivalences

 

and

 

the weighted low-rank approximation problem becomes equivalent to the parameter optimization problems

 

and

 

where   is the identity matrix of size  .

Alternating projections algorithm

edit

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 (  or  ) with the other one fixed. Although simultaneous minimization over both   and   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   (or  ) 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
endfor

Variable projections algorithm

edit

The alternating projections algorithm exploits the fact that the low rank approximation problem, parameterized in the image form, is bilinear in the variables   or  . The bilinear nature of the problem is effectively used in an alternative approach, called variable projections.[15]

Consider again the weighted low rank approximation problem, parameterized in the image form. Minimization with respect to the   variable (a linear least squares problem) leads to the closed form expression of the approximation error as a function of  

 

The original problem is therefore equivalent to the nonlinear least squares problem of minimizing   with respect to  . 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

edit

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   is linear, like we require all elements to be nonnegative, the problem is called structured low rank approximation.[16] 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  . 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,[17] 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.

See also

edit

References

edit
  1. ^ E. Schmidt, Zur Theorie der linearen und nichtlinearen Integralgleichungen, Math. Annalen 63 (1907), 433-476. doi:10.1007/BF01449770
  2. ^ C. Eckart, G. Young, The approximation of one matrix by another of lower rank. Psychometrika, Volume 1, 1936, Pages 211–8. doi:10.1007/BF02288367
  3. ^ L. Mirsky, Symmetric gauge functions and unitarily invariant norms, Q.J. Math. 11 (1960), 50-59. doi:10.1093/qmath/11.1.50
  4. ^ Srebro, Nathan; Jaakkola, Tommi (2003). Weighted Low-Rank Approximations (PDF). ICML'03.
  5. ^ Razenshteyn, Ilya; Song, Zhao; Woodruff, David P. (2016). Weighted Low Rank Approximations with Provable Guarantees. STOC '16 Proceedings of the forty-eighth annual ACM symposium on Theory of Computing.
  6. ^ Clarkson, Kenneth L.; Woodruff, David P. (2013). Low Rank Approximation and Regression in Input Sparsity Time. STOC '13 Proceedings of the forty-fifth annual ACM symposium on Theory of Computing. arXiv:1207.6365.
  7. ^ Nelson, Jelani; Nguyen, Huy L. (2013). OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings. FOCS '13. arXiv:1211.1002.
  8. ^ Sarlos, Tamas (2006). Improved approximation algorithms for large matrices via random projections. FOCS'06.
  9. ^ Song, Zhao; Woodruff, David P.; Zhong, Peilin (2017). Low Rank Approximation with Entrywise L1-Norm Error. STOC '17 Proceedings of the forty-ninth annual ACM symposium on Theory of Computing. arXiv:1611.00898.
  10. ^ Bringmann, Karl; Kolev, Pavel; Woodruff, David P. (2017). Approximation Algorithms for L0-Low Rank Approximation. NIPS'17. arXiv:1710.11253.
  11. ^ Chierichetti, Flavio; Gollapudi, Sreenivas; Kumar, Ravi; Lattanzi, Silvio; Panigrahy, Rina; Woodruff, David P. (2017). Algorithms for Lp Low-Rank Approximation. ICML'17. arXiv:1705.06730.
  12. ^ Bakshi, Ainesh L.; Woodruff, David P. (2018). Sublinear Time Low-Rank Approximation of Distance Matrices. NeurIPS. arXiv:1809.06986.
  13. ^ Indyk, Piotr; Vakilian, Ali; Wagner, Tal; Woodruff, David P. (2019). Sample-Optimal Low-Rank Approximation of Distance Matrices. COLT.
  14. ^ Boutsidis, Christos; Woodruff, David P.; Zhong, Peilin (2016). Optimal Principal Component Analysis in Distributed and Streaming Models. STOC. arXiv:1504.06729.
  15. ^ G. Golub and V. Pereyra, Separable nonlinear least squares: the variable projection method and its applications, Institute of Physics, Inverse Problems, Volume 19, 2003, Pages 1-26.
  16. ^ Chu, Moody T.; Funderlic, Robert E.; Plemmons, Robert J. (2003). "structured low-rank approximation". Linear Algebra and Its Applications. 366: 157–172. doi:10.1016/S0024-3795(02)00505-0.
  17. ^ "A General System for Heuristic Solution of Convex Problems over Nonconvex Sets" (PDF).
  • M. T. Chu, R. E. Funderlic, R. J. Plemmons, Structured low-rank approximation, Linear Algebra and its Applications, Volume 366, 1 June 2003, Pages 157–172 doi:10.1016/S0024-3795(02)00505-0
edit