Prodotto di matrici - algoritmo iterativo e ricorsivo a blocchi

Si lavora con matrici quadrate $n \times n$ a elementi complessi. Il caso rettangolare si tratta allo stesso modo; il caso quadrato non è una vera limitazione. Per le dimostrazioni sulla struttura a blocchi si assume $n = 2^k$ con $k \ge 0$, così ogni matrice si divide in quattro quadranti $\frac n2 \times \frac n2$ fino alle matrici $1 \times 1$.

Modello di costo. Costo unitario per ogni operazione aritmetica tra due scalari (somma, sottrazione, prodotto, divisione). Gli assegnamenti, le copie e il calcolo degli indici costano $0$. La complessità si esprime in funzione di $n$ (dimensione della matrice) oppure di $N = n^2$ (numero di elementi): sono la stessa quantità descritta in due modi, quindi va sempre dichiarata la metrica di taglia. Un algoritmo con $T(n) = \Theta(n^2)$ è lineare nella taglia $N$.

Notazione: le matrici sono lettere maiuscole, i loro elementi la lettera minuscola con due indici, $A = [a_{ij}]_{1 \le i,j \le n}$. Le precondizioni sull'input (qui: le matrici hanno la stessa dimensione) si chiamano vincoli statici. Gli algoritmi "basati sulla definizione" calcolano ogni elemento del risultato con la formula che lo definisce. Per la definizione di prodotto righe per colonne vedi 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 →.

Somma e sottrazione

$C = A + B$ è definita da $c_{ij} = a_{ij} + b_{ij}$.

SUM(A, B)
  n <- rows(A)
  for i <- 1 to n do
    for j <- 1 to n do
      c[i][j] <- a[i][j] + b[i][j]
  return C

Ogni for si traduce in una sommatoria con la stessa variabile, il cui termine generico è il lavoro di una iterazione. for annidati diventano sommatorie annidate:

$$T_{SUM}(n) = \sum_{i=1}^{n} \sum_{j=1}^{n} 1 = \sum_{i=1}^{n} n = n^2 = \Theta(n^2).$$

Qualunque algoritmo di somma deve eseguire, al caso peggiore, $n^2$ somme distinte (ogni $c_{ij}$ dipende da una coppia diversa di input, e non c'è modo di ottenerlo da altre somme): quindi $n^2$ non è migliorabile e l'algoritmo basato sulla definizione è ottimo. Lo stesso vale per la sottrazione SUB.

Prodotto righe per colonne, versione iterativa

$C = A \times B$ è definita da

$$c_{ij} = \sum_{k=1}^{n} a_{ik} b_{kj}, \qquad 1 \le i, j \le n.$$

MUL(A, B)
  n <- rows(A)
  for i <- 1 to n do
    for j <- 1 to n do
      c[i][j] <- a[i][1] * b[1][j]
      for k <- 2 to n do
        c[i][j] <- c[i][j] + a[i][k] * b[k][j]
  return C

Si inizializza con il primo prodotto (e non con $0$) perché sommare a $0$ sarebbe un'operazione in più. Conto le operazioni sostituendo i for con sommatorie. Per ogni coppia $(i,j)$ il lavoro è: $1$ prodotto iniziale, poi per ciascun $k = 2, \dots, n$ un prodotto e una somma, cioè $2$ operazioni:

$$T(n) = \sum_{i=1}^{n} \Bigl( \sum_{j=1}^{n} \bigl( 1 + \sum_{k=2}^{n} 2 \bigr) \Bigr) = \sum_{i=1}^{n} \sum_{j=1}^{n} \bigl(1 + 2(n-1)\bigr) = n^2 (2n - 1) = 2n^3 - n^2.$$

La sommatoria interna ha $n - 1$ termini (da $k=2$ a $k=n$), da cui il fattore $2(n-1)$. In totale i prodotti sono $n^3$ e le somme $n^3 - n^2$ (un prodotto per ogni terna $(i,j,k)$, e per ogni coppia $(i,j)$ una somma in meno dei prodotti). Quindi $T_{MUL}(n) = \Theta(n^3) = \Theta(N \sqrt N)$: non è lineare nel numero di elementi.

Esempio ($n = 2$).

$$\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} \begin{pmatrix} 5 & 6 \\ 7 & 8 \end{pmatrix} = \begin{pmatrix} 1\cdot 5 + 2 \cdot 7 & 1 \cdot 6 + 2 \cdot 8 \\ 3 \cdot 5 + 4 \cdot 7 & 3 \cdot 6 + 4 \cdot 8 \end{pmatrix} = \begin{pmatrix} 19 & 22 \\ 43 & 50 \end{pmatrix}.$$

Le operazioni sono $8$ prodotti e $4$ somme, cioè $12 = 2 \cdot 2^3 - 2^2$, come prevede la formula.

Prodotto ricorsivo a blocchi

