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!