0001 function checkhessian(problem, x, d)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061 if ~canGetCost(problem)
0062 error('It seems no cost was provided.');
0063 end
0064 if ~canGetGradient(problem)
0065 warning('manopt:checkhessian:nograd', ...
0066 'It seems no gradient was provided.');
0067 end
0068 if ~canGetHessian(problem)
0069 warning('manopt:checkhessian:nohess', ...
0070 'It seems no Hessian was provided.');
0071 end
0072
0073 x_isprovided = exist('x', 'var') && ~isempty(x);
0074 d_isprovided = exist('d', 'var') && ~isempty(d);
0075
0076 if ~x_isprovided && d_isprovided
0077 error('If d is provided, x must be too, since d is tangent at x.');
0078 end
0079
0080
0081 if ~x_isprovided
0082 x = problem.M.rand();
0083 end
0084 if ~d_isprovided
0085 d = problem.M.randvec(x);
0086 end
0087
0088
0089
0090
0091
0092
0093 storedb = StoreDB();
0094 xkey = storedb.getNewKey();
0095 f0 = getCost(problem, x, storedb, xkey);
0096 gradx = getGradient(problem, x, storedb, xkey);
0097 df0 = problem.M.inner(x, d, gradx);
0098 hessxd = getHessian(problem, x, d, storedb, xkey);
0099 d2f0 = problem.M.inner(x, d, hessxd);
0100
0101
0102
0103 if isfield(problem.M, 'exp')
0104 stepper = problem.M.exp;
0105 extra_message = '';
0106 else
0107 stepper = problem.M.retr;
0108 fprintf(['* M.exp() is not available: using M.retr() instead.\n' ...
0109 '* Please check the manifold documentation to see if\n' ...
0110 '* the retraction is second order. If not, the slope\n' ...
0111 '* test is allowed to fail at non-critical x.\n']);
0112 extra_message = ['(But do mind the message above: the slope may\n' ...
0113 'be allowed to be 2 at non-critical points x.)\n'];
0114 end
0115
0116
0117
0118
0119
0120 h = logspace(-8, 0, 51);
0121 value = zeros(size(h));
0122 for i = 1 : length(h)
0123 y = stepper(x, d, h(i));
0124 ykey = storedb.getNewKey();
0125 value(i) = getCost(problem, y, storedb, ykey);
0126 storedb.remove(ykey);
0127 end
0128
0129
0130
0131 model = polyval([.5*d2f0 df0 f0], h);
0132
0133
0134 err = abs(model - value);
0135
0136
0137 loglog(h, err);
0138 title(sprintf(['Hessian check.\nThe slope of the continuous line ' ...
0139 'should match that of the dashed\n(reference) line ' ...
0140 'over at least a few orders of magnitude for h.']));
0141 xlabel('h');
0142 ylabel('Approximation error');
0143
0144 line('xdata', [1e-8 1e0], 'ydata', [1e-16 1e8], ...
0145 'color', 'k', 'LineStyle', '--', ...
0146 'YLimInclude', 'off', 'XLimInclude', 'off');
0147
0148
0149 if ~all( err < 1e-12 )
0150
0151
0152
0153 isModelExact = false;
0154 window_len = 10;
0155 [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
0156 else
0157
0158
0159 isModelExact = true;
0160 range = 1:numel(h);
0161 poly = polyfit(log10(h), err, 1);
0162
0163 poly(end) = log10(poly(end));
0164
0165 title(sprintf(...
0166 ['Hessian check.\n'...
0167 'It seems the quadratic model is exact:\n'...
0168 'Model error is numerically zero for all h.']));
0169 end
0170 hold all;
0171 loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
0172 hold off;
0173
0174 if ~isModelExact
0175 fprintf('The slope should be 3. It appears to be: %g.\n', poly(1));
0176 fprintf(['If it is far from 3, then directional derivatives,\n' ...
0177 'the gradient or the Hessian might be erroneous.\n', ...
0178 extra_message]);
0179 else
0180 fprintf(['The quadratic model appears to be exact ' ...
0181 '(within numerical precision),\n'...
0182 'hence the slope computation is irrelevant.\n']);
0183 end
0184
0185
0186
0187 if isfield(problem.M, 'tangent')
0188 hess = getHessian(problem, x, d, storedb, xkey);
0189 phess = problem.M.tangent(x, hess);
0190 residual = problem.M.lincomb(x, 1, hess, -1, phess);
0191 err = problem.M.norm(x, residual);
0192 fprintf('Tangency residual should be zero, or very close; ');
0193 fprintf('residual: %g.\n', err);
0194 fprintf(['If it is far from 0, then the Hessian is not in the ' ...
0195 'tangent space.\n']);
0196 else
0197 fprintf(['Unfortunately, Manopt was unable to verify that the '...
0198 'output of the Hessian call is indeed a tangent ' ...
0199 'vector.\nPlease verify this manually.']);
0200 end
0201
0202
0203 d1 = problem.M.randvec(x);
0204 d2 = problem.M.randvec(x);
0205 h1 = getHessian(problem, x, d1, storedb, xkey);
0206 h2 = getHessian(problem, x, d2, storedb, xkey);
0207
0208
0209 a = randn(1);
0210 b = randn(1);
0211 ad1pbd2 = problem.M.lincomb(x, a, d1, b, d2);
0212 had1pbd2 = getHessian(problem, x, ad1pbd2, storedb, xkey);
0213 ahd1pbhd2 = problem.M.lincomb(x, a, h1, b, h2);
0214 errvec = problem.M.lincomb(x, 1, had1pbd2, -1, ahd1pbhd2);
0215 errvecnrm = problem.M.norm(x, errvec);
0216 had1pbd2nrm = problem.M.norm(x, had1pbd2);
0217 fprintf(['||a*H[d1] + b*H[d2] - H[a*d1+b*d2]|| should be zero, or ' ...
0218 'very close.\n\tValue: %g (norm of H[a*d1+b*d2]: %g)\n'], ...
0219 errvecnrm, had1pbd2nrm);
0220 fprintf('If it is far from 0, then the Hessian is not linear.\n');
0221
0222
0223 v1 = problem.M.inner(x, d1, h2);
0224 v2 = problem.M.inner(x, h1, d2);
0225 value = v1-v2;
0226 fprintf(['<d1, H[d2]> - <H[d1], d2> should be zero, or very close.' ...
0227 '\n\tValue: %g - %g = %g.\n'], v1, v2, value);
0228 fprintf('If it is far from 0, then the Hessian is not symmetric.\n');
0229
0230
0231
0232
0233
0234
0235
0236 if isfield(problem.M, 'tangent2ambient_is_identity') && ...
0237 ~problem.M.tangent2ambient_is_identity
0238
0239 fprintf('\n\n=== Special message ===\n');
0240 fprintf(['For this manifold, tangent vectors are represented\n' ...
0241 'differently from their ambient space representation.\n' ...
0242 'In practice, this means that when defining\n' ...
0243 'v = problem.ehess(x, u), one may need to call\n' ...
0244 'u = problem.M.tangent2ambient(x, u) first, so as to\n'...
0245 'transform u into an ambient vector, if this is more\n' ...
0246 'convenient. The output of ehess should be an ambient\n' ...
0247 'vector (it will be transformed to a tangent vector\n' ...
0248 'automatically).\n']);
0249
0250 end
0251
0252 if ~canGetHessian(problem)
0253 norm_grad = problem.M.norm(x, gradx);
0254 fprintf(['\nWhen using checkhessian with a finite difference ' ...
0255 'approximation, the norm of the residual\nshould be ' ...
0256 'compared against the norm of the gradient at the ' ...
0257 'point under consideration (%.2g).\nFurthermore, it ' ...
0258 'is expected that the FD operator is only approximately' ...
0259 ' symmetric.\nOf course, the slope can also be off.\n'], ...
0260 norm_grad);
0261 end
0262
0263 end