# Linear Interpolation

`⎕io=0` assumed throughout; works in 1-origin with the obvious modifications.

## Introduction

On Wednesday, a question arrived via Dyalog Support from an intern in Africa: If `M` is the matrix on the left, use linear interpolation to compute the result on the right.

``````   1 20         1 20
4 80         2 40
6 82         3 60
4 80
5 81
6 82
``````

## Linear Interpolation

Two points `(x0,y0)` and `(x1,y1)` specify a line; for any `x` there is a unique `y` on that line (assuming `x0≠x1`). The equation for the line derives as follows, starting from its slope `m`:

``````   m = (y1-y0) ÷ (x1-x0)
(y-y0) = m × (x-x0)
y = y0 + m × (x-x0)
``````

Therefore, if `⍺` is a 2-by-2 matrix of the two points and `⍵` are the x-values to be interpolated, then:

``````   g ← {(⊃⌽⍺)+(⍵-⊃⍺)÷÷/-⌿⍺}

⊢ M←1 4 6,⍪20 80 82
1 20
4 80
6 82

M[0 1;] g 2 3
40 60
M[1 2;] g 5
81
``````

## A New Twist, A New Solution

The problem as posed implicitly required that:

• The x-values are the positive integers bounded by `⊃⊖M`.
• Appropriate rows of the matrix are selected for a given x-value.
• The missing x-values and their interpolations are “slotted back” into the argument matrix.

These requirements are best met by `⍸`, interval index, a relatively new primitive function introduced in Dyalog APL version 16.0. The left argument `⍺` must be sorted and partitions the universe into disjoint contiguous intervals; `⍺⍸⍵` finds the index of the interval which contains an item of `⍵`. The result is `⎕io` dependent.

For the given matrix `M`, the partition (of the real numbers in this case) is depicted below. As in conventional mathematical notation, `[` denotes that the interval includes the left end-point and `)` denotes that the interval excludes the right end-point.

``````          1        4      6
─────────)[───────)[─────)[──────────
¯1       0       1       2

v←¯5 0 1 2.5 6 3 4 5 9 8 7

1 4 6 ⍸ v
¯1 ¯1 0 0 2 0 1 1 2 2 2

v ,[¯0.5] 1 4 6 ⍸ v
¯5  0 1 2.5 6 3 4 5 9 8 7
¯1 ¯1 0 0   2 0 1 1 2 2 2
``````

With `⍸` in hand, the problem can be solved as follows:

``````interpol←{
(x y)←↓⍉⍵
m←m,⊃⌽m←(2-/y)÷(2-/x)
j←0⌈x⍸i←1+⍳⊃⌽x
i,⍪y[j]+m[j]×i-x[j]
}

interpol M
1 20
2 40
3 60
4 80
5 81
6 82
``````

The problem of x-values less than the first end-point is finessed by applying `0⌈` to the interval indices, and that of x-values greater than or equal to the last end-point is finessed by repeating the last slope `m←m,⊃⌽m`.

It is possible to do the interpolation only on the missing indices (2 3 5 in this case) and insert them into the argument matrix. It seems neater to simply interpolate everything, and in so doing provide a check that the interpolated values equal the values given in the argument.

## An Alternative Interpolation

Interpolating according to two selected rows of a matrix of points treats the function as piecewise linear, with sharp inflection points where the lines join (different slopes between adjacent lines). A “holistic” alternative approach is possible: the matrix can be interpreted as specifying a single line and the interpolation is according to this single line. The `⌹` primitive function computes the coefficients of the line which best fits the points:

``````   ⎕rl←7*5  ⍝ for reproducible random numbers

⊢ M←t,⍪(?7⍴5)+¯17+3×t←?7⍴100
35  89
98 278
19  44
4  ¯5
62 170
49 133
25  59

M[;1] ⌹ 1,M[;,0]    ⍝ y-intercept and slope
¯15.3164 2.99731

interpola ← {(1,⍤0⊢⍵)+.×⍺[;1]⌹1,⍺[;,0]}

M[;1] ,[¯0.5] M interpola M[;0]
89      278    44      ¯5       170     133     59
89.5895 278.42 41.6325 ¯3.32713 170.517 131.552 59.6164

M interpola 33 35 37 39.7
83.5949 89.5895 95.5841 103.677
``````

## Finally

Our best wishes to the intern. Welcome to APL!