Полигармонический сплайн

В прикладной математике полигармонические сплайны используются для аппроксимации функций и интерполяции данных . Они очень полезны для интерполяции и подбора разрозненных данных во многих измерениях. Особые случаи включают тонкие пластинчатые шлицы. и натуральные кубические сплайны в одном измерении.

Определение

Полигармонический сплайн представляет собой линейную комбинацию полигармонических радиальных базисных функций (RBF), обозначаемых $\varphi$ плюс полиномиальный член:

$f(\mathbf {x} )=\sum _{i=1}^{N}w_{i}\varphi (|\mathbf {x} -\mathbf {c} _{i}|)+\mathbf {v} ^{\textrm {T}}{\begin{bmatrix}1\\\mathbf {x} \end{bmatrix}}$ ( 1 )

где

Полигармонические базисные функции

• $\mathbf {x} =[x_{1}\ x_{2}\ \cdots \ x_{d}]^{\textrm {T}}$ ( ${\textrm {T}}$ обозначает транспонирование матрицы, что означает $\mathbf {x}$ — вектор-столбец) — вектор с действительным знаком $d$ независимые переменные,
• $\mathbf {c} _{i}=[c_{i,1}\ c_{i,2}\ \cdots \ c_{i,d}]^{\textrm {T}}$ являются $N$ векторы того же размера, что и $\mathbf {x}$ (часто называемые центрами), которые должна интерполировать кривая или поверхность,
• $\mathbf {w} =[w_{1}\ w_{2}\ \cdots \ w_{N}]^{\textrm {T}}$ являются $N$ веса RBF,
• $\mathbf {v} =[v_{1}\ v_{2}\ \cdots \ v_{d+1}]^{\textrm {T}}$ являются $d+1$ веса полинома.

Полином с коэффициентами $\mathbf {v}$ повышает точность подгонки полигармонических сглаживающих сплайнов, а также улучшает экстраполяцию от центров $\mathbf {c} _{i}.$ На рисунке ниже показано сравнение сплайнов с полиномиальным членом и без полиномиального члена.

Полигармонические RBF имеют вид:

${\begin{aligned}\varphi (r)&={\begin{cases}r^{k}&{\text{with }}k=1,3,5,\ldots ,\\r^{k}\ln(r)&{\text{with }}k=2,4,6,\ldots \end{cases}}\\[5mm]r&=|\mathbf {x} -\mathbf {c} _{i}|={\sqrt {(\mathbf {x} -\mathbf {c} _{i})^{\mathrm {T} }\,(\mathbf {x} -\mathbf {c} _{i})}}.\end{aligned}}$

Другие значения показателя степени $k$ бесполезны (например, $\varphi (r)=r^{2}$), поскольку решения проблемы интерполяции может не существовать. Чтобы избежать проблем при $r=0$ (с $\log(0)=-\infty$), полигармонические RBF с натуральным логарифмом могут быть реализованы как:

$\varphi (r)={\begin{cases}r^{k-1}\ln(r^{r})&{\text{for }}r<1,\quad {\text{(this works because }}0^{0}{\text{ is defined)}}\\r^{k}\ln(r)&{\text{for }}r\geq 1.\end{cases}}$

или, проще говоря, добавление расширения непрерывности в $r=0$

$\varphi (r)={\begin{cases}0&{\text{for }}r<\epsilon ,\quad {\text{(for some very small value of }}\epsilon {\text{, e.g. if using floating point numbers in double precisions, }}\epsilon =10^{-200}{\text{)}}\\r^{k}\ln(r)&{\text{for }}r\geq \epsilon .\end{cases}}$

Веса $w_{i}$ и $v_{j}$ определяются так, что функция интерполирует $N$ заданные баллы $(\mathbf {c} _{i},f_{i})$ (для $i=1,2,\ldots ,N$) и выполняет $d+1$ условия ортогональности

$\sum _{i=1}^{N}w_{i}=0,\;\;\sum _{i=1}^{N}w_{i}\mathbf {c} _{i}=\mathbf {0} .$

В совокупности эти ограничения эквивалентны симметричной линейной системе уравнений

