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.

## CROSS-REFERENCE INFORMATION This function calls:
• randrot Generates uniformly random rotation matrices.
• randskew Generates random skew symmetric matrices with normal entries.
• solve_for_triu Solve the linear matrix equation AX + X'A' = H for X upper triangular.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
• multiskew Returns the skew-symmetric parts of the matrices in the 3D matrix X.
• multisym Returns the symmetric parts of the matrices in the 3D matrix X
• multitransp Transposing arrays of matrices.
• qr_unique Thin QR factorization ensuring diagonal of R is real, positive if possible.
• sylvester_nochecks Solve Sylvester equation without input checks.
This function is called by:

## 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 %
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)
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
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);
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