5.3. Introducción a NumPy#

Hasta ahora, dentro de Python, hemos trabajado con el módulo SymPy, que nos permite manejar cálculo simbólico, definir nuevas funciones, calcular límites o, como veremos más adelante, derivar e integrar funciones.

Sin embargo, las posibilidades de Python para el cálculo matemático no acaban, ni mucho menos, aquí. Ahora giraremos nuestra cabeza para buscar una herramienta que convierta a Python en una especie de super-calculadora, con enormes posibilidades a la hora de realizar grandes operaciones repetitivas con números, vectores o matrices. Nosotros le sacaremos el máximo jugo en la programación de los métodos numéricos que estudiaremos en esta materia pero, por supuesto, sus posibilidades son mucho mayores.

Yendo al grano, os vamos a presentar la librería NumPy (http://www.numpy.org/) que es la más utilizada para cálculos numéricos. Esta librería permite utilizar un enorme abanico de estructuras de datos y funciones para los cálculos numéricos. En esta sección os mostraremos algunas de estas funcionalidades.

NumPy es una librería enorme y tiene posibilidades muy extensas, por lo que sólo pretendemos realizar una introducción muy breve. Para descubrir nuevas funcionalidades y como utilizarla con muchos otros propósitos, la mejor fuente de información en línea son los buscadores, y, en particular, la comunidad http://stackoverflow.com/.

5.3.1. Objetivos#

  • Trabajar con objetos y métodos.

  • Introducción a los vectores unidimensionales de números (numpy.array).

  • Aplicar operacines numéricas elementales.

  • Manipulación de vectores numéricos (indexado, troceado, etc.).

  • Ejercicio: eficiencia de NumPy en funciones vectorizadas y no vectorizadas.

Importar el módulo NumPy

Antes de nada, para tener disponible NumPy en el código, debemos primero importarlo. Para hacer esto, es habitual importar NumPy utilizando el alias np:

import numpy as np

5.3.2. Programación orientada a objetos#

Como todos los módulos de Python, la librería NumPy está implementada siguiendo una estrategia de programación orientada a objetos. Por lo tanto, cualquier estructura de datos en Python, incluso la más simple, se debe entender como un objeto que pertenece a una clase, y sobre ella podemos realizar todas las operaciones implementadas sobre esta clase de objetos. Incluso un número es un objeto de una clase:

a = 3.3
print(a)
print(type(a))
isinstance(a,float)

3.3
<class 'float'>

True

Para comprobar los atributos y los métodos que podemos aplicar sobre un objeto en particular, podemos emplear la función dir(), que nos devuelve una lista con sus nombres. Como se puede ver en la lista, podemos distinguir dos tipos de atributos y métodos: aquéllos que van con prefijo y sufijo __*__ y aquéllos que no. Los del primer tipo se denominan privados y habitualmente hacen referencia a los cálculos más básicos que se pueden realizar dentro de la clase a la que pertenece el objeto.

dir(a)

['__abs__',
 '__add__',
 '__bool__',
 '__ceil__',
 '__class__',
 '__delattr__',
 '__dir__',
 '__divmod__',
 '__doc__',
 '__eq__',
 '__float__',
 '__floor__',
 '__floordiv__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getformat__',
 '__getnewargs__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__int__',
 '__le__',
 '__lt__',
 '__mod__',
 '__mul__',
 '__ne__',
 '__neg__',
 '__new__',
 '__pos__',
 '__pow__',
 '__radd__',
 '__rdivmod__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__rfloordiv__',
 '__rmod__',
 '__rmul__',
 '__round__',
 '__rpow__',
 '__rsub__',
 '__rtruediv__',
 '__setattr__',
 '__setformat__',
 '__sizeof__',
 '__str__',
 '__sub__',
 '__subclasshook__',
 '__truediv__',
 '__trunc__',
 'as_integer_ratio',
 'conjugate',
 'fromhex',
 'hex',
 'imag',
 'is_integer',
 'real']

Por ejemplo, podemos comprobar si un número real es mayor que otro de dos maneras diferentes: utilizando el operador lógico > o bien usando su método privado __ge__:

a=6.4; b=5.3
print(a.__gt__(b))
print(a > b)

True
True

En cualquier caso, en Python tanto los métodos privados como los públicos se utilizan de la misma manera. Por ejemplo:

a=3.4
print(a.is_integer())
print(a.__int__())
print(np.int(a))

False
3
3

NOTA: En este curso no vamos a trabajar ni será necesaria la implementación de código utilizando una programación orientada a objetos. Pero lo que sí será necesario es tomar conciencia de cuándo se utilizan objetos de diferentes clases y qué atributos y métodos tiene definido cada caso.

5.3.3. Vectores de números#

En Python existen multitud de formas de guardar datos numéricos (o no) como pueden ser la estrutura lista o tupla. En particular, las listas pueden contener una secuencia finita de números ordenados (y utilizando un índice para acceder a cada uno de los elementos de la lista). Además, son lo suficientemente flexibles como para contener datos de diferente naturaleza (combinación de números enteros, reales, listas de listas, etc.).

Pero esta flexibilidad de las listas en Python hace que su rendimiento computacional sea muy limitado. En la mayoría de las aplicaciones científicas en matemáticas (con sus aplicaciones a las diferentes partes de la ingeniería informática, inteligencia artificial, ciencia de datos,…), los problemas reales necesitan operaciones sobre enormes conjuntos de datos y, por tanto, la velocidad computacional es muy importante para estos grandes problemas. Para trabajar de forma eficiente con estos problemas, Numpy proporciona funciones especializadas y estructuras de datos adecuadas. En particular, para el caso de conjuntos de números de un mismo tipo (perdiendo parte de la flexibilidad de las listas, pero ganando eficacia computacional).

Vectores unidimensionales: Un vector unidimensional es una colección ordenada de números a los que se puede acceder mediante un índice (con lo que se preserva el orden). Por defecto, los vectores en NumPy son vectores fila.

Creación de vectores e indexado: Para crear un vector Numpy de longitud 10 e inicializado con ceros, se utilizaría la función np.zeros():

u = np.zeros(10)
print(u)
print(type(u))

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
<class 'numpy.ndarray'>

El tipo por defecto de los números que contienen los vectores en NumPy es float64 (que es el tipo guardado en np.float). Si queremos usar otros tipos, tendríamos que utilizar el argumento opcional dtype. El tipo de los números que contiene un vector puede comprobarse accediendo al atributo dtype de los vectores NumPy:

print(u.dtype)
w = np.zeros(5, dtype=np.int)
print(w)
print(type(w))
print(w.dtype)

float64
[0 0 0 0 0]
<class 'numpy.ndarray'>
int32

Lo que no es posible, por ejemplo, es añadir un valor del tipo cadena de texto (es decir, de tipo string) a un objeto numpy.ndarray, ya que todos los elementos del vector deben ser del mismo tipo (o, al menos, de un tipo que admita una conversión) y deben tener también el mismo tamaño. Para comprobar el tamaño de un vector, se puede usar la función len:

print(len(u))
v = np.zeros(10, dtype=np.int)
print(u + v) # Implicitamente se hace una conversión de tipo de int64 a float64
print(u + w) # ERROR: ¡Los vectores no tienen el mismo tamaño!

10
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-6dfd7306e4fa> in <module>
      2 v = np.zeros(10, dtype=np.int)
      3 print(u + v) # Implicitamente se hace una conversión de tipo de int64 a float64
----> 4 print(u + w) # ERROR: ¡Los vectores no tienen el mismo tamaño!

ValueError: operands could not be broadcast together with shapes (10,) (5,)

Una forma más específica de comprobar la dimensión de un vector es usar u.shape, que nos devuelve una tupla con las dimensiones del vector:

print(u.shape)

(10,)

shape nos informa sobre el tamaño del vector en cada dirección. En el caso de los vectores, solamente hay una única dirección, mentres que en conjuntos de datos con múltiples índices (matrices, o tensores \(n\)-dimensionales), shape nos informaría del tamaño de estas estructuras de datos en cada dirección. Por ejemplo, para crear una matriz de ceros de tipo entero de tamaño \(2\times 3\):

A =  np.zeros((2,3), dtype=np.int)
print(A)
print(A.shape)

[[0 0 0]
 [0 0 0]]
(2, 3)

Podemos cambiar las entradas de un vector utilizando la indexación,

print(u)
u[0] = 10.0
u[3] = -4.3
u[9] = 1.0
print(u)

[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 10.    0.    0.   -4.3   0.    0.    0.    0.    0.    1. ]

NOTA IMPORTANTE: Hay que recordar que en Python los valores de los índices empiezan en 0.

Evidentemente, existen otras maneras de crear vectores, como, por ejemplo, el uso de la función ones para crear un vector que contenga solamente unos:

w = np.ones(5)
print(w)
print(w.dtype)

[ 1.  1.  1.  1.  1.]
float64

o un vector de valores aleatorios:

w = np.random.rand(6)
print(w)

[ 0.02920593  0.14178457  0.82590731  0.27716738  0.25108964  0.59544243]

o también un vector de números de tipo numpy.array a partir de una lista Python de números:

u = [4.0, 8.0, 9.0, 11.0, -2.0]
v = np.array(u)
print(v)

[  4.   8.   9.  11.  -2.]

Existen otros dos métodos para crear vectores de números que nos serán de utilidad a lo largo del curso (y, en particular, cuando tengamos que pintar funciones con Matplotlib):

  • numpy.arange

  • numpy.linspace

La función arange crea un vector con valores enteros consecutivos (de forma semejante a la función de Python range). Por ejemplo, para crear el vector fila \(\vec{u}=(0, 1, 2, 3, 4, 5)\) usando arange,

u = np.arange(6)
print(u)
print(type(u))
print(u.dtype)

[0 1 2 3 4 5]
<type 'numpy.ndarray'>
int64

Podemos comprobar que el número \(6\) no está incluido ya que el rango de valores comienza en \(0\) y el vector solamente posee seis elementos. Para cambiar el valor numérico en el que comienza el vector:

u = np.arange(2, 6)
print(u)

[2 3 4 5]

La función linspace crea un vector de números equiespaciados con un valor inicial y final (ambos incluidos) y con un tamaño determinado:

w = np.linspace(0., 100., 6)
print(w)
print(w.dtype)

[   0.   20.   40.   60.   80.  100.]
float64

Esta función linspace, junto con la función meshgrid, se usará de forma habitual para pintar funciones con Matplotlib a lo largo de este curso.

5.3.4. Funciones y aritmética sobre vectores#

Los vectores en NumPy soportan las operaciones aritméticas básicas (producto, sumas, restas, etc.):

a = np.array([1.0, 0.2, 1.2])
b = np.array([2.0, 0.1, 2.1])
print(a)
print(b)

# Suma de a e b
c = a + b
print(c)

[ 1.   0.2  1.2]
[ 2.   0.1  2.1]
[ 3.   0.3  3.3]

y el producto de todos sus elementos por un valor escalar,

c = 10.0*a
print(c)

[ 10.   2.  12.]

y elevar todas sus componentes a una potencia:

a = np.array([2, 3, 4])
print(a**2)

[ 4  9 16]

También se pueden aplicar las funciones de cálculo usual a cada una de sus componentes:

# Crear un vector [0, π/2, π, 3π/2]
a = np.array([0.0, np.pi/2, np.pi, 3*np.pi/2])
print(a)

# Calcular el seno de cada componente del vector
b = np.sin(a)
print(b)

[0.         1.57079633 3.14159265 4.71238898]
[ 0.0000000e+00  1.0000000e+00  1.2246468e-16 -1.0000000e+00]

El código anterior calcula el seno de cada coeficiente del vector a. Debemos remarcar que la función que se está usando es np.sin, que depende directamente del módulo NumPy. El uso de cualquier otra implementación de la función en otros módulos (por ejemplo, en el módulo Math), podría dar lugar a error.

Evidentemente, también podríamos calcular el seno de cada coeficiente del vector, accediendo a cada uno de los elementos mediante su índice y haciendo los cálculos en el interior de un bucle for:

b = np.zeros(len(a))
for i in range(len(a)):
    b[i] = np.sin(a[i])

print(b)

[  0.00000000e+00   1.00000000e+00   1.22464680e-16  -1.00000000e+00]

En este caso el programa es más largo y difícil de leer. Además, en muchos casos será más lento. La manipulación de vectores y cualquiera de los cálculos realizados entre ellos sin acceder a los índices suelen conocerse como vectorización. Cuando sea posible emplearla, la vectorización incrementará el rendimiento y la velocidad de los códigos de cálculo. En el ejercicio final de este guión, se analizará el rendimiento de este tipo de técnicas.

5.3.5. Troceado de vectores#

Cuando se trabaja con vectores de números, es habitual tener que extraer un subconjunto de estos para crear un nuevo vector. Por ejemplo, obtener los tres primeros coeficientes de un vector o, en el caso de matrices, restringir los cálculos a su segunda columna. Este tipo de operaciones es el que se denomina troceado de vectores (o, en inglés, array slicing).

Vamos a explorar esto mediante varios ejemplos, trabajando con un vector de valores aleatorios:

a = np.random.rand(5)
print(a)

# Usar ':' implica el conjunto entero en el rango de los índices, es decir, desde 0 hasta (longitud-1)
b = a[:]
print("Troceado usando '[:]' {}".format(b))

# Usar '1:3' implica los índices 1 -> 3 (sin incluir a 3)
b = a[1:3]
print("Troceado usando '[1:3]': {}".format(b))

# Usar '2:-1' implica los índices 2 -> el segundo desde el final (sin incluirlo)
b = a[2:-1]
print("Troceado usando '[2:-1]': {}".format(b))

# Usar '2:-2' implica los índices 2 -> el tercero desde el final (sin incluirlo)
b = a[2:-2]
print("Troceado usando '[2:-2]': {}".format(b))

[0.63076568 0.47536652 0.04092676 0.277812   0.82745402]
Troceado usando '[:]' [0.63076568 0.47536652 0.04092676 0.277812   0.82745402]
Troceado usando '[1:3]': [0.47536652 0.04092676]
Troceado usando '[2:-1]': [0.04092676 0.277812  ]
Troceado usando '[2:-2]': [0.04092676]

NOTA: el uso del índice -1, se corresponde con el último elemento del vector. Del mismo modo, el índice -2 corresponde al penúltimo elemento, etc.. Este convenio de referenciar índices desde el final de un vector es muy útil, ya que con el uso de estos índices negativos se puede hacer referencia a los últimos coeficientes de un vector sin tener que hacer referencia expresa a su tamaño (o, incluso, sin conocerlo explícitamente).

Si lo que se quiere es trocear un vector desde el principio o desde el final del mismo hay que utilizar la sintaxis de índices con ‘:’

# Usar ':3' implica usar índices desde el principio hasta 3 (sin incluir el índice 3)
b = a[:3]
print("Troceado usando '[:3]': {}".format(b))

# Usar '4:' implica los índices desde 4 -> hasta el final
b = a[4:]
print("Troceado usando '[4:]': {}".format(b))

# Usar ':' implica todos los índices desde el comienzo hasta el final
b = a[:]
print("Troceado usando '[:]': {}".format(b))

Troceado usando '[:3]': [0.65319139 0.7371771  0.85358493]
Troceado usando '[4:]': [0.78166875]
Troceado usando '[:]': [0.65319139 0.7371771  0.85358493 0.24875806 0.78166875]

El troceado también se puede aplicar a matrices:

B = np.array([[1.3, 0], [0, 2.0]])
print(B)

# Extraer la segunda fila
row = B[1, :]
print(row)

# Extraer la primera columna (almacenada en un vector fila)
col = B[:, 0] 
print(col)

[[1.3 0. ]
 [0.  2. ]]
[0. 2.]
[1.3 0. ]

Existen muchas otras estrategias y sintaxis relacionada con el troceado de vectores, que quedan fuera del alcance de esta introdución breve a NumPy. Para una información más detallada, se puede consultar: https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html

5.3.6. Ejercicio: aceleración con NumPy en el cálculo de normas de vectores#

La norma euclídea (o módulo) de un vector \(\vec{x}=(x_{0},\ldots,x_{n-1})\in\mathbb{R}^{n}\) viene dada por

\[ \| x \| = \sqrt{\sum_{i=0}^{n-1} x_{i} x_{i}}, \]

donde \(x_{i}\) es el \(i\)-ésimo coeficiente de \(\vec{x}\). En resumen, su norma no es más que el cálculo del producto interior del vector por sí mismo, seguido por el cálculo de una raíz cuadrada. Para calcular el valor de la norma se pueden seguir diferentes estrategias: la primera de ellas puede ser usar un bucle para recorrer todos los coeficientes del vector y sumar el cuadrado de cada coeficiente, luego tomar la suma de todas estas cantidades y calcular su raíz cadrada.

Para evaluar el rendimiento computacional de esta implementación, vamos a coger un vector moderadamente grande (10 millones) y vamos a calcular el tiempo de cáculo:

# Crear un vector NumPy con 10 millones de valores aleatorios
x = np.random.rand(10000000)
print(type(x))

<class 'numpy.ndarray'>

def compute_norm(x):
    norm = 0.0
    for xi in x:
        norm += xi*xi
    return np.sqrt(norm)

%time norm = compute_norm(x)
print(norm)

Wall time: 3.84 s
1825.756466436551

El tiempo de cálculo que nos interesa es el que aparece bajo la etiqueta ‘Wall time’.

NOTA: Los detalle de como la etiqueta %time trabaja en este guión no son de importancia en este curso. Debemos indicar que esta forma de proceder es únicamente válida para consolas o ficheiros que se ejecuten bajo [I]Python y entornos de servidores Jupyter. Esta forma de trabajar es muy compacta y conveniente para medir el tiempo que se requiere para completar la ejecución de una línea de código.

Ejercicio 1: Usando la misma implementación basada en un bucle, accede a cada elemento del vector por su índice, comenzando desde el primero hasta el último. Haz lo mismo con un bucle que acceda en orde inverso, desde el último hasta el primero.

# ESCRIBE AQUÍ TU CÓDIGO

Ejercicio 2: Trata de emplear funciones NumPy para definir una función que evite el bucle sobre los coeficientes del vector y que no acceda elemento a elemento al mismo (es decir, haz una versión vectorizada del cálculo de la norma euclídea).

# ESCRIBE AQUÍ TU CÓDIGO

Ejercicio 3: Compara los tiempos de cálculo de las cuatro implementaciones para diferentes dimensiones del vector \(\vec{x}\), de tamaños \(10\), \(10^3\), \(10^5\), \(10^7\). A partir de estos tiempos de cálculo: ¿qué se puede deducir como conclusión?

# ESCRIBE AQUÍ TU CÓDIGO