Home > manopt > manifolds > rotations > rotationsfactory.m

rotationsfactory

PURPOSE ^

Returns a manifold structure to optimize over rotation matrices.

SYNOPSIS ^

function M = rotationsfactory(n, k)

DESCRIPTION ^

 Returns a manifold structure to optimize over rotation matrices.
 
 function M = rotationsfactory(n)
 function M = rotationsfactory(n, k)

 Special orthogonal group (the manifold of rotations): deals with matrices
 R of size n x n x k (or n x n if k = 1, which is the default) such that
 each n x n matrix is orthogonal, with determinant 1, i.e., X'*X = eye(n)
 if k = 1, or X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1.

 This is a description of SO(n)^k with the induced metric from the
 embedding space (R^nxn)^k, i.e., this manifold is a Riemannian
 submanifold of (R^nxn)^k endowed with the usual trace inner product.

 Tangent vectors are represented in the Lie algebra, i.e., as skew
 symmetric matrices. Use the function M.tangent2ambient(X, H) to switch
 from the Lie algebra representation to the embedding space
 representation. This is often necessary when defining
 problem.ehess(X, H).

 By default, the retraction is only a first-order approximation of the
 exponential. To force the use of a second-order approximation, call
 M.retr = M.retr2 after creating M. This switches from a QR-based
 computation to an SVD-based computation.

 By default, k = 1.

 See also: stiefelfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = rotationsfactory(n, k)
