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]: