Polyharmonisk spline

I tillämpad matematik används polyharmoniska splines för funktionsapproximation och datainterpolation . De är mycket användbara för att interpolera och anpassa spridda data i många dimensioner. Specialfall inkluderar tunna plattsplines och naturliga kubiska splines i en dimension.

Definition

En polyharmonisk spline är en linjär kombination av polyharmoniska radiella basfunktioner (RBF) betecknade med $\varphi$ plus en polynomterm:

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

 

 

 

 

()

var

Polyharmoniska basfunktioner

• $\mathbf {x} =[x_{1}\ x_{2}\ \cdots \ x_{d}]^{\textrm {T}}$ ( ${\textrm {T}}$ betecknar matristransponera, vilket betyder att $\mathbf {x}$ är en kolumnvektor) är en realvärderad vektor av $d$ oberoende variabler,
• $\mathbf {c} _{i}=[c_{i,1}\ c_{i,2}\ \cdots \ c_{i,d}]^{\textrm {T}}$ är $N$ vektorer av samma storlek som $\mathbf {x}$ (ofta kallade centra) som kurvan eller ytan måste interpolera,
• $\mathbf {w} =[w_{1}\ w_{2}\ \cdots \ w_{N}]^{\textrm {T}}$ är de $N$ vikterna för RBF:erna,
• $\mathbf {v} =[v_{1}\ v_{2}\ \cdots \ v_{d+1}]^{\textrm {T }}$ är vikterna $d+1$ för polynomet.

Polynomet med koefficienterna $\mathbf {v}$ förbättrar anpassningsnoggrannheten för polyharmoniska utjämningssplines och förbättrar även extrapoleringen bort från mitten $\mathbf {c} _{i}.$ Se figuren nedan för jämförelse av splines med polynomterm och utan polynomterm.

De polyharmoniska RBF:erna har formen:

${\begin{aligned}\varphi (r)&={\begin{cases}r^{k}&{\text{with }}k=1,3,5,\ldots ,\\r^ {k}\ln(r)&{\text{med }}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}}$

Andra värden för exponenten $k$ är inte användbara (som $\varphi (r)=r^{2}$ ), eftersom en lösning av interpolationsproblemet kanske inte existera. För att undvika problem vid $r=0$ (eftersom $\log(0)=-\infty$ ), kan de polyharmoniska RBF:erna med den naturliga logaritmen implementeras som:

$\varphi (r)={ \begin{cases}r^{k-1}\ln(r^{r})&{\text{för }}r<1,\quad {\text{(detta fungerar eftersom }}0^{0} {\text{ definieras)}}\\r^{k}\ln(r)&{\text{för }}r\geq 1.\end{cases}}$

eller, mer enkelt att lägga till en kontinuitetstillägg i $r=0$

$\varphi (r)={\begin{cases}0&{\text{för }}r<\epsilon ,\quad {\text{(för ett mycket litet värde av }}\epsilon {\text{, t.ex. om du använder flyttal med dubbel precision, }}\epsilon =10^{-200}{\text{)}}\\r^{k}\ln(r)&{\text{för }}r\ geq \epsilon .\end{cases}}$

Vikterna $w_{i}$ och $v_{j}$ bestäms så att funktionen interpolerar $N$ givna punkter $(\ mathbf {c} _{i},f_{i})$ (för $i=1,2,\ldots ,N$ ) och uppfyller $d+1$ ortogonalitetsförhållanden

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

Sammantaget är dessa begränsningar ekvivalenta med det symmetriska linjära ekvationssystemet

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

 

 

 

 

()

var

$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}}.$

