Álgebra Lineal en Python con NumPy (I): Operaciones básicas

Introducción

En esta entrada vamos a ver una introducción al Álgebra Lineal en Python con NumPy. En la mayoría de los artículos que hemos escrito hasta ahora en Pybonacci hemos tocado sin mencionarlos conceptos relativos al Álgebra Lineal: sin ir más lejos, el propio manejo de matrices o la norma de vectores. El lenguaje matricial es el punto de partida para una enorme variedad de desarrollos físicos y matemáticos, y por eso le vamos a dedicar un apartado especial para repasar las posibilidades que ofrece el paquete NumPy.

  1. Operaciones básicas
  2. Sistemas, autovalores y descomposiciones

En esta entrada se ha usado python 2.7.3 y numpy 1.6.1 y es compatible con python 3.2.3

Arrays y matrices

Como ya comentamos hace tiempo en nuestra introducción a Python, el paquete NumPy introdujo los arrays N-dimensionales, que no son más que colecciones homogéneas de elementos indexados usando N elementos. Los hemos utilizado constantemente usando las funciones de creación de arrays:

In [1]: import numpy as np
In [2]: np.array([1, 2, 3])  # Array de una lista
Out[2]: array([1, 2, 3])
In [3]: np.arange(5)  # Array de 5 enteros contando el 0
Out[3]: array([0, 1, 2, 3, 4])
In [4]: np.zeros((2, 3))  # Array de ceros de 2x3
Out[4]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
In [5]: np.linspace(0.0, 1.0, 5)  # Discretización de [0, 1] con 5 puntos
Out[5]: array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ])

Además de los arrays, con NumPy también podemos manejar matrices. Aunque parecen lo mismo, se utilizan de manera distinta; si alguien quiere investigar las diferencias, puede consultar la página NumPy para usuarios de MATLAB (en inglés).

Nota: Aunque la traducción al castellano de array sería precisamente vector, para evitar confusión voy a seguir utilizando la palabra inglesa, como he venido haciendo desde que se creó el blog.

Si queremos matrices, entonces hemos de usar la función matrix:

In [17]: v2 = np.matrix([0, 1, 2, 3])
In [18]: v2
Out[18]: matrix([[0, 1, 2, 3]])
In [19]: np.transpose(v2)
Out[19]:
matrix([[0],
        [1],
        [2],
        [3]])
In [20]: v2.T
Out[20]:
matrix([[0],
        [1],
        [2],
        [3]])
In [27]: np.matrix("""  # Para los nostálgicos
   ....: 1, 2;
   ....: 3, 4
   ....: """)
Out[27]:
matrix([[1, 2],
        [3, 4]])

Las ventajas de usar matrices en el fondo son muy pocas y además la mayoría de funciones de NumPy maneja arrays, así que tendrías que convertir entre ambos tipos constantemente. Hoy mostraremos brevemente el uso de las matrices también, pero mejor utilizar arrays y olvidarse 😉

Operaciones básicas, expansión

Suma

La suma de arrays y de matrices es igual: se realiza elemento a elemento. El producto por un escalar tampoco tiene misterio:

In [1]: import numpy as np
In [2]: a = np.array([[1, 2], [3, 4]])
In [3]: a + a
Out[3]:
array([[2, 4],
       [6, 8]])
In [4]: m = np.matrix([[1, 2], [3, 4]])
In [5]: m + m
Out[5]:
matrix([[2, 4],
        [6, 8]])
In [6]: a + m  # Podemos sumar arrays y matrices
Out[6]:
matrix([[2, 4],
        [6, 8]])
In [7]: 2 * a  # Producto por un escalar
Out[7]:
array([[2, 4],
       [6, 8]])
In [8]: -m
Out[8]:
matrix([[-1, -2],
        [-3, -4]])

Expansión o «broadcasting»

