Spline poliarmónico

Método de aproximación de funciones e interpolación de datos

En matemáticas aplicadas , los splines poliarmónicos se utilizan para la aproximación de funciones y la interpolación de datos . Son muy útiles para interpolar y ajustar datos dispersos en muchas dimensiones. Los casos especiales incluyen splines de placa delgada ^{[1] [2]} y splines cúbicos naturales en una dimensión. ^[3]

Definición

Un spline poliarmónico es una combinación lineal de funciones de base radiales (RBF) poliarmónicas denotadas por más un término polinomial: ${\displaystyle \varphi }$

dónde

Funciones base poliarmónicas

• ${\displaystyle \mathbf {x} =[x_{1}\ x_{2}\ \cdots \ x_{d}]^{\textrm {T}}}$( denota transposición de matriz, es decir , un vector columna) es un vector de valores reales de variables independientes, ${\displaystyle {\textrm {T}}}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle d}$
• ${\displaystyle \mathbf {c} _{i}=[c_{i,1}\ c_{i,2}\ \cdots \ c_{i,d}]^{\textrm {T}}}$son vectores del mismo tamaño que (a menudo llamados centros) que la curva o superficie debe interpolar, ${\displaystyle N}$ ${\displaystyle \mathbf {x} }$
• ${\displaystyle \mathbf {w} =[w_{1}\ w_{2}\ \cdots \ w_{N}]^{\textrm {T}}}$son los pesos de los RBF, ${\displaystyle N}$
• ${\displaystyle \mathbf {v} =[v_{1}\ v_{2}\ \cdots \ v_{d+1}]^{\textrm {T}}}$son los pesos del polinomio. ${\displaystyle d+1}$

El polinomio con los coeficientes mejora la precisión del ajuste para splines de suavizado poliarmónico y también mejora la extrapolación desde los centros. Consulte la figura a continuación para ver una comparación de splines con término polinomial y sin término polinomial. ${\displaystyle \mathbf {v} }$ ${\displaystyle \mathbf {c} _{i}.}$

Los RBF poliarmónicos tienen la forma:

${\displaystyle {\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}}}$

Otros valores del exponente no son útiles (como ), porque podría no existir una solución del problema de interpolación. Para evitar problemas en (ya que ), los RBF poliarmónicos con el logaritmo natural podrían implementarse como: ${\displaystyle k}$ ${\displaystyle \varphi (r)=r^{2}}$ ${\displaystyle r=0}$ ${\displaystyle \log(0)=-\infty }$

${\displaystyle \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}}}$

o, más simplemente, añadiendo una extensión de continuidad en ${\displaystyle r=0}$

${\displaystyle \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}}}$

Los pesos y se determinan de manera que la función interpole puntos dados (para ) y cumpla las condiciones de ortogonalidad. ${\displaystyle w_{i}}$ ${\displaystyle v_{j}}$ ${\displaystyle N}$ ${\displaystyle (\mathbf {c} _{i},f_{i})}$ ${\displaystyle i=1,2,\ldots ,N}$ ${\displaystyle d+1}$

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

En conjunto, estas restricciones son equivalentes al sistema de ecuaciones lineales simétricas.

dónde

${\displaystyle 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}}.}$

Para que este sistema de ecuaciones tenga una solución única, debe ser de rango completo. es de rango completo para condiciones muy suaves en los datos de entrada. Por ejemplo, en dos dimensiones, tres centros que forman un triángulo no degenerado aseguran que es de rango completo, y en tres dimensiones, cuatro centros que forman un tetraedro no degenerado aseguran que B es de rango completo. Como se explica más adelante, la transformación lineal resultante de la restricción del dominio de la transformación lineal al espacio nulo de es definida positiva. Esto significa que si es de rango completo, el sistema de ecuaciones ( 2 ) siempre tiene una solución única y se puede resolver utilizando un solucionador lineal especializado para matrices simétricas. Los pesos calculados permiten la evaluación del spline para cualquier ecuación que utilice ( 1 ). Muchos detalles prácticos de la implementación y el uso de splines poliarmónicos se explican en Fasshauer. ^[4] En Iske ^[5] los splines poliarmónicos se tratan como casos especiales de otros métodos de resolución múltiple en el modelado de datos dispersos. ${\displaystyle B}$ ${\displaystyle B}$ ${\displaystyle B}$ ${\displaystyle A}$ ${\textstyle B^{\textrm {T}}}$ ${\displaystyle B}$ ${\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}$