För att detta ekvationssystem ska ha en unik lösning måste $B$ vara full rang. $B$ är full rang för mycket milda tillstånd på indata. Till exempel, i två dimensioner säkerställer tre centra som bildar en icke-degenererad triangel att $B$ är full rang, och i tre dimensioner säkerställer fyra centra som bildar en icke-degenererad tetraeder att B är full rang. Såsom förklaras senare, är den linjära transformationen som är ett resultat av begränsningen av domänen för den linjära transformationen $A$ till nollutrymmet för ${\textstyle B^{\textrm {T}}}$ positiv definitiv. Detta betyder att om ${\displaystyle B} är full rang$ , har ekvationssystemet ( 2 ) alltid en unik lösning och den kan lösas med en linjär lösare specialiserad för symmetriska matriser. De beräknade vikterna tillåter utvärdering av spline för alla $\mathbf {x} \in \mathbb {R} ^{d}$ med hjälp av ekvation ( 1 ). Många praktiska detaljer för att implementera och använda polyharmoniska splines förklaras i Fasshauer. I Iske behandlas polyharmoniska splines som specialfall av andra multiupplösningsmetoder i spridd datamodellering.

Anledning till namnet "polyharmonisk"

En polyharmonisk ekvation är en partiell differentialekvation av formen $\Delta ^{m}f=0$ för alla naturliga tal $m$ , där $\Delta$ är Laplace operatör . Till exempel är den biharmoniska ekvationen $\Delta ^{2}f=0$ och den triharmoniska ekvationen är $\Delta ^{3}f=0$ . Alla de polyharmoniska radiella basfunktionerna är lösningar av en polyharmonisk ekvation (eller mer exakt, en modifierad polyharmonisk ekvation med en Dirac deltafunktion på höger sida istället för 0). Till exempel är den tunna plattans radiella basfunktionen en lösning av den modifierade 2-dimensionella biharmoniska ekvationen. Användning av 2D Laplace-operatorn ( $\Delta =\partial _{xx}+\partial _{yy}$ ) på den tunna plattans radiella basfunktionen ${\textstyle f_{\text{tp}}(x,y)=(x^{2}+y^{2})\log {\sqrt {x^{2}+y^{2}}}}$ antingen för hand eller med ett datoralgebrasystem visar att $\Delta f_{\text{tp}}= 4+4\log r$ . Att tillämpa Laplace-operatorn på $\Delta f_{\text{tp}}$ (detta är ${\displaystyle \Delta ^{2}f_{\text{tp}}} )$ ger 0 Men 0 är inte riktigt korrekt. För att se detta, ersätt $r^{2}=x^{2}+y^{2}$ med $\ rho ^{2}=x^{2}+y^{2}+h^{2}$ (där $h$ är ett litet tal som tenderar mot 0). Laplace-operatorn som tillämpas på $4\log \rho$ ger $\Delta ^{2}f_{\text{tp}}=8h^ {2}/\rho ^{4}$ . För $(x,y)=(0,0),$ närmar sig den högra sidan av denna ekvation oändligheten när $h$ närmar sig 0. För alla andra $(x,y)$ , den högra sidan närmar sig 0 när $h$ närmar sig 0. Detta indikerar att den högra sidan är en Dirac deltafunktion. Ett datoralgebrasystem kommer att visa det

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

Så den tunna plattans radiella basfunktionen är en lösning av ekvationen $\Delta ^{2}f_{\text{tp}}=8\pi \delta ( x,y)$ .

Tillämpning av 3D Laplacian ( ${\displaystyle \Delta =\partial _{xx}+\partial _{yy}+\partial _{zz}})$ på den biharmoniska RBF ${\textstil f_{\text{bi}}(x,y,z)={\sqrt {x^{2}+y^{ 2}+z^{2}}}}$ ger $\Delta f_{\text{bi}}=2/r$ och tillämpar 3D $\Delta ^{ 2}$ operator till den triharmoniska RBF $f_{\text{tri}}(x,y,z)= (x^{2}+y^{2}+z^{2})^{3/2}$ ger $\Delta ^{2}f_{\text{tri }}=24/r$ . Låter $\rho ^{2}=x^{2}+y^{2}+z^{2}+h^{2}$ och beräkning av $\Delta (1/\rho )=-3h^{2}/\rho ^{5}$ visar återigen att den högra sidan av PDE för de biharmoniska och triharmoniska RBF:erna är Dirac deltafunktioner. Eftersom

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

