Home > manopt > manifolds > fixedrank > fixedrankfactory_2factors_preconditioned.m

# fixedrankfactory_2factors_preconditioned

## PURPOSE

Manifold of m-by-n matrices of rank k with two factor quotient geometry.

## SYNOPSIS

function M = fixedrankfactory_2factors_preconditioned(m, n, k)

## DESCRIPTION

Manifold of m-by-n matrices of rank k with two factor quotient geometry.

function M = fixedrankfactory_2factors_preconditioned(m, n, k)

This geometry is tuned to least-squares problems such as low-rank matrix
completion with ell-2 loss.

A point X on the manifold is represented as a structure with two
fields: L and R. The matrices L (m-by-k) and R (n-by-k) are
full column-rank matrices such that X = L*R'.

Tangent vectors are represented as a structure with two fields: L, R.

Please cite the Manopt paper as well as the research paper:
@Techreport{mishra2012optimized,
Title   = {A {R}iemannian geometry for low-rank matrix completion},
Author  = {Mishra, B. and Adithya Apuroop, K. and Sepulchre, R.},
Journal = {Arxiv preprint arXiv:1211.1550},
Year    = {2012}
}

## CROSS-REFERENCE INFORMATION

This function calls:
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• hashmd5 Computes the MD5 hash of input data.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
This function is called by:

## SOURCE CODE

0001 function M = fixedrankfactory_2factors_preconditioned(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with two factor quotient geometry.
0003 %
0004 % function M = fixedrankfactory_2factors_preconditioned(m, n, k)
0005 %
0006 % This geometry is tuned to least-squares problems such as low-rank matrix
0007 % completion with ell-2 loss.
0008 %
0009 % A point X on the manifold is represented as a structure with two
0010 % fields: L and R. The matrices L (m-by-k) and R (n-by-k) are
0011 % full column-rank matrices such that X = L*R'.
0012 %
0013 % Tangent vectors are represented as a structure with two fields: L, R.
0014 %
0015 % Please cite the Manopt paper as well as the research paper:
0016 %     @Techreport{mishra2012optimized,
0017 %       Title   = {A {R}iemannian geometry for low-rank matrix completion},
0018 %       Author  = {Mishra, B. and Adithya Apuroop, K. and Sepulchre, R.},
0019 %       Journal = {Arxiv preprint arXiv:1211.1550},
0020 %       Year    = {2012}
0021 %     }
0022 %
0023 %
0025
0026 % This file is part of Manopt: www.manopt.org.
0027 % Original author: Bamdev Mishra, Dec. 30, 2012.
0028 % Contributors:
0029 % Change log:
0030 %
0031 %    Apr. 4, 2015 (BM):
0032 %        Cosmetic changes including avoiding storing the inverse of a
0033 %        k-by-k matrix.
0034 %
0035 %    Sep. 6, 2018 (NB):
0036 %        Removed M.exp() as it was not implemented.
0037
0038
0039     M.name = @() sprintf('LR''(tuned to least square problems) quotient manifold of %dx%d matrices of rank %d', m, n, k);
0040
0041     M.dim = @() (m+n-k)*k;
0042
0043
0044     % Some precomputations at the point X to be used in the inner product
0045     % (and pretty much everywhere else).
0046     function X = prepare(X)
0047         if ~all(isfield(X,{'LtL','RtR'}))
0048             L = X.L;
0049             R = X.R;
0050             X.LtL = L'*L;
0051             X.RtR = R'*R;
0052         end
0053     end
0054
0055
0056     % The choice of metric is motivated by symmetry and
0057     % tuned to least-squares cost function.
0058     M.inner = @iproduct;
0059     function ip = iproduct(X, eta, zeta)
0060         X = prepare(X);
0061         ip = trace(X.RtR*(eta.L'*zeta.L)) + trace(X.LtL*(eta.R'*zeta.R)); % Scaled metric
0062     end
0063
0064     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0065
0066     M.dist = @(x, y) error('fixedrankfactory_2factors_preconditioned.dist not implemented yet.');
0067
0068     M.typicaldist = @() 10*k;
0069
0072         X = prepare(X);
0073
0077     end
0078
0079     M.ehess2rhess = @ehess2rhess;
0080     function Hess = ehess2rhess(X, egrad, ehess, eta)
0081         X = prepare(X);
0082
0085
0086         % Directional derivative of the Riemannian gradient.
0087         Hess.L = ehess.L/X.RtR - 2*egrad.L*(X.RtR \ (symm(eta.R'*X.R) / X.RtR));
0088         Hess.R = ehess.R/X.LtL - 2*egrad.R*(X.LtL \ (symm(eta.L'*X.L) / X.LtL));
0089
0090         % We still need a correction factor for the non-constant metric.
0093
0094         % Project on the horizontal space.
0095         Hess = M.proj(X, Hess);
0096     end
0097
0098     M.proj = @projection;
0099     function etaproj = projection(X, eta)
0100         X = prepare(X);
0101
0102         % Projection onto the horizontal space.
0103         Lambda = 0.5*((eta.R'*X.R)/X.RtR  -   X.LtL\(X.L'*eta.L));
0104         etaproj.L = eta.L + X.L*Lambda;
0105         etaproj.R = eta.R - X.R*Lambda';
0106     end
0107
0108     M.tangent = M.proj;
0109
0110     M.tangent2ambient = @(X, eta) eta;
0111
0112     M.retr = @retraction;
0113     function Y = retraction(X, eta, t)
0114         if nargin < 3
0115             t = 1.0;
0116         end
0117         Y.L = X.L + t*eta.L;
0118         Y.R = X.R + t*eta.R;
0119
0120         % Numerical conditioning step: a simpler version.
0121         % We need to ensure that L and R are do not have very relative
0122         % skewed norms.
0123
0124         scaling = norm(X.L, 'fro')/norm(X.R, 'fro');
0125         scaling = sqrt(scaling);
0126         Y.L = Y.L / scaling;
0127         Y.R = Y.R * scaling;
0128
0129         % These are reused in the computations of gradient and Hessian.
0130         Y = prepare(Y);
0131     end
0132
0133
0134     M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];
0135
0136     M.rand = @random;
0137
0138     function X = random()
0139         X.L = randn(m, k);
0140         X.R = randn(n, k);
0141     end
0142
0143     M.randvec = @randomvec;
0144     function eta = randomvec(X)
0145         eta.L = randn(m, k);
0146         eta.R = randn(n, k);
0147         eta = projection(X, eta);
0148         nrm = M.norm(X, eta);
0149         eta.L = eta.L / nrm;
0150         eta.R = eta.R / nrm;
0151     end
0152
0153     M.lincomb = @lincomb;
0154
0155     M.zerovec = @(X) struct('L', zeros(m, k),'R', zeros(n, k));
0156
0157     M.transp = @(x1, x2, d) projection(x2, d);
0158
0159     % vec and mat are not isometries, because of the scaled inner metric.
0160     M.vec = @(X, U) [U.L(:) ; U.R(:)];
0161
0162     M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0163         'R', reshape(u((m*k+1):end), n, k));
0164
0165     M.vecmatareisometries = @() false;
0166
0167     % Auxiliary functions
0168     symm = @(M) .5*(M+M');
0169 end
0170
0171 % Linear combination of tangent vectors.
0172 function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
0173
0174     if nargin == 3
0175         d.L = a1*d1.L;
0176         d.R = a1*d1.R;
0177     elseif nargin == 5
0178         d.L = a1*d1.L + a2*d2.L;
0179         d.R = a1*d1.R + a2*d2.R;
0180     else