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.

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