Esercizio - prodotto di matrici circolanti e convoluzione ciclica con la DFT (esercizi su supercircolanti e circolanti, es. CR 10)

Esercizio supplementare 10 (con il 14) sulle matrici circolanti, dalle note ed esercizi del corso di Algoritmi per l'Ingegneria (Ing. Informatica UniPD) su supercircolanti e circolanti. Teoria usata: Matrici circolanti e supercircolanti - algoritmi divide and conquerle matrici supercircolanti (quadranti alternati ricorsivi) hanno solo n valori distinti: somma in n operazioni, prodotto in 4n log2 n + n con il trucco (A1+A2)(B1+B2) e (A1-A2)(B1-B2) invece delle 4 chiamate (Theta(n^2)); gli autovalori sono la trasformata di Hadamard della prima riga (FHT, n log2 n) e il prodotto passa dalla base diagonalizzante; le matrici circolanti sono polinomi nella matrice di rotazione R, il loro prodotto è la convoluzione ciclica delle prime righe e si diagonalizzano con la matrice di Fourier: prodotto in O(n log n) con la FFT.Matrici circolanti e supercircolanti - algoritmi divide and conquer →, Polinomi - rappresentazione per coefficienti e per valori e convoluzioneun polinomio di grado limitato da n si rappresenta con n coefficienti o con n valori su n punti distinti (teorema di interpolazione, matrice di Vandermonde). Per coefficienti la somma costa n e il prodotto, la convoluzione lineare, n^2 prodotti e (n-1)^2 somme; per valori il prodotto costa 2n-1 prodotti ma serve una base estesa di almeno 2n-1 punti. Le conversioni generali costano troppo: valutazione n volte con Horner 2n(n-1), interpolazione con Gauss-Jordan n^3 o con Lagrange n^2. Serve una base speciale: la DFT.Polinomi - rappresentazione per coefficienti e per valori e convoluzione →, 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à → e 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 →.

Testo

Una matrice $n \times n$ è circolante (CR) se ogni riga è la rotazione di una posizione verso destra della precedente; si rappresenta con la prima riga $r(A) = [A_{00}, \dots, A_{0,n-1}]$. Sia $R$ la CR con prima riga $[0, 1, 0, \dots, 0]$ (rotazione circolare a sinistra, $R^n = I$). Date due sequenze $a, b$ di lunghezza $n$, la convoluzione ciclica $c = a \circledast b$ è la sequenza con $p_c(z) = p_a(z)\,p_b(z) \bmod (z^n - 1)$, equivalentemente $c_k = \sum_{i=0}^{n-1} a_i\,b_{(k-i) \bmod n}$.

