Polinomi - rappresentazione per coefficienti e per valori e convoluzione

Il capitolo sulla FFT parte da un problema algebrico: moltiplicare due polinomi. L'idea è avere due rappresentazioni dello stesso oggetto, ognuna con operazioni facili diverse, e convertire in fretta da una all'altra. Qui si fissano le due rappresentazioni, i costi esatti delle operazioni e delle conversioni generali. Le basi su induzione e ricorsioni sono in Dimostrazioni, induzione e invariantiTecniche di dimostrazione (esempio, controesempio, assurdo), induzione con casi base multipli, invarianti di ciclo (inizializzazione, conservazione, uso alla fine), schema generale per provare la correttezza; esempi svolti su arrayMax, sequenza di bit e numeri di Fibonacci.Dimostrazioni, induzione e invarianti → e Algoritmi ricorsiviAlgoritmi ricorsivi come induzione eseguita; albero della ricorsione e record di attivazione nello stack; esempi ReverseArray, LinearSum, Power in tempo logaritmico, Fibonacci ricorsivo (esponenziale) e con memoizzazione; analisi della complessità con l'albero della ricorsione e correttezza per induzione.Algoritmi ricorsivi →; le sommatorie in SommatorieIl simbolo di sommatoria, le sue proprietà (linearità, additività, cambio di indice) e le somme notevoli di Gauss e geometrica.Sommatorie →.

I polinomi sono funzioni $p:\mathbb C\to\mathbb C$, con coefficienti complessi (richiami sui complessi in Numeri complessi, formula di Eulero ed esponenziali complessiUn numero complesso si scrive in forma cartesiana a+jb o polare |x|e^{jφ}; il prodotto moltiplica i moduli e somma le fasi. La formula di Eulero e^{jα}=cos α + j sin α lega esponenziali e sinusoidi e permette di trattare tutti i segnali del corso come somme di esponenziali complessi e^{(σ+jω)t}.Numeri complessi, formula di Eulero ed esponenziali complessi →). Il costo di un algoritmo è, in tutto il capitolo, il numero di operazioni aritmetiche tra scalari complessi (somme, differenze, prodotti, divisioni, ciascuna per 1).

Rappresentazione per coefficienti

Definizione (polinomio di grado limitato da $n$). Il polinomio $$p(x)=\sum_{i=0}^{n-1}a_i x^i=a_0+a_1x+\dots+a_{n-1}x^{n-1}$$ è rappresentato con $n$ coefficienti dal vettore $\mathbf a=(a_0,\dots,a_{n-1})\in\mathbb C^n$ e si scrive $p(x)\equiv\mathbf a$. Il numero $n$ è il limite di grado (degree bound); il grado effettivo è $\max\{i:a_i\neq 0\}\le n-1$.

Il vettore si può estendere con zeri (padding) senza cambiare il polinomio: $p(x)=x^2\equiv(0,0,1)\in\mathbb C^3$ e anche $(0,0,1,0)\in\mathbb C^4$. Le rappresentazioni per coefficienti dello stesso polinomio sono quindi infinite. Il padding serve agli algoritmi che pretendono una taglia fissa, per esempio una potenza di due (la FFT).

Rappresentazione per valori e teorema di interpolazione

Teorema (interpolazione). Date $n$ coppie $(x_i,y_i)\in\mathbb C\times\mathbb C$ con $x_i\ne x_j$ per $i\ne j$, esiste uno e un solo polinomio $p$ di grado limitato da $n$ tale che $p(x_i)=y_i$ per $0\le i\le n-1$.

