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:
Puede que ya sepas que en la biblioteca SciPy se pueden encontrar también funciones de Álgebra Lineal. ¿Cuáles usar? Puedes encontrar la respuesta en este enlace: https://docs.scipy.org/doc/scipy-0.18.1/reference/tutorial/linalg.html#scipy-linalg-vs-numpy-linalg
Como de momento sólo hemos usado NumPy no importaremos las funciones de SciPy, aunque como ves, es recomendable hacerlo.
In [1]:
import numpy as np
In [2]:
help(np.linalg)
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:
In [3]:
from numpy.linalg import norm, det
norm
Out[3]:
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.
In [4]:
np.dot
Out[4]:
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.
In [5]:
M = np.array([
[1, 2],
[3, 4]
])
v = np.array([1, -1])
In [6]:
v.T
Out[6]:
In [7]:
u = np.dot(M, v)
u
Out[7]:
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.
In [8]:
u, v
Out[8]:
In [9]:
np.allclose(u, v)
Out[9]:
In [10]:
np.isclose(0.0, 1e-8, atol=1e-10)
Out[10]:
En la versión 3.5 de Python se incorporó un nuevo operador @ para poder calcular hacer multiplicaciones de matrices de una forma más legible
In [11]:
u = M @ v
u
Out[11]:
In [12]:
mat = np.array([[1, 5, 8, 5],
[0, 6, 4, 2],
[9, 3, 1, 6]])
vec1 = np.array([5, 6, 2])
vec1 @ mat
Out[12]:
Si quieres saber más, puedes leer este artículo en Pybonacci escrito por Álex Sáez.
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}$$
In [13]:
from numpy.linalg import det
In [14]:
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)
In [15]:
C = A @ B
C
Out[15]:
In [16]:
det(C)
Out[16]:
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} $$
In [17]:
M = (np.array([[2, 0, 0],
[-1, 1, 0],
[3, 2, -1]])
@
np.array([[1, 1, 1],
[0, 1, 2],
[0, 0, 1]]))
M
Out[17]:
In [18]:
x = np.linalg.solve(M, np.array([-1, 3, 0]))
x
Out[18]:
In [19]:
np.allclose(M @ x, np.array([-1, 3, 0]))
Out[19]:
3- Hallar la inversa de la matriz $H$ y comprobar que $H H^{-1} = I$ (recuerda la función np.eye)
In [25]:
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);
In [30]:
np.linalg.det(H)
Out[30]:
In [31]:
Hinv = np.linalg.inv(H)
In [32]:
np.isclose(np.dot(Hinv, H), np.eye(6))
Out[32]:
In [33]:
np.set_printoptions(precision=3)
print(np.dot(Hinv, H))
4- Comprobar el número de condición de la matriz $H$.
In [34]:
np.linalg.cond(H)
Out[34]:
La cosa no queda ahí y también se pueden resolver problemas de autovalores y autovectores:
In [35]:
A = np.array([
[1, 0, 0],
[2, 1, 1],
[-1, 0, 1]
])
np.linalg.eig(A)
Out[35]:
Ya hemos aprendido a efectuar algunas operaciones útiles con NumPy. Estamos en condiciones de empezar a escribir programas más interesantes, pero aún queda lo mejor.