0001 function M = fixedrankfactory_2factors_preconditioned(m, n, k)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
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
0045
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
0057
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));
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
0070 M.egrad2rgrad = @egrad2rgrad;
0071 function rgrad = egrad2rgrad(X, egrad)
0072 X = prepare(X);
0073
0074
0075 rgrad.L = egrad.L/X.RtR;
0076 rgrad.R = egrad.R/X.LtL;
0077 end
0078
0079 M.ehess2rhess = @ehess2rhess;
0080 function Hess = ehess2rhess(X, egrad, ehess, eta)
0081 X = prepare(X);
0082
0083
0084 rgrad = egrad2rgrad(X, egrad);
0085
0086
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
0091 Hess.L = Hess.L + rgrad.L*(symm(eta.R'*X.R)/X.RtR) + eta.L*(symm(rgrad.R'*X.R)/X.RtR) - X.L*(symm(eta.R'*rgrad.R)/X.RtR);
0092 Hess.R = Hess.R + rgrad.R*(symm(eta.L'*X.L)/X.LtL) + eta.R*(symm(rgrad.L'*X.L)/X.LtL) - X.R*(symm(eta.L'*rgrad.L)/X.LtL);
0093
0094
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
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
0121
0122
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
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
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
0168 symm = @(M) .5*(M+M');
0169 end
0170
0171
0172 function d = lincomb(x, a1, d1, a2, d2)
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
0181 error('Bad use of fixedrankfactory_2factors_preconditioned.lincomb.');
0182 end
0183
0184 end
0185
0186
0187
0188
0189