Здравствуйте! Вопрос касается трикубической интерполяции (на регулярной сетке). Возможно, кто-то уже разбирался с этим и подскажет мне оптимальное решение. Я так понимаю, что существует два варианта  - посчитать все константы полинома и запомнить их (Lekien and Marsden, Tricubic Interpolation in Three Dimensions, 2005) или хранить только значения в вершинах и интерполировать значения «на лету» (http://en.wikipedia.org/wiki/Tricubic_interpolation). Для случая, когда надо получить только значение (интерполированное) функции, второй вариант вроде более привлекателен в смысле быстродействия (затраты памяти очевидно у второго метода всегда в 64 раза меньше).
Вопрос первый, в быстродействии возникает, если надо не только значение, но и производные функции (опять же из интерполированного полинома). Для случая честных значений констант полинома это весит около 1.5К умножений/суммирований. А кто делал по (второму) варианту со скалярными произведениями? Целесообразно ли персчитать дискретные точки в явный полином (через численные градиенты)?
Второй вопрос в точности, если исходная функция аналитическая и доступны ее производные. Предполагается, что функция достаточно гладкая для выбранного шага решетки, т.е. практически не более 1 перемены знака на кубик решетки. Сильно ли пострадает точность интерполяции при замене всех аналитических производных на интерполяцию по вершинам? Меня, конкретно, интересует применимость посчитаны градиентов для квазиньютоновских методов безусловной оптимизации.
П.С. Скорость построения самой решетки, да и в принципе ее размер в памяти для моей задачи не существенны. Важна скорость вычисления значений и градиентов. Ну и точность тоже.
Спасибо!

• Tricubic interpolation,ghost_in_machine, 00:03 , 18-Окт-10
  • Tricubic interpolation,ghost_in_machine, 00:37 , 18-Окт-10
  • Tricubic interpolation,AntonV, 15:25 , 01-Июн-11
    • Tricubic interpolation,yakovenko_aukr.net, 21:03 , 01-Июн-11
• Tricubic interpolation,pavlinux, 14:40 , 02-Июн-11
  • Tricubic interpolation,yakovenko_aukr.net, 04:14 , 03-Июн-11

Самому разобраться оказалось проще. Оставлю здесь выводы, может кому пригодяться (сам к примеру буду искать через год). Вообщем по варианту вычисления из вики считать удобно как f=CINT(x)^T.Bz.CINT(y), где Bz матрица 4х4 сформированнаяиз 16 вызовов CINT(z). Потому df/dx=dCINT(x)^T.Bz.CINT(y) и df/dy=CINT(x)^T.Bz.dCINT(y); df/dz удобнее формировать с нуля: 16 вызовами dCINT(z) получить матрицу dBz и потом df/dz=CINT(x)^T. dBz.CINT(y). Всего для вычисления на лету интерполированного значения и производных надо два по 16 произведений 4-вектора плюс умножения вектор-матрица-вектор ну и сами вектора т.е. где-то в 2 раза быстрее чем вариант из статьи. Ну точность поменьше, но в 16 раз меньше памяти (для 3D grid!!!) и 2 раза быстрее считать, ИМХО приемлемая компенсация даже для случая аналитических производных табулированной функции.
З.Ы. шоб потом не думать совсем dCINT(z)^T={ (4-3*z)*z-1, z*(9*z-10), (8-9*z)*z+1, z*(3*z-2) }
З.З.Ы. для своего проекта все равно оставил вариант из статьи везде, где есть аналитические производные :)))))

простите, не туда посмотрел, разница в флопсах всего ~10% 460 против 520

> Самому разобраться оказалось проще. Оставлю здесь выводы, может кому пригодяться

Спасибо большое. Пригодилось.

>> Самому разобраться оказалось проще. Оставлю здесь выводы, может кому пригодяться
> Спасибо большое. Пригодилось.

Раз люди пользуют то добавлю небольшой update:
1. Carlson R.E. and Fritsch F.N. MONOTONE PIECEWISE BICUBIC INTERPOLATION SIAM J. NUMER. ANAL. 22, 386-400, 1985. Тут написано как сделать решетку с сохнением знака полинома на интервалах.
2. Hyman J.M. ACCURATE MONOTONICITY PRESERVING CUBIC INTERPOLATION, SIAM J.Sci.Stat.Comput., 4, 645-654, 1983. Тут написано как сделать полином непрерывно растущим/спадающим на интервале (зафризить знак производной).

Пользуйте на здоровье :)

>  Важна скорость вычисления значений и градиентов.

Нынче, всё что больше 2-х измерений пихают в Cuda/OpenCL

>>  Важна скорость вычисления значений и градиентов.
> Нынче, всё что больше 2-х измерений пихают в Cuda/OpenCL

А как можно распараллелить многомерную кубическую интерполяцию (если не формировать матрицы в явном виде и параллелить уже матричное умножение)? Ведь каждый полином следующего измерения строиться на 4-х интеполированных точках пердыдущей размерности, а сама процедура интерполирования занимает всего несколько умножений/сложений. Вроде как затраты на синхронизацию большие (количество задействованых процессов уменьшаеться в четире раза и надо дождатся завершения всех 4-х перед следующим вычислением), а на собственно вычисления - еруендовые. Ладно бы измерений было с 10 (но такие гриды даже в современные размеры памяти не влезут) тогда можно параллелить по измерения, но ведь их всего 3 (ну максимум 4) в практических задачах..?