${\begin{bmatrix}A&B\\B^{\textrm {T}}&0\end{bmatrix}}\;{\begin{bmatrix}\mathbf {w} \\\mathbf {v} \end{bmatrix}}\;=\;{\begin{bmatrix}\mathbf {f} \\\mathbf {0} \end{bmatrix}}\;\;\;\;$ ( 2 )

где

$A_{i,j}=\varphi (|\mathbf {c} _{i}-\mathbf {c} _{j}|),\quad B={\begin{bmatrix}1&1&\cdots &1\\\mathbf {c} _{1}&\mathbf {c} _{2}&\cdots &\mathbf {c} _{N}\end{bmatrix}}^{\textrm {T}},\quad \mathbf {f} =[f_{1},f_{2},\ldots ,f_{N}]^{\textrm {T}}.$

Чтобы эта система уравнений имела единственное решение, необходимо $B$ должен быть полного ранга. $B$ имеет полный ранг для очень мягких условий входных данных. Например, в двух измерениях три центра, образующие невырожденный треугольник, гарантируют, что $B$ имеет полный ранг, а в трех измерениях четыре центра, образующие невырожденный тетраэдр, гарантируют, что B имеет полный ранг. Как будет объяснено позже, линейное преобразование, возникающее в результате ограничения области линейного преобразования $A$ в нулевое пространство ${\textstyle B^{\textrm {T}}}$ является положительно определенным. Это означает, что если $B$ является полным рангом, система уравнений ( 2 ) всегда имеет единственное решение и ее можно решить с помощью линейного решателя, специализирующегося на симметричных матрицах. Вычисленные веса позволяют оценить сплайн для любого $\mathbf {x} \in \mathbb {R} ^{d}$ используя уравнение ( 1 ). Многие практические детали реализации и использования полигармонических сплайнов объяснены в Фассхауэре. В Иске полигармонические сплайны рассматриваются как частные случаи других методов мультиразрешения при моделировании разбросанных данных.

Обсуждение

Основное преимущество интерполяции полигармоническим сплайном состоит в том, что обычно очень хорошие результаты интерполяции получаются для разрозненных данных без выполнения какой-либо «настройки», поэтому возможна автоматическая интерполяция. Это не относится к другим радиальным базисным функциям. Например, функция Гаусса $e^{-k\cdot r^{2}}$ необходимо настроить так, чтобы $k$ выбирается в соответствии с базовой сеткой независимых переменных. Если эта сетка неоднородна, то правильный выбор $k$ добиться хорошего результата интерполяции сложно или невозможно.

Основные недостатки:

• Для определения весов необходимо решить плотную линейную систему уравнений. Решение плотной линейной системы становится непрактичным, если размерность $N$ большой, поскольку требуемая память $O(N^{2})$ а количество требуемых операций равно $O(N^{3}).$
• Оценка вычисленной функции полигармонического сплайна при $M$ точки данных требуют $O(MN)$ операции. Во многих приложениях (примером является обработка изображений) $M$ намного больше, чем $N,$ и если оба числа большие, это непрактично.

Методы быстрого построения и оценки

Одним из простых подходов к ускорению построения и оценки модели является использование подмножества $k$ ближайшие узлы интерполяции для построения локальной модели каждый раз, когда мы оцениваем сплайн. В результате общее время, необходимое для построения и оценки модели при $M$ баллы меняются от $O(N^{3}+MN)$ к $O(k^{3}*M)$.Это может привести к улучшению тайминга, если $k$ намного меньше, чем $N$. Такой подход поддерживается некоторыми программными библиотеками, наиболее известной из которых является scipy.interpolate.RBFInterpolator .Основным недостатком является то, что он вносит небольшие разрывы в сплайн и требует настройки под конкретную задачу: правильный выбор числа соседей, $k$.В последнее время были разработаны методы преодоления вышеупомянутых трудностей без ущерба для основных преимуществ полигармонических сплайнов.

Во-первых, куча методов для быстрого $O(\log N)$ были предложены оценки:

• Битсон и др. представить метод интерполяции полигармонических сплайнов с $r^{2k-1}$ быть базовой функцией в одной точке в трех измерениях или меньше
• Черри и др. представить метод интерполяции полигармонических сплайнов с $r^{2k}\log r$ как базисная функция в одной точке в 4 измерениях или меньше

