Interpolazione a tratti, spline cubiche e minimi quadrati

Prerequisiti: Interpolazione polinomiale - forma di Lagrange, differenze divise ed erroreDati n+1 punti con ascisse distinte esiste un unico polinomio p_n di grado ≤ n che passa per tutti. Forma di Lagrange p = Σ y_i L_i(x) (L_i(x_j) = δ_ij); forma di Newton p = Σ f[x_0..x_k] Π_{j<k}(x − x_j) con le differenze divise (si aggiunge un nodo senza rifare tutto). Errore f(x) − p_n(x) = f^(n+1)(ξ)/(n+1)! · ω(x), ω(x) = Π(x − x_i): dipende da f e dalla scelta dei nodi. Nodi equispaziati con n grande danno il fenomeno di Runge; i nodi di Chebyshev minimizzano max|ω| = 2^(−n) su [−1,1].Interpolazione polinomiale - forma di Lagrange, differenze divise ed errore → (errore di interpolazione, Runge), Complemento ortogonale e proiezioni ortogonaliL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali → (minimi quadrati), Eliminazione di GaussCon tre operazioni elementari sulle righe (scambio, moltiplicazione per uno scalare non nullo, somma di un multiplo di un'altra riga) ogni matrice si riduce a scala senza cambiare il rango; serve a calcolare ranghi, risolvere sistemi, trovare relazioni di dipendenza e matrici che riducono a scala.Eliminazione di Gauss → (sistemi lineari).

L'interpolazione con un unico polinomio di grado alto ha due difetti: il fenomeno di Runge e il cattivo condizionamento. Con molti dati, o con dati sperimentali affetti da errore, due strategie migliori sono dividere in tratti (interpolazione a tratti, spline) oppure non passare per tutti i punti (approssimazione ai minimi quadrati).

Interpolazione lineare a tratti

Dati i nodi $x_0 < \dots < x_n$ e i valori $y_i$, in ogni intervallo $[x_i, x_{i+1}]$ si prende il segmento che unisce i punti: $$s(x) = y_i + \frac{y_{i+1} - y_i}{x_{i+1} - x_i}(x - x_i), \qquad x \in [x_i, x_{i+1}].$$ È continua ma ha "spigoli" nei nodi (derivata discontinua).

Teorema (errore). Se $f \in C^2$ e $h = \max_i (x_{i+1} - x_i)$, $$\max_{[a,b]}|f(x) - s(x)| \le \frac{h^2}{8}\,\max_{[a,b]}|f''|.$$

Dimostrazione. In ogni intervallo è l'interpolazione con $n=1$: $f - s = \frac{f''(\xi)}{2}(x - x_i)(x - x_{i+1})$ e $|(x-x_i)(x-x_{i+1})|$ ha massimo $\frac{h_i^2}{4}$ nel punto medio, da cui $\frac{h^2}{8}\max|f''|$ ∎.

L'errore tende a $0$ come $h^2$ senza condizioni sui nodi (a differenza del polinomio unico): dimezzando $h$ si divide l'errore per $4$.

Esempio. $f = e^x$ su $[0,1]$, $f'' \le e$. Con $h = 0{,}1$ ($10$ intervalli): maggiorante $\frac{0{,}01}{8}e = 3{,}40\cdot 10^{-3}$; errore effettivo (massimo calcolato su una griglia fine) $3{,}23\cdot 10^{-3}$. Per garantire $10^{-4}$: $h^2 \le \frac{8\cdot 10^{-4}}{e} = 2{,}94\cdot 10^{-4}$, $h \le 0{,}01716$, quindi $N \ge 58{,}3$, cioè $N = 59$ intervalli (errore effettivo $9{,}7\cdot 10^{-5}$ ✓).

(Il rapporto tra errori successivi per $h$ dimezzato è $\approx 4$: ordine $2$.)

Spline cubiche

Per ottenere una curva liscia senza usare un polinomio di grado alto si usano, in ogni intervallo, polinomi di grado $3$ raccordati in modo da avere derivata prima e seconda continue.

Definizione (spline cubica interpolante). È una funzione $S: [x_0, x_n] \to \mathbb{R}$ tale che

1. su ogni $[x_i, x_{i+1}]$, $S = S_i$ è un polinomio di grado $\le 3$;
2. $S(x_i) = y_i$ per $i = 0,\dots,n$;
3. $S \in C^2([x_0,x_n])$: $S_{i-1}'(x_i) = S_i'(x_i)$ e $S_{i-1}''(x_i) = S_i''(x_i)$ per $i = 1,\dots,n-1$.

Conteggio. Incognite: $4n$ coefficienti. Condizioni: $2n$ (ogni $S_i$ passa per i suoi due estremi) $+\ (n-1)$ (derivata prima) $+\ (n-1)$ (derivata seconda) $= 4n - 2$. Mancano due condizioni, che si assegnano agli estremi:

• spline naturale: $S''(x_0) = S''(x_n) = 0$ (la curva diventa un segmento fuori dall'intervallo, come un'asta elastica a cui non si applicano momenti agli estremi);
• spline vincolata (clamped): si assegnano $S'(x_0)$ e $S'(x_n)$;
• not-a-knot: $S'''$ continua nei nodi $x_1$ e $x_{n-1}$ (equivale a un unico polinomio nei primi e negli ultimi due intervalli; è il default di MATLAB spline).

Costruzione con i momenti. Si chiama $M_i = S''(x_i)$. Poiché $S_i''$ è lineare, vale $S_i''(x) = M_i\frac{x_{i+1} - x}{h_i} + M_{i+1}\frac{x - x_i}{h_i}$ con $h_i = x_{i+1} - x_i$. Integrando due volte e imponendo $S_i(x_i) = y_i$, $S_i(x_{i+1}) = y_{i+1}$: $$S_i(x) = \frac{M_i(x_{i+1} - x)^3 + M_{i+1}(x - x_i)^3}{6h_i} + \left(\frac{y_i}{h_i} - \frac{M_ih_i}{6}\right)(x_{i+1} - x) + \left(\frac{y_{i+1}}{h_i} - \frac{M_{i+1}h_i}{6}\right)(x - x_i).$$ La continuità di $S'$ nel nodo $x_i$ dà un'equazione per i momenti: $$\boxed{\frac{h_{i-1}}{6}M_{i-1} + \frac{h_{i-1} + h_i}{3}M_i + \frac{h_i}{6}M_{i+1} = \frac{y_{i+1} - y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}}}, \quad i = 1,\dots,n-1.$$ È un sistema tridiagonale a diagonale dominante (e quindi invertibile) nelle $n+1$ incognite $M_i$, con le due condizioni agli estremi ($M_0 = M_n = 0$ per la naturale). Si risolve in $O(n)$ operazioni. Per nodi equispaziati ($h_i = h$), moltiplicando l'equazione per $\frac{6}{h}$, diventa $$M_{i-1} + 4M_i + M_{i+1} = \frac{6}{h^2}\left(y_{i-1} - 2y_i + y_{i+1}\right).$$

Esempio svolto. Nodi $0,1,2,3$ ($h=1$), valori $y = 0,\ 1,\ 0,\ 1$; spline naturale: $M_0 = M_3 = 0$. Equazioni per $i=1,2$: $$4M_1 + M_2 = 6(0 - 2 + 0) = -12, \qquad M_1 + 4M_2 = 6(1 - 0 + 1) = 12.$$ Dalla prima $M_2 = -12 - 4M_1$; nella seconda $M_1 - 48 - 16M_1 = 12 \Rightarrow -15M_1 = 60 \Rightarrow M_1 = -4$, $M_2 = 4$. Gli intervalli:

• $S_0(x) = \frac{M_1x^3}{6} + \left(y_1 - \frac{M_1}6\right)x + y_0(1-x) = -\frac23x^3 + \frac53x$ (con $M_0=0$, $y_0 = 0$, $y_1 - M_1/6 = 1 + \frac23$);
• $S_1(x) = -\frac23(2-x)^3 + \frac23(x-1)^3 + \frac53(2-x) - \frac23(x-1)$;
• $S_2(x) = \frac23(3-x)^3 - \frac23(3-x) + (x-2)$.

Controlli: $S_0(1) = -\frac23 + \frac53 = 1$ ✓. Derivata in $x=1$: $S_0'(1) = -2 + \frac53 = -\frac13$ e $S_1'(1) = 2 - \frac53 - \frac23 = -\frac13$ ✓ (continua). Derivata seconda: $S_0''(x) = -4x$ vale $-4 = M_1$ in $x=1$; $S_1''(x) = -4(2-x) + 4(x-1)$ vale $-4$ in $x=1$ e $4 = M_2$ in $x=2$ ✓ (continua). Valori: $S(0{,}5) = -\frac23\cdot 0{,}125 + \frac53\cdot 0{,}5 = 0{,}75$; $S(1{,}5) = 0{,}5$ (centro di simmetria: $S(x) + S(3 - x) = 1$); $S(2{,}5) = 0{,}25$.

Il polinomio di grado $3$ sugli stessi dati è $\frac23x^3 - 3x^2 + \frac{10}3x$ e vale $1$ in $x = 0{,}5$ e $0$ in $x = 2{,}5$, mentre la spline dà $0{,}75$ e $0{,}25$: il polinomio oscilla di più.

Accuratezza. Per $f \in C^4$ e nodi con $h$ uniforme, la spline vincolata (o not-a-knot) soddisfa $$\max|f - S| \le \frac{5}{384}\,h^4\,\max|f^{(4)}|,$$ cioè errore $O(h^4)$ (dimezzando $h$ si divide per $16$). La spline naturale impone $S'' = 0$ agli estremi anche se $f''(x_0) \ne 0$: perde accuratezza vicino agli estremi (errore $O(h^2)$ lì), a meno che $f''$ sia nulla agli estremi.

Per $f(x) = \frac{1}{1+25x^2}$ su $[-1,1]$, errore massimo con $n$ intervalli equispaziati:

L'errore della spline diminuisce all'aumentare di $n$, quello del polinomio unico esplode (Runge). In MATLAB: spline(x,y,t), interp1(x,y,t,'linear'), interp1(x,y,t,'spline'), pchip (a tratti che non oscilla). In Python: np.interp (lineare), scipy.interpolate.CubicSpline.

Approssimazione ai minimi quadrati

Con dati sperimentali ($m$ misure affette da errore) non ha senso far passare un polinomio per tutti i punti: se $m$ è grande il polinomio oscilla seguendo il rumore. Si cerca un modello con pochi parametri (ad esempio una retta) che si avvicini ai dati nel senso dei minimi quadrati.

Problema ai minimi quadrati. Dati $m$ punti $(x_k, y_k)$ e $n+1 \le m$ funzioni di base $\varphi_0,\dots,\varphi_n$ (ad esempio $1, x, \dots, x^n$), trovare i coefficienti $c = (c_0,\dots,c_n)$ che minimizzano $$\Phi(c) = \sum_{k=1}^{m}\bigl(y_k - p(x_k)\bigr)^2 = \|y - Ac\|_2^2, \qquad p = \sum_j c_j\varphi_j,\ A_{kj} = \varphi_j(x_k).$$

Equazioni normali. $Ac$ descrive i vettori della forma "valori del modello nei punti", cioè il sottospazio generato dalle colonne di $A$. Il vettore di quel sottospazio più vicino a $y$ è la proiezione ortogonale di $y$ (Complemento ortogonale e proiezioni ortogonaliL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali →): il residuo $r = y - Ac$ è ortogonale a tutte le colonne di $A$, cioè $A^T(y - Ac) = 0$: $$\boxed{A^TA\,c = A^Ty.}$$ (Lo stesso risultato si ottiene ponendo a zero il gradiente di $\Phi$.) Se le colonne di $A$ sono linearmente indipendenti (per i polinomi di grado $n$: almeno $n+1$ ascisse distinte) la matrice $A^TA$ è simmetrica definita positiva e la soluzione è unica. Il minimo è $\|r\|^2 = \|y\|^2 - \|Ac\|^2$.

Retta di regressione $p(x) = a + bx$. Le equazioni normali si scrivono $$\begin{pmatrix} m & \sum x_k \\ \sum x_k & \sum x_k^2\end{pmatrix}\begin{pmatrix}a\\ b\end{pmatrix} = \begin{pmatrix}\sum y_k\\ \sum x_ky_k\end{pmatrix}, \qquad b = \frac{m\sum x_ky_k - \sum x_k\sum y_k}{m\sum x_k^2 - \left(\sum x_k\right)^2},\quad a = \frac{\sum y_k - b\sum x_k}{m}.$$

Esempio svolto. Dati: $(0;\,1{,}1)$, $(1;\,1{,}9)$, $(2;\,3{,}2)$, $(3;\,3{,}8)$, $(4;\,5{,}1)$, $m = 5$.

• $\sum x = 10$, $\sum y = 15{,}1$, $\sum x^2 = 0+1+4+9+16 = 30$, $\sum xy = 0 + 1{,}9 + 6{,}4 + 11{,}4 + 20{,}4 = 40{,}1$.
• Equazioni normali: $5a + 10b = 15{,}1$ e $10a + 30b = 40{,}1$.
• $b = \dfrac{5\cdot 40{,}1 - 10\cdot 15{,}1}{5\cdot 30 - 100} = \dfrac{200{,}5 - 151}{50} = 0{,}99$; $a = \dfrac{15{,}1 - 0{,}99\cdot 10}{5} = \dfrac{5{,}2}{5} = 1{,}04$.
• Retta: $p(x) = 1{,}04 + 0{,}99\,x$. Residui $y_k - p(x_k)$: $0{,}06;\ -0{,}13;\ 0{,}18;\ -0{,}21;\ 0{,}10$. Somma $= 0$ ✓ (la prima equazione normale dice che i residui sommano a $0$) e $\sum x_kr_k = 0 - 0{,}13 + 0{,}36 - 0{,}63 + 0{,}40 = 0$ ✓ (seconda equazione: ortogonalità a $x$).
• Somma dei quadrati dei residui: $0{,}0036 + 0{,}0169 + 0{,}0324 + 0{,}0441 + 0{,}0100 = 0{,}1070$.

Passando a una parabola il residuo scende solo a $0{,}1006$ (coefficienti $0{,}0214x^2 + 0{,}9043x + 1{,}0829$): il modello lineare è sufficiente. Con grado $4$ (cinque coefficienti per cinque punti) il polinomio interpola e il residuo è nullo: l'interpolazione è il caso limite dei minimi quadrati con $m = n+1$.

Modelli linearizzabili. Il modello $y = \alpha e^{\beta x}$ non è lineare nei parametri, ma $\ln y = \ln\alpha + \beta x$ sì: si fa la regressione lineare di $\ln y_k$ su $x_k$. (Attenzione: così si minimizzano i residui di $\ln y$, non di $y$.)

Aspetti numerici. Formare $A^TA$ eleva al quadrato il numero di condizionamento: nell'esempio $\text{cond}(A) = 4{,}74$ e $\text{cond}(A^TA) = 22{,}5$. Per polinomi di grado alto questo è disastroso, perciò i programmi risolvono il sistema con la fattorizzazione QR di $A$ (A\y in MATLAB con $A$ rettangolare, np.linalg.lstsq, polyfit), senza formare $A^TA$.

Codice

import numpy as np

def spline_naturale(x, y):
    x = np.asarray(x, float); y = np.asarray(y, float); n = len(x) - 1
    h = np.diff(x)
    A = np.zeros((n + 1, n + 1)); r = np.zeros(n + 1)
    A[0, 0] = A[n, n] = 1                                   # M_0 = M_n = 0
    for i in range(1, n):
        A[i, i-1:i+2] = [h[i-1] / 6, (h[i-1] + h[i]) / 3, h[i] / 6]
        r[i] = (y[i+1] - y[i]) / h[i] - (y[i] - y[i-1]) / h[i-1]
    M = np.linalg.solve(A, r)                               # momenti
    def S(t):
        i = min(max(np.searchsorted(x, t) - 1, 0), n - 1)
        return (M[i] * (x[i+1] - t)**3 + M[i+1] * (t - x[i])**3) / (6 * h[i]) \
            + (y[i] / h[i] - M[i] * h[i] / 6) * (x[i+1] - t) \
            + (y[i+1] / h[i] - M[i+1] * h[i] / 6) * (t - x[i])
    return S, M

S, M = spline_naturale([0, 1, 2, 3], [0, 1, 0, 1])
print(M, S(0.5), S(1.5), S(2.5))

x = np.array([0, 1, 2, 3, 4.]); y = np.array([1.1, 1.9, 3.2, 3.8, 5.1])
A = np.vander(x, 2)                                         # colonne [x, 1]
c = np.linalg.lstsq(A, y, rcond=None)[0]                    # MATLAB: A\y
r = y - A @ c
print(c, r @ r)
print(A.T @ r)                                              # ortogonalita' del residuo

Output:

[ 0. -4.  4.  0.] 0.7499999999999999 0.49999999999999994 0.25
[0.99 1.04] 0.1070000000000001
[9.76996262e-15 3.10862447e-15]

(La stampa dell'ultima riga mostra che il residuo è ortogonale alle colonne di $A$ a meno di $10^{-14}$.)

Errori tipici

• Dimenticare le due condizioni agli estremi: senza di esse il sistema dei momenti ha $n-1$ equazioni per $n+1$ incognite.
• Usare la spline naturale quando $f''$ agli estremi non è nulla e aspettarsi $O(h^4)$ fino ai bordi.
• Confondere interpolazione (passa per tutti i punti) e minimi quadrati (minimizza la somma dei quadrati dei residui, con pochi parametri).
• Risolvere i minimi quadrati con $A^TA$ per gradi alti: il condizionamento si eleva al quadrato.
• Scrivere nelle equazioni normali $\sum x_k^2$ come $\left(\sum x_k\right)^2$.

Domande d'esame

1. Quante condizioni servono per determinare una spline cubica interpolante su $n+1$ nodi e come si assegnano quelle mancanti? Risposta: $4n$ coefficienti, $4n-2$ condizioni (interpolazione e $C^2$); due in più agli estremi: naturale ($S''=0$), vincolata ($S'$ assegnata) o not-a-knot.
2. Stimare l'errore dell'interpolazione lineare a tratti. Risposta: $\frac{h^2}{8}\max|f''|$.
3. Scrivere le equazioni normali e interpretarle geometricamente. Risposta: $A^TAc = A^Ty$; il residuo $y - Ac$ è ortogonale alle colonne di $A$ (proiezione ortogonale di $y$).
4. Perché si preferisce la fattorizzazione QR alle equazioni normali? Risposta: $\text{cond}(A^TA) = \text{cond}(A)^2$.