Returns a manifold struct to optimize over k symmetric matrices of size n function M = symmetricfactory(n) function M = symmetricfactory(n, k) Returns M, a structure describing the Euclidean space of n-by-n symmetric matrices equipped with the standard Frobenius distance and associated trace inner product, as a manifold for Manopt. By default, k = 1. If k > 1, points and vectors are stored in 3D matrices X of size nxnxk such that each slice X(:, :, i), for i = 1:k, is symmetric.
0001 function M = symmetricfactory(n, k) 0002 % Returns a manifold struct to optimize over k symmetric matrices of size n 0003 % 0004 % function M = symmetricfactory(n) 0005 % function M = symmetricfactory(n, k) 0006 % 0007 % Returns M, a structure describing the Euclidean space of n-by-n symmetric 0008 % matrices equipped with the standard Frobenius distance and associated 0009 % trace inner product, as a manifold for Manopt. 0010 % 0011 % By default, k = 1. If k > 1, points and vectors are stored in 3D matrices 0012 % X of size nxnxk such that each slice X(:, :, i), for i = 1:k, is 0013 % symmetric. 0014 0015 % This file is part of Manopt: www.manopt.org. 0016 % Original author: Nicolas Boumal, Jan. 22, 2014. 0017 % Contributors: 0018 % Change log: 0019 % 0020 % Jan. 25, 2017 (NB): 0021 % M.tangent = M.proj now, instead of being identity. This is notably 0022 % necessary so that checkgradient will pick up on gradients that do 0023 % not lie in the appropriate tangent space. 0024 0025 if ~exist('k', 'var') || isempty(k) 0026 k = 1; 0027 end 0028 0029 M.name = @() sprintf('(Symmetric matrices of size %d)^%d', n, k); 0030 0031 M.dim = @() k*n*(n+1)/2; 0032 0033 M.inner = @(x, d1, d2) d1(:).'*d2(:); 0034 0035 M.norm = @(x, d) norm(d(:), 'fro'); 0036 0037 M.dist = @(x, y) norm(x(:)-y(:), 'fro'); 0038 0039 M.typicaldist = @() sqrt(k)*n; 0040 0041 M.proj = @(x, d) multisym(d); 0042 0043 M.egrad2rgrad = M.proj; 0044 0045 M.ehess2rhess = @(x, eg, eh, d) M.proj(x, eh); 0046 0047 M.tangent = M.proj; 0048 0049 M.exp = @exp; 0050 function y = exp(x, d, t) 0051 if nargin == 3 0052 y = x + t*d; 0053 else 0054 y = x + d; 0055 end 0056 end 0057 0058 M.retr = M.exp; 0059 0060 M.log = @(x, y) y-x; 0061 0062 M.hash = @(x) ['z' hashmd5(x(:))]; 0063 0064 M.rand = @() multisym(randn(n, n, k)); 0065 0066 M.randvec = @randvec; 0067 function u = randvec(x) %#ok<INUSD> 0068 u = multisym(randn(n, n, k)); 0069 u = u / norm(u(:), 'fro'); 0070 end 0071 0072 M.lincomb = @matrixlincomb; 0073 0074 M.zerovec = @(x) zeros(n, n, k); 0075 0076 M.transp = @(x1, x2, d) d; 0077 0078 M.pairmean = @(x1, x2) .5*(x1+x2); 0079 0080 0081 % Elaborate list of indices of diagonal entries of an nxnxk matrix. 0082 single_diag_entries = (1:(n+1):n^2)'; 0083 all_diag_entries = bsxfun(@plus, single_diag_entries, n^2*(0:(k-1))); 0084 all_diag_entries = all_diag_entries(:); 0085 0086 % Likewise, elaborate list of indices of upper-triangular entries. 0087 single_upper_triangle = find(triu(ones(n), 1)); 0088 all_upper_triangle = bsxfun(@plus, single_upper_triangle, n^2*(0:(k-1))); 0089 all_upper_triangle = all_upper_triangle(:); 0090 0091 % To vectorize a matrix, we extract all diagonal entries, then all 0092 % upper-triangular entries, the latter being scaled by sqrt(2) to 0093 % ensure isometry, that is: given two tangent vectors U and V at a 0094 % point X, M.inner(X, U, V) is equal to u'*v, where u = M.vec(X, U) and 0095 % likewise for v. This construction has the advantage of providing a 0096 % vectorized representation of matrices that has the same length as the 0097 % intrinsic dimension of the space they live in. 0098 M.vec = @(x, u_mat) [u_mat(all_diag_entries) ; ... 0099 sqrt(2)*u_mat(all_upper_triangle)]; 0100 M.mat = @matricize; 0101 function u_mat = matricize(X, u_vec) %#ok<INUSL> 0102 u_mat = zeros(n, n, k); 0103 u_mat(all_upper_triangle) = u_vec((k*n+1):end) / sqrt(2); 0104 u_mat = u_mat + multitransp(u_mat); 0105 u_mat(all_diag_entries) = u_vec(1:(k*n)); 0106 end 0107 M.vecmatareisometries = @() true; 0108 0109 end 0110 0111 % Former, easier versions for vec / mat. They had the disadvantage of 0112 % giving vector representations of length k*n^2, instead of k*n*(n+1)/2. 0113 % M.vec = @(x, u_mat) u_mat(:); 0114 % M.mat = @(x, u_vec) reshape(u_vec, [m, n]); 0115 % M.vecmatareisometries = @() true;