bilinear extrapolation based on interp2

[X,Y] = meshgrid(0:10);
Z = X.^2 + Y.^2;
[Xq,Yq] = meshgrid(0:0.25,10);
V = interp2(X,Y,Z,Xq,Vq,'linear');

I want to use interp2 to (bi)linearly interpolate the function Z(x,y) in the interior of the grid.

Although I have to take extrapolating with a pinch of salt, I want to come up with a bilinear extrapolation in case I have to evaluate outside the grid.

Say is an evaluation point outside the predefined grid, and is its closest point on the boundary. Then, the extrapoland is given by

where is the value of Z at t, similarly for the partial derivatives.

I would compute the derivative of the extrapoland as

Is this gradient correct? I am uncertain about that since t, the closest point on the boundary, is also a function of (x,y). But I do not know how would one differentiate that.

[Torsten]

If Z is only available on a meshgrid, which kind of approximation of the gradient dZ/dx and dZ/dy do you want to use in (t_x,t_y) ?

I think it's more reasonable to extrapolate the value in (x,y) outside the grid as p(t_x,t_y).

[SA-W]

"If Z is only available on a meshgrid, which kind of approximation of the gradient dZ/dx and dZ/dy do you want to use in (t_x,t_y) ?"

Thats a different question. Probably I assume a bilinear form of Z between the grid points such that I can compute (discontinuous) derivatives of Z. But lets assume I can compute the gradient everywhere inside the grid...

"I think it's more reasonable to extrapolate the value in (x,y) outside the grid as p(t_x,t_y)."

What do you mean by that?

[Torsten]

But lets assume I can compute the gradient everywhere inside the grid...

You need the gradient in the boundary point (t_x,t_y).

"I think it's more reasonable to extrapolate the value in (x,y) outside the grid as p(t_x,t_y)."

What do you mean by that?

I mean it's better to use a constant instead of a linear extrapolation. If a point (x,y) outside the grid domain is far apart, each kind of extrapolation tends to give unphysical results. And constant extrapolation is the petty evil.

[SA-W]

I can also compute the gradient at the boundary.

Of course, I assume that Z is quite smooth and p is close to the boundary to get something useful from linear extrapolation.

Given the extrapoland as stated in my question, do you think my gradient is correct?

[Torsten]

Z_extrap(x,y) = Z(t_x,t_y) + dZ/dx (t_x,t_y) * (x - t_x) + dZ/dy (t_x,t_y) * (y - t_y)

That's why I asked about how you want to approximate the gradients in the boundary point (t_x,t_y).

But as already said: I wouldn't use linear extrapolation.

[SA-W]

I would really to give that a try. But I need both the function value at the gradient.

So how would you compute the gradient of

Z_extrap(x,y) = Z(t_x,t_y) + dZ/dx (t_x,t_y) * (x - t_x) + dZ/dy (t_x,t_y) * (y - t_y)

[Torsten]

The gradient of Z_extrap or the dZ/dx and dZ/dy terms within the expression for Z_extrap ?

[SA-W]

"The gradient of Z_extrap or the dZ/dx and dZ/dy terms within the expression for Z_extrap ?"

The gradient of Z_extrap, i.e, the derivatives dZ_extrap/dx and dZ_extrap/dy.

[Torsten]

Determine the nearest point (t_x,t_y) to (x,y) on the boundary.

Depending on whether (t_x,t_y) is a corner point, a point on the vertical or a point on the horizontal boundary line, choose one-sided or central differences to approximate dZ/dx(t_x,t_y) and dZ/dy(t_x,t_y). If, e.g., (t_x,t_y) is on the lower horizontal boundary (not a corner point), approximate dZ/dx(t_x,t_y) = (Z(t_x+1,t_y)-Z(t_x-1,t_y))/(2*deltax) and dZ/dy(t_x,t_y) = (Z(t_x,t_y)-Z(t_x,t_y-1))/deltay (or with a one-sided difference approximation of higher order, if you like).

[SA-W]

As I said, dZ/dx and dZ/dy is known, also at the boundary, i.e. I know dZ/dx(t_x,t_y) and dZ/dy(t_x,t_y). So there is no need to do finite differencing.

But I think the issue is that Z(t_x,t_y), dZ/dx(t_x,t_y), and dZ/dy(t_x,t_y) are itself all functions of x and y because t_x,t_y depend on x and y.

[Torsten]

But I think the issue is that Z(t_x,t_y), dZ/dx(t_x,t_y), and dZ/dy(t_x,t_y) are itself all functions of x and y because t_x,t_y depend on x and y.

