Algoritmo di Strassen e ibridazione

Il prodotto ricorsivo a blocchi (Prodotto di matrici - algoritmo iterativo e ricorsivo a blocchimodello di costo sulle matrici n x n (operazioni tra scalari, taglia n oppure N = n^2); somma in n^2 operazioni (ottima); prodotto righe per colonne iterativo con costo esatto 2n^3 - n^2; proprietà frattale del prodotto a blocchi dimostrata sugli indici; algoritmo ricorsivo RMUL con T(n) = 8T(n/2) + n^2, Theta(n^3) e soluzione esatta 2n^3 - n^2 per induzione e per srotolamento: fa le stesse operazioni dell'iterativo.Prodotto di matrici - algoritmo iterativo e ricorsivo a blocchi →) ha $T(n) = 8T(\frac n2) + n^2$: la funzione di soglia $n^3$ è molto più grande del lavoro $n^2$ dei nodi interni, siamo nel caso 1 del Master Theorem - enunciato, casi e applicazioniMaster Theorem per T(n) = aT(n/b) + w(n): si confronta w(n) con la funzione soglia n^(log_b a). Caso 1 (w minore di una potenza polinomiale): Theta(n^(log_b a)); caso 2 (w della stessa taglia): Theta(n^(log_b a) log n); caso 3 (w maggiore e regolare): Theta(w(n)). Prova completa dai livelli dell'albero, versione estesa (parti intere, n0 qualsiasi), esempi, e l'uso scorretto quando nessun caso si applica (w = n log n).Master Theorem - enunciato, casi e applicazioni → e dominano le foglie. L'arietà $8$ è troppo alta. Per abbassare la complessità bisogna ridurre le chiamate ricorsive senza far crescere troppo il lavoro di divide e conquer. Strassen ha mostrato che bastano 7 chiamate su matrici $\frac n2 \times \frac n2$, al prezzo di $18$ somme/sottrazioni di matrici $\frac n2 \times \frac n2$.

Si assume ancora $n = 2^k$ e il modello di costo del prodotto di matrici (operazioni tra scalari a costo $1$).

Le formule

Divide (14 fattori, $10$ somme/sottrazioni, perché $A_5, A_6, B_4, B_7$ sono quadranti singoli):

Recurse: $P_i = A_i \times B_i$ per $1 \le i \le 7$, con sette chiamate ricorsive.

Conquer ($8$ somme/sottrazioni):

$$C_{11} = P_1 + P_2 - P_4 + P_6, \quad C_{12} = P_4 + P_5, \quad C_{21} = P_6 + P_7, \quad C_{22} = P_2 - P_3 + P_5 - P_7.$$

Non è noto come Strassen sia arrivato alle espressioni, ma per la correttezza non serve: basta verificarle sviluppando i prodotti (le matrici non commutano, quindi l'ordine dei fattori in ogni termine si mantiene).

Verifica di $C_{12}$ e $C_{21}$:

$$P_4 + P_5 = (A_{11} + A_{12})B_{22} + A_{11}(B_{12} - B_{22}) = A_{11}B_{22} + A_{12}B_{22} + A_{11}B_{12} - A_{11}B_{22} = A_{11}B_{12} + A_{12}B_{22} = C_{12},$$

$$P_6 + P_7 = A_{22}(B_{21} - B_{11}) + (A_{21} + A_{22})B_{11} = A_{22}B_{21} - A_{22}B_{11} + A_{21}B_{11} + A_{22}B_{11} = A_{21}B_{11} + A_{22}B_{21} = C_{21}.$$

Verifica di $C_{11}$. Si sviluppano i quattro prodotti che compaiono:

$$\begin{aligned} P_1 &= A_{12}B_{21} + A_{12}B_{22} - A_{22}B_{21} - A_{22}B_{22},\\ P_2 &= A_{11}B_{11} + A_{11}B_{22} + A_{22}B_{11} + A_{22}B_{22},\\ -P_4 &= -A_{11}B_{22} - A_{12}B_{22},\\ P_6 &= A_{22}B_{21} - A_{22}B_{11}. \end{aligned}$$

Nella somma si cancellano a coppie: $A_{12}B_{22}$ con $-A_{12}B_{22}$; $-A_{22}B_{21}$ con $A_{22}B_{21}$; $-A_{22}B_{22}$ con $A_{22}B_{22}$; $A_{11}B_{22}$ con $-A_{11}B_{22}$; $A_{22}B_{11}$ con $-A_{22}B_{11}$. Restano $A_{11}B_{11} + A_{12}B_{21} = C_{11}$.

Verifica di $C_{22}$. $P_2 - P_3 + P_5 - P_7$ con

$$\begin{aligned} P_2 &= A_{11}B_{11} + A_{11}B_{22} + A_{22}B_{11} + A_{22}B_{22},\\ -P_3 &= -A_{11}B_{11} - A_{11}B_{12} + A_{21}B_{11} + A_{21}B_{12},\\ P_5 &= A_{11}B_{12} - A_{11}B_{22},\\ -P_7 &= -A_{21}B_{11} - A_{22}B_{11}. \end{aligned}$$

Si cancellano $A_{11}B_{11}$, $A_{11}B_{22}$, $A_{22}B_{11}$, $A_{11}B_{12}$, $A_{21}B_{11}$ e restano $A_{21}B_{12} + A_{22}B_{22} = C_{22}$.

Esempio numerico ($n = 2$, i quadranti sono scalari): $A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}$, $B = \begin{pmatrix} 5 & 6 \\ 7 & 8 \end{pmatrix}$.

• $P_1 = (2 - 4)(7 + 8) = -30$
• $P_2 = (1 + 4)(5 + 8) = 65$
• $P_3 = (1 - 3)(5 + 6) = -22$
• $P_4 = (1 + 2) \cdot 8 = 24$
• $P_5 = 1 \cdot (6 - 8) = -2$
• $P_6 = 4 \cdot (7 - 5) = 8$
• $P_7 = (3 + 4) \cdot 5 = 35$

$C_{11} = -30 + 65 - 24 + 8 = 19$, $C_{12} = 24 - 2 = 22$, $C_{21} = 8 + 35 = 43$, $C_{22} = 65 + 22 - 2 - 35 = 50$: lo stesso risultato $\begin{pmatrix} 19 & 22 \\ 43 & 50 \end{pmatrix}$ del prodotto per definizione. Qui Strassen usa $7$ prodotti e $18$ somme/sottrazioni, in tutto $25$ operazioni, contro le $12$ dell'iterativo: sui casi piccoli non conviene.

Algoritmo SMUL

SMUL(A, B)
  n <- rows(A)
  if n = 1 then return [a11 * b11]                   { BASE }
  * quadranti Aij, Bij *
  { DIVIDE: 10 somme/sottrazioni di matrici n/2 x n/2 }
  A1 <- A12 - A22;  B1 <- B21 + B22
  A2 <- A11 + A22;  B2 <- B11 + B22
  A3 <- A11 - A21;  B3 <- B11 + B12
  A4 <- A11 + A12;  B4 <- B22
  A5 <- A11;        B5 <- B12 - B22
  A6 <- A22;        B6 <- B21 - B11
  A7 <- A21 + A22;  B7 <- B11
  { RECURSE }
  for i <- 1 to 7 do Pi <- SMUL(Ai, Bi)
  { CONQUER: 8 somme/sottrazioni }
  C11 <- P1 + P2 - P4 + P6;   C12 <- P4 + P5
  C21 <- P6 + P7;             C22 <- P2 - P3 + P5 - P7
  return C

La correttezza segue dalle quattro identità verificate sopra, per induzione sulla taglia (la stessa struttura del prodotto ricorsivo; vedi 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 →).

Analisi

Costi per $n > 1$: il divide fa $10$ operazioni su matrici di $\frac n2 \times \frac n2$, cioè $10 \cdot \frac{n^2}{4} = \frac 52 n^2$; il conquer ne fa $8$, cioè $8 \cdot \frac{n^2}{4} = 2n^2$; la ricorsione costa $7\,T(\frac n2)$. Il lavoro dei nodi interni è $\frac 52 n^2 + 2n^2 = \frac 92 n^2$:

$$T(n) = \begin{cases} 1 & n = 1, \\ 7\,T\!\left(\frac n2\right) + \frac 92 n^2 & n > 1. \end{cases}$$

Ordine di grandezza. La soglia è $n^{\log_2 7} \approx n^{2.81}$, sopra $n^2$: siamo ancora nel caso 1 e $T(n) = \Theta(n^{\log_2 7})$, contro il $\Theta(n^3)$ dell'iterativo.

Soluzione esatta con la formula generale (Albero della ricorrenza e formula generale per le ricorrenzeL'albero della ricorrenza ha la forma dell'albero delle chiamate dell'istanza peggiore e T(n) è la somma dei lavori dei nodi. Per la famiglia T(n) = s(n)T(f(n)) + w(n) con base T0 per n <= n0, contando i nodi per livello con le iterate di f e il livello f*(n,n0) si ottiene la formula chiusa T(n) = somma_l [prodotto_j s(f^(j)(n))] w(f^(l)(n)) + [prodotto s] T0, verificata per induzione e applicata a esempi (merge sort, Hanoi, ricorrenze con radice).Albero della ricorrenza e formula generale per le ricorrenze →): $s(n) = 7$, $f(n) = \frac n2$, $w(n) = \frac 92 n^2$, $n_0 = 1$, $T_0 = 1$, profondità $f^*(n, n_0) = \log_2 n - 1 = k - 1$. Il livello $\ell$ ha $7^\ell$ nodi di taglia $n/2^\ell$; le foglie sono $7^k = n^{\log_2 7}$ (perché $7^{\log_2 n} = n^{\log_2 7}$):

$$T(n) = \sum_{\ell=0}^{k-1} 7^\ell \cdot \frac 92 \left(\frac{n}{2^\ell}\right)^2 + 7^k = \frac 92 n^2 \sum_{\ell=0}^{k-1} \left(\frac 74\right)^\ell + 7^k.$$

La somma geometrica vale $\dfrac{(7/4)^k - 1}{7/4 - 1} = \dfrac 43 \left( \left(\frac 74\right)^k - 1 \right)$, e $\left(\frac 74\right)^k = \dfrac{7^k}{4^k} = \dfrac{7^k}{n^2}$. Quindi

$$T(n) = \frac 92 n^2 \cdot \frac 43 \left( \frac{7^k}{n^2} - 1 \right) + 7^k = 6 \cdot 7^k - 6n^2 + 7^k = 7 n^{\log_2 7} - 6 n^2.$$

Controllo: $n = 1$ dà $7 - 6 = 1$; $n = 2$ dà $7 \cdot 7 - 6 \cdot 4 = 25$, il numero di operazioni dell'esempio numerico; $n = 4$ dà $7 \cdot 49 - 6 \cdot 16 = 247$.

Costanti a confronto. Il termine dominante di Strassen ha costante $7$, quello dell'iterativo costante $2$: più che tripla. Per quale $n$ Strassen conviene? Si cerca il primo $n$ con

$$7 n^{\log_2 7} - 6n^2 < 2n^3 - n^2.$$

Numericamente il passaggio avviene a $n \approx 655$. Per le sole potenze di $2$ (le taglie ammesse) la prima è $1024$: matrici con più di un milione di elementi.

E il modello di costo ignora la ricorsione, i pattern di accesso alla memoria (meno regolari, quindi peggiori per la cache) e le matrici temporanee per i fattori $A_i, B_i$ e i prodotti $P_i$, che l'iterativo non richiede. In pratica il punto di incrocio è molto più alto di quello teorico.

Evoluzioni. Strassen non è più l'algoritmo asintoticamente migliore: esponenti via via più piccoli (le dispense riportano $\Theta(n^{2.372})$, Le Gall 2014) sono stati ottenuti con tecniche molto più complesse, con costanti enormi e senza implementazioni usabili. Strassen, pur con costanti più alte dell'iterativo, è invece usato nella pratica per matrici di grande taglia, tipicamente in forma ibrida.

Ibridazione

Definizione (ibridazione). Dati due algoritmi $A_1$ e $A_2$ per lo stesso problema, con $T_{A_1}(n) = o(T_{A_2}(n))$ ma con costanti moltiplicative più alte per $A_1$ nell'espressione esatta, l'ibridazione li combina in un algoritmo che ha lo stesso andamento asintotico del migliore ma, se possibile, costanti più basse.

Ibrido ingenuo (STUPID_HYBRID): se $n$ è sotto una soglia (per esempio $1024$) si esegue MUL, altrimenti SMUL. Il costo è il minimo delle due funzioni, ma per $n \to \infty$ l'algoritmo coincide con SMUL, quindi la costante di $n^{\log_2 7}$ resta $7$: nessun miglioramento asintotico, e l'iterativo viene eseguito al più una volta.

Cambiare il caso base. Si fissa una soglia $n_0$ (potenza di $2$) e si usa l'iterativo a ogni foglia dell'albero di Strassen, che viene interrotto alla taglia $n_0$:

H_MUL(A, B)
  n <- rows(A)
  if n <= n0 then return MUL(A, B)
  { DIVIDE, RECURSE, CONQUER come in SMUL, con chiamate ricorsive a H_MUL }

(Le dispense scrivono n < n0 nel codice e n <= n0 nella ricorrenza: qui si usa n <= n0, coerente con la ricorrenza.) La correttezza è la composizione di due algoritmi corretti: l'impianto è quello ricorsivo, e il caso base è risolto da un algoritmo corretto. La ricorrenza diventa

$$T(n, n_0) = \begin{cases} 2n^3 - n^2 & n \le n_0, \\ 7\,T\!\left(\frac n2, n_0\right) + \frac 92 n^2 & n > n_0. \end{cases}$$

Analisi parametrica

Si pone $n = 2^k$, $n_0 = 2^{k_0}$ e $m = k - k_0 = \log_2 \frac{n}{n_0}$ livelli interni: $f^*(n, n_0) = m - 1$ ultimo livello interno, e $7^m$ foglie, ciascuna con costo $2n_0^3 - n_0^2$:

$$T(n, n_0) = \sum_{\ell=0}^{m-1} 7^\ell \cdot \frac 92 \left(\frac{n}{2^\ell}\right)^2 + 7^m (2n_0^3 - n_0^2) = \frac 92 n^2 \cdot \frac 43 \left( \left(\tfrac 74\right)^m - 1 \right) + 7^m (2n_0^3 - n_0^2).$$

Ora $\left(\frac 74\right)^m = \frac{7^m}{4^m} = \frac{7^m n_0^2}{n^2}$, perché $4^m = (n/n_0)^2$. Sostituendo:

$$T(n, n_0) = 6 \cdot 7^m n_0^2 - 6n^2 + 7^m (2n_0^3 - n_0^2) = 7^m (2n_0^3 + 5n_0^2) - 6n^2.$$

Infine $7^m = \dfrac{n^{\log_2 7}}{n_0^{\log_2 7}}$, da cui

$$T(n, n_0) = n^{\log_2 7} \cdot \frac{2n_0 + 5}{n_0^{\log_2 7 - 2}} - 6n^2.$$

Controllo: per $n_0 = 1$ si ottiene $7 n^{\log_2 7} - 6n^2$ (Strassen puro); per $n_0 = n$ si ottiene $n^2(2n + 5) - 6n^2 = 2n^3 - n^2$ (iterativo).

Per $n = 16$, $n_0 = 8$: $T = 7 \cdot (2 \cdot 512 - 64) + \frac 92 \cdot 256 = 7 \cdot 960 + 1152 = 7872$, contro $7936$ dell'iterativo e $15271$ di Strassen puro.

Scelta di $n_0$

La costante del termine dominante è $g(n_0) = \dfrac{2n_0 + 5}{n_0^{\log_2 7 - 2}}$, da minimizzare sulle potenze di $2$. Posto $s = \log_2 7 - 2 \approx 0.807$, $g(n_0) = (2n_0 + 5) n_0^{-s}$ ha derivata

$$g'(n_0) = 2 n_0^{-s} - s(2n_0 + 5) n_0^{-s-1} = n_0^{-s-1} \bigl( 2n_0 - s(2n_0 + 5) \bigr),$$

positiva se e solo se $2n_0 (1 - s) > 5s$, cioè

$$n_0 > \frac{5 s}{2(1 - s)} = \frac{5(\log_2 7 - 2)}{2(3 - \log_2 7)} \approx 10.48.$$

Quindi $g$ scende fino a circa $10.48$ e poi cresce: tra le potenze di $2$ il minimo è in $n_0 = 8$ o $n_0 = 16$.

Si sceglie $n_0 = 8$: la costante passa da $7$ a $3.92$, quasi dimezzata. Con $n_0 = 8$ l'ibrido supera l'iterativo già da $n = 16$ (non da $655$ come Strassen puro). Per $n = 1024$ i costi relativi all'iterativo sono: $n_0 = 1$, $0.918$; $n_0 = 2$, $0.674$; $n_0 = 4$, $0.556$; $n_0 = 8$, $0.513$; $n_0 = 16$, $0.516$; $n_0 = 32$, $0.550$. Come sempre, l'analisi teorica non conta ricorsione, accessi alla memoria e matrici temporanee: nelle implementazioni la soglia ottimale è più alta e si trova per via sperimentale.

Verifica con un programma

Un'implementazione di iterativo, ricorsivo, SMUL e H_MUL che conta le operazioni (somme/sottrazioni di matrici a $n^2$, prodotti iterativi a $2n^3 - n^2$) per $n = 1, 2, \dots, 64$ e ogni $n_0 \le n$ dà gli stessi risultati del prodotto di numpy, e conteggi uguali a $2n^3 - n^2$, a $7n^{\log_2 7} - 6n^2$ e alla formula parametrica, per tutti i valori provati.

Esercizi collegati

• Esercizio - matrici ARB, prodotto con tre chiamate ricorsive (compito del 19 giugno 2024, es. 1): stessa idea (meno chiamate, più somme) su una famiglia di matrici strutturate, con $3$ prodotti invece di $4$.
• Esercizio - algoritmo di Strassen sulle matrici supercircolanti (esercizi su supercircolanti e circolanti, es. SC 6): come si semplifica Strassen sulle supercircolanti e perché non basta.
• Esercizio - domande di prima parte su moltiplicazione di matrici (simulazioni d'esame): conteggi con Strassen e con l'ibrido.

Domande d'esame

1. Quante chiamate ricorsive fa l'algoritmo di Strassen e di quante somme/sottrazioni di matrici $\frac n2 \times \frac n2$ ha bisogno? Risposta: $7$ chiamate, $10$ somme/sottrazioni nel divide e $8$ nel conquer, in tutto $18$; $T(n) = 7T(\frac n2) + \frac 92 n^2$, $\Theta(n^{\log_2 7})$.
2. Strassen puro batte l'algoritmo iterativo? Da quale taglia? Traccia: asintoticamente sì ($n^{2.81}$ contro $n^3$), ma $7 n^{\log_2 7} - 6n^2 < 2n^3 - n^2$ solo da $n \approx 655$, quindi per le potenze di $2$ da $1024$.
3. Perché STUPID_HYBRID non migliora la costante e H_MUL sì? Traccia: il primo sceglie solo in base alla taglia di partenza, per $n$ grande coincide con SMUL, costante $7$; il secondo usa l'iterativo alle foglie di tutto l'albero, costante $g(n_0) = (2n_0 + 5)/n_0^{\log_2 7 - 2}$, minima ($\approx 3.92$) con $n_0 = 8$.