de exakta PDE:erna som uppfylls av de biharmoniska och triharmoniska RBF:erna är $\Delta ^{2}f_{\text{bi}}=-8\pi \ delta (x,y,z)$ och $\Delta ^{3}f_{\text{tri}}=-96\pi \delta (x,y,z)$ .

Polyharmoniska utjämnande splines

Polyharmoniska splines minimerar

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

 

 

 

 

()

där ${\mathcal {B}}$ är någon ruta i $\mathbb {R} ^{d}$ som innehåller ett område av alla centra, $\lambda$ är något positivt konstant, och $\nabla ^{m}f$ är vektorn för alla ${\displaystyle m}:$ te ordningens partiella derivator av $f.$ Till exempel, i 2D $\nabla ^{1}f=(f_{x}\ f_{y})$ och $\nabla ^{2}f=(f_{xx}\ f_{xy}\ f_{yx}\ f_{yy})$ och i 3D $\nabla ^{2}f=(f_{xx}\ f_ {xy}\ f_{xz}\ f_{yx}\ f_{yy}\ f_{yz}\ f_{zx}\ f_{zy}\ f_{zz})$ . I 2D $|\nabla ^{2}f|^{2}=f_{xx}^{2}+2f_{xy}^{ 2}+f_{yy}^{2},$ gör integralen till den förenklade tunnplåtsenergin funktionell .

För att visa att polyharmoniska splines minimerar ekvationen ( 3 ), måste den passande termen omvandlas till en integral med hjälp av definitionen av Dirac delta-funktionen:

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

Så ekvation ( 3 ) kan skrivas som den funktionella

$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} .$

där $\alpha _{i}$ är ett multiindex som sträcker sig över alla partiella derivator av ordningen $m$ för $\mathbb {R} ^{d}.$ För att tillämpa Euler–Lagrange-ekvationen för en enda funktion av flera variabler och högre ordningens derivator,

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

och

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

är behövda. Att infoga dessa storheter i E−L-ekvationen visar det

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

 

 

 

 

()

En svag lösning $f(\mathbf {x} )$ av ( 4 ) uppfyller

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

 

 

 

 

()

för alla smidiga testfunktioner $g$ som försvinner utanför ${\mathcal {B}}.$ En svag lösning av ekvation ( 4 ) kommer fortfarande att minimera ( 3 ) samtidigt som man blir av med deltafunktionen genom integration.

Låt $f$ vara en polyharmonisk spline som definieras av ekvation ( 1 ). Följande beräkningar kommer att visa att $f$ uppfyller ( 5 ). Användning av $\Delta ^{m}$ på ekvation ( 1 ) ger

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

där $C_{2,2}=8\pi ,$$C_{2,3}=-8\pi ,$ och $C_{3,3}=-96\pi .$ Så ( 5 ) motsvarar

${\begin{aligned}\int _{\mathcal {B}}&\ summa _{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} \\&=\summa _{i=1}^{N}(f (\mathbf {c} _{i})-f_{i}+(-1)^{m}\lambda C_{m,d}w_{i})g(\mathbf {c} _{i}) \\&=0.\end{aligned}}$

 

 

 

 

()

Den enda möjliga lösningen på ( 6 ) för alla testfunktioner $g$ är

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

 

 

 

 

()

(vilket innebär interpolation om $\lambda =0$ ). Att kombinera definitionen av $f$ i ekvation ( 1 ) med ekvation ( 7 ) resulterar i nästan samma linjära system som ekvation ( 2 ) förutom att matrisen $A$ ersätts med $A+(-1)^{m}C_{m,d}\lambda I$ där $I$ är $N\ gånger N$ identitetsmatris. Till exempel, för 3D-triharmoniska RBFs, ersätts ${\displaystyle A} med$$A+96\pi \lambda I.$

Förklaring av ytterligare begränsningar

I ( 2 ) ges den nedre halvan av ekvationssystemet ( ${\displaystyle B^{\textrm {T}}\mathbf {w} =0} ) utan förklaring.$ Förklaringen kräver först att man härleder en förenklad form av ${\textstyle \int _{\mathcal {B}}|\nabla ^{m}f|^{2}\,d\mathbf {x} } när$ B $\displaystyle {\mathcal {B}} }$ är hela $\mathbb {R} ^{d}.$

