1A - Enoncé 3 novembre 2021¶

Correction de l’examen du 3 novembre 2021.

Exercice 1 : multiplication de matrices¶

On a besoin d’une fonction qui mesure le temps d’exécution d’une fonction.

[2]:

import time


def mesure_temps_fonction(fct, N=100):
    begin = time.perf_counter()
    for i in range(N):
        fct()
    return (time.perf_counter() - begin) / N


mesure_temps_fonction(lambda: time.sleep(0.1), N=10)

[2]:

0.10280692000000005

Q1 : Pourquoi (m1 @ m2) @ m3 est-il plus lent que m1 @ (m2 @ m3) ? (2 points)¶

Il y a deux options possible. Il suffit de compter le nombre d’opérations dans chaque option. Le coût d’une multiplication $M_{ab} \times m_{bc}$ est de l’ordre de $O(abc)$ . Donc :

cout((m1 @ m2) @ m3) ~ O(997 * 93 * 1003 + 997 * 1003 * 97) = 189998290
cout(m1 @ (m2 @ m3)) ~ O(93 * 1003 * 97 + 997 * 93 * 97) = 18042000

La seconde option est dix fois plus rapide.

[3]:

print(997 * 93 * 1003 + 997 * 1003 * 97, 93 * 1003 * 97 + 997 * 93 * 97)

189998290 18042000

[4]:

import numpy

m1 = numpy.random.randn(997, 93)
m2 = numpy.random.randn(93, 1003)
m3 = numpy.random.randn(1003, 97)

mesure_temps_fonction(lambda: m1 @ m2 @ m3)

[4]:

0.009560690999999987

[5]:

mesure_temps_fonction(lambda: (m1 @ m2) @ m3)

[5]:

0.009846866999999992

[6]:

mesure_temps_fonction(lambda: m1 @ (m2 @ m3))

[6]:

0.001004321000000008

Q2 : Ecrire une fonction qui calcule le nombre d’operations dans une multiplication de deux matrices (2 points)¶

[7]:

def n_ops(m1_shape, m2_shape):
    return m1_shape[0] * m2_shape[1] * m1_shape[1] * 2


n_ops(m1.shape, m2.shape)

[7]:

185998326

Q3 : Ecrire une fonction qui retourne le meilleur coût d’une multiplication de deux matrices et la meilleure option (2 points)¶

[8]:

def n_ops_3(sh1, sh2, sh3):
    m1_m2m3 = n_ops(sh1, (sh2[0], sh3[1])) + n_ops(sh2, sh3)
    m1m2_m3 = n_ops(sh1, sh2) + n_ops((sh1[0], sh2[1]), sh3)
    if m1m2_m3 < m1_m2m3:
        return m1m2_m3, 2
    else:
        return m1_m2m3, 1


n_ops_3(m1.shape, m2.shape, m3.shape)

[8]:

(36084000, 1)

Q4 : Ecrire une fonction qui effectue le produit de trois matrices le plus rapidement possible (2 points)¶

[9]:

from numpy.testing import assert_almost_equal


def produit3(m1, m2, m3):
    cout, meilleur = n_ops_3(m1.shape, m2.shape, m3.shape)
    if meilleur == 2:
        return (m1 @ m2) @ m3
    else:
        return m1 @ (m2 @ m3)


assert_almost_equal(produit3(m1, m2, m3), m1 @ (m2 @ m3))

Q5 : Vérifiez que vous retrouvez les mêmes résultats avec la fonction `mesure_temps` (2 points)¶

[10]:

mesure_temps_fonction(lambda: produit3(m1, m2, m3))

[10]:

0.0011657070000000048

On vérifie que c’est égal à :

[11]:

mesure_temps_fonction(lambda: m1 @ (m2 @ m3))

[11]:

0.0011010209999999887

Ici, vous avez le choix entre faire les questions 6 à 9 ou les questions 9 et 10.

Q6 : Ecrire une fonction qui retourne le meilleur coût d’une multiplication de 4 matrices et la meilleure option (3 points)¶

[12]:

m4 = numpy.random.randn(97, 20)

[13]:

