I tried on my Laptop and got:

• Matlab: 0.179527 s
• Octave: 0.687583 s

So there might be some work to do, but I admit that I haven't tried any compilation optimizations. So, the bottleneck could also be partially elsewhere. Don't know.

Thanks for the corrections.

The for loop has got only 4 iterations with big operations.  I tried to vectorize it, but it was not easy to do.  I would have the shifts in the third axis, but when selecting rows (or columns) of Z the shape is not preserved.  Reshaping would require transposing in one of the cases, if I did not overlook something.  My conclusion was, that this is unlikely to be faster, but complicates the code quite a lot.

Another optimization possibility could be the padding.  Padding is not always necessary.  It is only required if there are points near the boundary.  So this is not so much about performance, but more of the memory footprint, since it might avoid one copy of `Z`.  Thinking further, the full copy could always be avoided, when dealing with the near-boundary-points separately.  Then only the three first or last rows or columns had to be copied and one row or column padded, which would be cheaper.  The special case of copying the last row or column could also be caught, instead of using the algorithm with double padding there.  However, these are a lot of special cases for the code and the improvement might be rather low.

Regarding your suggestion with conv2() or filter2():  No, only if the interpolation is a pure translation (and optional crop), it is a real convolution.  In this case `h` would contain the same value, since it is the distance to the point before.  Usually it is not a convolution, or, it may be considered as a convolution with a different kernel everywhere.  So it is just similar to a convolution and conv2() or filter2() cannot be applied.

Feel free to open a new report if there isn't one already.

Sorry. I lost track over which change was tracked where, and what was already applied or not...

I pushed your patch with some minor changes here:
https://hg.savannah.gnu.org/hgweb/octave/rev/98484425bd1b

It's already a very big help that those are "real" hg patches. Thank you for that.
It's a convention for GNU projects to summarize (very briefly) in the commit message what has changed in which file. It would be even more helpful if you could add that, too. (Not a requirement. It'd just help us out a bit when reviewing/pushing. 😊)

Regarding the "xtest"s: Afaict, those are checking for the results that Matlab produces for these test cases. I left them in as is in this change. But it would be nice if there was a bug report that would track those. When there is one, the "xtest"s should be replaced by "test" with that bug id.

Did the latest patch get reviewed?  I know there was a series of other patches from Christof dealing with _pchip_deriv_ and this bug report may have been forgotten.

Rebased on new patch of bug #61863.

(file #52707)

OK, I submitted a new bug report (bug #61863) for the complex matrix issue. While doing this, I found out, that "pchip" needed some extra care. I also added a patch there.

Then I rebased the bicubic interpolation fix on the patch of bug #61863 with some minor improvements compared to the submission in comment #2. The new patch is attached.

(file #52694)

I implemented bicubic interpolation in interp2.m, see attached file. To make the implementation easier to understand, I also attached a PDF (LaTeX source) to explain the approach. I also added tests, not only for the bicubic interpolation, but for the sake of consistency also for the existing methods. I found a few bugs:

• Interpolating complex matrixes gave an NA + 0i instead of NA + NAi. I could fix this.
• Spline interpolation does only support meshgrid format, but not index pairs. Matlab supports that. There is already a note in _splinen_.m. I added a failing test using xtest.
• Spline interpolation returns a row, when interpolating a single column, though with the correct values. I added a failing test using xtest.

I have not opened new bug reports nor checked for existing ones regarding the two spline interpolation issues.

Where to go from here?

• I have not checked interp1.m, interp3.m, interpn.m. The cubic interpolation might also be Matlab incompatible.
• Bilinear interpolation could also be implemented by convolution easily. This might be faster, maybe only for specific situations. When switching from interp2 to the convolutional approach in imresize.m to support anti-aliasing (submitted in bug #61760, not yet pushed), the speed-up was quite big for some tests.

Please let me know if the implementation needs improvements.

(file #52682, file #52683)

Using "cubic" for the method argument of interp2 is just an alias for "pchip" (piecewise cubic Hermite interpolating polynomial). However in Matlab a cubic convolution is used, as stated in their documentation of interp2. This method is also described on Wikipedia's Bicubic Interpolation page.

A tiny example to show the difference:
Octave:

> interp2(magic(4), 1, 1.5, 'cubic')
ans = 8.1875

Matlab (R2021a):

> interp2(magic(4), 1, 1.5, 'cubic')
ans = 8.6250

The topic has already been discussed in bug #60322 (regarding imresize) and there is even a patch for imresize available to reach Matlab compatibility. Maybe this patch can be adapted for interp2 as well? What is your opinion on this?