# Alternative to finite differences for numerical computation of the Hessian of noisy function

I begin with the questions. I hope that they are not to vague.

My first question is: Is there any better way than least squares to get an estimate of the Hessian of a function $f$ which is affected by noise and is expensive to evaluate?

Second: If the least squares approach is suitable, what sampling strategies are there?

Third: Any pointers to further reading. I have tried to google, but I lack the terminology, so I might have had poor luck due to using the wrong keywords.

The problem in more detail:

I'm trying to compute a Hessian matrix of a multi-variable function $f:R^n\rightarrow R$ numerically. Analytic treatment is not within reach. The point $x$ at which I want to compute the hessian is an extremum with $f(x)=0$, and $\frac{\partial f}{\partial x} = 0$. Also it is known that the Hessian is symmetric and positive definite $0\le H^T=H=\frac{\partial^2 f}{\partial x^2}.$

The function $f$ itself is rather complex and takes several minutes to evaluate. Due to the definition of the function, it is known that it (theoretically) should have several continuous derivatives such that the Hessian is well defined.

Initially, I wanted to use finite differences to compute an estimate of the Hessian. However, due to finite precision arithmetics and other tolerances involved in the evaluation of the function $f$, a certain amount of numerical noise limits the precision of the evaluation, and the end result turns out to be drowned in noise. An attempt at an simple error analysis for an element in the Hessian matrix goes as follows. The finite difference approximation is:

$H_{ij}=\frac{f(x_i+h,x_j+k)-f(x_i+h,x_j-k)-f(x_i-h,x_j+k)+f(x_i-h,x_j-k) + O(\epsilon)}{4hk} + O(h^2) + O(hk) + O(k^2)$

so the error becomes

$\frac{O(\epsilon)}{4hk} + O(h^2) + O(hk) + O(k^2)$

Unfortunately, the $O(\epsilon)$ evaluation error is scaled by the inverse of the step-size and it is large enough that either the finite difference approximation becomes too poor or the scaled evaluation error becomes too large such that essentially no accuracy remains.

Next I have tried using least squares. From a Taylor expansion of $f$ we have

$f(x+h) \approx \underbrace{f(x)+h^T\frac{\partial f(x)}{\partial x}}_{=0} + \frac{1}{2}h^T\frac{\partial^2 f(x)}{\partial x^2}h$

My thinking is that we then could sample $y_k=f(h_k)$ locally for a number of local perturbations $h_k$. We then end up with a problem of the form

$H = \operatorname*{argmin}_H \sum_k (\frac{1}{2}h_k^THh_k - y_k)^2$, subject to: $H=H^T>0$

This problem is convex and easily solved with standard solvers.

However, here I have a problem with the sampling. Due to the cost of sampling the function, I don't want to use too many samples, but I see large variations in the estimated Hessian when too few samples are used. At the moment, I try to use latin hypercube sampling (LHS) to sample evenly but sparsely over the local domain (the dimension is about 10, but could go up to 20 in certain cases).

Thanks a lot for your thoughts on this problem!

• Ok, here's a suggestion: First do a wavelet transform of $f$, then denoise it with RiskShrink, and transform back. Next, since it's expensive to evaluate, apply a multivariate Chebyshev series. You'll be able to differentiate it in this representation quickly. Matlab has a wavelet toolbox that might allow you to succeed against the first step, Chebfun will help you with the second. – user14717 Jan 31 '18 at 15:33