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.
See also: stiefelfactory0001 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 % 0046 % See also: stiefelfactory 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) 0053 % Added egrad2rgrad and ehess2rhess 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 0103 M.egrad2rgrad = M.proj; 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); 0110 XtEgrad = multiprod(Xt, Egrad); 0111 symXtEgrad = multisym(XtEgrad); 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. 0290 % https://groups.google.com/g/manopttoolbox/c/KwdpyLiPUBw/m/aS-Yjq-pAwAJ 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