bug #56665: Incorrect results with ppval and interp1 when using complex query points

 Submitted by: Chad Oldfield Submitted on: Wed 24 Jul 2019 08:50:41 PM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Feature Request Status: Confirmed Assigned to: None Originator Name: Chad Oldfield Open/Closed: Open Release: 5.1.0 Operating System: Any

You are not logged in

Discussion

Thu 25 Jul 2019 06:14:10 PM UTC, comment #4:

The result of interp1 is not a differentiable function. You should not ask for the derivative between consecutive intervals. A part from those points, it is a polynomial which is safe to evaluate at complex points, with the note that the interval of reference is chosen according to the real part of the query point.

I think that the hints in the original submission should fix the problems there observed.

Marco Caliari <caliari>
Thu 25 Jul 2019 04:58:03 PM UTC, comment #3:

Agreed that it is debatable whether this is an "incorrect result" or a desired/undesired extension.

If this were to be done, the key to enabling the complex step method to work is that the program should follow the same path (i.e. through if statements, selection of indices, etc.) with the complex value x + 1i * step as it would with just x.  The concept scales from small problems like interpolation to large problems like computational fluid dynamics.  For the interpolation and piecewise polynomial evaluation, it means choosing the interval based on the real part of xbar.  It's possible this may impact the lookup function.

Would this be mathematically unsound for other uses?

My humble opinion in favour of adding this extension is that, as an environment primarily for numerical analysis, it is helpful to support accurate numerical methods.  The complex step method is great for numerical accuracy (as compared to finite differences), because it eliminations the cancellation error, so very small steps can be used to limit truncation errors.  As an example (not involving interpolation), the subtractive cancellation errors cause an inaccurate finite difference in the following derivative estimate with too small a step size, but the analytic version and complex step agree:

>> step = 1e-14;
>> (sin(pi/6 + step) - sin(pi/6 - step)) / (2 * step)
ans =  0.86320
>> cos(pi/6)
ans =  0.86603
>> imag(sin(pi/6 + 1i * step)) / step
ans =  0.86603

So this would allow better numerical accuracy in applications that include interpolation or piecewise polynomials as a part of a larger function that is to be differentiated numerically.

Thu 25 Jul 2019 04:28:37 PM UTC, comment #2:

I tried the sample code within Matlab and it throws nasty errors.  Octave would be no worse than Matlab if it were to detect and disallow complex query values.  Documentation is at https://www.mathworks.com/help/matlab/ref/interp1.html#btwp6lt-1-xq.

Specifically,

Query points, specified as a scalar, vector, matrix, or array of real numbers.

On the other hand, we also view Octave as a superset of Matlab.  If there were a way to add this capability that was mathematically sound I would be in favor of that.

Rik <rik5>
Thu 25 Jul 2019 03:00:41 PM UTC, comment #1:

pp = interp1([1 2], [2 4], 'pp')

is the polynomial p = @(x) 2*x and it is true that p(1.5+0.0i) should be 3+0.02i. But in a more general situation, pp is a piecewise polynomial function. When you want to evaluate it at a target point xbar, then ppval has to find the interval in which xbar lies and evaluate the corresponding polynomial. How do you identify the interval if xbar in complex? The complex step method works for analytical functions.

Matlab does not allow to evaluate pp at a complex point. I think octave should not allow, too.

Marco Caliari <caliari>
Wed 24 Jul 2019 08:50:41 PM UTC, original submission:

The following is a copy of the command window, obtained right after freshly starting Octave 5.1.0:

>> interp1([1 2], [2 4], 1.5 + 0.01i)

ans =  3.000000 - 0.020000i

>> pp = interp1([1 2], [2 4], 'pp');
>> ppval(pp, 1.5 + 0.01i)

ans =  3.000000 - 0.020000i

>> interp1([-2 2], [-4 4], 1.5 + 0.01i)

ans =   NA

The expected output for all of these has the opposite sign for the imaginary part:
ans =  3.000000 + 0.020000i

This can be correctly obtained, for example, using:

>> p = polyfit([1 2], [2 4], 1);
>> polyval(p, 1.5 + 0.01i)

The background behind this example is that it is convenient to use the complex step method to numerically calculate the derivative of a complicated function.  Such a complicated function may include interpolation (among various other steps and algorithms).  The method is described in the link below.  Much of the utility of this is that it avoids some of the numerical errors involved with finite differences.

Essentially, the derivative of the real-valued function f can be approximated by:

dfdx = imag(f(x + 1i * step)) / step

where step is a small value (eg. 1e-20).  The expected complex part of the results of all of the tests above is therefore 0.01 times the slope of the line, or 0.02.

In the examples below, the first two (interp1 and ppval) may be due to the use of ppval by interp1, and the use of the transpose operator at line 71 of ppval.m:

dx = (xi - x(idx))(:)';

It is possible that applying the array transpose would help:

dx = (xi - x(idx))(:).';

In the third example, it is possible that the NA result is due to the use of the complex magnitude in the comparison operators in line 380 of interp1.m:

outliers = (xi < minx) | ! (xi <= maxx);  # this even catches NaNs

When using the complex step method, the real-valued function needs to be written such that the same path is followed with a real vs. a complex input.  This could be achieved, for example, with:

outliers = (real(xi) < minx) | ! (real(xi) <= maxx);

Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

No files currently attached

Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by rik5 (Posted a comment)
• -email is unavailable- added by caliari (Posted a comment)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.