Home > manopt > manifolds > euclidean > symmetricfactory.m

symmetricfactory

PURPOSE ^

Returns a manifold struct to optimize over k symmetric matrices of size n

SYNOPSIS ^

function M = symmetricfactory(n, k)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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     if ~exist('k', 'var') || isempty(k)
0021         k = 1;
0022     end
0023 
0024     M.name = @() sprintf('(Symmetric matrices of size %d)^%d', n, k);
0025     
0026     M.dim = @() k*n*(n+1)/2;
0027     
0028     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0029     
0030     M.norm = @(x, d) norm(d(:), 'fro');
0031     
0032     M.dist = @(x, y) norm(x(:)-y(:), 'fro');
0033     
0034     M.typicaldist = @() sqrt(k)*n;
0035     
0036     M.proj = @(x, d) multisym(d);
0037     
0038     M.egrad2rgrad = M.proj;
0039     
0040     M.ehess2rhess = @(x, eg, eh, d) M.proj(x, eh);
0041     
0042     M.tangent = @(x, d) d;
0043     
0044     M.exp = @exp;
0045     function y = exp(x, d, t)
0046         if nargin == 3
0047             y = x + t*d;
0048         else
0049             y = x + d;
0050         end
0051     end
0052     
0053     M.retr = M.exp;
0054     
0055     M.log = @(x, y) y-x;
0056 
0057     M.hash = @(x) ['z' hashmd5(x(:))];
0058     
0059     M.rand = @() multisym(randn(n, n, k));
0060     
0061     M.randvec = @randvec;
0062     function u = randvec(x) %#ok<INUSD>
0063         u = multisym(randn(n, n, k));
0064         u = u / norm(u(:), 'fro');
0065     end
0066     
0067     M.lincomb = @matrixlincomb;
0068     
0069     M.zerovec = @(x) zeros(n, n, k);
0070     
0071     M.transp = @(x1, x2, d) d;
0072     
0073     M.pairmean = @(x1, x2) .5*(x1+x2);
0074     
0075     
0076     % Elaborate list of indices of diagonal entries of an nxnxk matrix.
0077     single_diag_entries = (1:(n+1):n^2)';
0078     all_diag_entries = bsxfun(@plus, single_diag_entries, n^2*(0:(k-1)));
0079     all_diag_entries = all_diag_entries(:);
0080     
0081     % Likewise, elaborate list of indices of upper-triangular entries.
0082     single_upper_triangle = find(triu(ones(n), 1));
0083     all_upper_triangle = bsxfun(@plus, single_upper_triangle, n^2*(0:(k-1)));
0084     all_upper_triangle = all_upper_triangle(:);
0085     
0086     % To vectorize a matrix, we extract all diagonal entries, then all
0087     % upper-triangular entries, the latter being scaled by sqrt(2) to
0088     % ensure isometry, that is: given two tangent vectors U and V at a
0089     % point X, M.inner(X, U, V) is equal to u'*v, where u = M.vec(X, U) and
0090     % likewise for v. This construction has the advantage of providing a
0091     % vectorized representation of matrices that has the same length as the
0092     % intrinsic dimension of the space they live in.
0093     M.vec = @(x, u_mat) [u_mat(all_diag_entries) ; ...
0094                          sqrt(2)*u_mat(all_upper_triangle)];
0095     M.mat = @matricize;
0096     function u_mat = matricize(X, u_vec) %#ok<INUSL>
0097         u_mat = zeros(n, n, k);
0098         u_mat(all_upper_triangle) = u_vec((k*n+1):end) / sqrt(2);
0099         u_mat = u_mat + multitransp(u_mat);
0100         u_mat(all_diag_entries) = u_vec(1:(k*n));
0101     end
0102     M.vecmatareisometries = @() true;
0103 
0104 end
0105 
0106 % Former, easier versions for vec / mat. They had the disadvantage of
0107 % giving vector representations of length k*n^2, instead of k*n*(n+1).
0108 % M.vec = @(x, u_mat) u_mat(:);
0109 % M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
0110 % M.vecmatareisometries = @() true;

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005