Una vez hemos visto el manejo básico de arrays en Python con NumPy, es hora de pasar a operaciones más interesantes como son las propias del Álgebra Lineal.
Los productos escalares y las inversiones de matrices están por todas partes en los programas científicos e ingenieriles, así que vamos a estudiar cómo se realizan en Python.
Como sabemos, las operaciones del álgebra lineal aparecen con mucha frecuencia a la hora de resolver sistemas de ecuaciones en derivadas parciales y en general al linealizar problemas de todo tipo, y suele ser necesario resolver sistemas con un número enorme de ecuaciones e incógnitas. Gracias a los arrays de NumPy podemos abordar este tipo de cálculos en Python, ya que todas las funciones están escritas en C o Fortran y tenemos la opción de usar bibliotecas optimizadas al límite.
El paquete de álgebra lineal en NumPy se llama linalg
, así que importando NumPy con la convención habitual podemos acceder a él escribiendo np.linalg
. Si imprimimos la ayuda del paquete vemos que tenemos funciones para:
import numpy as np
help(np.linalg)
Help on package numpy.linalg in numpy: NAME numpy.linalg DESCRIPTION Core Linear Algebra Tools ------------------------- Linear algebra basics: - norm Vector or matrix norm - inv Inverse of a square matrix - solve Solve a linear system of equations - det Determinant of a square matrix - lstsq Solve linear least-squares problem - pinv Pseudo-inverse (Moore-Penrose) calculated using a singular value decomposition - matrix_power Integer power of a square matrix Eigenvalues and decompositions: - eig Eigenvalues and vectors of a square matrix - eigh Eigenvalues and eigenvectors of a Hermitian matrix - eigvals Eigenvalues of a square matrix - eigvalsh Eigenvalues of a Hermitian matrix - qr QR decomposition of a matrix - svd Singular value decomposition of a matrix - cholesky Cholesky decomposition of a matrix Tensor operations: - tensorsolve Solve a linear tensor equation - tensorinv Calculate an inverse of a tensor Exceptions: - LinAlgError Indicates a failed linear algebra operation PACKAGE CONTENTS _umath_linalg info lapack_lite linalg setup DATA absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0... division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192... print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)... FILE /home/juanlu/.local/lib/python3.4/site-packages/numpy/linalg/__init__.py
Recordemos que si queremos usar una función de un paquete pero no queremos escribir la "ruta" completa cada vez, podemos usar la sintaxis from package import func
:
from numpy.linalg import norm, det
norm
<function numpy.linalg.linalg.norm>
El producto matricial usual (no el que se hace elemento a elemento, sino el del álgebra lineal) se calcula con la misma función que el producto matriz-vector y el producto escalar vector-vector: con la función dot
, que no está en el paquete linalg
sino directamente en numpy
y no hace falta importarlo.
np.dot
<function numpy.core.multiarray.dot>
Una consideración importante a tener en cuenta es que en NumPy no hace falta ser estricto a la hora de manejar vectores como si fueran matrices columna, siempre que la operación sea consistente. Un vector es una matriz con una sola dimensión: por eso si calculamos su traspuesta no funciona.
M = np.array([
[1, 2],
[3, 4]
])
v = np.array([1, -1])
v.T
array([ 1, -1])
u = np.dot(M, v)
u
array([-1, -1])
Para hacer comparaciones entre arrays de punto flotante se pueden usar las funciones np.allclose
y np.isclose
. La primera comprueba si todos los elementos de los arrays son iguales dentro de una tolerancia, y la segunda compara elemento a elemento y devuelve un array de valores True
y False
.
u, v
(array([-1, -1]), array([ 1, -1]))
np.allclose(u, v)
False
np.isclose(0.0, 1e-8, atol=1e-10)
False
1- Hallar el producto de estas dos matrices y su determinante:
$$\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 1 \\ -1 & 0 & 1 \end{pmatrix} \begin{pmatrix} 2 & 3 & -1 \\ 0 & -2 & 1 \\ 0 & 0 & 3 \end{pmatrix}$$from numpy.linalg import det
A = np.array([
[1, 0, 0],
[2, 1, 1],
[-1, 0, 1]
])
B = np.array([
[2, 3, -1],
[0, -2, 1],
[0, 0, 3]
])
print(A)
print(B)
[[ 1 0 0] [ 2 1 1] [-1 0 1]] [[ 2 3 -1] [ 0 -2 1] [ 0 0 3]]
C = np.dot(A, B)
C
array([[ 2, 3, -1], [ 4, 4, 2], [-2, -3, 4]])
det(C)
-12.0
2- Resolver el siguiente sistema:
$$ \begin{pmatrix} 2 & 0 & 0 \\ -1 & 1 & 0 \\ 3 & 2 & -1 \end{pmatrix} \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} -1 \\ 3 \\ 0 \end{pmatrix} $$M = np.dot(
np.array([
[2, 0, 0],
[-1, 1, 0],
[3, 2, -1]
]),
np.array([
[1, 1, 1],
[0, 1, 2],
[0, 0, 1]
]))
M
array([[ 2, 2, 2], [-1, 0, 1], [ 3, 5, 6]])
x = np.linalg.solve(M, np.array([-1, 3, 0]))
x
array([ 0.5, -4.5, 3.5])
np.allclose(np.dot(M, x), np.array([-1, 3, 0]))
True
3- Hallar la inversa de la matriz $H$ y comprobar que $H H^{-1} = I$ (recuerda la función np.eye
)
A = np.arange(1, 37).reshape(6,6)
A[1, 1::2] = 0
A[3, ::2] = 1
A[4, :] += 30
B = (2 ** np.arange(36)).reshape((6,6))
H = A + B
print(H)
[[ 2 4 7 12 21 38] [ 71 128 265 512 1035 2048] [ 4109 8206 16399 32784 65553 131090] [ 262145 524308 1048577 2097174 4194305 8388632] [ 16777271 33554488 67108921 134217786 268435515 536870972] [ 1073741855 2147483680 4294967329 8589934626 17179869219 34359738404]]
np.linalg.det(H)
456269.59150316147
Hinv = np.linalg.inv(H)
np.isclose(np.dot(Hinv, H), np.eye(6))
array([[False, False, False, False, False, False], [False, False, False, False, False, False], [False, False, False, False, False, False], [False, False, False, False, False, False], [False, True, False, False, True, False], [False, True, False, True, True, False]], dtype=bool)
np.set_printoptions(precision=3)
print(np.dot(Hinv, H))
[[ 2.188e+00 2.375e+00 5.000e+00 9.000e+00 1.900e+01 4.000e+01] [ 1.875e-01 1.438e+00 8.750e-01 1.750e+00 3.500e+00 6.000e+00] [ -5.000e-01 -5.000e-01 -5.000e-01 -2.000e+00 -6.000e+00 -1.600e+01] [ 1.250e-01 2.500e-01 5.000e-01 1.500e+00 2.000e+00 6.000e+00] [ 6.250e-02 0.000e+00 2.500e-01 -5.000e-01 1.000e+00 2.000e+00] [ -3.125e-02 0.000e+00 -1.250e-01 0.000e+00 0.000e+00 0.000e+00]]
4- Comprobar el número de condición de la matriz $H$.
np.linalg.cond(H)
1.0682988820367766e+17
Hasta ahora hemos visto cómo acceder a elementos aislados del array, pero la potencia de NumPy está en poder acceder a secciones enteras. Para ello se usa la sintaxis inicio:final:paso
: si alguno de estos valores no se pone toma un valor por defecto. Veamos ejemplos:
M = np.arange(36, dtype=float).reshape(4, 9)
M
array([[ 0., 1., 2., 3., 4., 5., 6., 7., 8.], [ 9., 10., 11., 12., 13., 14., 15., 16., 17.], [ 18., 19., 20., 21., 22., 23., 24., 25., 26.], [ 27., 28., 29., 30., 31., 32., 33., 34., 35.]])
# De la segunda a la tercera fila, incluida
M[1:3]
array([[ 9., 10., 11., 12., 13., 14., 15., 16., 17.], [ 18., 19., 20., 21., 22., 23., 24., 25., 26.]])
# Hasta la tercera fila sin incluir y de la segunda a la quinta columnas saltando dos
M[:2, 1:5:2]
#M[1:2:1, 1:5:2] # Equivalente
array([[ 1., 3.], [ 10., 12.]])
Pintar un tablero de ajedrez usando la función plt.matshow
.
tablero = np.zeros([8, 8], dtype=int)
tablero[0::2, 1::2] = 1
tablero[1::2, 0::2] = 1
tablero
array([[0, 1, 0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 1, 0, 1, 0], [0, 1, 0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 1, 0, 1, 0], [0, 1, 0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 1, 0, 1, 0], [0, 1, 0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 1, 0, 1, 0]])
%matplotlib inline
import matplotlib.pyplot as plt
plt.matshow(tablero, cmap=plt.cm.gray_r)
<matplotlib.image.AxesImage at 0x7fa7aa64a908>
Con E/S (I/O en inglés) entendemos leer y escribir datos archivos. Es algo que necesitaremos hacer con relativa frecuencia, y en NumPy es muy sencillo de hacer. Para el caso de la lectura se usa la función np.loadtxt
.
Para practicar, vamos a leer el archivo temperaturas.csv
que contiene datos diarios de temperaturas en Nueva York entre el 1 de enero de 2013 y el 1 de enero de 2014, obtenidos gratuitamente de http://ncdc.noaa.gov/. Como los hemos descargado en formato CSV habrá que tener algunas precauciones a la hora de leer el archivo.
!head ../static/temperaturas.csv # Esta línea no funciona en Windows
STATION,DATE,TMAX,TMIN GHCND:USW00094728,20130101,44,-33 GHCND:USW00094728,20130102,6,-56 GHCND:USW00094728,20130103,0,-44 GHCND:USW00094728,20130104,28,-11 GHCND:USW00094728,20130105,56,0 GHCND:USW00094728,20130106,78,11 GHCND:USW00094728,20130107,72,28 GHCND:USW00094728,20130108,89,17 GHCND:USW00094728,20130109,94,39
datos = np.loadtxt("../static/temperaturas.csv",
skiprows=1, # Saltamos una línea
usecols=(1, 2, 3), # Solo columnas 2, 3 y 4
delimiter=',') # Separados por comas
datos[:9]
array([[ 2.013e+07, 4.400e+01, -3.300e+01], [ 2.013e+07, 6.000e+00, -5.600e+01], [ 2.013e+07, 0.000e+00, -4.400e+01], [ 2.013e+07, 2.800e+01, -1.100e+01], [ 2.013e+07, 5.600e+01, 0.000e+00], [ 2.013e+07, 7.800e+01, 1.100e+01], [ 2.013e+07, 7.200e+01, 2.800e+01], [ 2.013e+07, 8.900e+01, 1.700e+01], [ 2.013e+07, 9.400e+01, 3.900e+01]])
La primera columna es un entero con formato "AAAAMMDD" que vamos a ignorar. Las temperaturas están medidas en décimas de grado Celsius, así que hay que pasarlas a grados Celsius. Vamos a calcular también la temperatura media.
Tmax = datos[:, 1] / 10
Tmin = datos[:, 2] / 10
Tavg = (Tmax + Tmin) / 2
Como vamos a ignorar la columna de las fechas tenemos que crear un dominio para el eje x. Simplemente construiremos un array de enteros desde 0 hasta 365.
x = np.arange(366)
Y ahora representamos la evolución de la temperatura media (por ejemplo de color negro), indicando "Daily summaries" en el título, "Days" en el eje x y "Temperature (C)" en el eje y, usando la interfaz orientada a objetos de matplotlib (función plt.subplots
). Podemos crear una zona rellena entre la máxima y la mínima con la función fill_between(x, max, min)
(por ejemplo de color #4f88b1). Si los límites del eje x no quedan como queremos podemos usar la función set_xlim(xmin, xmax)
.
fig, ax = plt.subplots()
ax.plot(x, Tavg, 'k')
ax.set_xlim(0, 366)
ax.fill_between(x, Tmin, Tmax, facecolor='#4f88b1', edgecolor='none')
ax.set_title("Resúmenes diarios")
ax.set_xlabel("Días")
ax.set_ylabel("Temperatura (C)")
<matplotlib.text.Text at 0x7fa7aa5a0240>
Ya hemos aprendido a efectuar algunas operaciones útiles con NumPy e incluso hemos hecho nuestro primer ejercicio de lectura de datos. Estamos en condiciones de empezar a escribir programas más interesantes, pero aún queda lo mejor.
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
%%html
<a href="https://twitter.com/Pybonacci" class="twitter-follow-button" data-show-count="false">Follow @Pybonacci</a>
<script>!function(d,s,id){var js,fjs=d.getElementsByTagName(s)[0],p=/^http:/.test(d.location)?'http':'https';if(!d.getElementById(id)){js=d.createElement(s);js.id=id;js.src=p+'://platform.twitter.com/widgets.js';fjs.parentNode.insertBefore(js,fjs);}}(document, 'script', 'twitter-wjs');</script>
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../static/styles/style.css'
HTML(open(css_file, "r").read())