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, i.e., X'*X = eye(n) if k = 1, or
 X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1. Furthermore,
 all these matrices have determinant +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.

 This is important:
 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), as the input H will then be a skew-symmetric matrix
 (but the output must not be, as the output is the Hessian in the
 embedding Euclidean space.)

 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. This is notably useful for the
 checkhessian tool.

 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, i.e., X'*X = eye(n) if k = 1, or
0010 % X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1. Furthermore,
0011 % all these matrices have determinant +1.
0012 %
0013 % This is a description of SO(n)^k with the induced metric from the
0014 % embedding space (R^nxn)^k, i.e., this manifold is a Riemannian
0015 % submanifold of (R^nxn)^k endowed with the usual trace inner product.
0016 %
0017 % This is important:
0018 % Tangent vectors are represented in the Lie algebra, i.e., as skew
0019 % symmetric matrices. Use the function M.tangent2ambient(X, H) to switch
0020 % from the Lie algebra representation to the embedding space
0021 % representation. This is often necessary when defining
0022 % problem.ehess(X, H), as the input H will then be a skew-symmetric matrix
0023 % (but the output must not be, as the output is the Hessian in the
0024 % embedding Euclidean space.)
0025 %
0026 % By default, the retraction is only a first-order approximation of the
0027 % exponential. To force the use of a second-order approximation, call
0028 % M.retr = M.retr2 after creating M. This switches from a QR-based
0029 % computation to an SVD-based computation. This is notably useful for the
0030 % checkhessian tool.
0031 %
0032 % By default, k = 1.
0033 %
0034 % See also: stiefelfactory
0035 
0036 % This file is part of Manopt: www.manopt.org.
0037 % Original author: Nicolas Boumal, Dec. 30, 2012.
0038 % Contributors:
0039 % Change log:
0040 %   Jan. 31, 2013 (NB)
0041 %       Added egrad2rgrad and ehess2rhess
0042 %   Oct. 21, 2016 (NB)
0043 %       Added M.retr2: a second-order retraction based on SVD.
0044 %   July 18, 2018 (NB)
0045 %       Fixed a bug in M.retr2 (only relevant for k > 1).
0046 %       Added inverse retraction as M.invretr.
0047 %       Retraction and inverse also available as M.retr_qr, M.invretr_qr.
0048 %       Renamed M.retr2 to M.retr_polar, and implemented M.invretr_polar.
0049 %       For backward compatibility, M.retr2 is still defined and is now
0050 %       equal to M.retr_polar. There is no M.invretr2.
0051 %   Sep.  06, 2018 (NB)
0052 %       Added M.isotransp, which is equal to M.transp since it is
0053 %       isometric.
0054 
0055     
0056     if ~exist('k', 'var') || isempty(k)
0057         k = 1;
0058     end
0059     
0060     if k == 1
0061         M.name = @() sprintf('Rotations manifold SO(%d)', n);
0062     elseif k > 1
0063         M.name = @() sprintf('Product rotations manifold SO(%d)^%d', n, k);
0064     else
0065         error('k must be an integer no less than 1.');
0066     end
0067     
0068     M.dim = @() k*nchoosek(n, 2);
0069     
0070     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0071     
0072     M.norm = @(x, d) norm(d(:));
0073     
0074     M.typicaldist = @() pi*sqrt(n*k);
0075     
0076     M.proj = @(X, H) multiskew(multiprod(multitransp(X), H));
0077     
0078     M.tangent = @(X, H) multiskew(H);
0079     
0080     M.tangent2ambient_is_identity = false;
0081     M.tangent2ambient = @(X, U) multiprod(X, U);
0082     
0083     M.egrad2rgrad = M.proj;
0084     
0085     M.ehess2rhess = @ehess2rhess;
0086     function Rhess = ehess2rhess(X, Egrad, Ehess, H)
0087         % Reminder : H contains skew-symmeric matrices. The actual
0088         % direction that the point X is moved along is X*H.
0089         Xt = multitransp(X);
0090         XtEgrad = multiprod(Xt, Egrad);
0091         symXtEgrad = multisym(XtEgrad);
0092         XtEhess = multiprod(Xt, Ehess);
0093         Rhess = multiskew( XtEhess - multiprod(H, symXtEgrad) );
0094     end
0095     
0096     M.retr_qr = @retraction_qr;
0097     function Y = retraction_qr(X, U, t)
0098         if nargin == 3
0099             tU = t*U;
0100         else
0101             tU = U;
0102         end
0103         Y = X + multiprod(X, tU);
0104         for kk = 1 : k
0105             % This QR-based retraction is only a first-order approximation
0106             % of the exponential map, not a second-order one.
0107             [Q, R] = qr(Y(:, :, kk));
0108             % The instruction with R ensures we are not flipping signs
0109             % of some columns, which should never happen in modern Matlab
0110             % versions but may be an issue with older versions.
0111             Y(:, :, kk) = Q * diag(sign(sign(diag(R))+.5));
0112             % This is guaranteed to always yield orthogonal matrices with
0113             % determinant +1. Simply look at the eigenvalues of a skew
0114             % symmetric matrix, then at those of identity plus that matrix,
0115             % and compute their product for the determinant: it's stricly
0116             % positive in all cases.
0117         end
0118     end
0119 
0120     M.invretr_qr = @inverse_retraction_qr;
0121     function S = inverse_retraction_qr(X, Y)
0122         
0123         % Assume k = 1 in this explanation:
0124         % If Y = Retr_X(XS) where S is a skew-symmetric matrix, then
0125         %  X(I+S) = YR
0126         % for some upper triangular matrix R. Multiply with X' on the left:
0127         %  I + S = (X'Y) R    (*)
0128         % Since S is skew-symmetric, add the transpose of this equation:
0129         %  2I + 0 = (X'Y) R + R' (X'Y)'
0130         % We can solve for R by calling solve_for_triu, then solve for S
0131         % using equation (*).
0132         R = zeros(n, n, k);
0133         XtY = multiprod(multitransp(X), Y);
0134         H = 2*eye(n);
0135         for kk = 1 : k
0136             R(:, :, kk) = solve_for_triu(XtY(:, :, kk), H);
0137         end
0138         % In exact arithmetic, taking the skew-symmetric part has the
0139         % effect of subtracting the identity from each slice; in inexact
0140         % arithmetic, taking the skew-symmetric part is beneficial to
0141         % further enforce tangency.
0142         S = multiskew(multiprod(XtY, R));
0143         
0144     end
0145     
0146     % A second order retraction is implemented here. To force its use,
0147     % after creating the factory M, execute M.retr = M.retr2.
0148     M.retr_polar = @retraction_polar;
0149     function Y = retraction_polar(X, U, t)
0150         if nargin == 3
0151             tU = t*U;
0152         else
0153             tU = U;
0154         end
0155         Y = X + multiprod(X, tU);
0156         for kk = 1 : k
0157             [Uk, ~, Vk] = svd(Y(:, :, kk));
0158             Y(:, :, kk) = Uk*Vk';
0159         end
0160     end
0161 
0162     M.invretr_polar = @inverse_retraction_polar;
0163     function S = inverse_retraction_polar(X, Y)
0164         
0165         % Assume k = 1 in this explanation:
0166         % If Y = Retr_X(XS) where S is a skew-symmetric matrix, then
0167         %  X(I+S) = YM
0168         % for some symmetric matrix M. Multiply with X' on the left:
0169         %  I + S = (X'Y) M    (*)
0170         % Since S is skew-symmetric, add the transpose of this equation:
0171         %  2I + 0 = (X'Y) M + M (X'Y)'
0172         % We can solve for M by calling sylvester, then solve for S
0173         % using equation (*).
0174         MM = zeros(n, n, k);
0175         XtY = multiprod(multitransp(X), Y);
0176         H = 2*eye(n);
0177         for kk = 1 : k
0178             MM(:, :, kk) = sylvester_nochecks(XtY(:, :, kk), XtY(:, :, kk)', H);
0179         end
0180         % In exact arithmetic, taking the skew-symmetric part has the
0181         % effect of subtracting the identity from each slice; in inexact
0182         % arithmetic, taking the skew-symmetric part is beneficial to
0183         % further enforce tangency.
0184         S = multiskew(multiprod(XtY, MM));
0185         
0186     end
0187 
0188     % By default, use QR retraction
0189     M.retr = M.retr_qr;
0190     M.invretr = M.invretr_qr;
0191     
0192     % For backward compatibility:
0193     M.retr2 = M.retr_polar;
0194     
0195     M.exp = @exponential;
0196     function Y = exponential(X, U, t)
0197         if nargin == 3
0198             exptU = t*U;
0199         else
0200             exptU = U;
0201         end
0202         for kk = 1 : k
0203             exptU(:, :, kk) = expm(exptU(:, :, kk));
0204         end
0205         Y = multiprod(X, exptU);
0206     end
0207     
0208     M.log = @logarithm;
0209     function U = logarithm(X, Y)
0210         U = multiprod(multitransp(X), Y);
0211         for kk = 1 : k
0212             % The result of logm should be real in theory, but it is
0213             % numerically useful to force it.
0214             U(:, :, kk) = real(logm(U(:, :, kk)));
0215         end
0216         % Ensure the tangent vector is in the Lie algebra.
0217         U = multiskew(U);
0218     end
0219 
0220     M.hash = @(X) ['z' hashmd5(X(:))];
0221     
0222     M.rand = @() randrot(n, k);
0223     
0224     M.randvec = @randomvec;
0225     function U = randomvec(X) %#ok<INUSD>
0226         U = randskew(n, k);
0227         nrmU = sqrt(U(:).'*U(:));
0228         U = U / nrmU;
0229     end
0230     
0231     M.lincomb = @matrixlincomb;
0232     
0233     M.zerovec = @(x) zeros(n, n, k);
0234     
0235     M.transp = @(x1, x2, d) d;
0236     M.isotransp = M.transp; % the transport is isometric
0237     
0238     M.pairmean = @pairmean;
0239     function Y = pairmean(X1, X2)
0240         V = M.log(X1, X2);
0241         Y = M.exp(X1, .5*V);
0242     end
0243     
0244     M.dist = @(x, y) M.norm(x, M.log(x, y));
0245     
0246     M.vec = @(x, u_mat) u_mat(:);
0247     M.mat = @(x, u_vec) reshape(u_vec, [n, n, k]);
0248     M.vecmatareisometries = @() true;
0249 
0250 end

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005