Kräv först att ${\textstyle \sum _{i=1}^{N}w_{i}=0.} Detta$ säkerställer att alla derivator av ordningen $m$ och högre av ${\textstyle f(\mathbf {x} )=\summa _{i=1}^{N}w_{ i}\varphi (|\mathbf {x} -\mathbf {c} _{i}|)}$ försvinner i det oändliga. Låt till exempel $m=3$ och $d=3$ och $\varphi$ vara den triharmoniska RBF. Då $\varphi _{zzy}=3y(x^{2}+y^{ 2})/(x^{2}+y^{2}+z^{2})^{3/2}$ (med tanke på $\varphi$ som en mappning från $\mathbb {R} ^{3}$ till $\mathbb {R}$ ). För ett givet centrum $\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}}}.$

På en linje $\mathbf {x} =\mathbf {a} +t\mathbf {b}$ för godtycklig punkt $\mathbf {a}$ och enhetsvektor ${ \displaystyle \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}}}.$

Att dividera både täljare och nämnare av detta med $t^{3}$ visar att ${\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},}$ en kvantitet oberoende av mitten $\mathbf {P} .$ Så på den givna raden,

$\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.$

Det räcker inte att kräva att ${\textstyle \sum _{i=1}^{N}w_{i}=0,}$ eftersom det i det följande är nödvändigt för $f_{\alpha }g_{\beta }$ försvinner i oändligheten, där $\alpha$ och $\beta$ är multiindex som $|\alpha |+|\beta |=2m-1.$ För triharmonisk $\varphi ,$$w_{i}u_{j}\varphi _{\alpha }(\mathbf {x} -\mathbf {c} _{i})\varphi _{\beta }(\mathbf {x} -\mathbf {d} _{j})$ (där $u_{j}$ och $\mathbf {d} _{j}$ är vikterna och mitten av $g$ ) är alltid summan av totala grad 5 polynom i $x,$$y,$ och $z$ dividerat med kvadratroten av ett totalt grad 8 polynom . Betrakta beteendet hos dessa termer på linjen $\mathbf {x} =\mathbf {a} +t\mathbf {b}$ när $t$ närmar sig oändligheten. Täljaren är ett polynom på grad 5 i $t.$ Att dividera täljare och nämnare med $t^{4}$ lämnar grad 4 och 5 termer i täljaren och en funktion av $\mathbf {b}$ endast i nämnaren . En grad 5 term dividerad med $t^{4}$ är en produkt av fem $b$ -koordinater och $t.$ Restriktionen ${\textstyle \sum w=0}$ (och ${\textstyle \sum u=0}$ ) gör att detta försvinner överallt på linjen. En grad 4 term dividerad med $t^{4}$ är antingen en produkt av fyra $b$ -koordinater och en $a$ -koordinat eller en produkt av fyra $b$ koordinater och en enda $c_{i}$ eller $d_{j}$ koordinat. ∑ ${\textstyle \sum w=0} gör att den första typen av$ försvinner överallt på linjen. De ytterligare begränsningarna ${\textstyle \sum _{i=1}^{N}w_{i}\mathbf {c} _{i}=0}$ kommer att göra den andra typen av termen försvinna.

Definiera nu den inre produkten av två funktioner $f,g:\mathbb {R} ^{d}\to \mathbb {R}$ definierade som en linjär kombination av polyharmoniska RBFs $\varphi _{m,d}$ med ${\textstyle \sum w=0}$ och ${\textstyle \sum w\mathbf {c} =0}$ som

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

Integrering av delar visar det

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

 

 

 

 

()

Låt till exempel $m=2$ och $d=2.$ Sedan

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

 

 

 

 

()

Att integrera den första termen av detta med delar en gång ger

${$

eftersom $f_{x}g_{xx}$ försvinner i det oändliga. Integrering med delar igen resulterar i ${\textstyle \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }fg_{xxxx}\,dx\,dy.}$