Во-вторых, Браун и др. предложили ускоренное построение модели путем применения итеративного решателя к линейной системе, предварительно подготовленной ACBF. Такой подход сокращает время работы с $O(N^{3})$ к $O(N^{2})$,и далее к $O(N\log N)$ в сочетании с методами ускоренной оценки.

Описанные выше подходы часто используются коммерческими библиотеками анализа геопространственных данных и некоторыми реализациями с открытым исходным кодом (например, ALGLIB ).Иногда методы декомпозиции области используются для улучшения асимптотического поведения.снижение требований к памяти с $O(N^{2})$ к $O(N)$, что делает полигармонические сплайны подходящими для наборов данных с более чем 1 000 000 точек.

Причина названия «полигармонический».

Полигармоническое уравнение – это уравнение в частных производных вида $\Delta ^{m}f=0$ для любого натурального числа $m$, где $\Delta$ — оператор Лапласа . Например, бигармоническое уравнение имеет вид $\Delta ^{2}f=0$ и тригармоническое уравнение $\Delta ^{3}f=0$. Все полигармонические радиальные базисные функции являются решениями полигармонического уравнения (или, точнее, модифицированного полигармонического уравнения с дельта-функцией Дирака в правой части вместо 0). Например, радиальная базисная функция тонкой пластины является решением модифицированного двумерного бигармонического уравнения. Применение 2D-оператора Лапласа ( $\Delta =\partial _{xx}+\partial _{yy}$) к радиальной базисной функции тонкой пластины ${\textstyle f_{\text{tp}}(x,y)=(x^{2}+y^{2})\log {\sqrt {x^{2}+y^{2}}}}$ либо вручную, либо с использованием системы компьютерной алгебры, показывает, что $\Delta f_{\text{tp}}=4+4\log r$. Применяя оператор Лапласа к $\Delta f_{\text{tp}}$ (Это $\Delta ^{2}f_{\text{tp}}$) дает 0. Но 0 — это не совсем правильно. Чтобы увидеть это, замените $r^{2}=x^{2}+y^{2}$ с $\rho ^{2}=x^{2}+y^{2}+h^{2}$ (где $h$ — некоторое небольшое число, стремящееся к 0). Оператор Лапласа применяется к $4\log \rho$ урожайность $\Delta ^{2}f_{\text{tp}}=8h^{2}/\rho ^{4}$. Для $(x,y)=(0,0),$ правая часть этого уравнения стремится к бесконечности, поскольку $h$ приближается к 0. Для любого другого $(x,y)$, правая часть приближается к 0, поскольку $h$ приближается к 0. Это указывает на то, что правая часть представляет собой дельта-функцию Дирака. Система компьютерной алгебры покажет, что

$\lim _{h\to 0}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }8h^{2}/(x^{2}+y^{2}+h^{2})^{2}\,dx\,dy=8\pi .$

Таким образом, радиальная базисная функция тонкой пластины является решением уравнения $\Delta ^{2}f_{\text{tp}}=8\pi \delta (x,y)$.

Применение 3D лапласиана ( $\Delta =\partial _{xx}+\partial _{yy}+\partial _{zz}$) к бигармоническому RBF ${\textstyle f_{\text{bi}}(x,y,z)={\sqrt {x^{2}+y^{2}+z^{2}}}}$ урожайность $\Delta f_{\text{bi}}=2/r$ и применение 3D $\Delta ^{2}$ оператор тригармонического RBF $f_{\text{tri}}(x,y,z)=(x^{2}+y^{2}+z^{2})^{3/2}$ урожайность $\Delta ^{2}f_{\text{tri}}=24/r$. Сдача в аренду $\rho ^{2}=x^{2}+y^{2}+z^{2}+h^{2}$ и вычисления $\Delta (1/\rho )=-3h^{2}/\rho ^{5}$ снова указывает на то, что правая часть УЧП для бигармонических и тригармонических RBF представляет собой дельта-функцию Дирака. С

$\lim _{h\to 0}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }-3h^{2}/(x^{2}+y^{2}+z^{2}+h^{2})^{5/2}\,dx\,dy\,dz=-4\pi ,$