Discusión

La principal ventaja de la interpolación spline poliarmónica es que normalmente se obtienen muy buenos resultados de interpolación para datos dispersos sin realizar ningún "ajuste", por lo que la interpolación automática es factible. Este no es el caso de otras funciones de base radial. Por ejemplo, la función gaussiana necesita ajustarse, de modo que se seleccione de acuerdo con la cuadrícula subyacente de las variables independientes. Si esta cuadrícula no es uniforme, una selección adecuada para lograr un buen resultado de interpolación es difícil o imposible. ${\displaystyle e^{-k\cdot r^{2}}}$ ${\displaystyle k}$ ${\displaystyle k}$

Las principales desventajas son:

• Para determinar los pesos, se debe resolver un sistema lineal denso de ecuaciones. Resolver un sistema lineal denso se vuelve impráctico si la dimensión es grande, ya que la memoria requerida es y el número de operaciones requeridas es ${\displaystyle N}$ ${\displaystyle O(N^{2})}$ ${\displaystyle O(N^{3}).}$
• La evaluación de la función spline poliarmónica calculada en los puntos de datos requiere operaciones. En muchas aplicaciones (el procesamiento de imágenes es un ejemplo), es mucho mayor que y si ambos números son grandes, esto no es práctico. ${\displaystyle M}$ ${\displaystyle O(MN)}$ ${\displaystyle M}$ ${\displaystyle N,}$

Métodos rápidos de construcción y evaluación

Un enfoque sencillo para acelerar la construcción y evaluación de modelos es utilizar un subconjunto de nodos de interpolación más cercanos para construir un modelo local cada vez que evaluamos el spline. Como resultado, el tiempo total necesario para la construcción y evaluación de modelos en los puntos cambia de a . Esto puede producir mejores tiempos si es mucho menor que . Este enfoque es recomendado por algunas bibliotecas de software, la más notable es scipy.interpolate.RBFInterpolator. El principal inconveniente es que introduce pequeñas discontinuidades en el spline y requiere un ajuste específico del problema: una elección adecuada del recuento de vecinos, . Recientemente, se han desarrollado métodos para superar las dificultades mencionadas anteriormente sin sacrificar las principales ventajas de los splines poliarmónicos. ${\displaystyle k}$ ${\displaystyle M}$ ${\displaystyle O(N^{3}+MN)}$ ${\displaystyle O(k^{3}*M)}$ ${\displaystyle k}$ ${\displaystyle N}$ ${\displaystyle k}$

En primer lugar, se propusieron una serie de métodos para una evaluación rápida : ${\displaystyle O(\log N)}$

• Beatson et al. ^[6] presentan un método para interpolar splines poliarmónicos con una función base en un punto en 3 dimensiones o menos. ${\displaystyle r^{2k-1}}$
• Cherrie et al. ^[7] presentan un método para interpolar splines poliarmónicos con una función base en un punto en 4 dimensiones o menos. ${\displaystyle r^{2k}\log r}$

En segundo lugar, Brown et al. propusieron una construcción de modelo acelerada mediante la aplicación de un solucionador iterativo a un sistema lineal preacondicionado por ACBF. ^[8] Este enfoque reduce el tiempo de ejecución de a , y además de a cuando se combina con técnicas de evaluación acelerada. ${\displaystyle O(N^{3})}$ ${\displaystyle O(N^{2})}$ ${\displaystyle O(N\log N)}$

Los enfoques anteriores suelen emplearse en bibliotecas comerciales de análisis de datos geoespaciales y en algunas implementaciones de código abierto (por ejemplo, ALGLIB ). A veces, se utilizan métodos de descomposición de dominios para mejorar el comportamiento asintótico, lo que reduce los requisitos de memoria de a , lo que hace que los splines poliarmónicos sean adecuados para conjuntos de datos con más de 1.000.000 de puntos. ${\displaystyle O(N^{2})}$ ${\displaystyle O(N)}$

Razón del nombre "poliarmónico"