Si cerca un algoritmo divide & conquer (vedi Divide and conquer - proprietà di sottostruttura e albero delle chiamateUn algoritmo divide and conquer nasce da una proprietà di sottostruttura (funzioni Divide e Conquer, base n0) con sottoistanze di taglia strettamente minore; l'albero delle chiamate ne descrive l'esecuzione, e il D&C non ha memoria (ripete le sottoistanze uguali). Esempio completo: MAX di un array, con correttezza per induzione, ricorrenza e conteggio esatto n-1 tramite il teorema sugli alberi binari pieni.Divide and conquer - proprietà di sottostruttura e albero delle chiamate →). Per $n = 2^k$ si scrivono $A$ e $B$ a quadranti:

$$A = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix}, \qquad B = \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix}, \qquad A_{hk}, B_{hk} \in \mathbb C^{\frac n2 \times \frac n2}.$$

Proprietà di sottostruttura (proprietà frattale)

Teorema (proprietà frattale del prodotto). I quadranti di $C = A \times B$ si ottengono con la stessa formula del prodotto righe per colonne di matrici $2 \times 2$, dove al posto degli scalari ci sono i quadranti: per $1 \le i, j \le 2$

$$C_{ij} = \sum_{k=1}^{2} A_{ik} B_{kj} = A_{i1} B_{1j} + A_{i2} B_{2j}.$$

Dimostrazione per $C_{11}$ (gli altri tre quadranti sono identici). Con $(M)_{ij}$ si indica l'elemento di posto $(i,j)$ della matrice $M$; per $1 \le i, j \le \frac n2$ vale $(C_{11})_{ij} = c_{ij}$ perché $C_{11}$ occupa le prime $\frac n2$ righe e colonne.

