# How do I numerically calculate a function from its noisy gradient?

I have the model $\ s(x,y)=x^2+y^2, 0 \leq x \leq 1, 0 \leq y \leq 1$.

Instead of observing the model directly I am observing the derivatives of the model + some noise (e):

$\ p(x,y)=s_x+e, q(x,y)=s_y+e$

From measurements of p(x,y and q(x,y) I want to estimate s(x). Say I know that s(0,0)=0.

According to the gradient theorem: $\ s(x,y) = \int_{(0,0)}^{(x,y)} [s_x,s_y]d\vec{r}$

regardless along which path we integrate.

As a small experiment (in Matlab) I added normal distributed noise, N(0,1), to p=2x and q=2y. Then I integrated first along x followed by along y: SXY. Next I integrated first along y followed by along x: SYX.

The results show that the gradient theorem do not hold in this case (because of the noise):   The Root Mean Square Errors relative to the model are:

ErmsXY =
0.1125
ErmsYX =
0.0920


How can I find a better (less RMS error and smoother) estimate of s from p and q?

EDIT:

From what I read; using the curve integral is referred to as local integration. There are also global integration methods where one instead tries to choose a S(x,y) that minimizes:

$\ \int_0^1\int_0^1 [|S_x - P|^2 + |S_y - Q|^2] dxdy$

Global integration methods are supposed to give better results when the gradient is noisy, but how do I do this in practice?

EDIT 2:

One approach that I have used is this:

first we introduce linear derivation operators: $\ s_x = D_x*s,s_y = D_y*s$.

The result is the following linear equation system:

$\ D_x*s=p, D_y*s=q$

Next find a Least Square Error solution to these equations. A LSE solution to these equations is supposed to be equivalent to minimizing the integral from above. How can this be shown?

The results are good: The RMS error is about 1/5 of that of SXY and SYX and the solution is also smoother.

However there are some drawbacks to this approach:

1. it is difficult to implement; must use central differences and "flatten" 2D s matrix into vector etc.

2. The derivation matrices are very large and sparse so they may consume a lot of RAM.

Another approach that seems potentially both simpler to code, less RAM consuming and faster is to use FFT. In the Fourier space these pdes become an algebraic equation. This is known as the Frankot-Chellappa algorithm, but unfortunately I have not got it to work on my example data.

You can filter either the gradient itself or the result, $s$. You would need to know the characteristics of the true gradients fairly well to know what the frequency bandwidth is. At that point you could design a low-pass filter that would preserve the signal but get rid of higher frequency noise.

• Thanks Jim. So I may e.g. take SXY and replace each value SXY(xi,yj) by a weighted sum over the value and its neighbours, where the weights may be e.g. a 2D gaussian?
– Andy
Jul 11 '12 at 12:48
• Sorry Jim. I had forgot to stress that I also want small RMS error relative to the model. I edited my question to take into account this. Smoothing gives a smoother result, but not a smaller RMS error?
– Andy
Jul 11 '12 at 13:21
• @Andy Yes, "a weighted sum over the value and its neighbors" is a pretty succinct description of filtering, and a 2D gaussian is one form of low-pass filter. Jul 11 '12 at 13:22
• @Andy For smaller error I would estimate the bandwidth by FFT'ing multiple "clean" (no noise added) $s$ results, and seeing where the highest frequency rolloff is (I'm assuming that they're not all the same). Design a LPF with that same rolloff- the Matlab "fdatool" can help with this- and then use that filter. It should improve your RMS. There will still be error, of course, but it should be reduced. Jul 11 '12 at 13:25
• Thanks Jim. But is there no way to combine the results from SXY and SYX in order to get a smaller RMS error?
– Andy
Jul 11 '12 at 13:58