Una ecuación poliarmónica es una ecuación diferencial parcial de la forma para cualquier número natural , donde es el operador de Laplace . Por ejemplo, la ecuación biarmónica es y la ecuación triarmónica es . Todas las funciones de base radial poliarmónicas son soluciones de una ecuación poliarmónica (o más precisamente, una ecuación poliarmónica modificada con una función delta de Dirac en el lado derecho en lugar de 0). Por ejemplo, la función de base radial de placa delgada es una solución de la ecuación biarmónica bidimensional modificada. ^[9] La aplicación del operador de Laplace 2D ( ) a la función de base radial de placa delgada ya sea a mano o usando un sistema de álgebra computacional muestra que . La aplicación del operador de Laplace a (esto es ) produce 0. Pero 0 no es exactamente correcto. Para ver esto, reemplace con (donde es un número pequeño que tiende a 0). El operador de Laplace aplicado a produce . Para el lado derecho de esta ecuación tiende a infinito cuando se acerca a 0. Para cualquier otro , el lado derecho se acerca a 0 cuando se acerca a 0. Esto indica que el lado derecho es una función delta de Dirac. Un sistema de álgebra computacional mostrará que ${\displaystyle \Delta ^{m}f=0}$ ${\displaystyle m}$ ${\displaystyle \Delta }$ ${\displaystyle \Delta ^{2}f=0}$ ${\displaystyle \Delta ^{3}f=0}$ ${\displaystyle \Delta =\partial _{xx}+\partial _{yy}}$ ${\textstyle f_{\text{tp}}(x,y)=(x^{2}+y^{2})\log {\sqrt {x^{2}+y^{2}}}}$ ${\displaystyle \Delta f_{\text{tp}}=4+4\log r}$ ${\displaystyle \Delta f_{\text{tp}}}$ ${\displaystyle \Delta ^{2}f_{\text{tp}}}$ ${\displaystyle r^{2}=x^{2}+y^{2}}$ ${\displaystyle \rho ^{2}=x^{2}+y^{2}+h^{2}}$ ${\displaystyle h}$ ${\displaystyle 4\log \rho }$ ${\displaystyle \Delta ^{2}f_{\text{tp}}=8h^{2}/\rho ^{4}}$ ${\displaystyle (x,y)=(0,0),}$ ${\displaystyle h}$ ${\displaystyle (x,y)}$ ${\displaystyle h}$

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

Entonces, la función de base radial de la placa delgada es una solución de la ecuación . ${\displaystyle \Delta ^{2}f_{\text{tp}}=8\pi \delta (x,y)}$

La aplicación del laplaciano 3D ( ) al RBF biarmónico da como resultado y la aplicación del operador 3D al RBF triarmónico da como resultado . Si se deja y se calcula nuevamente, se indica que el lado derecho de las EDP para los RBF biarmónicos y triarmónicos son funciones delta de Dirac. Dado que ${\displaystyle \Delta =\partial _{xx}+\partial _{yy}+\partial _{zz}}$ ${\textstyle f_{\text{bi}}(x,y,z)={\sqrt {x^{2}+y^{2}+z^{2}}}}$ ${\displaystyle \Delta f_{\text{bi}}=2/r}$ ${\displaystyle \Delta ^{2}}$ ${\displaystyle f_{\text{tri}}(x,y,z)=(x^{2}+y^{2}+z^{2})^{3/2}}$ ${\displaystyle \Delta ^{2}f_{\text{tri}}=24/r}$ ${\displaystyle \rho ^{2}=x^{2}+y^{2}+z^{2}+h^{2}}$ ${\displaystyle \Delta (1/\rho )=-3h^{2}/\rho ^{5}}$

${\displaystyle \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 ,}$

Las ecuaciones diferenciales parciales exactas satisfechas por los RBF biarmónicos y triarmónicos son y . ${\displaystyle \Delta ^{2}f_{\text{bi}}=-8\pi \delta (x,y,z)}$ ${\displaystyle \Delta ^{3}f_{\text{tri}}=-96\pi \delta (x,y,z)}$

Splines de suavizado poliarmónico

Los splines poliarmónicos minimizan

