Home > manopt > manifolds > fixedrank > fixedrankfactory_3factors.m

# fixedrankfactory_3factors

## PURPOSE

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

## SYNOPSIS

function M = fixedrankfactory_3factors(m, n, k)

## DESCRIPTION

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

function M = fixedrankfactory_3factors(m, n, k)

The first-order geometry follows the balanced quotient geometry described
in the paper,
"Linear regression under fixed-rank constraints: a Riemannian approach",
G. Meyer, S. Bonnabel and R. Sepulchre, ICML 2011.

The second-order geometry follows from the paper
"Fixed-rank matrix factorizations and Riemannian low-rank optimization",
B. Mishra, R. Meyer, S. Bonnabel and R. Sepulchre,
Computational Statistics, 29(3 - 4), pp. 591 - 621, 2014.

A point X on the manifold is represented as a structure with three
fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal,
while the matrix S (kxk) is a symmetric positive definite full rank
matrix.

Tangent vectors are represented as a structure with three fields: L, S
and R.

For first-order geometry, please cite the Manopt paper as well as the research paper:
@InProceedings{meyer2011linear,
Title        = {Linear regression under fixed-rank constraints: a {R}iemannian approach},
Author       = {Meyer, G. and Bonnabel, S. and Sepulchre, R.},
Booktitle    = {{28th International Conference on Machine Learning}},
Year         = {2011},
Organization = {{ICML}}
}
For second-order geometry, please cite the Manopt paper as well as the research paper:
@Article{mishra2014fixedrank,
Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
Journal = {Computational Statistics},
Year    = {2014},
Number  = {3-4},
Pages   = {591--621},
Volume  = {29},
Doi     = {10.1007/s00180-013-0464-z}
}

See also fixedrankembeddedfactory fixedrankfactory_2factors fixedrankfactory_3factors_preconditioned

## CROSS-REFERENCE INFORMATION

This function calls:
• stiefelfactory Returns a manifold structure to optimize over orthonormal matrices.
• hashmd5 Computes the MD5 hash of input data.
• lincomb Computes a linear combination of tangent vectors in the Manopt framework.
• lyapunov_symmetric Solves AX + XA = C when A = A', as a pseudo-inverse.
This function is called by:

## SOURCE CODE

