Description of checkmanifold

0001 function checkmanifold(M)
0002 % Run a collection of tests on a manifold obtained from a manopt factory
0003 %
0004 % function checkmanifold(M)
0005 %
0006 % M should be a manifold structure obtained from a Manopt factory. This
0007 % tool runs a collection of tests on M to verify (to some extent) that M is
0008 % indeed a valid description of a Riemannian manifold.
0009 %
0010 % This tool is work in progress: your suggestions for additional tests are
0011 % welcome on our forum or as pull requests on GitHub.
0012 %
0013 % See also: checkretraction
0014 
0015 % This file is part of Manopt: www.manopt.org.
0016 % Original author: Nicolas Boumal, Aug. 31, 2018.
0017 % Contributors:
0018 % Change log:
0019 %   April 12, 2020 (NB):
0020 %       Now checking M.dist(x, M.exp(x, v, t)) for several values of t
0021 %       because this test is only valid for norm(x, tv) <= inj(x).
0022 %   May 19, 2020 (NB):
0023 %       Now checking M.dim().
0024 %   Jan  8, 2021 (NB):
0025 %       Added partial checks of M.inner, M.tangent2ambient, M.proj, ...
0026 
0027     assert(isstruct(M), 'M must be a structure.');
0028     
0029     %% List required fields that must be function handles here
0030     list_of_functions = {'name', 'dim', 'inner', 'norm', 'typicaldist', ...
0031                          'proj', 'tangent', 'egrad2rgrad', 'retr', ...
0032                          'rand', 'randvec', 'zerovec', 'lincomb'};
0033     for k = 1 : numel(list_of_functions)
0034         field = list_of_functions{k};
0035         if ~(isfield(M, field) && isa(M.(field), 'function_handle'))
0036             fprintf('M.%s must be a function handle.\n', field);
0037         end
0038     end
0039     
0040     %% List recommended fields that should be function handles here
0041     list_of_functions = {'dist', 'ehess2rhess', 'exp', 'log', 'hash', ...
0042                          'transp', 'pairmean', 'vec', 'mat', ...
0043                          'vecmatareisometries'};
0044     for k = 1 : numel(list_of_functions)
0045         field = list_of_functions{k};
0046         if ~(isfield(M, field) && isa(M.(field), 'function_handle'))
0047             fprintf(['M.%s should ideally (but does not have to) be ' ...
0048                      'a function handle.\n'], field);
0049         end
0050     end
0051     
0052     %% Check random generators
0053     try
0054         x = M.rand();
0055         v = M.randvec(x);
0056         fprintf('Random tangent vector norm: %g (should be 1).\n', ...
0057                 M.norm(x, v));
0058         z = M.lincomb(x, 1, v, -1, v);
0059         fprintf('norm(v - v)_x = %g (should be 0).\n', M.norm(x, z));
0060     catch up %#ok<NASGU>
0061         fprintf('Couldn''t check rand, randvec, lincomb.\n');
0062     end
0063     
0064     %% Check inner product
0065     try
0066         x = M.rand();
0067         
0068         % Check symmetry
0069         u = M.randvec(x);
0070         v = M.randvec(x);
0071         uv = M.inner(x, u, v);
0072         vu = M.inner(x, v, u);
0073         fprintf('<u, v>_x = %g, <v, u>_x = %g, difference = %g (should be 0).\n', uv, vu, uv-vu);
0074         
0075         % Check linearity (and owing to symmetry: bilinearity)
0076         a = randn();
0077         b = randn();
0078         w = M.lincomb(x, a, u, b, v);
0079         z = M.randvec(x);
0080         wz = M.inner(x, w, z);
0081         wzbis = a*M.inner(x, u, z) + b*M.inner(x, v, z);
0082         fprintf('<au+bv, z>_x = %g, a<u, z>_x + b<v, z>_x = %g, difference = %g (should be 0).\n', wz, wzbis, wz-wzbis);
0083         
0084         % Should check positive definiteness too: it's somehow part of the
0085         % check for M.dim() below.
0086         
0087     catch up %#ok<NASGU>
0088         fprintf('Couldn''t check inner.\n');
0089     end
0090     
0091     %% Check tangent2ambient, proj, norm
0092     try
0093         x = M.rand();
0094         v = M.randvec(x);
0095         va = M.tangent2ambient(x, v);
0096         vp = M.proj(x, va);
0097         v_min_vp = M.lincomb(x, 1, v, -1, vp);
0098         df = M.norm(x, v_min_vp);
0099         fprintf('Norm of tangent vector minus its projection to tangent space: %g (should be zero).\n', df);
0100         
0101         % Should check that proj is linear, self-adjoint, idempotent.
0102         % The issue for generic code is that manifolds do not provide means
0103         % to generate random vectors in the ambient space.
0104         
0105     catch up %#ok<NASGU>
0106         fprintf('Couldn''t check tangent2ambient, proj, norm\n');
0107     end    
0108     
0109     %% Checking exp and dist
0110     try
0111         x = M.rand();
0112         v = M.randvec(x);
0113         for t = logspace(-8, 1, 10)
0114             y = M.exp(x, v, t);
0115             d = M.dist(x, y);
0116             err = d - abs(t)*M.norm(x, v);
0117             fprintf(['dist(x, M.exp(x, v, t)) - abs(t)*M.norm(x, v) = ' ...
0118                      '%g (t = %.1e; should be zero for small enough t).\n'], ...
0119                      err, t);
0120         end
0121     catch up %#ok<NASGU>
0122         fprintf('Couldn''t check exp and dist.\n');
0123         % Perhaps we want to rethrow(up) ?
0124         % Alternatively, we could check if exp and dist are available and
0125         % silently pass this test if not, but this way is more informative.
0126     end
0127     
0128     %% Checking mat, vec, vecmatareisometries
0129     try
0130         x = M.rand();
0131         u = M.randvec(x);
0132         v = M.randvec(x);
0133         U = M.vec(x, u);
0134         V = M.vec(x, v);
0135         if ~iscolumn(U) || ~iscolumn(V)
0136             fprintf('M.vec should return column vectors: they are not.\n');
0137         end
0138         if ~isreal(U) || ~isreal(V)
0139             fprintf('M.vec should return real vectors: they are not real.\n');
0140         end
0141         fprintf(['Unless otherwise stated, M.vec seems to return real ' ...
0142                  'column vectors, as intended.\n']);
0143         ru = M.norm(x, M.lincomb(x, 1, M.mat(x, U), -1, u));
0144         rv = M.norm(x, M.lincomb(x, 1, M.mat(x, V), -1, v));
0145         fprintf(['Checking mat/vec are inverse pairs: ' ...
0146                  '%g, %g (should be two zeros).\n'], ru, rv);
0147         a = randn(1);
0148         b = randn(1);
0149         fprintf('Checking if vec is linear: %g (should be zero).\n', ...
0150                 norm(M.vec(x, M.lincomb(x, a, u, b, v)) - (a*U + b*V)));
0151         if M.vecmatareisometries()
0152             fprintf('M.vecmatareisometries says true.\n');
0153         else
0154             fprintf('M.vecmatareisometries says false.\n');
0155         end
0156         fprintf('If true, this should be zero: %g.\n', ...
0157                     U(:).'*V(:) - M.inner(x, u, v));
0158     catch up %#ok<NASGU>
0159         fprintf('Couldn''t check mat, vec, vecmatareisometries.\n');
0160     end
0161     
0162     %% Checking dim
0163     dim_threshold = 200;
0164     if M.dim() <= dim_threshold
0165         x = M.rand();
0166         n = M.dim() + 1;
0167         B = cell(n, 1);
0168         for k = 1 : n
0169             B{k} = M.randvec(x);
0170         end
0171         G = grammatrix(M, x, B);
0172         eigG = sort(real(eig(G)), 'descend');
0173         fprintf('Testing M.dim() (works best when dimension is small):\n');
0174         fprintf('\tIf this number is machine-precision zero, then M.dim() may be too large: %g\n', eigG(n-1));
0175         fprintf('\tIf this number is not machine-precision zero, then M.dim() may be too small: %g\n', eigG(n));
0176     else
0177         fprintf('M.dim() not tested because it is > %d.\n', dim_threshold);
0178     end
0179     
0180     %% Recommend calling checkretraction
0181     fprintf('It is recommended also to call checkretraction.\n');
0182 
0183 end
0184
checkmanifold

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
