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;
M.invretr = M.invretr_polar;

after creating M. This switches from a QR-based computation to an
SVD-based computation.

For n = 3, the code uses Rodrigues formulas for exp and log. Since most
optimization algorithms use retractions by default, you can force those
algorithms to use the Rodrigues formulas as follows:

M.retr = M.exp;
M.invretr = M.log;

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.
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Matrix multiply 2-D slices of 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 a 3D array
• multitransp Transpose the matrix slices of an N-D array (no complex conjugate)
• 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 %
0029 %     M.retr = M.retr_polar;
0030 %     M.invretr = M.invretr_polar;
0031 %
0032 % after creating M. This switches from a QR-based computation to an
0033 % SVD-based computation.
0034 %
0035 % For n = 3, the code uses Rodrigues formulas for exp and log. Since most
0036 % optimization algorithms use retractions by default, you can force those
0037 % algorithms to use the Rodrigues formulas as follows:
0038 %
0039 %     M.retr = M.exp;
0040 %     M.invretr = M.log;
0041 %
0042 %
0043 % By default, k = 1.
0044 %
0045 %
0047
0048 % This file is part of Manopt: www.manopt.org.
0049 % Original author: Nicolas Boumal, Dec. 30, 2012.
0050 % Contributors: Spencer Kraisler
0051 % Change log:
0052 %   Jan. 31, 2013 (NB)
0054 %   Oct. 21, 2016 (NB)
0055 %       Added M.retr2: a second-order retraction based on SVD.
0056 %   July 18, 2018 (NB)
0057 %       Fixed a bug in M.retr2 (only relevant for k > 1).
0058 %       Added inverse retraction as M.invretr.
0059 %       Retraction and inverse also available as M.retr_qr, M.invretr_qr.
0060 %       Renamed M.retr2 to M.retr_polar, and implemented M.invretr_polar.
0061 %       For backward compatibility, M.retr2 is still defined and is now
0062 %       equal to M.retr_polar. There is no M.invretr2.
0063 %   Sep.  06, 2018 (NB)
0064 %       Added M.isotransp, which is equal to M.transp since it is
0065 %       isometric.
0066 %   June 18, 2019 (NB)
0067 %       Using qr_unique for the QR-based retraction.
0068 %   Nov. 19, 2019 (NB)
0069 %       Clarified that the isometric transport is not parallel transport
0070 %       along geodesics.
0071 %   June 10, 2022 (NB)
0072 %       Following input from Spencer Kraisler on Manopt forum, added code
0073 %       for Rodrigues formulas, now used instead of expm / logm for n = 3.
0074
0075
0076     if ~exist('k', 'var') || isempty(k)
0077         k = 1;
0078     end
0079
0080     if k == 1
0081         M.name = @() sprintf('Rotations manifold SO(%d)', n);
0082     elseif k > 1
0083         M.name = @() sprintf('Product rotations manifold SO(%d)^%d', n, k);
0084     else
0085         error('k must be an integer no less than 1.');
0086     end
0087
0088     M.dim = @() k*nchoosek(n, 2);
0089
0090     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0091
0092     M.norm = @(x, d) norm(d(:));
0093
0094     M.typicaldist = @() pi*sqrt(n*k);
0095
0096     M.proj = @(X, H) multiskew(multiprod(multitransp(X), H));
0097
0098     M.tangent = @(X, H) multiskew(H);
0099
0100     M.tangent2ambient_is_identity = false;
0101     M.tangent2ambient = @(X, U) multiprod(X, U);
0102
0104
0105     M.ehess2rhess = @ehess2rhess;
0106     function Rhess = ehess2rhess(X, Egrad, Ehess, H)
0107         % Reminder : H contains skew-symmeric matrices. The actual
0108         % direction that the point X is moved along is X*H.
0109         Xt = multitransp(X);
0112         XtEhess = multiprod(Xt, Ehess);
0113         Rhess = multiskew( XtEhess - multiprod(H, symXtEgrad) );
0114     end
0115
0116     % This QR-based retraction is only a first-order approximation
0117     % of the exponential map, not a second-order one.
0118     M.retr_qr = @retraction_qr;
0119     function Y = retraction_qr(X, U, t)
0120         % It is necessary to call qr_unique rather than simply qr to ensure
0121         % this is a retraction, to avoid spurious column sign flips.
0122         XU = multiprod(X, U);
0123         if nargin < 3
0124             Y = qr_unique(X + XU);
0125         else
0126             Y = qr_unique(X + t*XU);
0127         end
0128         % This is guaranteed to always yield orthogonal matrices with
0129         % determinant +1. Indeed: look at the eigenvalues of a skew
0130         % symmetric matrix, then at those of "identity plus that matrix",
0131         % and compute their product for the determinant: it is strictly
0132         % positive in all cases.
0133     end
0134
0135     M.invretr_qr = @inverse_retraction_qr;
0136     function S = inverse_retraction_qr(X, Y)
0137
0138         % Assume k = 1 in this explanation:
0139         % If Y = Retr_X(XS) where S is a skew-symmetric matrix, then
0140         %  X(I+S) = YR
0141         % for some upper triangular matrix R. Multiply with X' on the left:
0142         %  I + S = (X'Y) R    (*)
0143         % Since S is skew-symmetric, add the transpose of this equation:
0144         %  2I + 0 = (X'Y) R + R' (X'Y)'
0145         % We can solve for R by calling solve_for_triu, then solve for S
0146         % using equation (*).
0147         R = zeros(n, n, k);
0148         XtY = multiprod(multitransp(X), Y);
0149         H = 2*eye(n);
0150         for kk = 1 : k
0151             R(:, :, kk) = solve_for_triu(XtY(:, :, kk), H);
0152         end
0153         % In exact arithmetic, taking the skew-symmetric part has the
0154         % effect of subtracting the identity from each slice; in inexact
0155         % arithmetic, taking the skew-symmetric part is beneficial to
0156         % further enforce tangency.
0157         S = multiskew(multiprod(XtY, R));
0158
0159     end
0160
0161     % A second-order retraction is implemented here. To force its use,
0162     % after creating the factory M, execute M.retr = M.retr_polar.
0163     M.retr_polar = @retraction_polar;
0164     function Y = retraction_polar(X, U, t)
0165         if nargin == 3
0166             tU = t*U;
0167         else
0168             tU = U;
0169         end
0170         Y = X + multiprod(X, tU);
0171         for kk = 1 : k
0172             [Uk, ~, Vk] = svd(Y(:, :, kk));
0173             Y(:, :, kk) = Uk*Vk';
0174             % One can check that det(Uk*Vk') = det(X) for all skew-sym U.
0175             % That's because X + XU = X(I+U), and U is normal (UU' = U'U)
0176             % hence it can be unitarily diagonalized as U = WDW' with W
0177             % unitary and D diagonal; the eigenvalues of U are purely
0178             % imaginary because U+U' = 0, and they come in conjugate pairs
0179             % because U is real (there's a zero eigenvalue if n is odd).
0180             % This way, it follows that
0181             %   det(X+XU) = det(X)det(I+U) = det(X)det(I+D)
0182             % and det(I+D) is real, strictly positive provided n >= 2.
0183         end
0184     end
0185
0186     M.invretr_polar = @inverse_retraction_polar;
0187     function S = inverse_retraction_polar(X, Y)
0188
0189         % Assume k = 1 in this explanation:
0190         % If Y = Retr_X(XS) where S is a skew-symmetric matrix, then
0191         %  X(I+S) = YM
0192         % for some symmetric matrix M. Multiply with X' on the left:
0193         %  I + S = (X'Y) M    (*)
0194         % Since S is skew-symmetric, add the transpose of this equation:
0195         %  2I + 0 = (X'Y) M + M (X'Y)'
0196         % We can solve for M by calling sylvester, then solve for S
0197         % using equation (*).
0198         MM = zeros(n, n, k);
0199         XtY = multiprod(multitransp(X), Y);
0200         H = 2*eye(n);
0201         for kk = 1 : k
0202             MM(:, :, kk) = sylvester_nochecks(XtY(:, :, kk), XtY(:, :, kk)', H);
0203         end
0204         % In exact arithmetic, taking the skew-symmetric part has the
0205         % effect of subtracting the identity from each slice; in inexact
0206         % arithmetic, taking the skew-symmetric part is beneficial to
0207         % further enforce tangency.
0208         S = multiskew(multiprod(XtY, MM));
0209
0210     end
0211
0212     % By default, use QR retraction
0213     M.retr = M.retr_qr;
0214     M.invretr = M.invretr_qr;
0215
0216     % For backward compatibility:
0217     M.retr2 = M.retr_polar;
0218
0219     % Special case: for n = 3, we use Rodrigues formulas for expm / logm.
0220     if n == 3
0221         myexpm = @expm_SO3;
0222         mylogm = @logm_SO3;
0223     else
0224         myexpm = @expm;
0225         mylogm = @logm;
0226     end
0227
0228     M.exp = @exponential;
0229     function Y = exponential(X, U, t)
0230         if nargin == 3
0231             exptU = t*U;
0232         else
0233             exptU = U;
0234         end
0235         for kk = 1 : k
0236             exptU(:, :, kk) = myexpm(exptU(:, :, kk));
0237         end
0238         Y = multiprod(X, exptU);
0239     end
0240
0241     M.log = @logarithm;
0242     function U = logarithm(X, Y)
0243         U = multiprod(multitransp(X), Y);
0244         for kk = 1 : k
0245             % The result of logm should be real in theory, but it is
0246             % numerically useful to force it.
0247             U(:, :, kk) = real(mylogm(U(:, :, kk)));
0248         end
0249         % Ensure the tangent vector is in the Lie algebra.
0250         U = multiskew(U);
0251     end
0252
0253     M.hash = @(X) ['z' hashmd5(X(:))];
0254
0255     M.rand = @() randrot(n, k);
0256
0257     M.randvec = @randomvec;
0258     function U = randomvec(X) %#ok<INUSD>
0259         U = randskew(n, k);
0260         nrmU = sqrt(U(:).'*U(:));
0261         U = U / nrmU;
0262     end
0263
0264     M.lincomb = @matrixlincomb;
0265
0266     M.zerovec = @(x) zeros(n, n, k);
0267
0268     % Cheap vector transport
0269     M.transp = @(x1, x2, d) d;
0270     % This transporter is isometric (but it is /not/ parallel transport
0271     % along geodesics.)
0272     M.isotransp = M.transp;
0273
0274     M.pairmean = @pairmean;
0275     function Y = pairmean(X1, X2)
0276         V = M.log(X1, X2);
0277         Y = M.exp(X1, .5*V);
0278     end
0279
0280     M.dist = @(x, y) M.norm(x, M.log(x, y));
0281
0282     M.vec = @(x, u_mat) u_mat(:);
0283     M.mat = @(x, u_vec) reshape(u_vec, [n, n, k]);
0284     M.vecmatareisometries = @() true;
0285     M.lie_identity = @() repmat(eye(n), [1, 1, k]);
0286
0287 end
0288
0289 % The following code is based on a proposal by Spencer Kraisler, June 2022.
0291 %
0292 % Rodrigues formula for matrix exponential on SO(3): R = expm(phi)
0293 function R = expm_SO3(phi)
0294     phi_vee = [-phi(2, 3); phi(1, 3); -phi(1, 2)];
0295     norm_phi_vee = norm(phi_vee);
0296     if norm_phi_vee > 0
0297         q1 = sin(norm_phi_vee)/norm_phi_vee;
0298         r = norm_phi_vee / 2;
0299         q2 = (sin(r)/r).^2 / 2;
0300         R = eye(3) + q1*phi + q2*phi^2;
0301     else
0302         R = eye(3);
0303     end
0304 end
0305 % Rodrigues formula for inverting matrix exp on SO(3): phi = logm(R)
0306 function phi = logm_SO3(R)
0307     t = trace(R);
0308     norm_t = real(acos((t - 1)/2));
0309     if norm_t > 0 % could fail even when trace(R) < 3, because sensitive
0310         q = .5*norm_t/sin(norm_t);
0311     else
0312         q = .5;   % even with this, phi (below) could be nonzero
0313     end
0314     phi = q * [R(3, 2) - R(2, 3); R(1, 3) - R(3, 1); R(2, 1) - R(1, 2)];
0315     phi = [0 -phi(3) phi(2); phi(3) 0 -phi(1); -phi(2) phi(1) 0];
0316 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005