0002 % Returns a manifold structure to optimize over rotation matrices.
0003 %
0004 % function M = rotationsfactory(n)
0005 % function M = rotationsfactory(n, k)
0006 %
0007 % Special orthogonal group (the manifold of rotations): deals with matrices
0008 % R of size n x n x k (or n x n if k = 1, which is the default) such that
0009 % each n x n matrix is orthogonal, with determinant 1, i.e., X'*X = eye(n)
0010 % if k = 1, or X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1.
0011 %
0012 % This is a description of SO(n)^k with the induced metric from the
0013 % embedding space (R^nxn)^k, i.e., this manifold is a Riemannian
0014 % submanifold of (R^nxn)^k endowed with the usual trace inner product.
0015 %
0016 % Tangent vectors are represented in the Lie algebra, i.e., as skew
0017 % symmetric matrices. Use the function M.tangent2ambient(X, H) to switch
0018 % from the Lie algebra representation to the embedding space
0019 % representation. This is often necessary when defining
0020 % problem.ehess(X, H).
0021 %
0022 % By default, the retraction is only a first-order approximation of the
0023 % exponential. To force the use of a second-order approximation, call
0024 % M.retr = M.retr2 after creating M. This switches from a QR-based
0025 % computation to an SVD-based computation.
0026 %
0027 % By default, k = 1.
0028 %
0029 % See also: stiefelfactory
0030 
0031 % This file is part of Manopt: www.manopt.org.
0032 % Original author: Nicolas Boumal, Dec. 30, 2012.
0033 % Contributors:
0034 % Change log:
0035 %   Jan. 31, 2013 (NB)
0036 %       Added egrad2rgrad and ehess2rhess
0037 %   Oct. 21, 2016 (NB)
0038 %       Added M.retr2: a second-order retraction based on SVD.
0039 
0040     
0041     if ~exist('k', 'var') || isempty(k)
0042         k = 1;
0043     end
0044     
0045     if k == 1
0046         M.name = @() sprintf('Rotations manifold SO(%d)', n);
0047     elseif k > 1
0048         M.name = @() sprintf('Product rotations manifold SO(%d)^%d', n, k);
0049     else
0050         error('k must be an integer no less than 1.');
0051     end
0052     
0053     M.dim = @() k*nchoosek(n, 2);
0054     
0055     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0056     
0057     M.norm = @(x, d) norm(d(:));
0058     
0059     M.typicaldist = @() pi*sqrt(n*k);
0060     
0061     M.proj = @(X, H) multiskew(multiprod(multitransp(X), H));
0062     
0063     M.tangent = @(X, H) multiskew(H);
0064     
0065     M.tangent2ambient = @(X, U) multiprod(X, U);
0066     
0067     M.egrad2rgrad = M.proj;
0068     
0069     M.ehess2rhess = @ehess2rhess;
0070     function Rhess = ehess2rhess(X, Egrad, Ehess, H)
0071         % Reminder : H contains skew-symmeric matrices. The actual
0072         % direction that the point X is moved along is X*H.
0073         Xt = multitransp(X);
0074         XtEgrad = multiprod(Xt, Egrad);
0075         symXtEgrad = multisym(XtEgrad);
0076         XtEhess = multiprod(Xt, Ehess);
0077         Rhess = multiskew( XtEhess - multiprod(H, symXtEgrad) );
0078     end
0079     
0080     M.retr = @retraction;
0081     function Y = retraction(X, U, t)
0082         if nargin == 3
0083             tU = t*U;
0084         else
0085             tU = U;
0086         end
0087         Y = X + multiprod(X, tU);
0088         for i = 1 : k
0089             % This QR-based retraction is only a first-order approximation
0090             % of the exponential map, not a second-order one.
0091             [Q, R] = qr(Y(:, :, i));
0092             % The instruction with R ensures we are not flipping signs
0093             % of some columns, which should never happen in modern Matlab
0094             % versions but may be an issue with older versions.
0095             Y(:, :, i) = Q * diag(sign(sign(diag(R))+.5));
0096             % This is guaranteed to always yield orthogonal matrices with
0097             % determinant +1. Simply look at the eigenvalues of a skew
0098             % symmetric matrix, than at those of identity plus that matrix,
0099             % and compute their product for the determinant: it's stricly
0100             % positive in all cases.
0101         end
0102     end
0103     
0104     % A second order retraction is implemented here. To force its use,
0105     % after creating the factory M, execute M.retr = M.retr2.
0106     M.retr2 = @retraction2;
0107     function Y = retraction2(X, U, t)
0108         if nargin == 3
0109             tU = t*U;
0110         else
0111             tU = U;
0112         end
0113         Y = X + multiprod(X, tU);
0114         for i = 1 : k
0115             [Uk, ~, Vk] = svd(Y(:, :, k));
0116             Y(:, :, k) = Uk*Vk';
0117         end
0118     end
0119     
0120     M.exp = @exponential;
0121     function Y = exponential(X, U, t)
0122         if nargin == 3
0123             exptU = t*U;
0124         else
0125             exptU = U;
0126         end
0127         for i = 1 : k
0128             exptU(:, :, i) = expm(exptU(:, :, i));
0129         end
0130         Y = multiprod(X, exptU);
0131     end
0132     
0133     M.log = @logarithm;
0134     function U = logarithm(X, Y)
0135         U = multiprod(multitransp(X), Y);
0136         for i = 1 : k
0137             % The result of logm should be real in theory, but it is
0138             % numerically useful to force it.
0139             U(:, :, i) = real(logm(U(:, :, i)));
0140         end
0141         % Ensure the tangent vector is in the Lie algebra.
0142         U = multiskew(U);
0143     end
0144 
0145     M.hash = @(X) ['z' hashmd5(X(:))];
0146     
0147     M.rand = @() randrot(n, k);
0148     
0149     M.randvec = @randomvec;
0150     function U = randomvec(X) %#ok<INUSD>
0151         U = randskew(n, k);
0152         nrmU = sqrt(U(:).'*U(:));
0153         U = U / nrmU;
0154     end
0155     
0156     M.lincomb = @matrixlincomb;
0157     
0158     M.zerovec = @(x) zeros(n, n, k);
0159     
0160     M.transp = @(x1, x2, d) d;
0161     
0162     M.pairmean = @pairmean;
0163     function Y = pairmean(X1, X2)
0164         V = M.log(X1, X2);
0165         Y = M.exp(X1, .5*V);
0166     end
0167     
0168     M.dist = @(x, y) M.norm(x, M.log(x, y));
0169     
0170     M.vec = @(x, u_mat) u_mat(:);
0171     M.mat = @(x, u_vec) reshape(u_vec, [n, n, k]);
0172     M.vecmatareisometries = @() true;
0173 
0174 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005