точные УЧП, которым удовлетворяют бигармонические и тригармонические RBF, равны $\Delta ^{2}f_{\text{bi}}=-8\pi \delta (x,y,z)$ и $\Delta ^{3}f_{\text{tri}}=-96\pi \delta (x,y,z)$.

Полигармонические сглаживающие сплайны

Полигармонические сплайны минимизируют

$\sum _{i=1}^{N}(f(\mathbf {c} _{i})-f_{i})^{2}+\lambda \int _{{\mathcal {B}}\subset \mathbb {R} ^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x}$ ( 3 )

где ${\mathcal {B}}$ какая-то коробка внутри $\mathbb {R} ^{d}$ содержащая окрестность всех центров, $\lambda$ — некоторая положительная константа, и $\nabla ^{m}f$ это вектор всех $m$частные производные го порядка $f.$ Например, в 2D $\nabla ^{1}f=(f_{x}\ f_{y})$ и $\nabla ^{2}f=(f_{xx}\ f_{xy}\ f_{yx}\ f_{yy})$ и в 3D $\nabla ^{2}f=(f_{xx}\ f_{xy}\ f_{xz}\ f_{yx}\ f_{yy}\ f_{yz}\ f_{zx}\ f_{zy}\ f_{zz})$. В 2D $|\nabla ^{2}f|^{2}=f_{xx}^{2}+2f_{xy}^{2}+f_{yy}^{2},$ делая интеграл упрощенным функционалом энергии тонкой пластины .

Чтобы показать, что полигармонические сплайны минимизируют уравнение ( 3 ), подгоночный член необходимо преобразовать в интеграл, используя определение дельта-функции Дирака:

$\sum _{i=1}^{N}(f(\mathbf {c} _{i})-f_{i})^{2}=\int _{\mathcal {B}}\sum _{i=1}^{N}(f(\mathbf {x} )-f_{i})^{2}\delta (\mathbf {x} -\mathbf {c} _{i})\,d\mathbf {x} .$

Таким образом, уравнение ( 3 ) можно записать в виде функционала

$J[f]=\int _{\mathcal {B}}F(\mathbf {x} ,f,\partial ^{\alpha _{1}}f,\partial ^{\alpha _{2}}f,\ldots ,\partial ^{\alpha _{n}}f)\,d\mathbf {x} =\int _{\mathcal {B}}\left[\sum _{i=1}^{N}(f(\mathbf {x} )-f_{i})^{2}\delta (\mathbf {x} -\mathbf {c} _{i})+\lambda |\nabla ^{m}f|^{2}\right]\,d\mathbf {x} .$

где $\alpha _{i}$ - это мультииндекс , который пробегает все частные производные порядка $m$ для $\mathbb {R} ^{d}.$ Чтобы применить уравнение Эйлера – Лагранжа для одной функции многих переменных и производных более высокого порядка, величины

${\partial F \over \partial f}=2\sum _{i=1}^{N}(f(\mathbf {x} )-f_{i})\delta (\mathbf {x} -x_{i})$

и

$\sum _{i=1}^{n}\partial ^{\alpha _{i}}{\partial F \over \partial (\partial ^{\alpha _{i}}f)}=2\lambda \Delta ^{m}f$

необходимы. Подстановка этих величин в уравнение E−L показывает, что

$\left(\sum _{i=1}^{N}(f(\mathbf {x} )-f_{i})\delta (\mathbf {x} -\mathbf {c} _{i})\right)+(-1)^{m}\lambda \Delta ^{m}f=0.$ ( 4 )

Слабое решение $f(\mathbf {x} )$ из ( 4 ) удовлетворяет

$\int _{\mathcal {B}}\left(\sum _{i=1}^{N}(f(\mathbf {x} )-f_{i})\delta (\mathbf {x} -\mathbf {c} _{i})\right)g(\mathbf {x} )+(-1)^{m}\lambda (\Delta ^{m}f)g(\mathbf {x} )\,d\mathbf {x} =0$

( 5 )

