Esercizio - matrice di Fourier per una matrice generica a sinistra e a destra e per aI + bU (compito del 27 gennaio 2025, es. 2)

---

Teoria usata: 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à → ($\mathrm{DFT}_n(\mathbf a)=F_n\mathbf a$), 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 → (costo $2n\log_2n$), Trasformate notevoli ed esponenziazione di polinomi con la FFTalcune famiglie di vettori hanno DFT nota o più economica della FFT generale: il vettore costante (a,...,a) ha DFT (na,0,...,0), il vettore (a,0,...,0) ha DFT (a,...,a); un vettore (n,k)-sparso ha DFT formata da k copie della DFT{n/k} delle sue componenti di indice multiplo di k, e SPARSE_FFT costa 2 (n/k) log2(n/k) con ricorrenza T(n,k) = T(n/2,k/2). Per la potenza p(x)^k di un polinomio con n coefficienti: o si eleva a potenza per quadrature successive con LIN_CONV, O(nk log(nk)), oppure si valuta p con la FFT su h = 2^ceil(log2((n-1)k+1)) punti, si eleva alla k ogni valore e si antitrasforma, Theta(nk log(nk)).Trasformate notevoli ed esponenziazione di polinomi con la FFT → (DFT del vettore costante), Operazioni tra matriciLe matrici m×n formano uno spazio vettoriale (somma e prodotto per scalare elemento per elemento); il prodotto righe per colonne corrisponde alla composizione di funzioni lineari, è associativo ma non commutativo; la trasposta scambia righe e colonne e (AB)^T = B^T A^T.Operazioni tra matrici → (prodotto righe per colonne, trasposta di un prodotto).

Il prodotto di due matrici $n\times n$ qualsiasi richiede $\Theta(n^3)$ operazioni con l'algoritmo per definizione e circa $n^{2{,}8}$ con Algoritmo di Strassen e ibridazioneStrassen calcola il prodotto n x n con 7 chiamate ricorsive invece di 8, usando 10 somme/sottrazioni nel divide e 8 nel conquer: T(n) = 7T(n/2) + (9/2)n^2, Theta(n^(log2 7)) con soluzione esatta 7n^(log2 7) - 6n^2; la costante 7 è più che tripla di quella dell'iterativo (2) e Strassen batte l'iterativo solo da n = 655 (per le potenze di 2, da 1024). L'ibridazione cambia il caso base (iterativo per n <= n0): costo n^(log2 7)(2n0+5)/n0^(log2 7 - 2) - 6n^2, minimo per n0 = 8 con costante ~3.92.Algoritmo di Strassen e ibridazione →. Qui un fattore è la matrice di Fourier, che non è generica: moltiplicare $F$ per un vettore è una DFT e costa solo $\Theta(n\log n)$.

(a) $B=FA$ con $A$ qualsiasi

La colonna $j$ di $B=FA$ è $F$ per la colonna $j$ di $A$: $$B_{\bullet j}=F\,A_{\bullet j}=\mathrm{DFT}_n\bigl(A_{\bullet j}\bigr),\qquad j=0,\dots,n-1.$$ Quindi basta applicare la FFT a ciascuna delle $n$ colonne di $A$. Ogni FFT costa $2n\log_2n$ operazioni (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 →), perciò $$T(n)=n\cdot2n\log_2n=2n^2\log_2n=O(n^2\log n).$$ Estrarre le colonne e comporle in $B$ sono copie, che non si contano. Il costo è inferiore a quello di qualunque algoritmo noto per il prodotto di matrici generiche ($n^{2{,}37}$ nel miglior risultato asintotico noto, $n^3$ per definizione). Verificato con un'implementazione che conta le operazioni, per $n=2,4,\dots,32$ con $A$ complessa casuale: risultato uguale a $F\cdot A$ calcolato per definizione e numero di operazioni $2n^2\log_2n$ esatto.

(b) $B=AF$ con $A$ qualsiasi