donde es una caja en que contiene un entorno de todos los centros, es una constante positiva y es el vector de todas las derivadas parciales de orden th de Por ejemplo, en 2D y y en 3D . En 2D haciendo la integral la energía de placa delgada simplificada funcional . ${\displaystyle {\mathcal {B}}}$ ${\displaystyle \mathbb {R} ^{d}}$ ${\displaystyle \lambda }$ ${\displaystyle \nabla ^{m}f}$ ${\displaystyle m}$ ${\displaystyle f.}$ ${\displaystyle \nabla ^{1}f=(f_{x}\ f_{y})}$ ${\displaystyle \nabla ^{2}f=(f_{xx}\ f_{xy}\ f_{yx}\ f_{yy})}$ ${\displaystyle \nabla ^{2}f=(f_{xx}\ f_{xy}\ f_{xz}\ f_{yx}\ f_{yy}\ f_{yz}\ f_{zx}\ f_{zy}\ f_{zz})}$ ${\displaystyle |\nabla ^{2}f|^{2}=f_{xx}^{2}+2f_{xy}^{2}+f_{yy}^{2},}$

Para demostrar que los splines poliarmónicos minimizan la ecuación ( 3 ), el término de ajuste debe transformarse en una integral utilizando la definición de la función delta de Dirac:

${\displaystyle \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} .}$

Por lo tanto, la ecuación ( 3 ) puede escribirse como la función

${\displaystyle 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} .}$

donde es un multiíndice que abarca todas las derivadas parciales de orden para Para aplicar la ecuación de Euler-Lagrange para una sola función de múltiples variables y derivadas de orden superior, las cantidades ${\displaystyle \alpha _{i}}$ ${\displaystyle m}$ ${\displaystyle \mathbb {R} ^{d}.}$

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

y

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

son necesarios. Insertando estas cantidades en la ecuación E−L se muestra que

Una solución débil de ( 4 ) satisface ${\displaystyle f(\mathbf {x} )}$

para todas las funciones de prueba suaves que se desvanecen fuera de Una solución débil de la ecuación ( 4 ) aún minimizará ( 3 ) mientras se deshace de la función delta a través de la integración. ^[10] ${\displaystyle g}$ ${\displaystyle {\mathcal {B}}.}$

Sea un spline poliarmónico como se define en la ecuación ( 1 ). Los siguientes cálculos mostrarán que satisface ( 5 ). Al aplicar el operador a la ecuación ( 1 ), se obtiene ${\displaystyle f}$ ${\displaystyle f}$ ${\displaystyle \Delta ^{m}}$

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

donde y Entonces ( 5 ) es equivalente a ${\displaystyle C_{2,2}=8\pi ,}$ ${\displaystyle C_{2,3}=-8\pi ,}$ ${\displaystyle C_{3,3}=-96\pi .}$

La única solución posible para ( 6 ) para todas las funciones de prueba es ${\displaystyle g}$

(lo que implica interpolación si ). La combinación de la definición de en la ecuación ( 1 ) con la ecuación ( 7 ) da como resultado casi el mismo sistema lineal que la ecuación ( 2 ), excepto que la matriz se reemplaza por donde es la matriz identidad. Por ejemplo, para los RBF triarmónicos 3D, se reemplaza por ${\displaystyle \lambda =0}$ ${\displaystyle f}$ ${\displaystyle A}$ ${\displaystyle A+(-1)^{m}C_{m,d}\lambda I}$ ${\displaystyle I}$ ${\displaystyle N\times N}$ ${\displaystyle A}$ ${\displaystyle A+96\pi \lambda I.}$

Explicación de restricciones adicionales

En ( 2 ), la mitad inferior del sistema de ecuaciones ( ) se da sin explicación. La explicación requiere primero derivar una forma simplificada de cuando es todo ${\displaystyle B^{\textrm {T}}\mathbf {w} =0}$ ${\textstyle \int _{\mathcal {B}}|\nabla ^{m}f|^{2}\,d\mathbf {x} }$ ${\displaystyle {\mathcal {B}}}$ ${\displaystyle \mathbb {R} ^{d}.}$

