Low-rank approximation

From Wikipedia, the free encyclopedia
Jump to: navigation, search

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:



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


Basic low-rank approximation problem[edit]

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.[4] Let

be the singular value decomposition of and 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 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 denote the best rank approximation to and , respectively. Then, for any

Since , when and we conclude that for


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.

Image and kernel representations of the rank constraints[edit]

Using the equivalences


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


where is the identity matrix of size .

Alternating projections algorithm[edit]

The image representation of the rank constraint suggests a parameter optimization methods, 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

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.[5]

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 can find tons of real applications, including to recover a good solution from a 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.[6] And the more general form is named as convex-restricted low rank approximation.[7]

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,[8] 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]


  1. ^ I. Markovsky, Structured low-rank approximation and its applications, Automatica, Volume 44, Issue 4, April 2008, Pages 891–909. doi:10.1016/j.automatica.2007.09.011
  2. ^ I. Markovsky, J. C. Willems, S. Van Huffel, B. De Moor, and R. Pintelon, Application of structured total least squares for system identification and model reduction. IEEE Transactions on Automatic Control, Volume 50, Number 10, 2005, pages 1490–1500.
  3. ^ I. Markovsky, Low-Rank Approximation: Algorithms, Implementation, Applications, Springer, 2012, ISBN 978-1-4471-2226-5
  4. ^ 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
  5. ^ 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.
  6. ^ "structured low-rank approximation". Linear Algebra and its Applications. 366: 157–172. doi:10.1016/S0024-3795(02)00505-0. 
  7. ^ Zhang, Ying (2015). "restricted low-rank approximation". arXiv:1512.01748Freely accessible [cs.NA]. 
  8. ^ "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

External links[edit]