Dimostrazione. Imporre $p(x_i)=y_i$ significa risolvere il sistema lineare $\sum_{j=0}^{n-1}a_jx_i^j=y_i$ ($i=0,\dots,n-1$), cioè $X\mathbf a=\mathbf y$ con $X_{ij}=x_i^{\,j}$. La matrice $X$ è una matrice di Vandermonde $V(x_0,\dots,x_{n-1})$ (riga $i$ = potenze consecutive di $x_i$). Se $\det V\ne0$ il sistema ha una e una sola soluzione (vedi Matrice inversaL'inversa di una matrice quadrata A è la matrice A⁻¹ con A A⁻¹ = A⁻¹ A = I; esiste se e solo se rango(A) = n e si calcola con Gauss-Jordan riducendo (A | I) fino a (I | A⁻¹).Matrice inversa → e Sistemi lineari e teorema di Rouché-CapelliUn sistema lineare si scrive AX = B; ha soluzioni se e solo se B sta nell'immagine di A, cioè se rango(A) = rango(A|B) (Rouché-Capelli); le soluzioni sono una soluzione particolare più il nucleo e dipendono da n − r parametri.Sistemi lineari e teorema di Rouché-Capelli →). Il determinante vale $$\det V(x_0,\dots,x_{n-1})=\prod_{0\le i<j\le n-1}(x_j-x_i),$$ che è diverso da zero perché i punti sono distinti. $\square$

Dimostrazione del determinante, per induzione su $n$. Per $n=1$, $V=(1)$ e il prodotto è vuoto: vale $1$. Passo: in $V(x_0,\dots,x_{n-1})$ si sottrae, partendo dall'ultima colonna e andando verso sinistra, $x_0$ volte la colonna precedente: $C_j\leftarrow C_j-x_0C_{j-1}$ per $j=n-1,\dots,1$. Questa operazione non cambia il determinante (vedi Proprietà del determinante, Binet e LaplaceIl determinante è lineare in ogni riga, cambia segno scambiando due righe, non cambia sommando a una riga un multiplo di un'altra: così si calcola con Gauss riducendo a triangolare. Binet: det(AB) = det A · det B. Laplace: sviluppo lungo una riga o colonna con i complementi algebrici.Proprietà del determinante, Binet e Laplace →). Il nuovo elemento $(i,j)$ è $x_i^j-x_0x_i^{j-1}=x_i^{j-1}(x_i-x_0)$. La riga $0$ diventa $(1,0,\dots,0)$, quindi si sviluppa secondo la prima riga e resta il determinante della matrice $(n-1)\times(n-1)$ di elementi $x_i^{j-1}(x_i-x_0)$ con $i,j=1,\dots,n-1$. Dalla riga $i$ si porta fuori il fattore $(x_i-x_0)$: $$\det V(x_0,\dots,x_{n-1})=\prod_{i=1}^{n-1}(x_i-x_0)\cdot\det V(x_1,\dots,x_{n-1}),$$ e per ipotesi induttiva l'ultimo determinante è $\prod_{1\le i<j}(x_j-x_i)$. $\square$

Definizione (rappresentazione per valori). Dato $\mathbf x\in\mathbb C^n$ a componenti distinte (la base), il polinomio $p$ di grado limitato da $n$ è rappresentato dalla coppia $(\mathbf x,\mathbf y)$ con $y_i=p(x_i)$. Basi diverse danno rappresentazioni diverse. L'unicità è garantita dal teorema.

Esempio. $p(x)=x^2$ nella base $\mathbf x=(0,1,2)$ è $\mathbf y=(0,1,4)$. Anche questa rappresentazione si può estendere: aggiungendo il punto $x_3=3$ si aggiunge $y_3=9$, ma il padding è più elaborato di quello per coefficienti (bisogna calcolare un nuovo valore).

Somma e prodotto per coefficienti

Somma. $s(x)=p(x)+q(x)=\sum_i(a_i+b_i)x^i$, quindi $\mathbf s=\mathbf a+\mathbf b$ e il costo è $T_S(n)=n$ somme. Con limiti di grado diversi $m\le n$ bastano $m$ somme (il resto si copia): $\Theta(\min\{m,n\})$.