def n_ops_4(sh1, sh2, sh3, sh4):
    m1_m2m3m4 = n_ops(sh1, (sh2[0], sh4[1])) + n_ops_3(sh2, sh3, sh4)[0]
    m1m2_m3m4 = (
        n_ops(sh1, sh2) + n_ops((sh1[0], sh2[1]), (sh3[0], sh4[1])) + n_ops(sh3, sh4)
    )
    m1m2m3_m4 = n_ops_3(sh1, sh2, sh3)[0] + n_ops((sh1[0], sh3[1]), sh4)
    m = min(m1_m2m3m4, m1m2_m3m4, m1m2m3_m4)
    if m == m1_m2m3m4:
        return m, 1
    if m == m1m2_m3m4:
        return m, 2
    return m, 3


n_ops_4(m1.shape, m2.shape, m3.shape, m4.shape)

[13]:

(11331640, 1)

Q7 : Ecrire une fonction qui effectue le produit de 4 matrices le plus rapidement possible (3 points)¶

[14]:

def produit4(m1, m2, m3, m4):
    cout, meilleur = n_ops_4(m1.shape, m2.shape, m3.shape, m4.shape)
    if meilleur == 1:
        return m1 @ produit3(m2, m3, m4)
    if meilleur == 2:
        return (m1 @ m2) @ (m3 @ m4)
    return produit3(m1, m2, m3) @ m4


mesure_temps_fonction(lambda: produit4(m1, m2, m3, m4))

[14]:

0.000988687999999982

Q8 : Vérifiez que vous retrouvez les mêmes résultats avec la fonction mesure_temps et la matrice m4. (2 points)¶

[15]:

mesure_temps_fonction(lambda: ((m1 @ m2) @ m3) @ m4)

[15]:

0.010446371000000027

[16]:

mesure_temps_fonction(lambda: (m1 @ m2) @ (m3 @ m4))

[16]:

0.008082993999999956

[17]:

mesure_temps_fonction(lambda: m1 @ (m2 @ (m3 @ m4)))

[17]:

0.0008713240000000155

[18]:

mesure_temps_fonction(lambda: produit4(m1, m2, m3, m4))

[18]:

0.0009054390000000013

Q9 : On se penche sur le cas à une multiplication de N matrices, combien y a-t-il de multiplications de 2 matrices ? (2 points)¶