0001 function M = fixedrankfactory_3factors(m, n, k)
0002 % Manifold of m-by-n matrices of rank k with polar quotient geometry.
0003 %
0004 % function M = fixedrankfactory_3factors(m, n, k)
0005 %
0006 % The first-order geometry follows the balanced quotient geometry described
0007 % in the paper,
0008 % "Linear regression under fixed-rank constraints: a Riemannian approach",
0009 % G. Meyer, S. Bonnabel and R. Sepulchre, ICML 2011.
0010 %
0012 %
0013 % The second-order geometry follows from the paper
0014 % "Fixed-rank matrix factorizations and Riemannian low-rank optimization",
0015 % B. Mishra, R. Meyer, S. Bonnabel and R. Sepulchre,
0016 % Computational Statistics, 29(3 - 4), pp. 591 - 621, 2014.
0017 %
0018 % A point X on the manifold is represented as a structure with three
0019 % fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal,
0020 % while the matrix S (kxk) is a symmetric positive definite full rank
0021 % matrix.
0022 %
0023 % Tangent vectors are represented as a structure with three fields: L, S
0024 % and R.
0025 %
0026 %
0027 % For first-order geometry, please cite the Manopt paper as well as the research paper:
0028 %     @InProceedings{meyer2011linear,
0029 %       Title        = {Linear regression under fixed-rank constraints: a {R}iemannian approach},
0030 %       Author       = {Meyer, G. and Bonnabel, S. and Sepulchre, R.},
0031 %       Booktitle    = {{28th International Conference on Machine Learning}},
0032 %       Year         = {2011},
0033 %       Organization = {{ICML}}
0034 %     }
0035 % For second-order geometry, please cite the Manopt paper as well as the research paper:
0036 %     @Article{mishra2014fixedrank,
0037 %       Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
0038 %       Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
0039 %       Journal = {Computational Statistics},
0040 %       Year    = {2014},
0041 %       Number  = {3-4},
0042 %       Pages   = {591--621},
0043 %       Volume  = {29},
0044 %       Doi     = {10.1007/s00180-013-0464-z}
0045 %     }
0046 %
0047 %
0049
0050 % This file is part of Manopt: www.manopt.org.
0051 % Original author: Bamdev Mishra, Dec. 30, 2012.
0052 % Contributors:
0053 % Change log:
0054 %
0055 %    Apr. 18, 2018 (NB):
0056 %        Removed dependency on lyap.
0057 %
0058 %    Sep.  6, 2018 (NB):
0059 %        Removed M.exp() as it was not implemented.
0060
0061     M.name = @() sprintf('LSR'' quotient manifold of %dx%d matrices of rank %d', m, n, k);
0062
0063     M.dim = @() (m+n-k)*k;
0064
0065     % Choice of the metric on the orthnormal space is motivated by the symmetry present in the
0066     % space. The metric on the positive definite space is its natural metric.
0067     M.inner = @(X, eta, zeta) eta.L(:).'*zeta.L(:) + eta.R(:).'*zeta.R(:) ...
0068         + trace( (X.S\eta.S) * (X.S\zeta.S) );
0069
0070     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0071
0072     M.dist = @(x, y) error('fixedrankfactory_3factors.dist not implemented yet.');
0073
0074     M.typicaldist = @() 10*k;
0075
0076     skew = @(X) .5*(X-X');
0077     symm = @(X) .5*(X+X');
0078     stiefel_proj = @(L, H) H - L*symm(L'*H);
0079
0085     end
0086
0087
0088     M.ehess2rhess = @ehess2rhess;
0089     function Hess = ehess2rhess(X, egrad, ehess, eta)
0090
0091         % Riemannian gradient for the factor S.
0093
0094         % Directional derivatives of the Riemannian gradient.
0095         Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L);
0096         Hess.L = stiefel_proj(X.L, Hess.L);
0097
0098         Hess.R = ehess.R - eta.R*symm(X.R'*egrad.R);
0099         Hess.R = stiefel_proj(X.R, Hess.R);
0100
0101         Hess.S = X.S*symm(ehess.S)*X.S +  2*symm(eta.S*symm(egrad.S)*X.S);
0102
0103         % Correction factor for the non-constant metric on the factor S.
0104         Hess.S = Hess.S - symm(eta.S*(X.S\rgrad.S));
0105
0106         % Projection onto the horizontal space.
0107         Hess = M.proj(X, Hess);
0108     end
0109
0110
0111     M.proj = @projection;
0112     function etaproj = projection(X, eta)
0113         % First, projection onto the tangent space of the total space.
0114         eta.L = stiefel_proj(X.L, eta.L);
0115         eta.R = stiefel_proj(X.R, eta.R);
0116         eta.S = symm(eta.S);
0117
0118         % Then, projection onto the horizontal space.
0119         SS = X.S*X.S;
0120         AS = X.S*(skew(X.L'*eta.L) + skew(X.R'*eta.R) - 2*skew(X.S\eta.S))*X.S;
0121         omega = lyapunov_symmetric(SS, AS);
0122
0123         etaproj.L = eta.L - X.L*omega;
0124         etaproj.S = eta.S - (X.S*omega - omega*X.S);
0125         etaproj.R = eta.R - X.R*omega;
0126     end
0127
0128     M.tangent = M.proj;
0129     M.tangent2ambient = @(X, eta) eta;
0130
0131     M.retr = @retraction;
0132     function Y = retraction(X, eta, t)
0133         if nargin < 3
0134             t = 1.0;
0135         end
0136
0137         L = chol(X.S);
0138         Y.S = L'*expm(L'\(t*eta.S)/L)*L;
0139         Y.L = uf(X.L + t*eta.L);
0140         Y.R = uf(X.R + t*eta.R);
0141     end
0142
0143
0144     M.hash = @(X) ['z' hashmd5([X.L(:) ; X.S(:) ; X.R(:)])];
0145
0146     M.rand = @random;
0147     % Factors L and R are on Stiefel manifolds, hence we reuse
0148     % their random generators.
0149     stiefelm = stiefelfactory(m, k);
0150     stiefeln = stiefelfactory(n, k);
0151     function X = random()
0152         X.L = stiefelm.rand();
0153         X.R = stiefeln.rand();
0154         X.S = diag(1+rand(k, 1));
0155     end
0156
0157     M.randvec = @randomvec;
0158     function eta = randomvec(X)
0159         % A random vector on the horizontal space.
0160         eta.L = randn(m, k);
0161         eta.R = randn(n, k);
0162         eta.S = randn(k, k);
0163         eta = projection(X, eta);
0164         nrm = M.norm(X, eta);
0165         eta.L = eta.L / nrm;
0166         eta.R = eta.R / nrm;
0167         eta.S = eta.S / nrm;
0168     end
0169
0170     M.lincomb = @lincomb;
0171
0172     M.zerovec = @(X) struct('L', zeros(m, k), 'S', zeros(k, k), ...
0173         'R', zeros(n, k));
0174
0175     M.transp = @(x1, x2, d) projection(x2, d);
0176
0177     % vec and mat are not isometries, because of the scaled inner metric.
0178     M.vec = @(X, U) [U.L(:) ; U.S(:); U.R(:)];
0179     M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0180         'S', reshape(u((m*k+1): m*k + k*k), k, k), ...
0181         'R', reshape(u((m*k+ k*k + 1):end), n, k));
0182     M.vecmatareisometries = @() false;
0183
0184 end
0185
0186 % Linear combination of tangent vectors.
0187 function d = lincomb(x, a1, d1, a2, d2) %#ok<INLSL>
0188
0189     if nargin == 3
0190         d.L = a1*d1.L;
0191         d.R = a1*d1.R;
0192         d.S = a1*d1.S;
0193     elseif nargin == 5
0194         d.L = a1*d1.L + a2*d2.L;
0195         d.R = a1*d1.R + a2*d2.R;
0196         d.S = a1*d1.S + a2*d2.S;
0197     else
0199     end
0200
0201 end
0202
0203 function A = uf(A)
0204     [L, unused, R] = svd(A, 0); %#ok
0205     A = L*R';
0206 end

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005