Prodotto. Se $p\equiv\mathbf a$ e $q\equiv\mathbf b$ hanno limite di grado $n$, il prodotto $c(x)=p(x)q(x)$ ha grado effettivo al più $(n-1)+(n-1)=2n-2$: il suo limite di grado è $2n-1$. Il coefficiente $c_j$ raccoglie tutti i prodotti $a_kb_h$ con $k+h=j$, cioè $a_kb_{j-k}$ con entrambi gli indici in $[0,n-1]$: $$\begin{cases}0\le k\le n-1\\0\le j-k\le n-1\end{cases}\iff\max\{0,j-n+1\}\le k\le\min\{j,n-1\}.$$

Definizione (convoluzione lineare). Per $\mathbf a,\mathbf b\in\mathbb C^n$, il vettore $\mathbf c=\mathbf a\star\mathbf b\in\mathbb C^{2n-1}$ ha componenti $$c_j=\sum_{k=\max\{0,\,j-n+1\}}^{\min\{j,\,n-1\}}a_k\,b_{j-k},\qquad 0\le j\le 2n-2.$$ Se $p\equiv\mathbf a$ e $q\equiv\mathbf b$ allora $p(x)\,q(x)\equiv\mathbf a\star\mathbf b$. Con vettori di lunghezza $m$ e $n$ il risultato ha lunghezza $m+n-1$.

Esempio. $\mathbf a=(1,2,3)$, $\mathbf b=(4,5,6)$ ($n=3$, $\mathbf c\in\mathbb C^5$). Si ordinano i prodotti $a_kb_h$ per $k+h$:

Quindi $(1+2x+3x^2)(4+5x+6x^2)=4+13x+28x^2+27x^3+18x^4$ (controllato con numpy.convolve). Il numero di termini nella somma di $c_j$ è $\min\{j,n-1\}-\max\{0,j-n+1\}+1$: cresce da $1$ a $n$ e poi scende a $1$.

Teorema (costo della convoluzione per definizione). Calcolare $\mathbf a\star\mathbf b$ con $\mathbf a,\mathbf b\in\mathbb C^n$ richiede $n^2$ prodotti e $(n-1)^2$ somme, in tutto $2n^2-2n+1=\Theta(n^2)$ operazioni.

Dimostrazione. Ogni coppia $(k,h)$ con $0\le k,h\le n-1$ contribuisce a esattamente una componente, $c_{k+h}$: i prodotti sono $n^2$. Le componenti sono $2n-1$ e ogni componente con $t$ addendi richiede $t-1$ somme; la somma totale degli addendi è $n^2$, quindi le somme sono $n^2-(2n-1)=(n-1)^2$. $\square$

