Two-dimensional linear interpolation with SOS2 variables

Suppose we want to implement an interpolation scheme for some 2D function \(z=f(x,y)\). We assume we have some data points \((x_i,y_j,z_{i,j})\) where the \((x_i,y_j)\) values form a grid. The grid points do not have to be equidistant, but the rows and columns must align:

All these three examples are valid for this approach.

A 2D SOS2 based interpolation scheme can be expressed as:

We force that up to four neighboring elements can be non-zero in the matrix \(\lambda_{i,j}\). This is accomplished by calculating row- and column-sums and then saying up to two consecutive elements of these vectors are nonzero. This scheme is rather simple but it actually works.

A problem

As indicated in the comments below, this formulation will not always give you unique solutions.

Consider the really small grid \(x=[0,1]\), \(y=[0,1]\) with the following data (using \(z=x\cdot y\)):

\[\begin{matrix}x&y&z\\ \hline 0&0&0\\0&1&0\\1&0&0\\1&1&1\end{matrix}\]

Our grid of \((x,y,z)\) values now looks like:

\[\left(\begin{matrix}(0,1,0)&(1,1,1)\\  (0,0,0)&(1,0,0)\end{matrix}\right)\]

Assume we want to find an approximation for \(f(0.5, 0.5)\). Using the above formulation we can see the following values:

\[\lambda=\left(\begin{matrix}0&0.5\\0.5&0\end{matrix}\right)\]	\(z=0.5\)
\[\lambda=\left(\begin{matrix}0.25&0.25\\0.25&0.25\end{matrix}\right)\]	\(z=0.25\)
\[\lambda=\left(\begin{matrix}0.5&0\\0&0.5\end{matrix}\right)\]	\(z=0\)

So depending on the situation we can find different interpolated values for \(z\). Obviously this is not very satisfactory. In the next paragraphs I will discuss some countermeasures we can take.

Triangles

Instead of interpolating between four points we can also interpolate between three points. Schematically our grid becomes:

Note that we also could have chosen diagonals that slope down (I discuss this in more detail later). If we have \(m\) \(x\)-values and \(n\) \(y\)-values, we end up with \(n+m-1\) diagonals. We can see this easily from:

Basically the idea is to select also two neighboring diagonals. This is again easily expressed using a SOS2 set. The \(\lambda\)’s along a diagonal are described by \(\lambda_{i,i+t}\) where \(t\) is a fixed shift or offset (note that \(t\) will be negative as well as zero and positive). Obviously, the diagonal that passes through \(\lambda_{1,1}\), \(\lambda_{2,2}\), etc., has \(t=0\). I.e. the condition to add is:

\[\boxed{\begin{align}&\textbf{sos2 variable }\beta_t\\&\beta_t = \sum_i \lambda_{i,i+t}\end{align}}\]

(1)

Does this help us with the non-uniqueness? The solutions with

\[\lambda=\left(\begin{matrix}0.25&0.25\\0.25&0.25\end{matrix}\right)\]	\(z=0.25\)
\[\lambda=\left(\begin{matrix}0.5&0\\0&0.5\end{matrix}\right)\]	\(z=0\)

are no longer feasible. We are left with only:

\[\lambda=\left(\begin{matrix}0&0.5\\0.5&0\end{matrix}\right)\]

\(z=0.5\)

Unfortunately this is not the best approximation. Note that the point \((x,y)=(0.5,0.5)\) is actually part of two possible triangles so we can see two possible configurations in cases like this.

\[\lambda=\left(\begin{matrix}0&0.5\\0.5&\end{matrix}\right)\]	\(z=0.5\)
\[\lambda=\left(\begin{matrix}&0.5\\0.5&0\end{matrix}\right)\]	\(z=0.5\)

These configurations have the same approximation. If a point is precisely on a diagonal it is part of two triangles, but as we have seen this is not an issue as only two \(\lambda_{i,j}\)’s are nonzero, making the interpolation unique again.

Modeling in GAMS

Equation (1) is not directly implementable in GAMS. If variable \(\beta_t\) is indexed by \(t\) and variable \(\lambda_{i,j}\) by \((i,j)\), we need to map between \(t\) and \((i,j)\). One possible way to do this is as follows:

sets  t  /t1*t3/  map(i,j,t) ; abort$(card(t)<>card(i)+card(j)-1) "set t has wrong size"; map(i,j,t)$(ord(j)=ord(i)+ord(t)-card(j)) = yes; display map; sos2 variable beta(t); equation defbeta(t); defbeta(t).. beta(t) =e= sum(map(i,j,t), lambda(i,j));

The advantage of using an auxiliary set map is that we can develop this independently of the equation and debug this before running actual the model. The resulting equation is now very simple.

For our small 2x2 example, the mapping set looks like:

----     50 SET map                 t1          t2          t3 i1.j1                     YES i1.j2                                 YES i2.j1         YES i2.j2                     YES

This corresponds to:

Downward sloping diagonals

Instead of upward sloping diagonals, we can use downward sloping ones instead.

 

The gams code to populate our mapping set ‘map’ for this case is not very complicated:

map(i,j,t)$(ord(i)+ord(j)=ord(t)-1) = yes;

Alternative

If we can make the function \(f(x,y)\) separable, we may be able to use a different scheme. E.g. assume \(f(x,y)=x\cdot y\) with \(x,y>0\). We can make this separable by applying a log transformation: \(\log(z)=\log(x)+\log(y)\). We can implement each of the functions \(z’=\log(z),x’=\log(x),y’=\log(y)\) using standard 1D piecewise linear formulations. Then we just add the constraint: \(z’=x’+y’\).

Of course another alternative is using a nonlinear solver.

References

1. H.P.Williams, “Model Building in Mathematical Programming”, 5th edition, Wiley, 2013.
2. Misener, R. & Floudas, C.A., “Piecewise-Linear Approximations of Multidimensional Functions”, J. Optim. Theory Appl. (2010), pp.145-120.
   In this paper SOS1 sets are used for the \(\delta\) and \(\gamma\) variables and SOS2 sets for the \(\beta\) variables. Sadly, Chris Floudas passed away last month while vacationing in Greece.
3. GAMS: Piecewise linear functions with SOS2 variables. Basic GAMS syntax and rules for SOS2 sets.
4. 2D Interpolation With SOS2 variables. Example GAMS model.