Primero, se requiere que Esto asegura que todas las derivadas de orden y superior a se anulen en el infinito. Por ejemplo, sea y y el RBF triarmónico. Entonces (considerando como una aplicación de a ). Para un centro dado ${\textstyle \sum _{i=1}^{N}w_{i}=0.}$ ${\displaystyle m}$ ${\textstyle f(\mathbf {x} )=\sum _{i=1}^{N}w_{i}\varphi (|\mathbf {x} -\mathbf {c} _{i}|)}$ ${\displaystyle m=3}$ ${\displaystyle d=3}$ ${\displaystyle \varphi }$ ${\displaystyle \varphi _{zzy}=3y(x^{2}+y^{2})/(x^{2}+y^{2}+z^{2})^{3/2}}$ ${\displaystyle \varphi }$ ${\displaystyle \mathbb {R} ^{3}}$ ${\displaystyle \mathbb {R} }$ ${\displaystyle \mathbf {P} =(P_{1},P_{2},P_{3}),}$

${\displaystyle \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}}}.}$

En una línea para un punto arbitrario y un vector unitario ${\displaystyle \mathbf {x} =\mathbf {a} +t\mathbf {b} }$ ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {b} ,}$

${\displaystyle \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}}}.}$

Dividiendo tanto el numerador como el denominador de esto por se muestra que una cantidad independiente del centro . Entonces, en la línea dada, ${\displaystyle 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},}$ ${\displaystyle \mathbf {P} .}$

${\displaystyle \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.}$

No es suficiente exigir eso porque en lo que sigue es necesario que se desvanezca en el infinito, donde y son multiíndices tales que Para triarmónico (donde y son los pesos y centros de ) es siempre una suma de polinomios totales de grado 5 en y dividido por la raíz cuadrada de un polinomio total de grado 8. Considere el comportamiento de estos términos en la línea cuando se acerca al infinito. El numerador es un polinomio de grado 5 en Dividir numerador y denominador por deja los términos de grado 4 y 5 en el numerador y una función de solo en el denominador. Un término de grado 5 dividido por es un producto de cinco coordenadas y La restricción (y ) hace que esto se desvanezca en todas partes de la línea. Un término de grado 4 dividido por es un producto de cuatro coordenadas y una coordenada o un producto de cuatro coordenadas y una sola coordenada o . La restricción hace que el primer tipo de término se desvanezca en todas partes de la línea. Las restricciones adicionales harán que el segundo tipo de término se desvanezca. ${\textstyle \sum _{i=1}^{N}w_{i}=0,}$ ${\displaystyle f_{\alpha }g_{\beta }}$ ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle |\alpha |+|\beta |=2m-1.}$ ${\displaystyle \varphi ,}$ ${\displaystyle w_{i}u_{j}\varphi _{\alpha }(\mathbf {x} -\mathbf {c} _{i})\varphi _{\beta }(\mathbf {x} -\mathbf {d} _{j})}$ ${\displaystyle u_{j}}$ ${\displaystyle \mathbf {d} _{j}}$ ${\displaystyle g}$ ${\displaystyle x,}$ ${\displaystyle y,}$ ${\displaystyle z}$ ${\displaystyle \mathbf {x} =\mathbf {a} +t\mathbf {b} }$ ${\displaystyle t}$ ${\displaystyle t.}$ ${\displaystyle t^{4}}$ ${\displaystyle \mathbf {b} }$ ${\displaystyle t^{4}}$ ${\displaystyle b}$ ${\displaystyle t.}$ ${\textstyle \sum w=0}$ ${\textstyle \sum u=0}$ ${\displaystyle t^{4}}$ ${\displaystyle b}$ ${\displaystyle a}$ ${\displaystyle b}$ ${\displaystyle c_{i}}$ ${\displaystyle d_{j}}$ ${\textstyle \sum w=0}$ ${\textstyle \sum _{i=1}^{N}w_{i}\mathbf {c} _{i}=0}$

Ahora defina el producto interno de dos funciones definidas como una combinación lineal de RBF poliarmónicos con y como ${\displaystyle f,g:\mathbb {R} ^{d}\to \mathbb {R} }$ ${\displaystyle \varphi _{m,d}}$ ${\textstyle \sum w=0}$ ${\textstyle \sum w\mathbf {c} =0}$

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

La integración por partes muestra que

Por ejemplo, sea y Entonces ${\displaystyle m=2}$ ${\displaystyle d=2.}$

Integrando el primer término de esto por partes una vez obtenemos

${\displaystyle \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}$