Si los objetos que estamos operando no tienen las mismas dimensiones, NumPy puede adaptar algunas de ellas para completar la operación. Esto se denomina broadcasting en inglés, y a falta de una traducción mejor y a menos que alguien tenga algo que objetar lo voy a denominar «expansión». Porque yo lo valgo.

NumPy alinea los arrays a la derecha y empieza a comprobar las dimensiones por el final. Si son todas compatibles realiza la operación, y si no lanza un error. Dos dimensiones son compatibles si

  1. Son iguales, o
  2. Una de ellas es 1.

Si los objetos no tienen el mismo número de dimensiones, se asume que las restantes son 1. Si dos dimensiones son distintas pero compatibles, la menor se expande hasta tener el tamaño de la mayor. Vamos a ver un par de ejemplos:

In [29]: c = np.zeros((3, 3))
In [30]: c
Out[30]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
In [31]: c.shape
Out[31]: (3, 3)
In [32]: n = np.array([1, 2, 3])
In [33]: n
Out[33]: array([1, 2, 3])
In [34]: n.shape
Out[34]: (3,)
In [35]: n + c  # Se expande la primera dimensión de n
Out[35]:
array([[ 1.,  2.,  3.],
       [ 1.,  2.,  3.],
       [ 1.,  2.,  3.]])

 

In [61]: xx = np.arange(4)
In [62]: yy = np.arange(4).reshape((4, 1))
In [63]: xx
Out[63]: array([0, 1, 2, 3])
In [64]: xx.shape
Out[64]: (4,)
In [65]: yy
Out[65]:
array([[0],
       [1],
       [2],
       [3]])
In [66]: yy.shape
Out[66]: (4, 1)
In [67]: xx + yy * 1j  # Se expanden la segunda dimensión de y y la primera de x
Out[67]:
array([[ 0.+0.j,  1.+0.j,  2.+0.j,  3.+0.j],
       [ 0.+1.j,  1.+1.j,  2.+1.j,  3.+1.j],
       [ 0.+2.j,  1.+2.j,  2.+2.j,  3.+2.j],
       [ 0.+3.j,  1.+3.j,  2.+3.j,  3.+3.j]])

Como podéis ver, esto abre un mundo de posibilidades. Si queréis ampliarlas, podéis consultar el material de donde he sacado esto en la página Broadcasting de la guía del usuario de NumPy.

Producto

En el producto es donde surgen las diferencias fundamentales entre arrays y matrices. Mientras que con los arrays el producto se hace elemento a elemento (y si las dimensiones no coinciden se expanden como hemos visto antes), el producto de matrices tiene, como ya sabemos, una definición totalmente distinta.

Para el producto matricial:

In [69]: a = np.array([[1, 0], [2, -1]])
In [70]: a
Out[70]:
array([[ 1,  0],
       [ 2, -1]])
In [71]: np.dot(a, a)
Out[71]:
array([[1, 0],
       [0, 1]])
In [72]: m = np.matrix('1 0; 2 -1')
In [73]: m
Out[73]:
matrix([[ 1,  0],
        [ 2, -1]])
In [74]: m * m
Out[74]:
matrix([[1, 0],
        [0, 1]])

Y para el producto elemento a elemento

In [80]: a
Out[80]:
array([[ 1,  0],
       [ 2, -1]])
In [81]: a * a
Out[81]:
array([[1, 0],
       [4, 1]])
In [82]: m
Out[82]:
matrix([[ 1,  0],
        [ 2, -1]])
In [83]: np.multiply(m, m)
Out[83]:
matrix([[1, 0],
        [4, 1]])

Nótese que con las matrices el producto «por defecto» es el producto tal y como lo conocemos en Álgebra Lineal.

Ya hemos utilizado la función dot, que realiza un producto escalar entre sus dos argumentos, en este caso generalizado para matrices. NumPy también ofrece la función inner, que no es más que producto interior. Matemáticamente, el producto escalar es un caso particular del interior; en NumPy, ambos son idénticos para objetos de dimensión menor o igual que 2. Para arrays de dimensiones mayores que 2, se diferencian en las dimensiones sobre las que se suma.

