Probleme – La multiplication chaînée de matrices

1. Le coût dépend du parenthésage. Multiplier une matrice $a\times b$ par une $b\times c$ coûte $abc$ multiplications scalaires. Avec $d=(10,100,5,50)$ :

$$((A_1{A}_2)A_3):10\cdot 100\cdot 5+10\cdot 5\cdot 50=5000+2500=7500,$$

$$(A_1(A_2{A}_3)):100\cdot 5\cdot 50+10\cdot 100\cdot 50=25000+50000=75000.$$

Un facteur dix, pour trois matrices. Sur six matrices, l'écart se compte en ordres de grandeur.

2. Le sous-problème. $C[i][j]$ = coût minimal pour calculer $A_i\times \cdots \times A_j$. Le dernier produit effectué coupe la chaîne en un point $k$ :

$$C[i][i]=0,C[i][j]=\underset{i\le k<j}{\min }(C[i][k]+C[k+1][j]+d_{i-1}{d}_k{d}_j).$$

La sous-structure optimale se voit ici très bien : si le parenthésage optimal coupe en $k$, alors les deux moitiés doivent elles-mêmes être parenthésées de façon optimale — sinon on améliorerait le tout en les remplaçant, sans toucher au reste.

L'ordre de remplissage n'est ni « par ligne » ni « par colonne » : $C[i][j]$ dépend de cases dont l'intervalle est plus court. On remplit donc par longueur d'intervalle croissante : d'abord tous les $C[i][i]$, puis tous les intervalles de longueur $2$, etc. C'est la première fois du chapitre où l'ordre de remplissage n'est pas l'ordre naturel des indices, et c'est là que se trompent les mises en œuvre.

3. Le programme.

/* Cout minimal du produit A_1 x ... x A_n, ou A_i est d[i-1] x d[i].
   Remplit aussi coupe[i][j] : le point de coupure optimal.
   Preconditions : n >= 1, d de taille n+1, d[i] >= 1.
   Complexite : Theta(n^3) en temps, Theta(n^2) en memoire. */
int chaine(const int d[], int n, int cout[][N_MAX + 1], int coupe[][N_MAX + 1]) {
    assert(n >= 1);
    for (int i = 1; i <= n; i = i + 1) { cout[i][i] = 0; }
    /* INVARIANT : au debut du tour l, tous les cout[i][j] avec j-i+1 < l
       sont definitifs. */
    for (int l = 2; l <= n; l = l + 1) {
        for (int i = 1; i + l - 1 <= n; i = i + 1) {
            int j = i + l - 1;
            cout[i][j] = INT_MAX;
            for (int k = i; k < j; k = k + 1) {
                int q = cout[i][k] + cout[k + 1][j] + d[i - 1] * d[k] * d[j];
                if (q < cout[i][j]) { cout[i][j] = q; coupe[i][j] = k; }
            }
        }
    }
    return cout[1][n];
}

/* Affiche le parenthesage optimal de A_i x ... x A_j.
   Precondition : coupe a ete rempli par chaine(). */
void afficher(int coupe[][N_MAX + 1], int i, int j) {
    if (i == j) { printf("A%d", i); return; }
    printf("(");
    afficher(coupe, i, coupe[i][j]);
    afficher(coupe, coupe[i][j] + 1, j);
    printf(")");
}

Mesure sur l'instance classique $d=(30,35,15,5,10,20,25)$, six matrices :

ligne i=1 : 0  15750  7875  9375  11875  15125
ligne i=2 : -      0  2625  4375   7125  10500
ligne i=3 : -      -     0   750   2500   5375
ligne i=4 : -      -     -     0   1000   3500
ligne i=5 : -      -     -     -      0   5000
ligne i=6 : -      -     -     -      -      0
optimum = 15125   parenthesage = ((A1(A2A3))((A4A5)A6))

Une énumération exhaustive des parenthésages confirme $15125$ et le même arbre.

4. Les complexités. Il y a $\Theta (n^2)$ sous-problèmes, chacun coûtant $\Theta (n)$ à résoudre (la boucle sur $k$) : $\Theta (n^3)$ en temps, $\Theta (n^2)$ en mémoire. On ne peut pas réduire la mémoire ici, puisque la reconstruction lit la table coupe.

Le nombre de parenthésages de $n$ matrices est le $(n-1)$-ième nombre de Catalan,

$$C_{n-1}=\frac{1}{n}\left(\genfrac{}{}{0ex}{}{2n-2}{n-1}\right)\sim \frac{4^{n-1}}{(n-1{)}^{3/2}\sqrt{\pi }},$$

soit $C_9=4862$ pour $n=10$ et $C_{29}=1002242216651368$, environ $10^{15}$, pour $n=30$. L'énumération est hors de portée dès $n=25$, la table est instantanée jusqu'à $n=500$ : c'est le gain de la programmation dynamique, chiffré.