ya que se anula en el infinito. Integrando por partes nuevamente resulta en ${\displaystyle f_{x}g_{xx}}$ ${\textstyle \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }fg_{xxxx}\,dx\,dy.}$

Así que, integrando por partes dos veces para cada término de ( 9 ), obtenemos

${\displaystyle \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.}$

Dado que ( 8 ) muestra que ${\textstyle (\Delta ^{m}f)(\mathbf {x} )=\sum _{i=1}^{N}w_{i}C_{m,d}\delta (\mathbf {x-\mathbf {c} _{i}} ),}$

${\displaystyle {\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}}}$

Así que si y ${\textstyle \sum w=0}$ ${\textstyle \sum w\mathbf {c} =0,}$

Ahora se puede explicar el origen de las restricciones . Aquí hay una generalización de lo definido anteriormente para incluir posiblemente monomios hasta grado En otras palabras, donde es un vector columna de todos los monomios de grado de las coordenadas de La mitad superior de ( 2 ) es equivalente a Entonces, para obtener una spline de suavizado, se debe minimizar el campo escalar definido por ${\displaystyle B^{\textrm {T}}\mathbf {w} =0}$ ${\displaystyle B}$ ${\displaystyle B}$ ${\displaystyle m-1.}$ $${\displaystyle 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}}}$$ ${\displaystyle \mathbf {c} _{i}^{j}}$ ${\displaystyle j}$ ${\displaystyle \mathbf {c} _{i}.}$ ${\displaystyle A\mathbf {w} +B\mathbf {v} -\mathbf {f} =0.}$ ${\displaystyle F:\mathbb {R} ^{N+d+1}\rightarrow \mathbb {R} }$

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

Las ecuaciones