Qui $F$ moltiplica a destra. Si usa che la trasposta di un prodotto è il prodotto delle trasposte in ordine inverso, $(AF)^T=F^TA^T$, e che $F$ è simmetrica: $F_{ij}=\omega_n^{ij}=\omega_n^{ji}=F_{ji}$, quindi $F^T=F$. Allora $$B^T=(AF)^T=F^TA^T=F\,A^T.$$ Il prodotto $FA^T$ si calcola come in (a), applicando la FFT alle colonne di $A^T$, cioè alle righe di $A$. Infine $B=(FA^T)^T$. Per capirlo senza trasposte: la riga $i$ di $B$ è la riga $i$ di $A$ per $F$, e per la simmetria di $F$ è la DFT della riga $i$ di $A$. Le trasposizioni sono solo copie ($a'_{ji}\leftarrow a_{ij}$), non contate, quindi $$T(n)=2n^2\log_2n=O(n^2\log n),$$ come in (a). (Verificato numericamente che $AF=(FA^T)^T$ per $n=2,\dots,32$.)

(c) $B=FA$ con $A=aI+bU$

La matrice $A$ ha $a+b$ sulla diagonale e $b$ fuori: si scrive $$A=a\,I_n+b\,U_n,\qquad U_n=\text{matrice con tutti gli elementi uguali a }1.$$ Per la linearità del prodotto $$FA=a\,F\,I_n+b\,F\,U_n=a\,F+b\,F\,U_n.$$ Ogni colonna di $U_n$ è il vettore $(1,1,\dots,1)^T$, e per la DFT del vettore costante (Trasformate notevoli ed esponenziazione di polinomi con la FFTalcune famiglie di vettori hanno DFT nota o più economica della FFT generale: il vettore costante (a,...,a) ha DFT (na,0,...,0), il vettore (a,0,...,0) ha DFT (a,...,a); un vettore (n,k)-sparso ha DFT formata da k copie della DFT_{n/k} delle sue componenti di indice multiplo di k, e SPARSE_FFT costa 2 (n/k) log2(n/k) con ricorrenza T(n,k) = T(n/2,k/2). Per la potenza p(x)^k di un polinomio con n coefficienti: o si eleva a potenza per quadrature successive con LIN_CONV, O(nk log(nk)), oppure si valuta p con la FFT su h = 2^ceil(log2((n-1)k+1)) punti, si eleva alla k ogni valore e si antitrasforma, Theta(nk log(nk)).Trasformate notevoli ed esponenziazione di polinomi con la FFT →) $F\,(1,\dots,1)^T=(n,0,\dots,0)^T$. Quindi $F\,U_n$ ha la prima riga fatta di $n$ e tutte le altre nulle: $F\,U_n=n\,\mathbf e_0\,\mathbf 1^T$ (con $\mathbf e_0=(1,0,\dots,0)^T$ e $\mathbf 1=(1,\dots,1)^T$). Perciò $$FA=a\,F+n\,b\,\mathbf e_0\mathbf 1^T,\qquad\text{cioè}\qquad(FA)_{ij}=\begin{cases}a+nb&i=0,\\a\,\omega_n^{ij}&i\ge1.\end{cases}$$ (per $i=0$: $a\cdot\omega_n^{0}+nb=a+nb$.)

Algoritmo e costo. Anche se $F$ ha $n^2$ elementi, essi appartengono all'insieme di soli $n$ valori $\{\omega_n^0,\dots,\omega_n^{n-1}\}$ (perché $\omega_n^{ij}=\omega_n^{ij\bmod n}$). Basta quindi calcolare:

• i prodotti $a\,\omega_n^k$ per $k=1,\dots,n-1$ ($a\cdot\omega_n^0=a$ non richiede operazioni): $n-1$ prodotti;
• $nb$: $1$ prodotto;
• $a+nb$: $1$ somma.

Tutti i $n^2$ elementi di $FA$ sono poi copie di questi $n+1$ valori. Quindi $$T(n)=(n-1)+1+1=n+1=\Theta(n),$$ molto meno dei $n^2$ elementi della matrice (e del $\Theta(n^2\log n)$ di (a)). La soluzione del corso riporta $T(n)=n$: conta $n-1$ prodotti per le potenze e un prodotto per $nb$, ma dimentica la somma $a+nb$ della prima riga; contandola si ottiene $n+1$ (stesso ordine, $\Theta(n)$).

Esempio ($n=4$). Con $\omega_4=i$: $$FA=\begin{pmatrix}a+4b&a+4b&a+4b&a+4b\\a&ia&-a&-ia\\a&-a&a&-a\\a&-ia&-a&ia\end{pmatrix}.$$ Prodotti necessari: $ia$, $-a$, $-ia$ (con $\omega_4^1,\omega_4^2,\omega_4^3$), più $4b$ e $a+4b$. Confrontato con $F\cdot A$ calcolato per definizione con $a,b$ complessi casuali, per $n=2,\dots,32$.