для всех функций плавного тестирования $g$ которые исчезают за пределами ${\mathcal {B}}.$ Слабое решение уравнения ( 4 ) все равно минимизирует ( 3 ), избавляясь при этом от дельта-функции путем интегрирования.

Позволять $f$ быть полигармоническим сплайном, как определено уравнением ( 1 ). Следующие расчеты покажут, что $f$ удовлетворяет ( 5 ). Применение $\Delta ^{m}$ оператор к уравнению ( 1 ) дает

$\Delta ^{m}f=\sum _{i=1}^{M}w_{i}C_{m,d}\delta (\mathbf {x} -\mathbf {c} _{i})$

где $C_{2,2}=8\pi ,$ $C_{2,3}=-8\pi ,$ и $C_{3,3}=-96\pi .$ Итак ( 5 ) эквивалентно

${\begin{aligned}\int _{\mathcal {B}}&\sum _{i=1}^{N}\delta (\mathbf {x} -\mathbf {c} _{i})(f(\mathbf {x} )-f_{i}+(-1)^{m}\lambda C_{m,d}w_{i})g(\mathbf {x} )\,d\mathbf {x} \\&=\sum _{i=1}^{N}(f(\mathbf {c} _{i})-f_{i}+(-1)^{m}\lambda C_{m,d}w_{i})g(\mathbf {c} _{i})\\&=0.\end{aligned}}$ ( 6 )

Единственное возможное решение ( 6 ) для всех тестовых функций $g$ является

$f(\mathbf {c} _{j})-f_{j}+(-1)^{m}\lambda C_{m,d}w_{j}=0\quad {\text{for}}\ j=1,2,\ldots ,N$ ( 7 )

(что подразумевает интерполяцию, если $\lambda =0$). Объединив определение $f$ в уравнении ( 1 ) с уравнением ( 7 ) дают почти ту же линейную систему, что и уравнение ( 2 ), за исключением того, что матрица $A$ заменяется на $A+(-1)^{m}C_{m,d}\lambda I$ где $I$ это $N\times N$ идентификационная матрица. Например, для трехмерных тригармонических RBF: $A$ заменяется на $A+96\pi \lambda I.$

Пояснение дополнительных ограничений

В ( 2 ) нижняя половина системы уравнений ( $B^{\textrm {T}}\mathbf {w} =0$) дается без пояснений. Объяснение сначала требует вывода упрощенной формы ${\textstyle \int _{\mathcal {B}}|\nabla ^{m}f|^{2}\,d\mathbf {x} }$ когда ${\mathcal {B}}$ это все из $\mathbb {R} ^{d}.$

Во-первых, требуйте, чтобы ${\textstyle \sum _{i=1}^{N}w_{i}=0.}$ Это гарантирует, что все производные порядка $m$ и выше ${\textstyle f(\mathbf {x} )=\sum _{i=1}^{N}w_{i}\varphi (|\mathbf {x} -\mathbf {c} _{i}|)}$ исчезнуть в бесконечности. Например, пусть $m=3$ и $d=3$ и $\varphi$ быть тригармоническим RBF. Затем $\varphi _{zzy}=3y(x^{2}+y^{2})/(x^{2}+y^{2}+z^{2})^{3/2}$ (учитывая $\varphi$ как отображение из $\mathbb {R} ^{3}$ к $\mathbb {R}$). Для данного центра $\mathbf {P} =(P_{1},P_{2},P_{3}),$

$\varphi _{zzy}(\mathbf {x} -\mathbf {P} )={\frac {3(y-P_{2})((y-P_{2})^{2}+(x-P_{1})^{2})}{((x-P_{1})^{2}+(y-P_{2})^{2}+(z-P_{3})^{2})^{3/2}}}.$

На линии $\mathbf {x} =\mathbf {a} +t\mathbf {b}$ для произвольной точки $\mathbf {a}$ и единичный вектор $\mathbf {b} ,$

$\varphi _{zzy}(\mathbf {x} -\mathbf {P} )={\frac {3(a_{2}+b_{2}t-P_{2})((a_{2}+b_{2}t-P_{2})^{2}+(a_{1}+b_{1}t-P_{1})^{2})}{((a_{1}+b_{1}t-P_{1})^{2}+(a_{2}+b_{2}t-P_{2})^{2}+(a_{3}+b_{3}t-P_{3})^{2})^{3/2}}}.$

