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.retr_polar after creating M. This switches from a QR-based
 computation to an SVD-based computation. If used, you may want to modify
 M.invretr accordingly.

 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.retr_polar after creating M. This switches from a QR-based
0029 % computation to an SVD-based computation. If used, you may want to modify
0030 % M.invretr accordingly.
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 %   June 18, 2019 (NB)
0055 %       Using qr_unique for the QR-based retraction.
0056 %   Nov. 19, 2019 (NB)
0057 %       Clarified that the isometric transport is not parallel transport
0058 %       along geodesics.
0059 
0060     
0061     if ~exist('k', 'var') || isempty(k)
0062         k = 1;
0063     end
0064     
0065     if k == 1
0066         M.name = @() sprintf('Rotations manifold SO(%d)', n);
0067     elseif k > 1
0068         M.name = @() sprintf('Product rotations manifold SO(%d)^%d', n, k);
0069     else
0070         error('k must be an integer no less than 1.');
0071     end
0072     
0073     M.dim = @() k*nchoosek(n, 2);
0074     
0075     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0076     
0077     M.norm = @(x, d) norm(d(:));
0078     
0079     M.typicaldist = @() pi*sqrt(n*k);
0080     
0081     M.proj = @(X, H) multiskew(multiprod(multitransp(X), H));
0082     
0083     M.tangent = @(X, H) multiskew(H);
0084     
0085     M.tangent2ambient_is_identity = false;
0086     M.tangent2ambient = @(X, U) multiprod(X, U);
0087     
0088     M.egrad2rgrad = M.proj;
0089     
0090     M.ehess2rhess = @ehess2rhess;
0091     function Rhess = ehess2rhess(X, Egrad, Ehess, H)
0092         % Reminder : H contains skew-symmeric matrices. The actual
0093         % direction that the point X is moved along is X*H.
0094         Xt = multitransp(X);
0095         XtEgrad = multiprod(Xt, Egrad);
0096         symXtEgrad = multisym(XtEgrad);
0097         XtEhess = multiprod(Xt, Ehess);
0098         Rhess = multiskew( XtEhess - multiprod(H, symXtEgrad) );
0099     end
0100     
0101     % This QR-based retraction is only a first-order approximation
0102     % of the exponential map, not a second-order one.
0103     M.retr_qr = @retraction_qr;
0104     function Y = retraction_qr(X, U, t)
0105         % It is necessary to call qr_unique rather than simply qr to ensure
0106         % this is a retraction, to avoid spurious column sign flips.
0107         XU = multiprod(X, U);
0108         if nargin < 3
0109             Y = qr_unique(X + XU);
0110         else
0111             Y = qr_unique(X + t*XU);
0112         end
0113         % This is guaranteed to always yield orthogonal matrices with
0114         % determinant +1. Indeed: look at the eigenvalues of a skew
0115         % symmetric matrix, then at those of "identity plus that matrix",
0116         % and compute their product for the determinant: it is strictly
0117         % positive in all cases.
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.retr_polar.
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     % Cheap vector transport
0236     M.transp = @(x1, x2, d) d;
0237     % This transporter is isometric (but it is /not/ parallel transport
0238     % along geodesics.)
0239     M.isotransp = M.transp;
0240     
0241     M.pairmean = @pairmean;
0242     function Y = pairmean(X1, X2)
0243         V = M.log(X1, X2);
0244         Y = M.exp(X1, .5*V);
0245     end
0246     
0247     M.dist = @(x, y) M.norm(x, M.log(x, y));
0248     
0249     M.vec = @(x, u_mat) u_mat(:);
0250     M.mat = @(x, u_vec) reshape(u_vec, [n, n, k]);
0251     M.vecmatareisometries = @() true;
0252 
0253 end

Generated on Tue 19-May-2020 18:46:12 by m2html © 2005