0001 function M = fixedrankfactory_2factors_subspace_projection(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
0040
0041
0042
0043
0044
0045
0046 M.name = @() sprintf('LR'' quotient manifold of %dx%d matrices of rank %d', m, n, k);
0047
0048 M.dim = @() (m+n-k)*k;
0049
0050
0051
0052 function X = prepare(X)
0053 if ~all(isfield(X,{'RtR'}) == 1)
0054 X.RtR = X.R'*X.R;
0055 end
0056 end
0057
0058
0059
0060 M.inner = @iproduct;
0061 function ip = iproduct(X, eta, zeta)
0062 X = prepare(X);
0063
0064 ip = eta.L(:).'*zeta.L(:) + trace(X.RtR\(eta.R'*zeta.R));
0065 end
0066
0067 M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0068
0069 M.dist = @(x, y) error('fixedrankfactory_2factors_subspace_projection.dist not implemented yet.');
0070
0071 M.typicaldist = @() 10*k;
0072
0073 skew = @(X) .5*(X-X');
0074 symm = @(X) .5*(X+X');
0075 stiefel_proj = @(L, H) H - L*symm(L'*H);
0076
0077 M.egrad2rgrad = @egrad2rgrad;
0078 function rgrad = egrad2rgrad(X, egrad)
0079 X = prepare(X);
0080
0081 rgrad.L = stiefel_proj(X.L, egrad.L);
0082 rgrad.R = egrad.R*X.RtR;
0083 end
0084
0085
0086 M.ehess2rhess = @ehess2rhess;
0087 function Hess = ehess2rhess(X, egrad, ehess, eta)
0088 X = prepare(X);
0089
0090
0091 rgrad = egrad2rgrad(X, egrad);
0092
0093
0094 Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L);
0095 Hess.L = stiefel_proj(X.L, Hess.L);
0096
0097 Hess.R = ehess.R*X.RtR + 2*egrad.R*symm(eta.R'*X.R);
0098
0099
0100 Hess.R = Hess.R - rgrad.R*(X.RtR\(symm(X.R'*eta.R))) - eta.R*(X.RtR\(symm(X.R'*rgrad.R))) + X.R*(X.RtR\(symm(eta.R'*rgrad.R)));
0101
0102
0103 Hess = M.proj(X, Hess);
0104 end
0105
0106
0107 M.proj = @projection;
0108 function etaproj = projection(X, eta)
0109 X = prepare(X);
0110
0111 eta.L = stiefel_proj(X.L, eta.L);
0112 SS = X.RtR;
0113 AS1 = 2*X.RtR*skew(X.L'*eta.L)*X.RtR;
0114 AS2 = 2*skew(X.RtR*(X.R'*eta.R));
0115 AS = skew(AS1 + AS2);
0116
0117 Omega = nested_sylvester(SS, AS);
0118 etaproj.L = eta.L - X.L*Omega;
0119 etaproj.R = eta.R - X.R*Omega;
0120 end
0121
0122 M.tangent = M.proj;
0123 M.tangent2ambient = @(X, eta) eta;
0124
0125 M.retr = @retraction;
0126 function Y = retraction(X, eta, t)
0127 if nargin < 3
0128 t = 1.0;
0129 end
0130 Y.L = uf(X.L + t*eta.L);
0131 Y.R = X.R + t*eta.R;
0132
0133
0134 Y = prepare(Y);
0135 end
0136
0137
0138 M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];
0139
0140 M.rand = @random;
0141
0142
0143 stiefelm = stiefelfactory(m, k);
0144 function X = random()
0145 X.L = stiefelm.rand();
0146 X.R = randn(n, k);
0147 end
0148
0149 M.randvec = @randomvec;
0150 function eta = randomvec(X)
0151 eta.L = randn(m, k);
0152 eta.R = randn(n, k);
0153 eta = projection(X, eta);
0154 nrm = M.norm(X, eta);
0155 eta.L = eta.L / nrm;
0156 eta.R = eta.R / nrm;
0157 end
0158
0159 M.lincomb = @lincomb;
0160
0161 M.zerovec = @(X) struct('L', zeros(m, k),...
0162 'R', zeros(n, k));
0163
0164 M.transp = @(x1, x2, d) projection(x2, d);
0165
0166
0167 M.vec = @(X, U) [U.L(:) ; U.R(:)];
0168 M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
0169 'R', reshape(u((m*k+1):end), n, k));
0170 M.vecmatareisometries = @() false;
0171
0172
0173 end
0174
0175
0176 function d = lincomb(x, a1, d1, a2, d2)
0177
0178 if nargin == 3
0179 d.L = a1*d1.L;
0180 d.R = a1*d1.R;
0181 elseif nargin == 5
0182 d.L = a1*d1.L + a2*d2.L;
0183 d.R = a1*d1.R + a2*d2.R;
0184 else
0185 error('Bad use of fixedrankfactory_2factors_subspace_projection.lincomb.');
0186 end
0187
0188 end
0189
0190 function A = uf(A)
0191 [L, unused, R] = svd(A, 0);
0192 A = L*R';
0193 end
0194
0195 function omega = nested_sylvester(sym_mat, asym_mat)
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207 [V, lambda] = eig(sym_mat, 'vector');
0208 X = lyapunov_symmetric_eig(V, lambda, asym_mat);
0209 omega = lyapunov_symmetric_eig(V, lambda, X);
0210
0211 end