Så att integrera med delar två gånger för varje term av ( 9 ) ger

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

Eftersom ${\textstyle (\Delta ^{m}f)(\mathbf {x} )=\ summa _{i=1}^{N}w_{i}C_{m,d}\delta (\mathbf {x-\mathbf {c} _{i}} ),} ( 8$ ) visar att

${\begin{aligned}\langle f,f\rangle &=(-1)^{m}\int _{\mathbb {R} ^{d}}f(\mathbf {x} )\summa _{i=1}^{N}w_{i}(-1)^{m}C_{m,d}\delta (\mathbf {x-\mathbf {c} _{i}} )\,d \mathbf {x} =(-1)^{m}C_{m,d}\summa _{i=1}^{N}w_{i}f(\mathbf {c} _{i})\\ &=(-1)^{m}C_{m,d}\summa _{i=1}^{N}\summa _{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}}$

Så om ${\textstyle \sum w=0}$ och ${\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} .$

 

 

 

 

()

Nu kan ursprunget för begränsningarna $B^{\textrm {T}}\mathbf {w} =0$ förklaras. Här $B$ en generalisering av $B$ som definierats ovan för att eventuellt inkludera monomer upp till grad $m-1.$ Med andra ord,

$$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}}$$

där $\mathbf {c} _{i}^{j}$ är en kolumnvektor av alla grad $j$ monomer av koordinaterna för $\mathbf {c} _{i}.$ Den övre halvan av ( 2 ) motsvarar $A\mathbf {w} +B\mathbf {v} -\mathbf {f} =0.$ Så för att få en utjämningsspline bör man minimera skalärfältet $F:\mathbb {R} ^{N+d+1 }\rightarrow \mathbb {R}$ definierad av

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

Ekvationerna

${\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 {för}}\ i=1,2,\ldots ,N$

och

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

(där $A_{i*}$ anger rad $i$ av $A$ ) är ekvivalenta med de två linjära ekvationssystemen $A(A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\mathbf {w} )=0$ och ${\displaystyle B^{\textrm {T}}(A\mathbf {w} +B\mathbf {v} -\mathbf {f} )=0.} Eftersom A {\$ displaystyle $}$ är inverterbart är det första systemet ekvivalent med $A\mathbf {w} +B\mathbf {v} -\mathbf {f} +\lambda C\ mathbf {w} =0.$ Så det första systemet innebär att det andra systemet är ekvivalent med $B^{\textrm {T}}\mathbf {w} =0.$ Precis som i föregående utjämningsspline-koefficientderivat, blir den övre halvan av ( 2 ) $(A+\lambda CI)\mathbf {w} +B\mathbf {v} =\mathbf {f} .$