• 10. Dimostrare che, se $A$ e $B$ sono matrici CR, allora $C = AB$ è la matrice CR con prima riga $r(C) = r(A) \circledast r(B)$. Dunque moltiplicare matrici CR equivale a convolverne ciclicamente le prime righe. Suggerimento: si scrivano $A$ e $B$ come polinomi nella matrice $R$, $A = \sum_{k=0}^{n-1} A_{0k}R^k$, e si sfrutti la relazione $R^n = I$.
• 14. Dimostrare che gli autovalori $(\alpha_0, \dots, \alpha_{n-1})$ di una matrice circolante $A$ formano la DFT del vettore $r(A)^T$, cioè $(\alpha_0, \dots, \alpha_{n-1})^T = F^{(n)}\,r(A)^T$ (l'argomento ricalca quello sviluppato per le matrici SC).

Per fissare le idee si svolge anche un esempio numerico completo con $n = 4$, $r(A) = (1, 2, 0, 3)$ e $r(B) = (2, 1, 1, 0)$: prodotto per definizione, convoluzione ciclica, e prodotto passando dalla DFT.

Preliminare: $A = \sum_k A_{0k} R^k$

La matrice $R$ ha $R_{h,j} = 1$ se $j = (h + 1) \bmod n$ e $0$ altrimenti (infatti la sua prima riga è $e_1$ e ogni riga è la precedente ruotata a destra). Per induzione $(R^k)_{h,j} = 1$ se $j = (h + k) \bmod n$ e $0$ altrimenti: $R^k$ sposta di $k$ posizioni. Allora

$$\Bigl( \sum_{k=0}^{n-1} A_{0k} R^k \Bigr)_{h,j} = A_{0,(j-h) \bmod n} = A_{h,j},$$

perché solo il termine con $k = (j - h) \bmod n$ ha un $1$ in posizione $(h, j)$, e $A_{h,j} = A_{0,(j-h) \bmod n}$ per la definizione di CR. Quindi ogni CR è un polinomio in $R$ di grado minore di $n$ con coefficienti $r(A)$.

Punto 10: $r(AB) = r(A) \circledast r(B)$

Siano $A = \sum_i a_i R^i$ e $B = \sum_j b_j R^j$, con $a = r(A)$, $b = r(B)$. Il prodotto di polinomi nella stessa matrice si espande come per i polinomi ordinari (le potenze di $R$ commutano con tutto), con $R^i R^j = R^{i+j}$:

$$AB = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_i b_j R^{i+j}.$$

Siccome $R^n = I$, vale $R^{m} = R^{m \bmod n}$ per ogni $m \ge 0$. Raggruppando i termini con lo stesso $k = (i + j) \bmod n$:

$$AB = \sum_{k=0}^{n-1} \Bigl( \sum_{\substack{i, j \\ i + j \equiv k \,(\mathrm{mod}\ n)}} a_i b_j \Bigr) R^k = \sum_{k=0}^{n-1} c_k R^k, \qquad c_k = \sum_{i=0}^{n-1} a_i\,b_{(k-i) \bmod n},$$

dove per ogni $i$ c'è un solo $j \in \{0, \dots, n-1\}$ con $i + j \equiv k$, cioè $j = (k - i) \bmod n$. Il membro destro è un polinomio in $R$, quindi una CR (le CR sono tutte e sole le combinazioni lineari delle potenze di $R$), con prima riga $(c_0, \dots, c_{n-1}) = a \circledast b$. Dunque $AB$ è CR e $r(AB) = r(A) \circledast r(B)$. $\blacksquare$

Equivalenza delle due definizioni di convoluzione ciclica. $p_a(z)p_b(z) = \sum_{i,j} a_ib_j z^{i+j}$, e ridurre modulo $z^n - 1$ significa sostituire $z^n$ con $1$, quindi $z^{i+j}$ con $z^{(i+j) \bmod n}$: il coefficiente di $z^k$ è la somma degli $a_ib_j$ con $i + j \equiv k$, cioè $c_k$ sopra.

Punto 14: gli autovalori sono la DFT

Si usano il polinomio in $R$ e la diagonalizzazione di $R$. Sia $\omega = \omega_n = e^{i2\pi/n}$ e, per $j = 0, \dots, n-1$, $w_j = (\omega^{0j}, \omega^{1j}, \dots, \omega^{(n-1)j})^T$.

1. $R w_j = \omega^j w_j$. La componente $k$ di $Rw_j$ è la componente $k + 1$ di $w_j$ (l'effetto di $R$ è $Rx = (x_1, \dots, x_{n-1}, x_0)$): $\omega^{(k+1)j} = \omega^j \omega^{kj}$ per $k < n - 1$, e per $k = n - 1$ si ha $w_{j,0} = 1 = \omega^{nj} = \omega^j\omega^{(n-1)j}$.
2. $A w_j = p_A(\omega^j) w_j$, con $p_A(z) = \sum_k A_{0k}z^k$: $A w_j = \sum_k A_{0k} R^k w_j = \sum_k A_{0k} \omega^{jk} w_j$.
3. Quindi $w_j$ è autovettore di ogni CR, con autovalore $\alpha_j = \sum_{k=0}^{n-1} A_{0k}\,\omega^{jk} = (F^{(n)} r(A)^T)_j$, dove $F^{(n)}$ è la matrice di Fourier, $F_{jk} = \omega^{jk}$ (simmetrica; le sue colonne sono i $w_j$). Gli $n$ autovettori $w_0, \dots, w_{n-1}$ sono indipendenti (matrice di Vandermonde sui punti distinti $\omega^0, \dots, \omega^{n-1}$), quindi $A = F\,\mathrm{diag}(\alpha)\,F^{-1}$ e questi sono tutti gli autovalori. $\blacksquare$

Seconda strada (come per le SC): se $v$ è autovettore di una matrice con prima riga $r(A)$ e $v_0 \ne 0$, dalla prima componente di $Av = \alpha v$ si ha $r(A) \cdot v = \alpha v_0$; con $v = w_j$, $v_0 = 1$, quindi $\alpha_j = r(A) \cdot w_j = \sum_k A_{0k}\omega^{jk}$.

Conseguenza (teorema di convoluzione ciclica in forma matriciale). Gli autovalori di $AB$ sono i prodotti $\alpha_j\beta_j$ (stessa base $F$), quindi

$$r(C)^T = \bigl(F^{(n)}\bigr)^{-1} \bigl[ \bigl(F^{(n)} r(A)^T\bigr) \otimes \bigl(F^{(n)} r(B)^T\bigr) \bigr],$$

con $\otimes$ prodotto componente per componente: la convoluzione ciclica si calcola con tre DFT e $n$ prodotti, cioè $O(n \log n)$ 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 →). (Le note del corso scrivono "FHT" in questo punto; serve la FFT, perché la trasformata di Hadamard diagonalizza solo le supercircolanti.)

Esempio numerico ($n = 4$)

$r(A) = (1, 2, 0, 3)$, $r(B) = (2, 1, 1, 0)$:

$$A = \begin{pmatrix} 1 & 2 & 0 & 3 \\ 3 & 1 & 2 & 0 \\ 0 & 3 & 1 & 2 \\ 2 & 0 & 3 & 1 \end{pmatrix}, \qquad B = \begin{pmatrix} 2 & 1 & 1 & 0 \\ 0 & 2 & 1 & 1 \\ 1 & 0 & 2 & 1 \\ 1 & 1 & 0 & 2 \end{pmatrix}, \qquad A = I + 2R + 3R^3, \; B = 2I + R + R^2.$$

1. Per definizione (prima riga di $AB$ = prima riga di $A$ per le colonne di $B$). Le colonne di $B$ sono $(2, 0, 1, 1)$, $(1, 2, 0, 1)$, $(1, 1, 2, 0)$, $(0, 1, 1, 2)$ e la prima riga di $A$ è $(1, 2, 0, 3)$:

$$c_0 = 1\cdot 2 + 2 \cdot 0 + 0 \cdot 1 + 3 \cdot 1 = 5, \quad c_1 = 1 + 4 + 0 + 3 = 8, \quad c_2 = 1 + 2 + 0 + 0 = 3, \quad c_3 = 0 + 2 + 0 + 6 = 8.$$

2. Convoluzione ciclica, $c_k = \sum_i a_i b_{(k-i) \bmod 4}$ con $a = (1, 2, 0, 3)$, $b = (2, 1, 1, 0)$:

• $c_0 = a_0b_0 + a_1b_3 + a_2b_2 + a_3b_1 = 2 + 0 + 0 + 3 = 5$;
• $c_1 = a_0b_1 + a_1b_0 + a_2b_3 + a_3b_2 = 1 + 4 + 0 + 3 = 8$;
• $c_2 = a_0b_2 + a_1b_1 + a_2b_0 + a_3b_3 = 1 + 2 + 0 + 0 = 3$;
• $c_3 = a_0b_3 + a_1b_2 + a_2b_1 + a_3b_0 = 0 + 2 + 0 + 6 = 8$.

Stesso risultato, $r(AB) = (5, 8, 3, 8)$ ✓.

3. Con la DFT ($\omega_4 = i$, $F_{jk} = i^{jk}$):

$$F^{(4)} = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{pmatrix}.$$

$F a = (1 + 2 + 0 + 3,\; 1 + 2i + 0 - 3i,\; 1 - 2 + 0 - 3,\; 1 - 2i + 0 + 3i) = (6, 1 - i, -4, 1 + i)$: questi sono gli autovalori di $A$. Controlli: la somma è $6 + (1 - i) - 4 + (1 + i) = 4 = \mathrm{tr}(A)$; il prodotto è $6 \cdot (-4) \cdot (1 - i)(1 + i) = -48 = \det A$. Con $w_1 = (1, i, -1, -i)$: la prima riga di $A$ dà $1 + 2i + 0 - 3i = 1 - i$ ✓ e la seconda riga $(3, 1, 2, 0)$ dà $3 + i - 2 + 0 = 1 + i = (1 - i)\,i$ ✓.

$F b = (2 + 1 + 1 + 0,\; 2 + i - 1 + 0,\; 2 - 1 + 1 - 0,\; 2 - i - 1 + 0) = (4, 1 + i, 2, 1 - i)$.

Prodotto componente per componente: $y = (6 \cdot 4,\; (1 - i)(1 + i),\; (-4) \cdot 2,\; (1 + i)(1 - i)) = (24, 2, -8, 2)$.

Antitrasformata: $c = \frac 14\,\overline{F}\,y$, cioè $c_k = \frac 14 \sum_j y_j\,i^{-jk}$:

$$c_0 = \tfrac{24 + 2 - 8 + 2}{4} = 5, \qquad c_1 = \tfrac{24 + 2(-i) - 8(-1) + 2(i)}{4} = \tfrac{32}{4} = 8,$$

$$c_2 = \tfrac{24 + 2(-1) - 8(1) + 2(-1)}{4} = 3, \qquad c_3 = \tfrac{24 + 2(i) - 8(-1) + 2(-i)}{4} = 8.$$

Ancora $(5, 8, 3, 8)$ ✓: i tre metodi coincidono.

Verifica con un programma. Per $n = 1, 2, 3, 4, 5, 8$ e prime righe casuali, il prodotto delle matrici circolanti piene coincide con $\mathrm{circ}(a \circledast b)$, i prodotti commutano, $R^k$ e la combinazione $\sum a_k R^k$ riproducono la matrice, gli $w_j$ sono autovettori con autovalori $F a$ e la formula con la DFT dà la convoluzione; i valori dell'esempio ($(5, 8, 3, 8)$, $(6, 1 - i, -4, 1 + i)$, $\det A = -48$) sono quelli calcolati da numpy.