Algoritmo de De Boor
# Algoritmo de De Boor

En el subcampo matemático del análisis numérico, el algoritmo de De Boor es un algoritmo de tiempo polinomial y numéricamente estable para evaluar curvas spline en forma B-spline. Es una generalización del algoritmo Casteljau para las curvas de Bézier. El algoritmo fue ideado por Carl R. De Boor. Se han creado variantes simplificadas y potencialmente más rápidas del algoritmo de De Boor, pero sufren una estabilidad comparativamente menor.

## Introducción

El algoritmo de De Boor es un esquema eficiente y numéricamente estable para evaluar una curva spline $\mathbf {S} (x)$ en posición $x$. La curva se construye a partir de una suma de funciones B-spline $B_{i,p}(x)$ multiplicada con valores vectoriales potencialmente constantes $\mathbf {c} _{i}$, llamados puntos de control.

    $\mathbf {S} (x)=\sum _{i}\mathbf {c} _{i}B_{i,p}(x).$

Las B-splines de orden $p+1$ son funciones polinómicas unitarias de grado $p$ definidas sobre una cuadrícula de nudos ${t_{0},\dots ,t_{i},\dots ,t_{m}}$ (se utilizan índices basados en cero en adelante). El algoritmo de De Boor utiliza operaciones $O(p^2)$ + O(p) para evaluar la curva de spline. Nota: El artículo principal sobre B-splines y las publicaciones clásicas utilizan una notación diferente: la B-spline es indexada como $B_{i,n}(x)$$n=p+1$.

## Soporte local

Las B-splines tienen soporte local, lo que significa que los polinomios son positivos solamente en un ámbito finito y cero en otros lugares. La fórmula de recursión Cox-De Boor muestra esto:

    $B_{i,0}(x):=\left\{{\begin{matrix}1&\mathrm {si} \quad t_{i}\leq x<t_{i+1}\\0&\mathrm {de\quad otra\quad forma} \end{matrix}}\right.$
    $B_{i,p}(x):={\frac {x-t_{i}}{t_{i+p}-t_{i}}}B_{i,p-1}(x)+{\frac {t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}}}B_{i+1,p-1}(x).$

Permitiendo que el índice $k$ defina el intervalo de nudos que contiene la posición, $x\in [t_{k},t_{k+1}]$. Podemos ver en la fórmula de recursión que solo B-splines con $i=k-p,\dots ,k$ no son ceros para este intervalo de nudos. Así, la suma se reduce a:

    $\mathbf {S} (x)=\sum _{i=k-p}^{k}\mathbf {c} _{i}B_{i,p}(x).$

De $i\geq 0$ se deduce que $k\geq p$. Del mismo modo, vemos en la recursividad que la ubicación del nudo más alto está en el índice $k+1+p$ . Esto significa que cualquier intervalo de nudos $[t_{k},t_{k+1})$ que se use realmente debe tener al menos $p$ nudos adicionales antes y después. En un programa de computadora, esto generalmente se logra repitiendo la primera y la última ubicación de nudo utilizada $p$ veces. Por ejemplo, para $p=3$ y ubicaciones de nudos reales $(0,1,2)$, uno podría rellenar el vector de nudos como $(0,0,0,0,1,2,2,2,2)$.

## El algoritmo

Con estas definiciones, ahora podemos describir el algoritmo de De Boor . El algoritmo no calcula las funciones de las B-splines $B_{i,p}(x)$ directamente. En cambio evalúa $\mathbf {S} (x)$ a través de una fórmula de recursión equivalente.

Deja a $\mathbf {d} _{i,r}$ ser un nuevo punto de control con $\mathbf {d} _{i,0}:=\mathbf {c} _{i}$ para $i=k-p,\dots ,k$. Para $r=1,\dots ,p$ la siguiente recursión es aplicada:

    $\mathbf {d} _{i,r}=(1-\alpha _{i,r})\mathbf {d} _{i-1,r-1}+\alpha _{i,r}\mathbf {d} _{i,r-1};\quad i=k-p+r,\dots ,k$
    $\alpha _{i,r}={\frac {x-t_{i}}{t_{i+1+p-r}-t_{i}}}.$

Una vez las iteraciones están completas, tenemos $\mathbf {S} (x)=\mathbf {d} _{k,p}$, esto significa que $\mathbf {d} _{k,p}$ es el resultado deseado.

El algoritmo De Boor es más eficaz que un cálculo explícito de B-splines $B_{i,p}(x)$ con la fórmula de recursión Cox-De Boor, porque no calcula términos que están garantizados para ser multiplicados por cero.

## Optimizaciones

El algoritmo anterior no está optimizado para la implementación en un ordenador. Requiere memoria para $(p+1)+p+\dots +1=(p+1)(p+2)/2$ puntos de control provisional $\mathbf {d} _{i,r}$. Cada punto de control provisional es escrito exactamente una vez y leídos dos veces. Al invertir la iteración sobre $i$ (contando hacia abajo, en vez de hacia arriba), podemos ejecutar el algoritmo con memoria para solo $p+1$ puntos de control provisionales, permitiendo a $\mathbf {d} _{i,r}$ reutilizar la memoria para $\mathbf {d} _{i,r-1}$. De modo parecido, hay solo un valor de $\alpha$ utilizado en cada paso, de manera que podemos reutilizar la memoria también.

Además, es más conveniente utilizar un índice basado en cero $j=0,\dots ,p$ para los puntos de control provisionales. La relación al índice anterior es $i=j+k-p$. Por ello obtenemos el algoritmo mejorado:

Sea $\mathbf {d} _{j}:=\mathbf {c} _{j+k-p}$ para $j=0,\dots ,p$. Iterar para $r=1,\dots ,p$:

    $\mathbf {d} _{j}:=(1-\alpha _{j})\mathbf {d} _{j-1}+\alpha _{j}\mathbf {d} _{j};\quad j=p,\dots ,r\quad$($j$ debe ser contado hacia abajo)

$\alpha _{j}:={\frac {x-t_{j+k-p}}{t_{j+1+k-r}-t_{j+k-p}}}.$

Después que las iteraciones se completan, el resultado es $\mathbf {S} (x)=\mathbf {d} _{p}$.

## Implementación de ejemplo

El código siguiente en el lenguaje de programación de Python es una implementación inocente ("naive") del algoritmo optimizado.

