0001 function M = spheresymmetricfactory(n)
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 M.name = @() sprintf('Sphere of symmetric matrices of size %d', n);
0031
0032 M.dim = @() n*(n+1)/2 - 1;
0033
0034 M.inner = @(x, d1, d2) d1(:).'*d2(:);
0035
0036 M.norm = @(x, d) norm(d, 'fro');
0037
0038 M.dist = @(x, y) real(2*asin(.5*norm(x - y, 'fro')));
0039
0040 M.typicaldist = @() pi;
0041
0042 M.proj = @proj;
0043 function xdot = proj(x, d)
0044 d = (d+d.')/2;
0045 xdot = d - x*(x(:).'*d(:));
0046 end
0047
0048 M.tangent = @proj;
0049
0050
0051
0052 M.egrad2rgrad = @proj;
0053
0054 M.ehess2rhess = @ehess2rhess;
0055 function rhess = ehess2rhess(x, egrad, ehess, u)
0056
0057
0058
0059 rhess = proj(x, ehess) - (x(:)'*egrad(:))*u;
0060 end
0061
0062 M.exp = @exponential;
0063
0064 M.retr = @retraction;
0065
0066 M.log = @logarithm;
0067 function v = logarithm(x1, x2)
0068 v = proj(x1, x2 - x1);
0069 di = M.dist(x1, x2);
0070
0071 if di > 1e-6
0072 nv = norm(v, 'fro');
0073 v = v * (di / nv);
0074 end
0075 end
0076
0077 M.hash = @(x) ['z' hashmd5(x(:))];
0078
0079 M.rand = @() random(n);
0080
0081 M.randvec = @(x) randomvec(n, x);
0082
0083 M.lincomb = @matrixlincomb;
0084
0085 M.zerovec = @(x) zeros(n);
0086
0087 M.transp = @(x1, x2, d) proj(x2, d);
0088
0089 M.pairmean = @pairmean;
0090 function y = pairmean(x1, x2)
0091 y = x1+x2;
0092 y = y / norm(y, 'fro');
0093 end
0094
0095
0096 M.vec = @(x, u_mat) u_mat(:);
0097 M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0098 M.vecmatareisometries = @() false;
0099
0100 end
0101
0102
0103 function y = exponential(x, d, t)
0104
0105 if nargin == 2
0106
0107 td = d;
0108 else
0109 td = t*d;
0110 end
0111
0112 nrm_td = norm(td, 'fro');
0113
0114 if nrm_td > 0
0115 y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
0116 else
0117 y = x;
0118 end
0119
0120 end
0121
0122
0123 function y = retraction(x, d, t)
0124
0125 if nargin == 2
0126 t = 1;
0127 end
0128
0129 y = x + t*d;
0130 y = y / norm(y, 'fro');
0131
0132 end
0133
0134
0135 function x = random(n)
0136
0137 x = randn(n);
0138 x = (x + x.')/2;
0139 x = x/norm(x, 'fro');
0140
0141 end
0142
0143
0144 function d = randomvec(n, x)
0145
0146 d = randn(n);
0147 d = (d + d.')/2;
0148 d = d - x*(x(:).'*d(:));
0149 d = d / norm(d, 'fro');
0150
0151 end