${\displaystyle {\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}$

y

${\displaystyle {\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}$

(donde denota la fila de ) son equivalentes a los dos sistemas de ecuaciones lineales y Dado que es invertible, el primer sistema es equivalente a Por lo que el primer sistema implica que el segundo sistema es equivalente a Al igual que en la derivación del coeficiente de spline de suavizado anterior, la mitad superior de ( 2 ) se convierte en ${\displaystyle A_{i*}}$ ${\displaystyle i}$ ${\displaystyle A}$ ${\displaystyle A(A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\mathbf {w} )=0}$ ${\displaystyle B^{\textrm {T}}(A\mathbf {w} +B\mathbf {v} -\mathbf {f} )=0.}$ ${\displaystyle A}$ ${\displaystyle A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\mathbf {w} =0.}$ ${\displaystyle B^{\textrm {T}}\mathbf {w} =0.}$ ${\displaystyle (A+\lambda CI)\mathbf {w} +B\mathbf {v} =\mathbf {f} .}$

Esta derivación del sistema de ecuaciones de spline de suavizado poliarmónico no asumió las restricciones necesarias para garantizar que Pero las restricciones necesarias para garantizar esto, y son un subconjunto de lo cual es verdadero para el punto crítico de Por lo tanto es verdadero para el formado a partir de la solución del sistema de ecuaciones de spline de suavizado poliarmónico. Porque la integral es positiva para toda la transformación lineal resultante de la restricción del dominio de la transformación lineal a tal que debe ser definida positiva. Este hecho permite transformar el sistema de ecuaciones de spline de suavizado poliarmónico en un sistema de ecuaciones simétrico definido positivo que se puede resolver el doble de rápido utilizando la descomposición de Cholesky. ^[9] ${\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,}$ ${\displaystyle B^{\textrm {T}}w=0}$ ${\displaystyle w}$ ${\displaystyle F.}$ ${\textstyle \int _{{\mathcal {\mathbb {R} }}^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =Cw^{\textrm {T}}Aw}$ ${\displaystyle f}$ ${\displaystyle w\neq 0,}$ ${\displaystyle A}$ ${\displaystyle w}$ ${\displaystyle B^{T}w=0}$

Ejemplos

La siguiente figura muestra la interpolación a través de cuatro puntos (marcados por "círculos") utilizando diferentes tipos de splines poliarmónicos. La "curvatura" de las curvas interpoladas crece con el orden del spline y la extrapolación en el límite izquierdo ( x  < 0) es razonable. La figura también incluye las funciones de base radial φ = exp(− r ² ) que también dan una buena interpolación. Finalmente, la figura también incluye el spline no poliarmónico phi = r ² para demostrar que esta función de base radial no puede pasar por los puntos predefinidos (la ecuación lineal no tiene solución y se resuelve en un sentido de mínimos cuadrados).

Interpolación con diferentes splines poliarmónicos que deben pasar los 4 puntos predefinidos marcados por un círculo (la interpolación con phi = r ² no es útil, ya que el sistema de ecuaciones lineales del problema de interpolación no tiene solución; se resuelve en un sentido de mínimos cuadrados, pero luego no pasa por los centros)

La siguiente figura muestra la misma interpolación que en la primera figura, con la única excepción de que los puntos a interpolar se escalan por un factor de 100 (y el caso phi = r ² ya no está incluido). Dado que φ = (escala· r ) ^k = (escala ^k )· r ^k , el factor (escala ^k ) se puede extraer de la matriz A del sistema de ecuaciones lineales y, por lo tanto, la solución no se ve influenciada por la escala. Esto es diferente para la forma logarítmica del spline, aunque la escala no tiene mucha influencia. Este análisis se refleja en la figura, donde la interpolación no muestra muchas diferencias. Nótese que, para otras funciones de base radial, como φ = exp(− kr ² ) con k  = 1, la interpolación ya no es razonable y sería necesario adaptar  k .

La misma interpolación que en la primera figura, pero los puntos a interpolar están escalados por 100

La siguiente figura muestra la misma interpolación que en la primera figura, con la única excepción de que no se tiene en cuenta el término polinómico de la función (y ya no se incluye el caso phi = r ² ). Como se puede ver en la figura, la extrapolación para x  < 0 ya no es tan "natural" como en la primera figura para algunas de las funciones base. Esto indica que el término polinómico es útil si se produce la extrapolación.

La misma interpolación que en la primera figura, pero sin el término polinomial

Véase también

• Ponderación de distancia inversa
• Función de base radial
• Superficie de subdivisión (alternativa emergente a las superficies basadas en splines)
• Ranura

Referencias

1. RL Harder y RN Desmarais: Interpolación mediante splines de superficie. Journal of Aircraft, 1972, número 2, págs. 189-191
2. J. Duchon: Splines que minimizan seminormas invariantes a la rotación en espacios de Sobolev. Teoría constructiva de funciones de varias variables, W. Schempp y K. Zeller (eds), Springer, Berlín, págs. 85-100
3. Wendland, Holger (2005). Aproximación de datos dispersos . Cambridge University Press. pág. 9. ISBN 0521843359.
4. GF Fasshauer GF: Métodos de aproximación sin malla con MATLAB. World Scientific Publishing Company, 2007, ISPN-10: 9812706348
5. A. Iske: Métodos multirresolución en modelado de datos dispersos, Notas de clase en ciencia computacional e ingeniería, 2004, vol. 37, ISBN 3-540-20479-2 , Springer-Verlag, Heidelberg. 
6. RK Beatson, MJD Powell y AM Tan: Evaluación rápida de splines poliarmónicos en tres dimensiones. IMA Journal of Numerical Analysis, 2007, 27, págs. 427–450.
7. JB Cherrie; RK Beatson; DL Ragozin (2000), Evaluación rápida de funciones de base radial: métodos para splines poliarmónicos de cuatro dimensiones
8. Damian Brown; Leevan Ling; Edward Kansa; Jeremy Levesley (2000), Sobre métodos aproximados de preacondicionamiento cardinal para resolver ecuaciones en derivadas parciales con funciones de base radial
9. Powell, MJD (1993). "Algunos algoritmos para la interpolación de splines de placas delgadas en funciones de dos variables" (PDF) . Informe técnico del Departamento de Matemáticas Aplicadas y Física Teórica de la Universidad de Cambridge . Archivado desde el original (PDF) el 25 de enero de 2016. Consultado el 7 de enero de 2016 .
10. Evans, Lawrence (1998). Ecuaciones diferenciales parciales . Providence: American Mathematical Society. pp. 450−452. ISBN 0-8218-0772-2.

Enlaces externos

Código de computadora

• Spline poliarmónico, un ejemplo interactivo con código fuente de Matlab/Octave