.. _td2acenoncesession2Arst: ======================================== 2A.data - Calcul Matriciel, Optimisation ======================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/td2a/td2a_cenonce_session_2A.ipynb|*` `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. .. code:: ipython3 %matplotlib inline .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Numpy arrays ~~~~~~~~~~~~ La convention d’import classique de numpy est la suivante: .. code:: ipython3 import numpy as np Creation d’un array: notion de datatype, et dimensions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 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… .. code:: ipython3 l = [1, 42, 18 ] a = np.array(l) print(a) print(a.dtype) print(a.ndim) print(a.shape) print(a.size) a .. parsed-literal:: [ 1 42 18] int32 1 (3,) 3 .. parsed-literal:: array([ 1, 42, 18]) On peut indiquer explicitement le `dtype `__ lors de la création de l’array. Sinon, Numpy sélectionne automatiquement le `dtype `__. Numpy ajoute un grand nombre de `dtype `__ à ceux de Python. Allez jeter un oeil à la `liste `__. .. code:: ipython3 b = np.array(l, dtype=float) print(b) print(b.dtype) .. parsed-literal:: [ 1. 42. 18.] float64 .. code:: ipython3 l[0] = 1.0 bb = np.array(l) print(bb) print(bb.dtype) .. parsed-literal:: [ 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. .. code:: ipython3 a[0] = 2.5 a .. parsed-literal:: array([ 2, 42, 18]) On peut forcer le casting dans un autre type avec `astype `__ : .. code:: ipython3 aa = a.astype(float) aa[0] = 2.5 aa .. parsed-literal:: 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 .. code:: ipython3 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)) .. parsed-literal:: [[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 Indexation, Slicing, Fancy indexing ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: ipython3 print(c) .. parsed-literal:: [[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. .. code:: ipython3 print(c[1,3]) print(c[1,:3]) print(c[:,4]) .. parsed-literal:: 8 [5 6 7] [4 9 4] Si on utilise pas un couple sur un array 2d on récupère un array 1d .. code:: ipython3 print(c[1], c[1].shape) print(c[1][:3]) .. parsed-literal:: [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. .. code:: ipython3 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) .. parsed-literal:: 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… .. code:: ipython3 list_python = range(10) list_python[[True, True, False, True]] # déclenche une exception :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in () 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 .. code:: ipython3 list_python[[2, 3, 2, 7]] # déclenche une exception :: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in () ----> 1 list_python[[2, 3, 2, 7]] # déclenche une exception TypeError: range indices must be integers or slices, not list View contre Copy ^^^^^^^^^^^^^^^^ Créons un array :math:`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) .. code:: ipython3 d = np.arange(1, 6, 0.5) d .. parsed-literal:: 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, :math:`e` est une view sur :math:`d`. Lorsqu’on modifie :math:`e`, :math:`d` aussi est modifié. (Remarquez au passage que numpy fournit quelques constantes bien pratiques….) .. code:: ipython3 e = d e[[0,2, 4]] = - np.pi e .. parsed-literal:: array([-3.14159265, 1.5 , -3.14159265, 2.5 , -3.14159265, 3.5 , 4. , 4.5 , 5. , 5.5 ]) .. code:: ipython3 d .. parsed-literal:: 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 :math:`d` indirectement, il faut travailler sur une copie de :math:`d` (`deep copy `__). .. code:: ipython3 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) .. parsed-literal:: [-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) .. code:: ipython3 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) .. parsed-literal:: 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] Manipulation de shape ^^^^^^^^^^^^^^^^^^^^^ 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. .. code:: ipython3 g = np.arange(12) print(g) g.reshape((4,3)) .. parsed-literal:: [ 0 1 2 3 4 5 6 7 8 9 10 11] .. parsed-literal:: 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 `__ .. code:: ipython3 g.reshape((4,3), order='F') .. parsed-literal:: 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. .. code:: ipython3 np.zeros_like(g) .. parsed-literal:: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) .. code:: ipython3 np.ones_like(g) .. parsed-literal:: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) On peut aussi concatener ou stacker `horizontalement `__/`verticalement `__ différents arrays. .. code:: ipython3 np.concatenate((g, np.zeros_like(g))) #Attention à la syntaxe: le type d'entrée est un tuple! .. parsed-literal:: 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]) .. code:: ipython3 gmat = g.reshape((1, len(g))) np.concatenate((gmat, np.ones_like(gmat)), axis=0) .. parsed-literal:: 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]]) .. code:: ipython3 np.concatenate((gmat, np.ones_like(gmat)), axis=1) .. parsed-literal:: 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]]) .. code:: ipython3 np.hstack((g, g)) .. parsed-literal:: 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]) .. code:: ipython3 np.vstack((g,g)) .. parsed-literal:: 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]]) Exercice 1: Echiquier et Crible d’Erathosthène ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Exercice 1-A Echiquier: Créer une matrice échiquier (des 1 et des 0 alternés) de taille 8x8, de deux façons différentes - en vous servant de slices - en vous servant de la fonction `tile `__ - Exercice 1-B Piège lors d’une extraction 2d: - Définir la matrice :math:`M = \left(\begin{array}{ccccc} 1 & 5 & 9 & 13 & 17 \\ 2 & 6 & 10 & 14 & 18 \\ 3 & 7 & 11 & 15 & 19 \\ 4 & 8 & 12 & 16 & 20 \\ \end{array}\right)` - En **extraire** la matrice :math:`\left(\begin{array}{ccc} 6 & 18 & 10 \\ 7 & 19 & 11 \\ 5 & 17 & 9 \\ \end{array}\right)` - Exercice 1-C Crible d’Erathosthène: On souhaite implémenter un `crible d’Erathosthène `__ pour trouver les nombres premiers inférieurs à :math:`N=1000`. - partir d’un array de booléens de taille N+1, tous égaux à True. - Mettre 0 et 1 à False car ils ne sont pas premiers - pour chaque entier :math:`k` entre 2 et :math:`\sqrt{N}`: - si :math:`k` est premier: on passe ses multiples (entre :math:`k^2` et :math:`N`) à False - on print la liste des entiers premiers .. code:: ipython3 #Exo1a-1: #Exo1a-2: .. code:: ipython3 #Exo1B: .. code:: ipython3 #Exo1C: Manipulation et Opérations sur les arrays ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 `__ Opérations élément par élément ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ On déclare :math:`a` et :math:`b` sur lesquelles nous allons illustrer quelques opérations .. code:: ipython3 a = np.ones((3,2)) b = np.arange(6).reshape(a.shape) print(a) b .. parsed-literal:: [[ 1. 1.] [ 1. 1.] [ 1. 1.]] .. parsed-literal:: 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 (:math:`a` contient des float, :math:`b` contient des int), numpy adopte le type le plus “grand” (au sens de l’inclusion). .. code:: ipython3 print( (a + b)**2 ) print( np.abs( 3*a - b ) ) f = lambda x: np.exp(x-1) print( f(b) ) .. parsed-literal:: [[ 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 `__ : .. code:: ipython3 b .. parsed-literal:: array([[0, 1], [2, 3], [4, 5]]) .. code:: ipython3 1/b .. parsed-literal:: 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. .. parsed-literal:: array([[ inf, 1. ], [ 0.5 , 0.33333333], [ 0.25 , 0.2 ]]) Broadcasting ^^^^^^^^^^^^ Que se passe-t-il si les dimensions sont différentes? .. code:: ipython3 c = np.ones(6) c .. parsed-literal:: array([ 1., 1., 1., 1., 1., 1.]) .. code:: ipython3 b+c # déclenche une exception :: --------------------------------------------------------------------------- ValueError Traceback (most recent call last) in () ----> 1 b+c # déclenche une exception ValueError: operands could not be broadcast together with shapes (3,2) (6,) .. code:: ipython3 c = np.arange(3).reshape((3,1)) print(b,c, sep='\n') b+c .. parsed-literal:: [[0 1] [2 3] [4 5]] [[0] [1] [2]] .. parsed-literal:: 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 `__ : .. code:: ipython3 a = np.zeros((3,3)) a[:,0] = -1 b = np.array(range(3)) print(a + b) .. parsed-literal:: [[-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 : .. code:: ipython3 print(b.shape) print(b[:,np.newaxis].shape) print(b[np.newaxis,:].shape) .. parsed-literal:: (3,) (3, 1) (1, 3) .. code:: ipython3 print( a + b[np.newaxis,:] ) print( a + b[:,np.newaxis] ) print(b[:,np.newaxis]+b[np.newaxis,:]) print(b + b) .. parsed-literal:: [[-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] Réductions ^^^^^^^^^^ 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. .. code:: ipython3 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)) .. parsed-literal:: [[0 1 2 3 4] [5 6 7 8 9]] 45 [ 5 7 9 11 13] [10 35] .. code:: ipython3 print(np.all(c[0] < c[1])) print(c.min(), c.max()) print(c.min(axis=1)) .. parsed-literal:: True 0 9 [0 5] Algèbre linéaire ~~~~~~~~~~~~~~~~ Vous avez un éventail de fonctions pour faire de l’algèbre linéaire dans `numpy `__ ou dans `scipy `__. Cela peut vous servir si vous cherchez à faire une décomposition matricielle particulière (`LU `__, `QR `__, `SVD `__,…), si vous vous intéressez aux valeurs propres d’une matrice, etc… Exemples simples ^^^^^^^^^^^^^^^^ Commençons par construire deux arrays 2d correspondant à une matrice triangulaire inférieure et une matrice diagonale : .. code:: ipython3 A = np.tril(np.ones((3,3))) A .. parsed-literal:: array([[ 1., 0., 0.], [ 1., 1., 0.], [ 1., 1., 1.]]) .. code:: ipython3 b = np.diag([1,2, 3]) b .. parsed-literal:: array([[1, 0, 0], [0, 2, 0], [0, 0, 3]]) On a vu que les multiplications entre array s’effectuaient élément par élement. Si l’on souhaite faire des multiplications matricielles, il faut utiliser la fonction `dot `__. La version 3.5 introduit un nouvel opérateur [@](https://docs.python.org/3.6/whatsnew/3.5.html#pep-465-a-dedicated-infix-operator-for-matrix-multiplication) qui désigne explicitement la multiplication matricielle. .. code:: ipython3 print(A.dot(b)) print(A*b) print(A.dot(A)) .. parsed-literal:: [[ 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 :math:`A` .. code:: ipython3 print(np.linalg.det(A)) inv_A = np.linalg.inv(A) print(inv_A) print(inv_A.dot(A)) .. parsed-literal:: 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 :math:`Ax = b`\ … .. code:: ipython3 x = np.linalg.solve(A, np.diag(b)) print(np.diag(b)) print(x) print(A.dot(x)) .. parsed-literal:: [1 2 3] [ 1. 1. 1.] [ 1. 2. 3.] … ou encore obtenir les valeurs propres de :math:`A`. .. code:: ipython3 np.linalg.eig(A) .. parsed-literal:: (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]])) .. code:: ipython3 np.linalg.eigvals(A) .. parsed-literal:: array([ 1., 1., 1.]) Numpy Matrix ^^^^^^^^^^^^ `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). .. code:: ipython3 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) .. parsed-literal:: [[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.* .. code:: ipython3 m * m .. parsed-literal:: matrix([[ 30, 36, 42], [ 66, 81, 96], [102, 126, 150]]) .. code:: ipython3 a * a .. parsed-literal:: array([[ 1, 4, 9], [16, 25, 36], [49, 64, 81]]) .. code:: ipython3 m * a # La priorité des matrix est plus importantes que celles des arrays .. parsed-literal:: matrix([[ 30, 36, 42], [ 66, 81, 96], [102, 126, 150]]) .. code:: ipython3 print(m**2) print(a**2) .. parsed-literal:: [[ 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 .. code:: ipython3 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))) .. parsed-literal:: 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]] Génération de nombres aléatoires et statistiques ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Le module `numpy.random `__ apporte à python la possibilité de générer un échantillon de taille :math:`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. Par exemple, on peut obtenir un array 4x3 de tirages gaussiens standard (soit en utilisant `randn `__ ou `normal `__): .. code:: ipython3 np.random.randn(4,3) .. parsed-literal:: 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. .. code:: ipython3 N = int(1e7) from random import normalvariate %timeit [normalvariate(0,1) for _ in range(N)] .. parsed-literal:: 9.04 s ± 149 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) .. code:: ipython3 %timeit np.random.randn(N) .. parsed-literal:: 301 ms ± 7.46 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) Exercice 2 : marches aléatoires ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Simulez (**en une seule fois!**) 10000 marches aléatoires de taille 1000, partant de 0 et de pas +1 ou -1 équiprobables - Faites un graphe représentant la racine de la moyenne des carrés des positions (=cumul des pas à un instant donné) en fonction du temps - Quels sont les amplitudes maximales et minimales atteintes parmi l’ensemble des marches aléatoires? - Combien de marches s’éloigne de plus de 50 de l’origine? - Parmi celles qui le font, quelle est la moyenne des temps de passage (i.e. le premier moment où ces marches dépassent +/-50)? 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 `__. Exercice 3 : retrouver la série aléatoire à partir des marches aléatoires ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 ? Optimisation avec scipy ~~~~~~~~~~~~~~~~~~~~~~~ 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): - Pour minimiser une fonction scalaire d’une ou plusieurs variables:\ `scipy.optimize.minimize `__ - Pour minimiser une fonction scalaire d’une variable uniquement:\ `scipy.optimize.minimize_scalar `__ 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* .. code:: ipython3 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* .. code:: ipython3 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) .. parsed-literal:: (-0.5, 0.5) .. image:: td2a_cenonce_session_2A_128_1.png On voit que le minimum se trouve près de :math:`[-\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. .. code:: ipython3 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") .. parsed-literal:: --- 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 à :math:`-0.4052` en :math:`[-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. .. code:: ipython3 for i in range(4): optim_res = optimize.minimize(fun, x0, method='BFGS') print("---\nMethod:{} - Test:{}\n".format(m,i),optim_res, "\n") .. parsed-literal:: --- 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. .. code:: ipython3 for m in methods: print("Method:{}:".format(m)) %timeit optim_res = optimize.minimize(fun, x0, method=m) print('############') .. parsed-literal:: 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. .. code:: ipython3 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) .. parsed-literal:: 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. Exercice 4: simulation, régression, estimation par maximisation de la vraisemblance ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - On commence par simuler la variable :math:`Y = 3 X_1 -2 X_2 +2 + \epsilon` où :math:`X_1,X_2,\epsilon \sim \mathcal{N}(0,1)` - On souhaite ensuite retrouver les coefficients dans la `régression linéaire `__ de :math:`Y` sur :math:`X_1` et :math:`X_2` dans un modèle avec constante, par la méthode des Moindres Carrés Ordinaires. On rappelle que la forme matricielle de l’estimateur des MCO est :math:`\hat{\beta} = (X'X)^{-1}X'Y` - Enfin, :math:`Y` étant normale, on souhaite estimer ses paramètres par maximisation de vraisemblance: - La densité s’écrit: :math:`f(x, \mu, \sigma) = \frac{1}{\sigma \sqrt{2\pi} } e^{ -\frac{(x-\mu)^2}{2\sigma^2} }` - La log-vraisemblance: :math:`\ln\mathcal{L}(\mu,\sigma^2) = \sum_{i=1}^n \ln f(x_i;\,\mu,\sigma^2) = -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln\sigma^2 - \frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2`. - L’écriture des conditions au premier ordre donne une formule fermée pour les estimateurs du maximum de vraisemblance: :math:`\hat{\mu} = \overline{x} \equiv \frac{1}{n}\sum_{i=1}^n x_i`, :math:`\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \overline{x})^2`. - Vérifiez en les implémentant directement que vous trouvez bien la même solution que le minimum obtenu en utilisant `scipy.optimize.minimize `__ pour minimiser l’opposé de la log-vraissemblance. Exercice 5 : Optimisation quadratique (sous contraintes) avec cvxopt ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Voir l’exercice 1 `ici `__ Références ~~~~~~~~~~ - `100 numpy exercises `__ - `Un tutoriel bien fait et très complet sur numpy `__. L’un des auteurs n’est autre que Gaël Varoquaux qui sera présent pour la séance 3. Voir aussi l’\ `ensemble du tutoriel `__ et notamment la `partie optimisation `__