Il y a en toujours N-1. On considère le produit $M_1 \times... \times M_n$ . La multiplication commence toujours par une multiplication de deux matrices consécutives quelles qu’elles soient. On les suppose aux positions $(i, i+1)$ . On note le résultat $MM_i$ . Après ce produit, il faudra faire : $(M_1 \times ... \times M_{i-1} \times MM_i \times M_{i+2} \times ... \times M_n$ , soit une multiplication de $N-2$ matrices. On obtient le résultat par récurrence.

[19]:

Ici s’arrête l’énoncé pour ceux qui ont choisit de répondre aux question 6 à 9.

Q10 : Résoudre l’optimisation de multiplication de N matrices.¶

On l’envisage de façon récursive. La première solution effectue plein de calculs en double mais nous verront comment la modifier.

[20]:

def n_ops_N(shapes):
    if len(shapes) <= 1:
        raise RuntimeError("Unexpected list of shapes: %r." % shapes)
    if len(shapes) == 2:
        return n_ops(*shapes), 1
    if len(shapes) == 3:
        return n_ops_3(*shapes)
    best_cost = None
    best_pos = None
    for i in range(1, len(shapes)):
        if i == 1:
            cost = (
                n_ops(shapes[0], (shapes[1][0], shapes[-1][1])) + n_ops_N(shapes[1:])[0]
            )
            best_cost = cost
            best_pos = i
        elif i == len(shapes) - 1:
            cost = n_ops_N(shapes[:-1])[0] + n_ops(
                (shapes[0][0], shapes[-2][1]), shapes[-1]
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i
        else:
            cost = (
                n_ops_N(shapes[:i])[0]
                + n_ops_N(shapes[i:])[0]
                + n_ops((shapes[0][0], shapes[i - 1][1]), (shapes[i][0], shapes[-1][1]))
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i

    if best_pos is None:
        raise RuntimeError(shapes)
    return best_cost, best_pos


n_ops_N([m1.shape, m2.shape, m3.shape, m4.shape])

[20]:

(11331640, 1)

[21]:

n_ops_4(m1.shape, m2.shape, m3.shape, m4.shape)

[21]:

(11331640, 1)

[22]:

def product_N(inputs):
    if len(inputs) <= 1:
        raise RuntimeError(
            "List inputs must contain at least two elements bot has %d." % len(inputs)
        )
    cost, pos = n_ops_N([i.shape for i in inputs])
    if len(inputs) == 2:
        return inputs[0] @ inputs[1]
    if pos == 1:
        right = product_N(inputs[1:])
        return inputs[0] @ right
    if pos == len(shapes) - 1:
        left = product_N(inputs[:-1])
        return left @ inputs[-1]
    else:
        left = product_N(inputs[: pos + 1])
        right = product_N(inputs[pos + 1 :])
        return left @ right


assert_almost_equal(m1 @ m2 @ m3 @ m4, product_N([m1, m2, m3, m4]))

[23]:

mesure_temps_fonction(lambda: produit4(m1, m2, m3, m4))

[23]:

0.0009734980000000349

[24]:

mesure_temps_fonction(lambda: product_N([m1, m2, m3, m4]))

[24]:

0.0009873739999999031

Ici s’arrête ce qui est attendu comme réponse à la question 10.

Les calculs en double…

On vérifie en ajoutant une ligne pour afficher tous les appels à n_ops_N.

[25]:

def n_ops_N(shapes, verbose=False):
    if verbose:
        print("n_ops_N(%r)" % shapes)
    if len(shapes) <= 1:
        raise RuntimeError("Unexpected list of shapes: %r." % shapes)
    if len(shapes) == 2:
        return n_ops(*shapes), 1
    if len(shapes) == 3:
        return n_ops_3(*shapes)
    best_cost = None
    best_pos = None
    for i in range(1, len(shapes)):
        if i == 1:
            cost = (
                n_ops(shapes[0], (shapes[1][0], shapes[-1][1]))
                + n_ops_N(shapes[1:], verbose=verbose)[0]
            )
            best_cost = cost
            best_pos = i
        elif i == len(shapes) - 1:
            cost = n_ops_N(shapes[:-1], verbose=verbose)[0] + n_ops(
                (shapes[0][0], shapes[-2][1]), shapes[-1]
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i
        else:
            cost = (
                n_ops_N(shapes[:i], verbose=verbose)[0]
                + n_ops_N(shapes[i:], verbose=verbose)[0]
                + n_ops((shapes[0][0], shapes[i - 1][1]), (shapes[i][0], shapes[-1][1]))
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i

    if best_pos is None:
        raise RuntimeError(shapes)
    return best_cost, best_pos


m5 = numpy.random.randn(20, 17)

n_ops_N([m1.shape, m2.shape, m3.shape, m4.shape, m5.shape], verbose=True)

n_ops_N([(997, 93), (93, 1003), (1003, 97), (97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97), (97, 20), (20, 17)])
n_ops_N([(1003, 97), (97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97)])
n_ops_N([(97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97), (97, 20)])
n_ops_N([(997, 93), (93, 1003)])
n_ops_N([(1003, 97), (97, 20), (20, 17)])
n_ops_N([(997, 93), (93, 1003), (1003, 97)])
n_ops_N([(97, 20), (20, 17)])
n_ops_N([(997, 93), (93, 1003), (1003, 97), (97, 20)])
n_ops_N([(93, 1003), (1003, 97), (97, 20)])
n_ops_N([(997, 93), (93, 1003)])
n_ops_N([(1003, 97), (97, 20)])
n_ops_N([(997, 93), (93, 1003), (1003, 97)])

[25]:

(9697854, 1)

On voit deux appels identiques n_ops_N([(97, 20), (20, 17)]) et n_ops_N([(93, 1003), (1003, 97), (97, 20)]). Ce n’est pas trop problématique pour un petit nombre de matrices mais cela pourrait le devenir si ce même algorithme était appliquée à autre chose.

Plutôt que de réécrire l’algorithme différemment, on se propose d’ajouter un paramètre pour garder la trace des résultats déjà retournés.

[26]:

def n_ops_N_opt(shapes, cache=None, verbose=False):
    if cache is None:
        cache = {}
    key = tuple(shapes)
    if key in cache:
        # On s'arrête, déjà calculé.
        return cache[key]

    if verbose:
        print("n_ops_N(%r)" % shapes)
    if len(shapes) <= 1:
        raise RuntimeError("Unexpected list of shapes: %r." % shapes)

    if len(shapes) == 2:
        res = n_ops(*shapes), 1
        cache[key] = res
        return res

    if len(shapes) == 3:
        res = n_ops_3(*shapes)
        cache[key] = res
        return res

    best_cost = None
    best_pos = None
    for i in range(1, len(shapes)):
        if i == 1:
            cost = (
                n_ops(shapes[0], (shapes[1][0], shapes[-1][1]))
                + n_ops_N_opt(shapes[1:], verbose=verbose, cache=cache)[0]
            )
            best_cost = cost
            best_pos = i
        elif i == len(shapes) - 1:
            cost = n_ops_N_opt(shapes[:-1], verbose=verbose, cache=cache)[0] + n_ops(
                (shapes[0][0], shapes[-2][1]), shapes[-1]
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i
        else:
            cost = (
                n_ops_N_opt(shapes[:i], verbose=verbose, cache=cache)[0]
                + n_ops_N_opt(shapes[i:], verbose=verbose, cache=cache)[0]
                + n_ops((shapes[0][0], shapes[i - 1][1]), (shapes[i][0], shapes[-1][1]))
            )
            if cost < best_cost:
                best_cost = cost
                best_pos = i

    if best_pos is None:
        raise RuntimeError(shapes)

    res = best_cost, best_pos
    cache[key] = res
    return res


n_ops_N_opt([m1.shape, m2.shape, m3.shape, m4.shape, m5.shape], verbose=True)

n_ops_N([(997, 93), (93, 1003), (1003, 97), (97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97), (97, 20), (20, 17)])
n_ops_N([(1003, 97), (97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97)])
n_ops_N([(97, 20), (20, 17)])
n_ops_N([(93, 1003), (1003, 97), (97, 20)])
n_ops_N([(997, 93), (93, 1003)])
n_ops_N([(997, 93), (93, 1003), (1003, 97)])
n_ops_N([(997, 93), (93, 1003), (1003, 97), (97, 20)])
n_ops_N([(1003, 97), (97, 20)])

[26]:

(9697854, 1)

La liste est moins longue et tous les appels sont uniques. On met à jour la fonction product_N.

[27]:

def product_N_opt(inputs, cache=None):
    if len(inputs) <= 1:
        raise RuntimeError(
            "List inputs must contain at least two elements bot has %d." % len(inputs)
        )
    cost, pos = n_ops_N_opt([i.shape for i in inputs], cache=cache)
    if len(inputs) == 2:
        return inputs[0] @ inputs[1]
    if pos == 1:
        right = product_N_opt(inputs[1:], cache=cache)
        return inputs[0] @ right
    if pos == len(shapes) - 1:
        left = product_N_opt(inputs[:-1], cache=cache)
        return left @ inputs[-1]
    else:
        left = product_N_opt(inputs[: pos + 1], cache=cache)
        right = product_N_opt(inputs[pos + 1 :], cache=cache)
        return left @ right


assert_almost_equal(m1 @ m2 @ m3 @ m4, product_N([m1, m2, m3, m4]))

[28]:

mesure_temps_fonction(lambda: product_N([m1, m2, m3, m4, m5]))

[28]:

0.0010903469999999516

[29]:

mesure_temps_fonction(lambda: product_N_opt([m1, m2, m3, m4, m5]))

[29]:

0.0009383259999999893

[30]:

mesure_temps_fonction(lambda: m1 @ m2 @ m3 @ m4 @ m5)

[30]:

0.01018160299999991

Tout fonctionne.

Remarques lors de la correction¶

Il y a eu peu d’erreurs lors des premières questions. Par la suite des erreurs fréquentes sont apparues.

Il ne fallait pas utiliser de produits matriciel dans les fonctions de coûts. L’intérêt est d’utiliser ces fonctions pour décider du calcul à faire, pour déterminer le calcul optimal. Et le calcu de ce coût doit être négligeable par rapport au coût matriciel lui-même sinon l’intérêt en est fortement réduit.

Le produit de 4 matrices ne pouvait pas faire intervenir m1 @ m2 @ m3 car cette notation ne précise pas explicitement l’ordre à suivre.

Enfin, les mesures de temps étaient destinées à repérer les erreurs de code éventuelles. Si la mesure donne l’inverse ce qui est attendu, c’est qu’il y a sans doute une erreur de code. De même, si la mesure de temps dure très longtemps, c’est aussi une indication que le code est probablement erroné.

[31]:

Notebook on github