No. It's a Taylor approximation of Z_extrap(x,y) of first-order with center point (t_x,t_y), the point nearest to (x,y).

[SA-W]

If I go from (x,y) to (x+ deltaX,y), then t_x and t_y may be different, hence also Z(t_x,t_y), dZ/dx(t_x,t_y), and dZ/dy(t_x,t_y) change with x in a discontinuous manner.

But I do not know. Let's see what MattJ can say about that.

Answer 1

[Matt J]

Although you've been advised against it, if you truly must extrapolate bilinearly then griddedInterpolant will do it for you,

[X,Y] = deal(0:10);
[Xq,Yq] = deal(0:0.25,10);

Z = X.^2 + (Y').^2;

F=griddedInterpolant({X,Y}, Z,'linear','linear');

V=F({Xq,Yq})

[SA-W]

As I stated in the comments, I also must be able to compute derivatives of the interpolated/extrapolated function. I think there is no matlab built-in that allows me also to query those derivatives from an extrapolation scheme.

I think the gradient of Z_extrap as defined in my question is discontinous as the closest point at the boundary also depends on p=(x,y). What do you mean?

[Matt J]

The derivative of,

F=griddedInterpolant({X,Y}, Z,'linear','linear');

would be,

F=griddedInterpolant({X,Y}, gradient(Z),'previous','nearest');

although I would probably use linear interpolation instead in most practical scenarios,

F=griddedInterpolant({X,Y}, gradient(Z),'linear','nearest');

[SA-W]

Makes sense. Anyway, I also want to explore

Z_extrap(x,y) = Z(t_x,t_y) + dZ/dx (t_x,t_y) * (x - t_x) + dZ/dy (t_x,t_y) * (y - t_y)

Are dZ_extrap/dx and dZ_extrap/dy continuous and is it correct how I computed them in my question?

[Matt J]

The derivative is not correct, for reasons you already mentioned in your post. But more importantly, the extrapolation scheme you've defined is not bilinear, which you said is what you wanted.

[Torsten]

@Matt J

If you define

Z_extrap(x,y) = Z(t_x,t_y) + dZ/dx (t_x,t_y) * (x - t_x) + dZ/dy (t_x,t_y) * (y - t_y)

why do you think

dZ_extrap/dx (t_x,t_y) = dZ/dx (t_x,t_y)
dZ_extrap/dy (t_x,t_y) = dZ/dy (t_x,t_y)

is not correct ?

[Matt J]

@Torsten Because the dependence of t_x and t_y on (x,y) is not accounted for in the derivative.

[SA-W]

@Matt J Say I would choose t = (t_x, t_y) to be always the same point regardless of x and y, would

dZ_extrap/dx (t_x,t_y) = dZ/dx (t_x,t_y)

dZ_extrap/dy (t_x,t_y) = dZ/dy (t_x,t_y)

be correct then?

And to make the extrapolation bilinear, I would have to add a mixed term x*y?

[SA-W]

@Matt J I have seen you placed a comment here, but deleted it again probably accidentally?

I think you wrote that

dZ_extrap/dx (t_x,t_y) = dZ/dx (t_x,t_y)

dZ_extrap/dy (t_x,t_y) = dZ/dy (t_x,t_y)

is the correct, although (t_x, t_y) depends on (x,y). Is that true?

[Matt J]

No, I didn't delete it accidentally. I withdrew it. I've changed my mind about everything, including my original answer.

I believe the computation of the derivatives for the extrapolation defined in your original post would be as below. That is, the form of the derivative will depend on which of 8 regions around the domain (x,y) is located in (I show only 6 of them). In particular, the derivatives over the northeast, northwest, southeast, and southwest regions will be constant. The appealing thing here is that the derivatives can be computed by a combination of nearest-neighbor extrapolation of the first derivatives and linear extrapolation of the second derivatives.

[SA-W]

Interesting. Thanks for this nice picture. Would you share your thoughts why the term shows up only in some regions and. The methodology of finding is the same in all regions. So I would have expected to have the additional terms in all regions.

[Matt J]

In your original post, you proposed formulas for the derivatives based on the assumption that (tx,ty) are constant independent of (x,y). This assumption, and therefore the formulas, is valid in those external corner regions south/north-east/west of the rectangular domain. The nearest point will be the nearest corner of the rectangle regardless of the position of (x,y) within those regions.

[SA-W]

Makes sense, also

in the two "middle regions" shown in your graphic since the nearest point in those regions does not change with x. I think you obtained

in the same regions by applying chain rule to , but why are the remaining two summands and independent of y?

[Matt J]

They aren't.

[SA-W]

Then, why are there the derivatives and not showing up in your formula?

[Matt J]

Ah. They should.

And in fact (y-t_y) is always zero in that region, so you should probably discard the other two terms I had.

[SA-W]

Ah. They should.

So, the graphic you created some comments ago is not correct in those regions?

And in fact (y-t_y) is always zero in that region, so you should probably discard the other two terms I had.

Why is (y-t_y = 0) there? This would mean that the closest point to (x,y) has the same y-coordinates, which is most likely not the case.

[Matt J]

So, the graphic you created some comments ago is not correct in those regions?

Yes.

Why is (y-t_y = 0) there? This would mean that the closest point to (x,y) has the same y-coordinates, which is most likely not the case.

It is the case in the regions to the right and left of the rectangle. Similarly x=t_x in the regions above and below.

[SA-W]

It is the case in the regions to the right and left of the rectangle

Why? Say either (t_x1, t_y1) or (t_x2, t_y2) is the closest point on the boundary of the grid. The y-coordinate of (x,y) can be anywhere in the interval [t_y1, t_y2]

[Matt J]

The nearest point (tx,ty) from a point (x,y) to the right (or left) of the rectangle is its orthogonal projection on the rectangle. An orthogonal line drawn from (x,y) to the rectangle will be perfectly horizontal, and therefore satisfy t_y=y.

[Matt J]

We can also see it from the explicit formulae for t_x and t_y

t_x=max(min(x,xmax), xmin)
t_y=max(min(y,ymax), ymin)

Clearly whever ymin<=y<=ymax we will have t_y=y.

[SA-W]

Oh, we talked at cross purposes. What I meant by closest point (but not precisely explained) is the closest vertex at the boundary.

Does that change a lot in the formulaes left and right to the rectangle? I think the equations for northeast/northwest/southeast/southwest are still correct.

[Matt J]

It changes it, and I don't really see the rationale for using the closest vertex. That will create a discontinuity at the boundaries of the rectangle. The first order Taylor expansion of Z taken at an arbitrary point just inside the boundaries of the rectangle will be very different from the first order Taylor expansion at the nearest corner.

[SA-W]

I will probably do the extrapolation outside matlab and finding the closest vertex is easier there.

The first order Taylor expansion of Z taken at an arbitrary point just inside the boundaries of the rectangle will be very different from the first order Taylor expansion at the nearest corner.

The first degree taylor expansion of Z at is

Of course, is different for the closest vertex than for the closest point. Consequently, the partial derivtives and Z in the above are evaluated at a different point. But that should be all if I look at the above equation. Do I miss something here?

[Matt J]

I will probably do the extrapolation outside matlab and finding the closest vertex is easier there.

That seems doubtful. The closest point has a very simple closed-form formula, which I gave above.

Consequently, the partial derivtives and Z in the above are evaluated at a different point. But that should be all if I look at the above equation. Do I miss something here?

Yes. You seem to have ignored,

(1) The affect of changing (tx,ty) on the multiplicative factors (x-tx) and (y-ty). They can become much larger for the closest vertex than the closest point.

(2) That changing the partial derivatives by a large amount is also a big deal.

[SA-W]

That seems doubtful. The closest point has a very simple closed-form formula, which I gave above.

This formula is true only for a rectangle, but I want to implement a try to define Z(x,y) (actually, it is the f(x,y) we are talking about in a different chat), on a convex hull of a set of points, i.e. the boundary of Z(x,y) is a polygon.

You seem to have ignored,

(1) The affect of changing (tx,ty) on the multiplicative factors (x-tx) and (y-ty). They can become much larger for the closest vertex than the closest point.

(2) That changing the partial derivatives by a large amount is also a big deal.

True. So given

where is the closest point at the boundary. We derived that

in southeast, southwest, northeast, northwest of the rectangle.

In the region between southeast and northeast, is also not a function of x, i.e.,

but is a function of y. So is

true?

[Matt J]

It's true, but the analysis is going to be moot on an arbitrary convex polygon. The breakdown of regions and the derivative formulas applicable to them become much more complicated in that case.

Also, if we are really talking about the convex function f(x,y) from your other posts, I should point out that none of the 2D extrapolation schemes we've discussed above, whether standard bilinear or your proposal, preserve convexity.

Finally, you said that in many cases f will have the separable form f(x,y)=h(x)+g(y). In that case, the extrapolation problem boils down to 1D extrapolations of h and g, which seems trivial, both in the computation of f and of the extrapolated derivatives.

[SA-W]

Also, if we are really talking about the convex function f(x,y) from your other posts, I should point out that none of the 2D extrapolation schemes we've discussed above, whether standard bilinear or your proposal, preserve convexity.

I know. But the extrapolation should be used very seldomly compared to in-domain evaluations. I just have to "offer" it, because if not, my pde solver must return NaN, hence my objective function. The question arises why not simple increasing the rectangle (the domain). Well, then I have the issue of insufficient parameter activation that we also addressed in the other post. But maybe I can try a bigger domain with the new basis that you suggested.

It's true, but the analysis is going to be moot on an arbitrary convex polygon.

Yes, lets forget about the polygon. But is it also complicated to derive the formulae between southeast and northeast in case we want to be the closest vertex?

[Matt J]

Yes, lets forget about the polygon. But is it also complicated to derive the formulae between southeast and northeast in case we want to be the closest vertex?

As I said in my last comment, if f(x,y)=h(x)+g(y) then this is a 1D extrapoltion problem in which case "nearest point" and "nearest vertex" mean the same thing.

[SA-W]

Yes, but if f(x,y) is not separable in x and y?

To provide a proof-of-concept, I can not always assume f(x,y)=g(x)+h(y). Although most analytical versions of f are of this form, f is in general not known. And since it is not known and we are looking for better models of f, the assumption f(x,y)=g(x)+h(y) is too strong in my opinion and can only be justified for pure re-identification problems where f(x,y) is known a priori.

So I would appreciate if we can derive the derivatives (using closest vertex instead of closest point) in case f is defined in a rectangular grid.

[Matt J]

You would have to draw the Voronoi diagram of the vertices. By definition of the Voronoi diagram, tx and ty will be constant within each Voronoi cell and you can take derivatives there while freely assuming tx and ty to be constant.

Note that taking the derivatives on the boundaries of the cells will in general not be possible. Since the extrapolated function is potentially discontinuous there, as we've discussed, it will likewise be non-differentiable.

[SA-W]

I see. The different Veronoi cells probably translate into a bunch of if-statements in my pde solver, hence too much overhead.

As for the region of the rectangle between northeast and northwest, are the derivatives given by

?

[Matt J]

I've lost track of what we're talking about. Are we still talking about the nearest-vertex scheme? If so, the "region between northeast and northwest" is not applicable anymore. As mentioned in my last comment, it's the Voronoi cells that are applicable.

[SA-W]

Sorry, my equations in the last comment pertain to the closest point strayegy.

[Matt J]

Yes, for closest point, the formulas look correct to me. It should be easy to test them in any case.

[SA-W]

How would you test them? With finite differences?

[Matt J]

yes, with finite differences.

Answer 2

[John D'Errico]

Ugh. This is just a bad idea, Extrapolation in general is a bad thing to do. A better word is probably excrapolation. You will get exactly what you should expect from the name I just made up. There will be difficulties in trying to extrapolate, because if you try to extrapolate one cell, it is not going to be consistent with the extrapolant from the cell immediately next to it.

It is FAR better to use a tool that can at least try to intelligently extrapolate, using the shape of the surface near the boundary. In this case, I'll suggest my inpaint_nans tool, as found on the file exchange.

A = NaN(100,100);
A(30:70,30:70) = sin((X+Y)/5);
Ahat = inpaint_nans(A);
surf(Ahat)
hold on
H = surf(A);
H.FaceColor = 'r';

I doubt you can do too much better than that. It is actually not too bad, considering the extent of the extrapolation. (Yes, one of the things on my lengthy list of round-tuits is a better version of inpaint_nans for problems like this. But inpaint_nans is nearly 20 years old, and I fleshed out the idea for that improvement 20 years ago. Sigh.)

https://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint_nans?s_tid=ta_fx_results

[SA-W]

Thank you!

I am aware of your nice FEX tool inpaint_nans. But as stated in the comments under the question, I need both function values and the first derivatives of the extrapoland which I can not easily do with your tool.

So I think I have to come up with an extrapolation scheme that gives me more control.

I think I provided one of the simplest extrapolation scheme in my question.

[John D'Errico]

The derivatives that you compute will not be very good, in the sense that they are themselves just based on local finite differences. And then you base the extrapolant on those local derivatives. Do you see this is a bad idea?

I'll strongly argue that you FIRST perform the extrapolation using a good tool like inpaint_nans. and THEN you can compute a derivative, from the results. This is better because inpaint_nans implicitly uses that same derivative information around the perimeter to infer the shape of the extrapolated surface. But it also uses all of that information at once, to infer a SMOOTHLY behaved extrapolant.

You are trying to solve the problem from the wrong direction. Hey, it is your choice of course. Of course, if you already know how to solve the problem, then why did you ask this question in the first place?