Разделив числитель и знаменатель этого результата на $t^{3}$ показывает, что ${\textstyle \lim _{t\to \infty }\varphi _{zyy}(\mathbf {x} -\mathbf {P} )=3b_{2}(b_{2}^{2}+b_{1}^{2})/(b_{1}^{2}+b_{2}^{2}+b_{3}^{2})^{3/2},}$ величина, не зависящая от центра $\mathbf {P} .$ Итак, в данной строке

$\lim _{t\to \infty }f_{zyy}(\mathbf {x} )=\lim _{t\to \infty }\sum _{i=1}^{N}w_{i}\varphi _{zyy}(\mathbf {x} -\mathbf {c} _{i})=\left(\sum _{i=1}^{N}w_{i}\right)3b_{2}(b_{2}^{2}+b_{1}^{2})/(b_{1}^{2}+b_{2}^{2}+b_{3}^{2})^{3/2}=0.$

Требовать этого недостаточно ${\textstyle \sum _{i=1}^{N}w_{i}=0,}$ потому что в дальнейшем необходимо $f_{\alpha }g_{\beta }$ исчезнуть в бесконечности, где $\alpha$ и $\beta$ являются мультииндексами такими, что $|\alpha |+|\beta |=2m-1.$ Для тригармоники $\varphi ,$ $w_{i}u_{j}\varphi _{\alpha }(\mathbf {x} -\mathbf {c} _{i})\varphi _{\beta }(\mathbf {x} -\mathbf {d} _{j})$ (где $u_{j}$ и $\mathbf {d} _{j}$ это веса и центры $g$) всегда является суммой полных полиномов 5-й степени от $x,$ $y,$ и $z$ разделить на квадратный корень из полинома 8-й степени. Рассмотрим поведение этих слагаемых на линии $\mathbf {x} =\mathbf {a} +t\mathbf {b}$ как $t$ приближается к бесконечности. Числитель представляет собой многочлен пятой степени. $t.$ Деление числителя и знаменателя на $t^{4}$ оставляет члены 4-й и 5-й степени в числителе и функцию $\mathbf {b}$ только в знаменателе. Член 5-й степени, разделенный на $t^{4}$ является произведением пяти $b$ координаты и $t.$ ${\textstyle \sum w=0}$ (и ${\textstyle \sum u=0}$) ограничение заставляет это исчезнуть повсюду на линии. Член 4-й степени, разделенный на $t^{4}$ либо является произведением четырех $b$ координаты и $a$ координата или произведение четырех $b$ координаты и один $c_{i}$ или $d_{j}$ координировать. ${\textstyle \sum w=0}$ ограничение приводит к тому, что члены первого типа исчезают повсюду на линии. Дополнительные ограничения ${\textstyle \sum _{i=1}^{N}w_{i}\mathbf {c} _{i}=0}$ приведет к исчезновению терминов второго типа.

Теперь определим внутренний продукт двух функций $f,g:\mathbb {R} ^{d}\to \mathbb {R}$ определяется как линейная комбинация полигармонических RBF $\varphi _{m,d}$ с ${\textstyle \sum w=0}$ и ${\textstyle \sum w\mathbf {c} =0}$ как

$\langle f,g\rangle =\int _{\mathbb {R} ^{d}}(\nabla ^{m}f)\cdot (\nabla ^{m}g)\,d\mathbf {x} .$

Интегрирование по частям показывает, что

$\langle f,g\rangle =(-1)^{m}\int _{\mathbb {R} ^{d}}f(\Delta ^{m}g)\,d\mathbf {x} .$

( 8 )

Например, пусть $m=2$ и $d=2.$ Затем

$\langle f,g\rangle =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }(f_{xx}g_{xx}+2f_{xy}g_{xy}+f_{yy}g_{yy})\,dx\,dy.$

( 9 )

Интегрирование первого члена этого по частям однажды дает

$\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{xx}g_{xx}\,dx\,dy=\int _{-\infty }^{\infty }f_{x}g_{xx}{\big |}_{-\infty }^{\infty }\,dy-\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{x}g_{xxx}\,dx\,dy=-\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{x}g_{xxx}\,dx\,dy$