1. Si parte dalla definizione e si spezza la somma su $k$ in due metà, la prima per $k \le \frac n2$ e la seconda per $k > \frac n2$: $$(C_{11})_{ij} = \sum_{k=1}^{n} a_{ik} b_{kj} = \sum_{k=1}^{n/2} a_{ik} b_{kj} + \sum_{k=n/2+1}^{n} a_{ik} b_{kj}.$$
2. Nella seconda somma si pone $k = h + \frac n2$ con $1 \le h \le \frac n2$ (traslazione dell'indice, non cambia i termini): $$= \sum_{k=1}^{n/2} a_{ik} b_{kj} + \sum_{h=1}^{n/2} a_{i,h+n/2}\, b_{h+n/2,\,j}.$$
3. Si riconoscono i quadranti. Nella prima somma $i, k \le \frac n2$ e $k, j \le \frac n2$: $a_{ik} = (A_{11})_{ik}$ e $b_{kj} = (B_{11})_{kj}$. Nella seconda $i \le \frac n2$ e $h + \frac n2 > \frac n2$: $a_{i,h+n/2} = (A_{12})_{ih}$ e $b_{h+n/2,\,j} = (B_{21})_{hj}$ (riga $h$ del blocco in basso, colonna $j$ del blocco a sinistra): $$= \sum_{k=1}^{n/2} (A_{11})_{ik} (B_{11})_{kj} + \sum_{h=1}^{n/2} (A_{12})_{ih} (B_{21})_{hj} = (A_{11} B_{11})_{ij} + (A_{12} B_{21})_{ij}.$$

Quindi $C_{11} = A_{11}B_{11} + A_{12}B_{21}$. $\blacksquare$

La proprietà dice che un prodotto $n \times n$ si riduce a 8 prodotti $\frac n2 \times \frac n2$ e a 4 somme di matrici $\frac n2 \times \frac n2$.

Algoritmo RMUL

RMUL(A, B)
  n <- rows(A)
  if n = 1 then return [a11 * b11]       { matrice 1x1 }
  * A = [A11 A12; A21 A22], B = [B11 B12; B21 B22] *
  for i <- 1 to 2 do
    for j <- 1 to 2 do
      Cij <- SUM(RMUL(Ai1, B1j), RMUL(Ai2, B2j))
  return C

La riga fra asterischi è un metacomando: non è un'istruzione. In un linguaggio reale la divisione in quadranti si realizza passando alla funzione ricorsiva i parametri che indicano la sottomatrice (riga e colonna di partenza e dimensione), senza copiare niente. Costi: nel caso base c'è un prodotto; nel caso ricorsivo, 8 chiamate su taglia $\frac n2$ più 4 somme di matrici $\frac n2 \times \frac n2$, ciascuna di costo $(\frac n2)^2$ per $T_{SUM}$, cioè $4 \cdot \frac{n^2}{4} = n^2$:

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

Analisi con il Master Theorem

Qui $a = 8$, $b = 2$, $w(n) = n^2$ e la funzione di soglia è $n^{\log_b a} = n^{\log_2 8} = n^3$. Poiché $w(n) = n^2 = O(n^{3 - \varepsilon})$ con $\varepsilon = 1$, 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 →: $T(n) = \Theta(n^3)$. A dominare è il lavoro alle foglie ($n^3$ foglie, una per terna $(i,j,k)$): l'arietà $8$ dell'albero fa crescere troppo in fretta il numero di foglie, mentre il lavoro nei nodi interni è piccolo. Non c'è nessun miglioramento rispetto all'iterativo.

Soluzione esatta: stesso conteggio dell'iterativo

Proposizione. $T(n) = 2n^3 - n^2$ per ogni $n = 2^k$: RMUL esegue esattamente le stesse operazioni di MUL, solo in ordine diverso.

Dimostrazione per induzione su $k$ (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 →).

• Base ($n = 1$): $T(1) = 1$ e $2 \cdot 1^3 - 1^2 = 1$.
• Ipotesi: $T\bigl(\frac n2\bigr) = 2\bigl(\frac n2\bigr)^3 - \bigl(\frac n2\bigr)^2$.
• Passo: $$T(n) = 8\,T\!\left(\tfrac n2\right) + n^2 = 8\left( 2\tfrac{n^3}{8} - \tfrac{n^2}{4} \right) + n^2 = 2n^3 - 2n^2 + n^2 = 2n^3 - n^2. \qquad \blacksquare$$

La stessa formula si ottiene srotolando l'albero (formula generale: vedi 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 →). Al livello $\ell$ ci sono $8^\ell$ nodi di taglia $n/2^\ell$, ciascuno con lavoro $\bigl(n/2^\ell\bigr)^2$; i nodi interni sono i livelli $\ell = 0, \dots, k-1$ e le $8^k = n^3$ foglie valgono $1$ ciascuna:

$$T(n) = \sum_{\ell=0}^{k-1} 8^\ell \frac{n^2}{4^\ell} + 8^k = n^2 \sum_{\ell=0}^{k-1} 2^\ell + n^3 = n^2 (2^k - 1) + n^3 = n^3 - n^2 + n^3 = 2n^3 - n^2.$$

La somma geometrica $\sum_{\ell=0}^{k-1} 2^\ell = 2^k - 1 = n - 1$.

Esempio ($n = 4$). Il livello $0$ ha una chiamata con $4$ somme $2 \times 2$ (cioè $4 \cdot 4 = 16$ operazioni); il livello $1$ ha $8$ chiamate di taglia $2$, ognuna con $4$ somme di matrici $1 \times 1$ ($4$ operazioni), cioè $32$; il livello $2$ ha $64$ foglie, $64$ prodotti. Totale $16 + 32 + 64 = 112 = 2 \cdot 4^3 - 4^2$.

Verifica con un programma

def mul(A, B):                      # iterativo, conta le operazioni
    n, ops = len(A), 0
    C = [[0]*n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            c = A[i][0]*B[0][j]; ops += 1
            for k in range(1, n):
                c += A[i][k]*B[k][j]; ops += 2
            C[i][j] = c
    return C, ops

# n = 1, 2, 4, ..., 64: ops == 2*n**3 - n**2

Per $n = 8$ il conteggio dà $2 \cdot 512 - 64 = 960$ operazioni, e la versione ricorsiva dà lo stesso valore.

Il Master Theorem non basta a battere il cubo: serve ridurre il numero di chiamate ricorsive da $8$ a meno, accettando di aumentare il lavoro dei nodi interni. È l'idea di 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 →.

Esercizi collegati

• Esercizio - prodotto di matrici complesse, ricorrenza e induzione parametrica (compito del 8 luglio 2024, es. 1): conteggio esatto in operazioni reali di questo stesso algoritmo.
• Esercizio - domande di prima parte su moltiplicazione di matrici (simulazioni d'esame): la domanda sul numero di operazioni per $n = 8$.

Domande d'esame

1. Quante operazioni tra scalari esegue il prodotto righe per colonne ricorsivo basato sulla definizione per $n = 8$? Risposta: $2n^3 - n^2 = 2 \cdot 512 - 64 = 960$ ($n^3 = 512$ prodotti e $n^3 - n^2 = 448$ somme, perché ogni prodotto si somma nel $c_{ij}$ tranne i $64$ primi). Traccia: l'algoritmo ricorsivo fa le stesse operazioni dell'iterativo.
2. Perché il prodotto ricorsivo a blocchi con $8$ chiamate non migliora l'iterativo? Traccia: $a = 8$, soglia $n^3$, lavoro $n^2$ polinomialmente minore, caso 1 del Master Theorem: dominano le foglie, $\Theta(n^3)$; per di più il conteggio esatto è identico, $2n^3 - n^2$.
3. Perché la somma di matrici è ottima? Traccia: ogni $c_{ij}$ richiede la propria somma, quindi servono $n^2$ operazioni e $\Theta(n^2)$ è anche un limite inferiore.