Denna härledning av det polyharmoniska utjämningsspline-ekvationssystemet antog inte de begränsningar som var nödvändiga för att garantera att ${\textstyle \int _{{\mathcal {\mathbb {R} }}^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =Cw^{\textrm {T}}Aw.}$ Men de begränsningar som krävs för att garantera detta, ${\textstyle \sum w=0}$ och ${\textstyle \sum w\mathbf {c} =0,}$ är en delmängd av $B^{\textrm {T}}w=0$ vilket är sant för den kritiska punkten $w$ i $F.$ Så ${\textstyle \int _{{\mathcal {\mathbb {R} }}^{d}}|\nabla ^{m}f|^{2}\,d\mathbf {x} =Cw^{\textrm {T}}Aw}$ är sant för $f$ som bildas från lösningen av det polyharmoniska utjämningsspline-ekvationssystemet. Eftersom integralen är positiv för alla $w\neq 0,$ den linjära transformationen som härrör från begränsningen av domänen för linjär transformation $A$ till $w$ så att $B^{T}w=0$ måste vara positivt definitivt. Detta faktum gör det möjligt att transformera det polyharmoniska utjämnande spline-ekvationssystemet till ett symmetriskt positivt definitivt system av ekvationer som kan lösas dubbelt så snabbt med Cholesky-nedbrytningen.

Exempel

Nästa figur visar interpolationen genom fyra punkter (markerade med "cirklar") med olika typer av polyharmoniska splines. "Krökningen" för de interpolerade kurvorna växer med splineordningen och extrapoleringen vid den vänstra gränsen ( x < 0) är rimlig. Figuren inkluderar även de radiella basfunktionerna φ = exp(− r ² ) vilket ger en bra interpolation också. Slutligen inkluderar figuren även den icke-polyharmoniska spline phi = r ² för att visa att denna radiella basfunktion inte kan passera genom de fördefinierade punkterna (den linjära ekvationen har ingen lösning och löses i minsta kvadraters mening).

Interpolation med olika polyharmoniska splines som ska passera de 4 fördefinierade punkterna markerade av en cirkel (interpolationen med phi = r ² är inte användbar, eftersom det linjära ekvationssystemet för interpolationsproblemet inte har någon lösning; det är löst i minsta kvadraters mening, men passerar sedan inte centren)

Nästa figur visar samma interpolation som i den första figuren, med det enda undantaget att punkterna som ska interpoleras skalas med en faktor 100 (och fallet phi = r ² ingår inte längre). Eftersom φ = (skala· r ) ^k = (skala ^k )· r ^k , kan faktorn (skala ^k ) extraheras från matrisen A i det linjära ekvationssystemet och därför påverkas inte lösningen av skalningen. Detta är annorlunda för den logaritmiska formen av spline, även om skalningen inte har mycket inflytande. Denna analys återspeglas i figuren, där interpolationen inte visar stora skillnader. Observera att för andra radiella basfunktioner, såsom φ = exp(− kr ² ) med k = 1, är interpolationen inte längre rimlig och det skulle vara nödvändigt att anpassa k .

Samma interpolation som i den första figuren, men punkterna som ska interpoleras skalas med 100

Nästa figur visar samma interpolation som i den första figuren, med det enda undantaget att polynomtermen för funktionen inte beaktas (och fallet phi = r ² ingår inte längre). Som framgår av figuren är extrapoleringen för x < 0 inte längre lika "naturlig" som i den första figuren för vissa av basfunktionerna. Detta indikerar att polynomtermen är användbar om extrapolering sker.

Samma interpolation som i den första figuren, men utan polynomtermen

Diskussion

Den största fördelen med polyharmonisk spline-interpolation är att vanligtvis mycket bra interpolationsresultat erhålls för spridda data utan att utföra någon "avstämning", så automatisk interpolering är möjlig. Detta är inte fallet för andra radiella basfunktioner. Till exempel måste den Gaussiska funktionen $e^{-k\cdot r^{2}}$ ställas in, så att $k$ väljs enligt det underliggande rutnätet för de oberoende variablerna. Om detta rutnät är olikformigt är ett korrekt val av $k$ för att uppnå ett bra interpolationsresultat svårt eller omöjligt.

De största nackdelarna är:

• För att bestämma vikterna måste ett tätt linjärt ekvationssystem lösas. Att lösa ett tätt linjärt system blir opraktiskt om dimensionen $N$ är stor, eftersom minnet som krävs är $O(N^{2})$ och antalet operationer som krävs är $O(N^{3}).$
• Att utvärdera den beräknade polyharmoniska splinefunktionen vid $M$ datapunkter kräver $O(MN)$ operationer. I många applikationer (bildbehandling är ett exempel) $M$ mycket större än $N,$ och om båda siffrorna är stora är detta inte praktiskt.

Nyligen har metoder utvecklats för att övervinna ovannämnda svårigheter. Till exempel Beatson et al. presentera en metod för att interpolera polyharmoniska splines vid en punkt i 3 dimensioner i $O(\log N)$ operationer istället för $) {\displaystyle O(N)}$ operationer

Se även

• Omvänd avståndsviktning
• Radiell basfunktion
• Indelningsyta (framväxande alternativ till splinebaserade ytor)
• Spline

externa länkar

Datorkod

• Polyharmonic Spline , Ett interaktivt exempel med Matlab/Octave-källkod