Con NumPy podemos calcular también el producto vectorial y el producto exterior de dos vectores, utilizando las funciones cross y outer respectivamente:

In [102]: u = np.arange(3)
In [103]: v = 1 + np.arange(3)  # ¡Expansión sobre el 1!
In [104]: u
Out[104]: array([0, 1, 2])
In [105]: v
Out[105]: array([1, 2, 3])
In [106]: np.cross(u, v)  # Producto vectorial
Out[106]: array([-1,  2, -1])
In [107]: np.outer(u, v)  # Producto exterior
Out[107]:
array([[0, 0, 0],
       [1, 2, 3],
       [2, 4, 6]])

Norma y determinante

NumPy ofrece también diversas funciones para calcular la norma y el determinante de un array. Para calcular el determinante de un array de dimensión 2 utilizamos la función linalg.det. También podemos conocer la traza con la función trace:

In [111]: a
Out[111]:
array([[ 1,  0],
       [ 2, -1]])
In [112]: np.linalg.det(a)
Out[112]: -1.0
In [113]: np.trace(a)
Out[113]: 0

Con NumPy podemos obtener también la norma de un array de cualquier dimensión y el número de condición de un array bidimensional. Para ello utilizamos las funciones linalg.norm y linalg.cond respectivamente, que aceptan un segundo argumento: el tipo de norma utilizado. Así, podemos conocer la norma-2 de un vector \(||v||_2\) o la norma infinito o del supremo de una matriz \(||A||_{infty}\). Para una lista completa, puedes consultar la documentación.

In [121]: a
Out[121]:
array([[ 1,  0],
       [ 2, -1]])
In [122]: np.linalg.norm(a)
Out[122]: 2.4494897427831779
In [123]: np.linalg.norm(a, np.inf)
Out[123]: 3
In [124]: np.linalg.norm(a, 1)
Out[124]: 3
In [125]: v
Out[125]: array([1, 2, 3])
In [126]: np.linalg.norm(v)
Out[126]: 3.7416573867739413
In [127]: np.linalg.norm(v, np.inf)
Out[127]: 3
In [128]: np.linalg.norm(v, 1)
Out[128]: 6
In [129]: np.linalg.norm(v, -np.inf)
Out[129]: 1
In [130]: np.linalg.cond(a)
Out[130]: 5.828427124746189
In [131]: np.linalg.cond(a, np.inf)
Out[131]: 9.0

Convenio de Einstein

Finalmente, y como regalo a los que hayan llegado hasta aquí abajo, os enseñamos que con NumPy podemos utilizar el convenio de sumación de Einstein, utilizando la función einsum. Para los que conozcan un poco de álgebra tensorial, esta función les va a encantar: con ella podemos hacer todo lo que hemos visto anteriormente y mucho más. En la documentación viene explicado detenidamente cómo utilizarla y hay unos cuantos ejemplos; vamos a ver algunos:

In [135]: a
Out[135]:
array([[ 1,  0],
       [ 2, -1]])
In [136]: np.trace(a)
Out[136]: 0
In [137]: np.einsum('ii', a)  # Traza
Out[137]: 0
In [138]: np.diag(a)
Out[138]: array([ 1, -1])
In [139]: np.einsum('ii -> i', a)  # Diagonal
Out[139]: array([ 1, -1])
In [140]: np.dot(a, a)
Out[140]:
array([[1, 0],
       [0, 1]])
In [142]: np.einsum('ij, jk', a, a)  # Producto matricial
Out[142]:
array([[1, 0],
       [0, 1]])
In [143]: u = np.array([1, -2])
In [144]: u
Out[144]: array([ 1, -2])
In [145]: np.dot(a, u)
Out[145]: array([1, 4])
In [146]: np.einsum('ij, j', a, u)  # Matriz por vector
Out[146]: array([1, 4])
In [147]: np.dot(u, u)
Out[147]: 5
In [151]: np.einsum('i, i', u, u)  # Producto escalar por sí mismo
Out[151]: 5
In [153]: np.outer(u, u)
Out[153]:
array([[ 1, -2],
       [-2,  4]])