Per $n=3$: $9$ prodotti e $4$ somme (verificato contando le operazioni di un'implementazione in Python per $n=1,\dots,11$). Per $n=8$: $64+49=113$.

L'operazione $\star$ è commutativa e associativa (perché lo è il prodotto di polinomi): si usa nell'Esponenziazione veloce per quadrature successivepotenza x^n con un operatore associativo con elemento neutro, costo = numero di applicazioni dell'operatore; STUPID_POWER costa n-1, la sottostruttura x^n = (x^floor(n/2))^2 (per n dispari per x) dà POWER con T(n) = floor(log2 n) + (numero di bit a 1 di n) - 1, Theta(log n) (caso peggiore 2 floor(log2 n), n = 2^t - 1); versioni iterative POWER2 e POWER3 dalla rappresentazione binaria di n; esempio con la matrice di Fibonacci; il modello a costo unitario vale solo se l'operatore non fa crescere gli operandi.Esponenziazione veloce per quadrature successive → applicata ai polinomi.

Somma e prodotto per valori

Due polinomi si possono sommare o moltiplicare per valori solo se sono rappresentati sulla stessa base $\mathbf x$. Allora $$s(x_i)=p(x_i)+q(x_i)\ \Rightarrow\ \mathbf y_s=\mathbf y_p+\mathbf y_q,\qquad c(x_i)=p(x_i)\,q(x_i)\ \Rightarrow\ \mathbf y_c=\mathbf y_p\odot\mathbf y_q,$$ dove $\odot$ è il prodotto componente per componente. Il costo è $\Theta(n)$ per entrambe: il prodotto per valori è lineare, contro il quadratico per coefficienti.

Attenzione al limite di grado. Il prodotto ha limite di grado $2n-1$ e quindi per essere determinato servono almeno $2n-1$ punti. Se i due fattori sono dati su soli $n$ punti, $\mathbf y_p\odot\mathbf y_q$ non rappresenta il prodotto. Controesempio: $p=q=1+x$ sulla base $(0,1)$ dà $\mathbf y_p=\mathbf y_q=(1,2)$ e $\mathbf y_p\odot\mathbf y_q=(1,4)$, ma il polinomio di grado limitato da $2$ per $(0,1)$ e $(1,4)$ è $1+3x$, mentre $(1+x)^2=1+2x+x^2$. Su una base estesa a $3$ punti, $(0,1,2)$, si ha $\mathbf y_p=(1,2,3)$ e $\mathbf y_p\odot\mathbf y_q=(1,4,9)$, che sono proprio i valori di $(1+x)^2$ in $0,1,2$.

Conversione coefficienti $\to$ valori: la valutazione

Per valutare $p\equiv\mathbf a$ in un punto $z$ si può applicare la definizione:

V(a, z)
n <- a.len
y <- a[0]; pow <- 1
for i <- 1 to n-1
    pow <- pow * z
    y <- y + a[i] * pow
return y

Dopo il primo comando dell'iterazione $i$, pow vale $z^i$, quindi y accumula $\sum_{j\le i}a_jz^j$. Ogni iterazione fa una somma e due prodotti: $T_V(n)=3(n-1)$.

Regola di Horner. Si riscrive il polinomio annidando le potenze: $$p(z)=a_0+z\bigl(a_1+z(a_2+\dots+z(a_{n-2}+z\,a_{n-1})\dots)\bigr).$$

HORNER(a, z)
n <- a.len
y <- a[n-1]
for j <- 2 to n
    y <- a[n-j] + z * y
return y

Correttezza (invariante). Alla fine dell'iterazione $j$ ($2\le j\le n$) vale $y=\sum_{i=n-j}^{n-1}a_i\,z^{\,i-n+j}$. Prima del ciclo, $y=a_{n-1}$, che è la formula con "$j=1$". Se vale per $j-1$, cioè $y=\sum_{i=n-j+1}^{n-1}a_iz^{\,i-n+j-1}$, allora il nuovo valore è $a_{n-j}+z\,y=a_{n-j}+\sum_{i=n-j+1}^{n-1}a_iz^{\,i-n+j}$, e il termine isolato $a_{n-j}$ è quello con $i=n-j$ ed esponente $0$: la formula vale per $j$. A $j=n$ si ha $y=\sum_{i=0}^{n-1}a_iz^i=p(z)$. $\square$

Ogni iterazione fa una somma e un prodotto: $T_H(n)=2(n-1)$, un terzo in meno di $V$.

Esempio. $p(x)=3+4x+5x^2+6x^3$ in $z=2$: $y=6$; $j=2$: $5+2\cdot6=17$; $j=3$: $4+2\cdot17=38$; $j=4$: $3+2\cdot38=79$. Controllo: $3+8+20+48=79$. Sono $6=2\cdot3$ operazioni contro le $9$ di V.

La valutazione su una base $\mathbf x\in\mathbb C^n$ ripete Horner $n$ volte:

VAL(a, x)
n <- a.len
for i <- 0 to n-1
    y[i] <- HORNER(a, x[i])
return y

con costo $T_{VAL}(n)=n\cdot2(n-1)=\Theta(n^2)$. Non esiste di meglio per una base generica: la FFT sfrutterà una base speciale.

Conversione valori $\to$ coefficienti: l'interpolazione

Sistema lineare. Si risolve $V(\mathbf x)\,\mathbf a=\mathbf y$ con il metodo di Gauss-Jordan (vedi 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 →): $n$ pivot, ciascuno richiede un numero costante di operazioni su $n^2$ elementi, quindi $T_I(n)=\Theta(n^3)$.

Esempio. Per $(0,0),(1,1),(2,4)$ il sistema è $a_0=0$, $a_0+a_1+a_2=1$, $a_0+2a_1+4a_2=4$, con unica soluzione $\mathbf a=(0,0,1)$ cioè $p(x)=x^2$.

Formula di Lagrange. Il polinomio interpolante è $$p(x)=\sum_{k=0}^{n-1}y_k\,\frac{\prod_{j\ne k}(x-x_j)}{\prod_{j\ne k}(x_k-x_j)}.$$ Ogni frazione vale $1$ in $x_k$ e $0$ negli altri $x_j$, quindi $p(x_j)=y_j$, e ha grado $n-1$: per l'unicità è il polinomio cercato. I coefficienti si ottengono in $\Theta(n^2)$ così:

1. si calcolano i coefficienti di $Q(x)=\prod_{j=0}^{n-1}(x-x_j)$, di grado $n$: si moltiplica per un fattore $(x-x_j)$ alla volta, e ogni moltiplicazione costa $O(n)$, in totale $O(n^2)$;
2. per ogni $k$ si calcolano i coefficienti di $Q_k(x)=Q(x)/(x-x_k)$ con la regola di Ruffini (divisione per un binomio, $O(n)$);
3. si calcola $Q_k(x_k)$ con Horner ($O(n)$) e si somma $\dfrac{y_k}{Q_k(x_k)}\,Q_k(x)$ ai coefficienti di $p$ ($O(n)$).

Le fasi 2 e 3 si ripetono per $n$ valori di $k$, ognuna di costo $O(n)$: $O(n^2)$ in tutto.

Esempio. Punti $(0,0),(1,1),(2,4)$. $Q=x(x-1)(x-2)=x^3-3x^2+2x$. Allora $Q_0=x^2-3x+2$ con $Q_0(0)=2$; $Q_1=x^2-2x$ con $Q_1(1)=-1$; $Q_2=x^2-x$ con $Q_2(2)=2$. Quindi $$p=0\cdot\tfrac{Q_0}{2}+1\cdot\tfrac{Q_1}{-1}+4\cdot\tfrac{Q_2}{2}=-(x^2-2x)+2(x^2-x)=x^2.$$ (Implementazione con Fraction provata su polinomi casuali di $1$-$8$ coefficienti.)

Perché serve una base speciale

Le due conversioni generali costano $\Theta(n^2)$ (valutazione e Lagrange) o più. Il percorso "coefficienti $\to$ valori su una base estesa $\to$ prodotto $\odot$ $\to$ coefficienti" costerebbe quindi $\Theta(n^2)$, come la convoluzione diretta: il prodotto lineare per valori non basta. Se però esistesse una base $\mathbf x$ per cui valutazione e interpolazione costano $\Theta(n\log n)$, il prodotto di polinomi costerebbe $\Theta(n\log n)$. Questa base esiste: le radici dell'unità, e la trasformazione si chiama DFT (vedi DFT, matrice di Fourier e radici n-esime dell'unitàla DFT di ordine n valuta un polinomio di grado limitato da n sulle n radici n-esime dell'unità, y_i = p(ω^i) con ω = e^{2πi/n}; è il prodotto per la matrice di Fourier F_n, una Vandermonde invertibile. Le radici formano un gruppo moltiplicativo e valgono tre lemmi: cancellazione (ω_{dn}^{dk} = ω_n^k, ω_n^{n/2} = -1), dimezzamento (i quadrati delle n radici sono le n/2 radici n/2-esime, ciascuna due volte) e somma (la somma delle potenze di ω^k vale n se k è multiplo di n, 0 altrimenti). Teorema di convoluzione lineare: (c|0) = DFT^{-1}{2n}(DFT{2n}(a|0) ⊙ DFT_{2n}(b|0)).DFT, matrice di Fourier e radici n-esime dell'unità →), calcolata con la FFT (Fast Fourier Transform - algoritmo ricorsivo e trasformata inversala FFT calcola DFT_n (n potenza di due) con divide et impera: si separano i coefficienti di indice pari e dispari, p(x) = p0(x^2) + x p1(x^2), e per il lemma di dimezzamento bastano due DFT di ordine n/2; y_i = y0_i + w^i y1_i e y_{i+n/2} = y0_i - w^i y1_i. La ricorrenza T(n) = 2T(n/2) + 2n con T(1) = 0 ha soluzione esatta 2 n log2 n. L'inversa ha matrice F^{-1}_{ij} = w^{-ij}/n: si ottiene con la FFT sulle potenze negative (poi si divide per n fuori dalla ricorsione) oppure con DFT_n(y) rovesciato. LIN_CONV moltiplica due polinomi con n coefficienti in 12 n log2(2n) + 4n operazioni.Fast Fourier Transform - algoritmo ricorsivo e trasformata inversa →).

Il punto di vista del segnale è lo stesso, a parole diverse: la convoluzione lineare di due sequenze finite è la risposta di un sistema FIR (Sistemi LTI e convoluzione discretaUn sistema lineare e tempo-invariante (LTI) è completamente descritto dalla sua risposta impulsiva h[n]: y[n] = Σ x[k] h[n-k] = x[n] * h[n] (somma di convoluzione). Si ricava scomponendo x in impulsi traslati e usando linearità e invarianza. La convoluzione è commutativa, associativa e distributiva; una cascata di LTI equivale a un solo filtro con h = h1 * h2 e l'ordine non conta. Se x ha lunghezza N e h lunghezza L, y ha lunghezza N+L-1. Si calcola col metodo della finestra scorrevole o a tabella. I FIR sono LTI; altri sistemi (n·x[n], x[-n], x²) non lo sono.Sistemi LTI e convoluzione discreta →, Calcolo della convoluzione e sue proprietàIl supporto della convoluzione è la somma dei supporti, $\operatorname{rect}*\operatorname{rect}=\Lambda$, e due esponenziali causali danno $(e^{-bt}-e^{-at})/(a-b)$. Si calcola con il metodo grafico a casi (ribaltare, traslare, individuare gli intervalli di sovrapposizione). Proprietà: lineare, commutativa, associativa, $\delta$ è l'elemento neutro, la traslazione si somma, l'area è il prodotto delle aree.Calcolo della convoluzione e sue proprietà →) e coincide con il prodotto dei polinomi nella variabile $z^{-1}$.

Domande d'esame

• Per vettori di $8$ componenti, quante moltiplicazioni e somme richiede la convoluzione lineare per definizione? Risposta: $n^2=64$ prodotti e $(n-1)^2=49$ somme, $113$ operazioni.
• Quanti punti servono per rappresentare per valori il prodotto di un polinomio con $5$ coefficienti per uno con $4$ coefficienti? Perché non bastano $5$? Risposta: il prodotto ha limite di grado $5+4-1=8$, servono almeno $8$ punti; con $5$ punti il prodotto componente per componente determina solo un polinomio di grado limitato da $5$, diverso dal prodotto.
• Quante operazioni richiede la regola di Horner per un polinomio con $n$ coefficienti? E la valutazione su $n$ punti? Risposta: $2(n-1)$ per un punto, $2n(n-1)=\Theta(n^2)$ per $n$ punti.

Esercizi svolti: Esercizio - matrice V per V trasposta con V di Vandermonde, formula degli elementi e calcolo in tempo quadratico (compito del 4 settembre 2024, es. 2); Esercizio - convoluzione con il vettore delle potenze di due in tempo lineare (compito del 5 settembre 2025, es. 2); Esercizio - domande di prima parte su polinomi, DFT e FFT (simulazioni d'esame).