numpy arrays sont la première chose à considérer pour accélérer un algorithme. Les matrices sont présentes dans la plupart des algorithmes et numpy optimise les opérations qui s'y rapporte. Ce notebook est un parcours en diagonal.
%matplotlib inline
from jyquickhelper import add_notebook_menu
add_notebook_menu()
La convention d'import classique de numpy est la suivante:
import numpy as np
On part d'une liste python contenant des entiers. On peut créer un array numpy à partir de cette liste. Cet array possède des attributs indiquant le data type, le nombre de dimensions de l'array, etc...
l = [1, 42, 18 ]
a = np.array(l)
print(a)
print(a.dtype)
print(a.ndim)
print(a.shape)
print(a.size)
a
[ 1 42 18] int32 1 (3,) 3
array([ 1, 42, 18])
b = np.array(l, dtype=float)
print(b)
print(b.dtype)
[ 1. 42. 18.] float64
l[0] = 1.0
bb = np.array(l)
print(bb)
print(bb.dtype)
[ 1. 42. 18.] float64
Assigner un float dans un array de type int va caster le float en int, et ne modifie pas le dtype de l'array.
a[0] = 2.5
a
array([ 2, 42, 18])
On peut forcer le casting dans un autre type avec astype :
aa = a.astype(float)
aa[0] = 2.5
aa
array([ 2.5, 42. , 18. ])
A partir d'une liste de listes, on obtient un array bi-dimmensionnel.
On peut le transposer ou encore l'aplatir en un array 1d
c = np.array([range(5), range(5,10), range(5)])
print(c)
print("ndim:{}".format(c.ndim))
print("shape:{}".format(c.shape))
print(c.transpose()) #same as c.T
print("shape transposed:{}".format(c.T.shape))
print(c.flatten())
print("ndim flattened:{}".format(c.flatten().ndim))
[[0 1 2 3 4] [5 6 7 8 9] [0 1 2 3 4]] ndim:2 shape:(3, 5) [[0 5 0] [1 6 1] [2 7 2] [3 8 3] [4 9 4]] shape transposed:(5, 3) [0 1 2 3 4 5 6 7 8 9 0 1 2 3 4] ndim flattened:1
print(c)
[[0 1 2 3 4] [5 6 7 8 9] [0 1 2 3 4]]
L'indexation des array multidimensionnels fonctionne avec des tuples.
La syntaxe ':'
permet d'obtenir tous les éléments de la dimension.
print(c[1,3])
print(c[1,:3])
print(c[:,4])
8 [5 6 7] [4 9 4]
Si on utilise pas un couple sur un array 2d on récupère un array 1d
print(c[1], c[1].shape)
print(c[1][:3])
[5 6 7 8 9] (5,) [5 6 7]
On peut aussi utiliser l'indexation par un array (ou une liste python) de booléens ou d'entiers (un mask). Cela s'appelle le fancy indexing. Un mask d'entiers permet de désigner les éléments que l'on souhaite extraire via la liste de leurs indices, on peut aussi répéter l'indice d'un élément pour répéter l'élement dans l'array que l'on extrait.
ar = np.arange(1,10) #arange est l'equivalent de range mais retourne un numpy array
print('ar = ',ar)
idx = np.array([1, 4, 3, 2, 1, 7, 3])
print('idx = ',idx)
print("ar[idx] =", ar[idx])
print('######')
idx_bool = np.ones(ar.shape, dtype=bool)
idx_bool[idx] = False
print('idx_bool = ', idx_bool)
print('ar[idx_bool] = ', ar[idx_bool])
print('######', 'Que se passe-t-il dans chacun des cas suivants?', '######' )
try:
print('ar[np.array([True, True, False, True])] = ', ar[np.array([True, True, False, True])])
except Exception as e:
# l'expression ar[[True, True, False, True]] déclenche une erreur depuis numpy 1.13
print("Erreur", e)
ar = [1 2 3 4 5 6 7 8 9] idx = [1 4 3 2 1 7 3] ar[idx] = [2 5 4 3 2 8 4] ###### idx_bool = [ True False False False False True True False True] ar[idx_bool] = [1 6 7 9] ###### Que se passe-t-il dans chacun des cas suivants? ###### Erreur boolean index did not match indexed array along dimension 0; dimension is 9 but corresponding boolean dimension is 4
Pourquoi parle-t-on de fancy indexing? Essayez d'indexer des listes python de la même manière...
list_python = range(10)
list_python[[True, True, False, True]] # déclenche une exception
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-15-d185468fd957> in <module>() 1 list_python = range(10) ----> 2 list_python[[True, True, False, True]] # déclenche une exception TypeError: range indices must be integers or slices, not list
list_python[[2, 3, 2, 7]] # déclenche une exception
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-16-769cbcb8bc01> in <module>() ----> 1 list_python[[2, 3, 2, 7]] # déclenche une exception TypeError: range indices must be integers or slices, not list
Créons un array $d$. En plus de renvoyer directement un array, la fonction arange permet aussi d'utiliser un step flottant. (Essayer avec le range de python pour voir)
d = np.arange(1, 6, 0.5)
d
array([ 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5])
Un point important est que l'on ne recopie pas un array lorsqu'on effectue une assignation ou un slicing d'un array. On travaille dans ce cas avec une View sur l'array d'origine (shallow copy). Toute modification sur la View affecte l'array d'origine.
Dans l'exemple qui suit, $e$ est une view sur $d$. Lorsqu'on modifie $e$, $d$ aussi est modifié. (Remarquez au passage que numpy fournit quelques constantes bien pratiques....)
e = d
e[[0,2, 4]] = - np.pi
e
array([-3.14159265, 1.5 , -3.14159265, 2.5 , -3.14159265, 3.5 , 4. , 4.5 , 5. , 5.5 ])
d
array([-3.14159265, 1.5 , -3.14159265, 2.5 , -3.14159265, 3.5 , 4. , 4.5 , 5. , 5.5 ])
Si on ne veut pas modifier $d$ indirectement, il faut travailler sur une copie de $d$ (deep copy).
d = np.linspace(1,5.5,10) #Question subsidiaire: en quoi est-ce différent de np.arange avec un step float?
f = d.copy()
f[:4] = -np.e #il s'agit du nombre d'euler, pas de l'array e ;)
print(f)
print(d)
[-2.71828183 -2.71828183 -2.71828183 -2.71828183 3. 3.5 4. 4.5 5. 5.5 ] [ 1. 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5]
Ce point est important car source classique d'erreurs silencieuses: les erreurs les plus vicieuses car l'output sera faux mais python ne râlera pas...
Il faut un peu de temps pour s'habituer mais on finit par savoir de manière naturelle quand on travaille sur une view, quand on a besoin de faire une copie explicitement, etc... En tout cas, vérifiez vos sorties, faites des tests de cohérence, cela ne nuit jamais.
Retenez par exemple que le slicing vous renvoie une view sur l'array, alors que le fancy indexing effectue une copie.
(Au passage, remarquez le NaN (=NotaNumber) déjà introduit lors de la séance 1 sur pandas qui est un module basé sur numpy)
print('d = ',d)
slice_of_d = d[2:5]
print('\nslice_of_d = ', slice_of_d)
slice_of_d[0] = np.nan
print('\nd = ', d)
mask = np.array([2, 3, 4])
fancy_indexed_subarray = d[mask]
print('\nfancy_indexed_subarray = ', fancy_indexed_subarray)
fancy_indexed_subarray[0] = -2
print('\nd = ', d)
d = [ 1. 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5] slice_of_d = [ 2. 2.5 3. ] d = [ 1. 1.5 nan 2.5 3. 3.5 4. 4.5 5. 5.5] fancy_indexed_subarray = [ nan 2.5 3. ] d = [ 1. 1.5 nan 2.5 3. 3.5 4. 4.5 5. 5.5]
La méthode reshape permet de changer la forme de l'array. Il existe de nombreuses manipulations possibles.
On précise à reshape la forme souhaitée: par un entier si on veut un array 1d de cette longueur, ou un couple pour un array 2d de cette forme.
g = np.arange(12)
print(g)
g.reshape((4,3))
[ 0 1 2 3 4 5 6 7 8 9 10 11]
array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11]])
Par défaut, reshape utilise l'énumération dans l'ordre du langage C (aussi appelé "row first" ), on peut préciser que l'on souhaite utiliser l'ordre de Fortran ("column first"). Ceux qui connaissent Matlab et R sont habitués à l'ordre "column-first". Voir l'article wikipedia
g.reshape((4,3), order='F')
array([[ 0, 4, 8], [ 1, 5, 9], [ 2, 6, 10], [ 3, 7, 11]])
On peut utiliser -1 sur une dimension, cela sert de joker: numpy infère la dimension nécessaire ! On peut créer directement des matrices de 0 et de 1 à la dimension d'un autre array.
np.zeros_like(g)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
np.ones_like(g)
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
On peut aussi concatener ou stacker horizontalement/verticalement différents arrays.
np.concatenate((g, np.zeros_like(g))) #Attention à la syntaxe: le type d'entrée est un tuple!
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
gmat = g.reshape((1, len(g)))
np.concatenate((gmat, np.ones_like(gmat)), axis=0)
array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
np.concatenate((gmat, np.ones_like(gmat)), axis=1)
array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
np.hstack((g, g))
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
np.vstack((g,g))
array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]])
#Exo1a-1:
#Exo1a-2:
#Exo1B:
#Exo1C:
Il existe un très grand nombre de routines pour manipuler les arrays numpy: Vous trouverez sans doute utiles les pages spécifiques aux routines de stats ou de maths
On déclare $a$ et $b$ sur lesquelles nous allons illustrer quelques opérations
a = np.ones((3,2))
b = np.arange(6).reshape(a.shape)
print(a)
b
[[ 1. 1.] [ 1. 1.] [ 1. 1.]]
array([[0, 1], [2, 3], [4, 5]])
Les opérations arithmétiques avec les scalaires, ou entre arrays s'effectuent élément par élément. Lorsque le dtype n'est pas le même ($a$ contient des float, $b$ contient des int), numpy adopte le type le plus "grand" (au sens de l'inclusion).
print( (a + b)**2 )
print( np.abs( 3*a - b ) )
f = lambda x: np.exp(x-1)
print( f(b) )
[[ 1. 4.] [ 9. 16.] [ 25. 36.]] [[ 3. 2.] [ 1. 0.] [ 1. 2.]] [[ 0.36787944 1. ] [ 2.71828183 7.3890561 ] [ 20.08553692 54.59815003]]
Remarquez que la division par zéro ne provoque pas d'erreur mais introduit la valeur inf :
b
array([[0, 1], [2, 3], [4, 5]])
1/b
c:\Python36_x64\lib\site-packages\ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in true_divide """Entry point for launching an IPython kernel.
array([[ inf, 1. ], [ 0.5 , 0.33333333], [ 0.25 , 0.2 ]])
Que se passe-t-il si les dimensions sont différentes?
c = np.ones(6)
c
array([ 1., 1., 1., 1., 1., 1.])
b+c # déclenche une exception
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-39-882b3e9536b7> in <module>() ----> 1 b+c # déclenche une exception ValueError: operands could not be broadcast together with shapes (3,2) (6,)
c = np.arange(3).reshape((3,1))
print(b,c, sep='\n')
b+c
[[0 1] [2 3] [4 5]] [[0] [1] [2]]
array([[0, 1], [3, 4], [6, 7]])
L'opération précédente fonctionne car numpy effectue ce qu'on appelle un broadcasting de c
: une dimension étant commune, tout se passe comme si on dupliquait c sur la dimension non-partagée avec b. Vous trouverez une explication visuelle simple ici :
a = np.zeros((3,3))
a[:,0] = -1
b = np.array(range(3))
print(a + b)
[[-1. 1. 2.] [-1. 1. 2.] [-1. 1. 2.]]
Par contre, il peut parfois être utile de préciser la dimension sur laquelle on souhaite broadcaster, on ajoute alors explicitement une dimension :
print(b.shape)
print(b[:,np.newaxis].shape)
print(b[np.newaxis,:].shape)
(3,) (3, 1) (1, 3)
print( a + b[np.newaxis,:] )
print( a + b[:,np.newaxis] )
print(b[:,np.newaxis]+b[np.newaxis,:])
print(b + b)
[[-1. 1. 2.] [-1. 1. 2.] [-1. 1. 2.]] [[-1. 0. 0.] [ 0. 1. 1.] [ 1. 2. 2.]] [[0 1 2] [1 2 3] [2 3 4]] [0 2 4]
On parle de réductions lorsque l'opération réduit la dimension de l'array. Il en existe un grand nombre. Elles existent souvent sous forme de fonction de numpy ou de méthodes d'un array numpy. On n'en présente que quelques unes, mais le principe est le même : par défaut elles opèrent sur toutes les dimensions, mais on peut via l'argument axis préciser la dimension selon laquelle on souhaite effectuer la réduction.
c = np.arange(10).reshape((2,-1)) #Note: -1 is a joker!
print(c)
print(c.sum())
print(c.sum(axis=0))
print(np.sum(c, axis=1))
[[0 1 2 3 4] [5 6 7 8 9]] 45 [ 5 7 9 11 13] [10 35]
print(np.all(c[0] < c[1]))
print(c.min(), c.max())
print(c.min(axis=1))
True 0 9 [0 5]
Commençons par construire deux arrays 2d correspondant à une matrice triangulaire inférieure et une matrice diagonale :
A = np.tril(np.ones((3,3)))
A
array([[ 1., 0., 0.], [ 1., 1., 0.], [ 1., 1., 1.]])
b = np.diag([1,2, 3])
b
array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
print(A.dot(b))
print(A*b)
print(A.dot(A))
[[ 1. 0. 0.] [ 1. 2. 0.] [ 1. 2. 3.]] [[ 1. 0. 0.] [ 0. 2. 0.] [ 0. 0. 3.]] [[ 1. 0. 0.] [ 2. 1. 0.] [ 3. 2. 1.]]
On peut calculer l'inverse ou le déterminant de $A$
print(np.linalg.det(A))
inv_A = np.linalg.inv(A)
print(inv_A)
print(inv_A.dot(A))
1.0 [[ 1. 0. 0.] [-1. 1. 0.] [ 0. -1. 1.]] [[ 1. 0. 0.] [ 0. 1. 0.] [ 0. 0. 1.]]
... résoudre des systèmes d'equations linéaires du type $Ax = b$...
x = np.linalg.solve(A, np.diag(b))
print(np.diag(b))
print(x)
print(A.dot(x))
[1 2 3] [ 1. 1. 1.] [ 1. 2. 3.]
... ou encore obtenir les valeurs propres de $A$.
np.linalg.eig(A)
(array([ 1., 1., 1.]), array([[ 0.00000000e+00, 0.00000000e+00, 4.93038066e-32], [ 0.00000000e+00, 2.22044605e-16, -2.22044605e-16], [ 1.00000000e+00, -1.00000000e+00, 1.00000000e+00]]))
np.linalg.eigvals(A)
array([ 1., 1., 1.])
Matrix est une sous classe spécialisée pour le calcul matriciel. Il s'agit d'un array numpy 2d qui conserve sa dimension 2d à travers les opérations. Pensez aux différences que cela implique... On peut les construire classiquement depuis les array ou les objets pythons, ou via une string à la Matlab ( où les points virgules indiquent les lignes).
m = np.matrix(' 1 2 3; 4 5 6; 7 8 9')
a = np.arange(1,10).reshape((3,3))
print(m)
print(a)
print(m[0], a[0])
print(m[0].shape, a[0].shape)
[[1 2 3] [4 5 6] [7 8 9]] [[1 2 3] [4 5 6] [7 8 9]] [[1 2 3]] [1 2 3] (1, 3) (3,)
Matrix surcharge par ailleurs les opérateurs * et ** pour remplacer les opérations élément par élément par les opérations matricielles. Enfin, une Matrix possède des attributs supplémentaires. Notamment, Matrix.I qui désigne l'inverse, Matrix.A l'array de base.
Il est probable que cela évolue puisque Python 3.5 a introduit le symbol @
pour la multiplication matricielle.
m * m
matrix([[ 30, 36, 42], [ 66, 81, 96], [102, 126, 150]])
a * a
array([[ 1, 4, 9], [16, 25, 36], [49, 64, 81]])
m * a # La priorité des matrix est plus importantes que celles des arrays
matrix([[ 30, 36, 42], [ 66, 81, 96], [102, 126, 150]])
print(m**2)
print(a**2)
[[ 30 36 42] [ 66 81 96] [102 126 150]] [[ 1 4 9] [16 25 36] [49 64 81]]
La syntaxe est plus légère pour effectuer du calcul matriciel
m[0,0]= -1
print("det", np.linalg.det(m), "rank",np.linalg.matrix_rank(m))
print(m.I*m)
a[0,0] = -1
print("det", np.linalg.det(a), "rank",np.linalg.matrix_rank(a))
print(a.dot(np.linalg.inv(a)))
det 6.0 rank 3 [[ 1.00000000e+00 3.99680289e-15 4.49640325e-15] [ 0.00000000e+00 1.00000000e+00 0.00000000e+00] [ -8.88178420e-16 0.00000000e+00 1.00000000e+00]] det 6.0 rank 3 [[ 1.00000000e+00 8.88178420e-16 -1.77635684e-15] [ -6.66133815e-16 1.00000000e+00 0.00000000e+00] [ -3.33066907e-16 2.66453526e-15 1.00000000e+00]]
Le module numpy.random apporte à python la possibilité de générer un échantillon de taille $n$ directement, alors que le module natif de python ne produit des tirages que un par un. Le module numpy.random est donc bien plus efficace si on veut tirer des échantillon conséquents. Par ailleurs, scipy.stats fournit des méthodes pour un très grand nombre de distributions et quelques fonctions classiques de statistiques.
np.random.randn(4,3)
array([[-0.53862576, 0.7316812 , -0.43393759], [-0.39077735, -1.48022294, 0.61423791], [ 1.29123337, -2.92158205, -2.33375479], [-0.63012998, 0.37943656, 0.33758665]])
Pour se convaincre que numpy.random est plus efficace que le module random de base de python. On effectue un grand nombre de tirages gaussiens standard, en python pur et via numpy.
N = int(1e7)
from random import normalvariate
%timeit [normalvariate(0,1) for _ in range(N)]
9.04 s ± 149 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.random.randn(N)
301 ms ± 7.46 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Simulez (en une seule fois!) 10000 marches aléatoires de taille 1000, partant de 0 et de pas +1 ou -1 équiprobables
Vous aurez peut-être besoin des fonctions suivantes: np.abs, np.mean, np.max, np.where, np.argmax, np.any, np.cumsum, np.random.randint.
L'exercice précédent montre comment générer une marche aléatoire à partir d'une série temporelle aléatoire. Comment retrouver la série initiale à partir de la marche aléatoire ?
Le module scipy.optimize fournit un panel de méthodes d'optimisation. En fonction du problème que vous souhaitez résoudre, il vous faut choisir la méthode adéquate. Je vous conseille vivement la lecture de ce tutoriel sur l'optimisation numérique, écrit par Gaël Varoquaux.
Récemment, l'ensemble des solvers ont été regroupés sous deux interfaces, même si on peut toujours faire appel à chaque solver directement, ce qui n'est pas conseillé car les entrées sorties ne sont pas normalisées (par contre vous devrez sans doute aller voir l'aide de chaque méthode pour vous en servir):
Vous obtiendrez en sortie un objet de type scipy.optimize.OptimizeResult.
Dans la suite, je développe un petit exemple inspiré du tutoriel de la toolbox d'optimisation de Matlab. Par ailleurs, la documentation de cette toolbox est plutôt claire et peut toujours vous servir lorsque que vous avez besoin de vous rafraichir la mémoire sur l'optimisation numérique.
On commence par définir la fonction bowl_peak
def bowl_peak(x,y):
return x*np.exp(-x**2-y**2)+(x**2+y**2)/20
On va ensuite chercher un exemple dans la gallerie matplotlib pour la représenter: contour3d_demo3. On modifie légèrement le code pour l'utiliser avec bowl_peak
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm #colormaps
min_val = -2
max_val = 2
fig = plt.figure()
ax = fig.gca(projection='3d')
x_axis = np.linspace(min_val,max_val,100)
y_axis = np.linspace(min_val,max_val,100)
X, Y = np.meshgrid(x_axis, y_axis, copy=False, indexing='xy')
Z = bowl_peak(X,Y)
#X, Y, Z = axes3d.get_test_data(0.05)
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, alpha=0.2)
cset = ax.contour(X, Y, Z, zdir='z', offset=-0.5, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=min_val, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='y', offset=max_val, cmap=cm.coolwarm)
ax.set_xlabel('X')
ax.set_xlim(min_val, max_val)
ax.set_ylabel('Y')
ax.set_ylim(min_val, max_val)
ax.set_zlabel('Z')
ax.set_zlim(-0.5, 0.5)
(-0.5, 0.5)
On voit que le minimum se trouve près de $[-\frac{1}{2}, 0]$. On va utiliser ce point pour initialiser l'optimisation. On va tester différentes méthodes et comparer les sorties obtenues.
from scipy import optimize
x0 = np.array([-0.5, 0])
fun = lambda x: bowl_peak(x[0],x[1])
methods = [ 'Nelder-Mead', 'CG', 'BFGS', 'Powell', 'COBYLA', 'L-BFGS-B' ]
for m in methods:
optim_res = optimize.minimize(fun, x0, method=m)
print("---\nMethod:{}\n".format(m),optim_res, "\n")
--- Method:Nelder-Mead final_simplex: (array([[ -6.69025421e-01, -1.44567490e-04], [ -6.69110179e-01, -1.81386054e-04], [ -6.68989849e-01, -2.01126337e-04]]), array([-0.40523686, -0.40523685, -0.40523684])) fun: -0.40523685823917283 message: 'Optimization terminated successfully.' nfev: 38 nit: 20 status: 0 success: True x: array([ -6.69025421e-01, -1.44567490e-04]) --- Method:CG fun: -0.4052368583334503 jac: array([ -2.12926418e-04, 3.72529030e-09]) message: 'Desired error not necessarily achieved due to precision loss.' nfev: 24 nit: 1 njev: 3 status: 2 success: False x: array([ -6.69183901e-01, -3.71395638e-09]) --- Method:BFGS fun: -0.40523687025688715 hess_inv: array([[ 0.52865446, 0. ], [ 0. , 1. ]]) jac: array([ -6.08339906e-06, 0.00000000e+00]) message: 'Optimization terminated successfully.' nfev: 28 nit: 6 njev: 7 status: 0 success: True x: array([ -6.69075034e-01, -7.45058060e-09]) --- Method:Powell direc: array([[ 0.00000000e+00, 1.00000000e+00], [ -6.85432298e-04, -4.67045589e-11]]) fun: -0.40523687026669025 message: 'Optimization terminated successfully.' nfev: 62 nit: 2 status: 0 success: True x: array([ -6.69071822e-01, -1.15386055e-08]) --- Method:COBYLA fun: -0.4052368678399868 maxcv: 0.0 message: 'Optimization terminated successfully.' nfev: 32 status: 1 success: True x: array([ -6.69108584e-01, -4.89154557e-05]) --- Method:L-BFGS-B fun: -0.40523687026621352 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64> jac: array([ 1.35447209e-06, 0.00000000e+00]) message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' nfev: 15 nit: 3 status: 0 success: True x: array([ -6.69071114e-01, -8.35621530e-09])
On trouve un minimum à $-0.4052$ en $[-0.669, 0.000]$ pour toutes les méthodes qui convergent. Notez le message de sortie de 'CG' qui signifie que le gradient ne varie plus assez. Personnellement, je ne trouve pas ce message de sortie très clair. Le point trouvé est bien l'optimum cherché pourtant. Notez aussi le nombre d'évaluations de la fonction (nfev) pour chaque méthode, et le nombre d'évaluation de gradient (njev) pour les méthodes qui reposent sur un calcul de gradient.
Remarquez aussi que si on relance Anneal plusieurs fois, on n'est pas assuré d'obtenir la même solution, puisqu'il s'agit d'une métaheuristique.
for i in range(4):
optim_res = optimize.minimize(fun, x0, method='BFGS')
print("---\nMethod:{} - Test:{}\n".format(m,i),optim_res, "\n")
--- Method:L-BFGS-B - Test:0 fun: -0.40523687025688715 hess_inv: array([[ 0.52865446, 0. ], [ 0. , 1. ]]) jac: array([ -6.08339906e-06, 0.00000000e+00]) message: 'Optimization terminated successfully.' nfev: 28 nit: 6 njev: 7 status: 0 success: True x: array([ -6.69075034e-01, -7.45058060e-09]) --- Method:L-BFGS-B - Test:1 fun: -0.40523687025688715 hess_inv: array([[ 0.52865446, 0. ], [ 0. , 1. ]]) jac: array([ -6.08339906e-06, 0.00000000e+00]) message: 'Optimization terminated successfully.' nfev: 28 nit: 6 njev: 7 status: 0 success: True x: array([ -6.69075034e-01, -7.45058060e-09]) --- Method:L-BFGS-B - Test:2 fun: -0.40523687025688715 hess_inv: array([[ 0.52865446, 0. ], [ 0. , 1. ]]) jac: array([ -6.08339906e-06, 0.00000000e+00]) message: 'Optimization terminated successfully.' nfev: 28 nit: 6 njev: 7 status: 0 success: True x: array([ -6.69075034e-01, -7.45058060e-09]) --- Method:L-BFGS-B - Test:3 fun: -0.40523687025688715 hess_inv: array([[ 0.52865446, 0. ], [ 0. , 1. ]]) jac: array([ -6.08339906e-06, 0.00000000e+00]) message: 'Optimization terminated successfully.' nfev: 28 nit: 6 njev: 7 status: 0 success: True x: array([ -6.69075034e-01, -7.45058060e-09])
On va évaluer le temps de calcul nécessaire à chaque méthode.
for m in methods:
print("Method:{}:".format(m))
%timeit optim_res = optimize.minimize(fun, x0, method=m)
print('############')
Method:Nelder-Mead: 894 µs ± 40.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############ Method:CG: 788 µs ± 49 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############ Method:BFGS: 592 µs ± 4.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############ Method:Powell: 1.01 ms ± 26.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############ Method:COBYLA: 193 µs ± 12.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############ Method:L-BFGS-B: 222 µs ± 18.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) ############
On peut aussi fournir des arguments supplémentaires à la fonction qu'on optimise. Par exemple, les données lorsque vous maximisez une log-vraissemblance. En voici un exemple: on considère une version rescaled de la fonction bowl_peak. Vous pourriez aussi utiliser une lambda fonction.
def shifted_scaled_bowlpeak(x,a,b,c):
return (x[0]-a)*np.exp(-((x[0]-a)**2+(x[1]-b)**2))+((x[0]-a)**2+(x[0]-b)**2)/c
a = 2
b = 3
c = 10
optim_res = optimize.minimize(shifted_scaled_bowlpeak, x0, args=(a,b,c), method='BFGS')
print(optim_res)
print('#######')
optim_res = optimize.minimize(lambda x:shifted_scaled_bowlpeak(x,a,b,c), x0, method='BFGS')
print(optim_res)
fun: 0.05000000675226609 hess_inv: array([[ 1.40782352e+00, -1.59338758e+02], [ -1.59338758e+02, 7.19318682e+05]]) jac: array([ -9.78726894e-06, 5.63450158e-08]) message: 'Optimization terminated successfully.' nfev: 96 nit: 23 njev: 24 status: 0 success: True x: array([ 2.49997551, -1.22943768]) ####### fun: 0.05000000675226609 hess_inv: array([[ 1.40782352e+00, -1.59338758e+02], [ -1.59338758e+02, 7.19318682e+05]]) jac: array([ -9.78726894e-06, 5.63450158e-08]) message: 'Optimization terminated successfully.' nfev: 96 nit: 23 njev: 24 status: 0 success: True x: array([ 2.49997551, -1.22943768])
Vous pouvez continuer ce petit benchmark en ajoutant le gradient et la hessienne... les calculs seront plus précis et plus rapides.
Voir l'exercice 1 ici