с $f_{x}g_{xx}$ исчезает в бесконечности. Повторное интегрирование по частям дает ${\textstyle \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }fg_{xxxx}\,dx\,dy.}$

Таким образом, интегрирование по частям дважды для каждого члена ( 9 ) дает

$\langle f,g\rangle =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f(g_{xxxx}+2g_{xxyy}+g_{yyyy})\,dx\,dy=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f(\Delta ^{2}g)\,dx\,dy.$

С ${\textstyle (\Delta ^{m}f)(\mathbf {x} )=\sum _{i=1}^{N}w_{i}C_{m,d}\delta (\mathbf {x-\mathbf {c} _{i}} ),}$ ( 8 ) показывает, что

${\begin{aligned}\langle f,f\rangle &=(-1)^{m}\int _{\mathbb {R} ^{d}}f(\mathbf {x} )\sum _{i=1}^{N}w_{i}(-1)^{m}C_{m,d}\delta (\mathbf {x-\mathbf {c} _{i}} )\,d\mathbf {x} =(-1)^{m}C_{m,d}\sum _{i=1}^{N}w_{i}f(\mathbf {c} _{i})\\&=(-1)^{m}C_{m,d}\sum _{i=1}^{N}\sum _{j=1}^{N}w_{i}w_{j}\varphi (\mathbf {c} _{i}-\mathbf {c} _{j})=(-1)^{m}C_{m,d}\mathbf {w} ^{\textrm {T}}A\mathbf {w} .\end{aligned}}$

Итак, если ${\textstyle \sum w=0}$ и ${\textstyle \sum w\mathbf {c} =0,}$

$\int _{\mathbb {R} ^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =(-1)^{m}C_{m,d}\mathbf {w} ^{\textrm {T}}A\mathbf {w} .$

( 10 )

Теперь о происхождении ограничений $B^{\textrm {T}}\mathbf {w} =0$ можно объяснить. Здесь $B$ представляет собой обобщение $B$ определено выше, чтобы возможно включать мономы до степени $m-1.$ Другими словами, $$B={\begin{bmatrix}1&1&\dots &1\\\mathbf {c} _{1}&\mathbf {c} _{2}&\dots &\mathbf {c} _{N}\\\vdots &\vdots &\dots &\vdots \\\mathbf {c} _{1}^{m-1}&\mathbf {c} _{2}^{m-1}&\dots &\mathbf {c} _{N}^{m-1}\end{bmatrix}}^{\textrm {T}}$$где $\mathbf {c} _{i}^{j}$ вектор-столбец любой степени $j$ мономы координат $\mathbf {c} _{i}.$ Верхняя половина ( 2 ) эквивалентна $A\mathbf {w} +B\mathbf {v} -\mathbf {f} =0.$ Таким образом, чтобы получить сглаживающий сплайн, необходимо минимизировать скалярное поле $F:\mathbb {R} ^{N+d+1}\rightarrow \mathbb {R}$ определяется

$F(\mathbf {w} ,\mathbf {v} )=|A\mathbf {w} +B\mathbf {v} -\mathbf {f} |^{2}+\lambda C\mathbf {w} ^{\textrm {T}}A\mathbf {w} .$

Уравнения

${\frac {\partial F}{\partial w_{i}}}=2A_{i*}(A\mathbf {w} +B\mathbf {v} -\mathbf {f} )+2\lambda CA_{i*}\mathbf {w} =0\quad {\textrm {for}}\ i=1,2,\ldots ,N$

и

${\frac {\partial F}{\partial v_{i}}}=2B_{i*}^{\textrm {T}}(A\mathbf {w} +B\mathbf {v} -\mathbf {f} )=0\quad {\textrm {for}}\ i=1,2,\ldots ,d+1$

(где $A_{i*}$ обозначает строку $i$ из $A$) эквивалентны двум системам линейных уравнений $A(A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\mathbf {w} )=0$ и $B^{\textrm {T}}(A\mathbf {w} +B\mathbf {v} -\mathbf {f} )=0.$ С $A$ обратима, первая система эквивалентна $A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\mathbf {w} =0.$ Таким образом, первая система подразумевает, что вторая система эквивалентна $B^{\textrm {T}}\mathbf {w} =0.$ Как и в предыдущем выводе коэффициента сглаживающего сплайна, верхняя половина ( 2 ) становится $(A+\lambda CI)\mathbf {w} +B\mathbf {v} =\mathbf {f} .$

Этот вывод системы уравнений полигармонического сглаживающего сплайна не предполагал ограничений, необходимых для того, чтобы гарантировать, что ${\textstyle \int _{{\mathcal {\mathbb {R} }}^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =Cw^{\textrm {T}}Aw.}$ Но ограничения, необходимые для того, чтобы гарантировать это, ${\textstyle \sum w=0}$ и ${\textstyle \sum w\mathbf {c} =0,}$ являются подмножеством $B^{\textrm {T}}w=0$ что справедливо для критической точки $w$ из $F.$ Так ${\textstyle \int _{{\mathcal {\mathbb {R} }}^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =Cw^{\textrm {T}}Aw}$ верно для $f$ формируется из решения системы уравнений полигармонического сглаживающего сплайна. Поскольку интеграл положителен для всех $w\neq 0,$ линейное преобразование, возникающее в результате ограничения области линейного преобразования $A$ к $w$ такой, что $B^{T}w=0$ должно быть положительно определенным. Это обстоятельство позволяет преобразовать систему уравнений полигармонического сглаживающего сплайна в симметричную положительно определенную систему уравнений, которая может быть решена в два раза быстрее с использованием разложения Холецкого.

Примеры

На следующем рисунке показана интерполяция по четырем точкам (отмеченным «кружками») с использованием различных типов полигармонических сплайнов. «Кривизна» интерполируемых кривых растет с увеличением порядка сплайна, и экстраполяция на левой границе ( x < 0) является разумной. На рисунке также указаны радиальные базисные функции φ = exp(− r ² ), что также дает хорошую интерполяцию. Наконец, на рисунке присутствует также неполигармонический сплайн phi = r ² чтобы продемонстрировать, что эта радиальная базисная функция не может проходить через заранее определенные точки (линейное уравнение не имеет решения и решается методом наименьших квадратов).

Интерполяция с различными полигармоническими сплайнами, которые должны проходить через 4 предопределенные точки, отмеченные кружком (интерполяция с phi = r ² бесполезно, так как система линейных уравнений интерполяционной задачи не имеет решения; решается методом наименьших квадратов, но тогда не проходит центры)

На следующем рисунке показана та же интерполяция, что и на первом рисунке, с той лишь разницей, что интерполируемые точки масштабируются в 100 раз (и случай phi = r ² больше не включается). Поскольку φ = (масштаб · r ) ^к = (масштаб ^к ) · р ^к , фактор (масштаб ^к ) можно извлечь из матрицы A системы линейных уравнений, и поэтому масштабирование не влияет на решение. Это отличается от логарифмической формы сплайна, хотя масштабирование не имеет большого влияния. Этот анализ отражен на рисунке, где интерполяция показывает не большие различия. Обратите внимание, что для других радиальных базисных функций, таких как φ = exp(− kr ² ) при k = 1 интерполяция больше не является разумной и необходимо будет адаптировать k .

Та же интерполяция, что и на первом рисунке, но точки, подлежащие интерполяции, масштабируются на 100.

На следующем рисунке показана та же интерполяция, что и на первом рисунке, с той лишь разницей, что полиномиальный член функции не учитывается (и случай phi = r ² больше не включается). Как видно из рисунка, экстраполяция при x для некоторых базисных функций < 0 уже не такая «естественная», как на первом рисунке. Это указывает на то, что полиномиальный член полезен в случае экстраполяции.

Та же интерполяция, что и на первом рисунке, но без полиномиального члена

См. также

• Обратное взвешивание расстояния
• Радиальная базисная функция
• Поверхность подразделения (новая альтернатива сплайновым поверхностям)
• Сплайн

Внешние ссылки

Компьютерный код

• Полигармонический сплайн . Интерактивный пример с исходным кодом Matlab/Octave.