In [154]: np.einsum('i, j', u, u)  # Producto exterior
Out[154]:
array([[ 1, -2],
       [-2,  4]])

Las posibilidades que esto ofrece son enormes. Puedes probar a jugar con tensores de mayor orden, como por ejemplo el símbolo de Levi-Civita.

https://gist.github.com/2689795

¡Y esto ha sido todo por hoy! Ha salido un artículo un poco largo porque he intentado mostrar el abanico de posibilidades que ofrece NumPy. Espero que te haya resultado interesante y que no olvides que puedes difundir el artículo, así como hacernos llegar tus comentarios.

¡Un saludo!

Referencias

Juan Luis Cano

Estudiante de ingeniería aeronáutica y con pasión por la programación y el software libre. Obsesionado con mejorar los pequeños detalles y con ganas de cambiar el mundo. Divulgando Python en español a través de Pybonacci y la asociación Python España.

More Posts - Website

Follow Me:
TwitterLinkedIn

29 thoughts on “Álgebra Lineal en Python con NumPy (I): Operaciones básicas

  1. Para ser sincero, después de unos años de usar numpy, creo que nunca he usado ‘matrix’ para nada 🙁

    Por otra parte, muchas gracias por este repaso. Después de muchos años de no usar algunas cosas no está mal poder recordarlas.

    Saludos.

    1. Ya, yo tampoco le veo mucho sentido a matrix, pero por lo menos había que mencionarlo por si acaso.

      Gracias a ti Kiko! Yo tengo más de lo que aprender que de lo que enseñar, me parece a mí 😛

      Un saludo!

  2. Ojo con traducir array como vector. Hacéis bien en no hacerlo; mejor dejarlo en inglés, o como “arreglo” o disposición.

  3. Muchas gracias, muy interesante las publicaciones. Por favor si me podrían ayudar a desarrollar el Algoritmo de Thomas en diferencias finitas.. les estaré muy agradecido

    1. Hola lincol, ¡gracias a ti por el comentario! Por si necesitas inspiración, buscando “thomas algorithm python” en Google he encontrado unos cuantos resultados 🙂 ¡Un saludo!

  4. Sabes cómo definir una función de error en scikit-learn? Llevo un montón de tiempo y no consigo que funcione me he leído toda la documentación. Y tengo la duda colgada en Stake… sin responder una sola vez.

    1. Hola Daniel, buscando en Google he encontrado tu pregunta de Stack Overflow así que he puesto una recompensa de 50 puntos, a ver si hay suerte 🙂 ¡Un saludo!

  5. Hola. Me han contestado y he solucionado el problema! Muchas gracias por lo de la recompensa. No sé cómo funciona, pero creo que yo no podía hacerlo. En cualquier caso, gracias!

    1. ¡Genial! Supongo que no porque la recompensa mínima son 50 puntos que se te descuentan de tu reputación y tú no tenías suficientes. Me alegro de que por fin te lo hayan aclarado, marca una de las dos respuestas como buena y a esa le doy la recompensa 🙂 ¡Saludos!

    1. Daniel, tienes que marcar una respuesta como buena haciendo clic en el “tick” verde que aparece – de esa forma indicas que una de las respuestas resolvió tu problema (lo cual no quita que votes positivo la otra también). ¡Un saludo!

  6. Hola Juan, tenes idea como hago para anexar listas a una variable numpy.ndarray()? hay algo parecido a x.append(y) de las listas comunes? dando vueltas por internet no me queda claro

    Saludos

    Gabriel

  7. Hola buenos días, estoy empezando en NumPy y querría saber como se hace un producto escalar de dos vectores, gracias

Comments are closed.