Álgebra Lineal con NumPy¶

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.

Álgebra lineal¶

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:

Funciones básicas (norma de un vector, inversa de una matriz, determinante, traza)
Resolución de sistemas
Autovalores y autovectores
Descomposiciones matriciales (QR, SVD)
Pseudoinversas

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)

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/.miniconda36/envs/aeropython36/lib/python3.6/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:

In [3]:

from numpy.linalg import norm, det
norm

Out[3]:

<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.

In [4]:

np.dot

Out[4]:

<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.

In [5]:

M = np.array([
    [1, 2],
    [3, 4]
])
v = np.array([1, -1])

In [6]:

v.T

Out[6]:

array([ 1, -1])

In [7]:

u = np.dot(M, v)
u

Out[7]:

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.

In [8]:

u, v

Out[8]:

(array([-1, -1]), array([ 1, -1]))

In [9]:

np.allclose(u, v)

Out[9]:

False

In [10]:

np.isclose(0.0, 1e-8, atol=1e-10)

Out[10]:

False

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]:

array([-1, -1])

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]:

array([23, 67, 66, 49])

Si quieres saber más, puedes leer este artículo en Pybonacci escrito por Álex Sáez.

Ejercicios¶

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)

[[ 1  0  0]
 [ 2  1  1]
 [-1  0  1]]
[[ 2  3 -1]
 [ 0 -2  1]
 [ 0  0  3]]

In [15]:

C = A @ B
C

Out[15]:

array([[ 2,  3, -1],
       [ 4,  4,  2],
       [-2, -3,  4]])

In [16]:

det(C)

Out[16]:

-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} $$

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]:

array([[ 2,  2,  2],
       [-1,  0,  1],
       [ 3,  5,  6]])

In [18]:

x = np.linalg.solve(M, np.array([-1, 3, 0]))
x

Out[18]:

array([ 0.5, -4.5,  3.5])

In [19]:

np.allclose(M @ x, np.array([-1, 3, 0]))

Out[19]:

True

3- Hallar la inversa de la matriz $H$ y comprobar que $H H^{-1} = I$ (recuerda la función np.eye)

In [20]:

# aeropython: preserve
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]]

In [21]:

np.linalg.det(H)

Out[21]:

456427.51147512451

In [22]:

Hinv = np.linalg.inv(H)

In [23]:

np.isclose(np.dot(Hinv, H), np.eye(6))

Out[23]:

array([[False, False, False, False, False, False],
       [False, False, False, False, False, False],
       [False, False, False, False, False, False],
       [False, False,  True,  True, False, False],
       [False, False, False, False, False, False],
       [False, False, False, False, False, False]], dtype=bool)

In [24]:

np.set_printoptions(precision=3)
print(np.dot(Hinv, H))

[[  0.562  -1.     -2.25   -4.5    -8.    -16.   ]
 [ -0.094   0.875  -0.25   -0.5    -1.     -3.   ]
 [  0.125   0.5     2.5     3.      4.      8.   ]
 [  0.062   0.125   0.      1.      1.      2.   ]
 [ -0.062  -0.25   -0.5    -1.     -1.     -2.   ]
 [  0.031   0.062   0.25    0.5     0.5     2.   ]]

¡No funciona! Y no solo eso sino que los resultados varían de un ordenador a otro.

4- Comprobar el número de condición de la matriz $H$.

In [25]:

np.linalg.cond(H)

Out[25]:

1.0682988875862586e+17

La matriz está mal condicionada.

Autovalores y autovectores¶

La cosa no queda ahí y también se pueden resolver problemas de autovalores y autovectores:

In [26]:

A = np.array([
    [1, 0, 0],
    [2, 1, 1],
    [-1, 0, 1]
])

np.linalg.eig(A)

Out[26]:

(array([ 1.,  1.,  1.]), array([[  0.000e+00,   0.000e+00,   4.930e-32],
        [  1.000e+00,  -1.000e+00,  -1.000e+00],
        [  0.000e+00,   2.220e-16,   2.220e-16]]))

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.

Álgebra lineal en Python con NumPy en Pybonacci

####

¡Síguenos en Twitter!
###### Follow @AeroPython
###### Este notebook ha sido realizado por: Juan Luis Cano, Mabel Delgado y Álex Sáez
#####
Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.

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

In [27]:

# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/aeropython.css'
HTML(open(css_file, "r").read())

Out[27]:

/* This template is inspired in the one used by Lorena Barba in the numerical-mooc repository: https://github.com/numerical-mooc/numerical-mooc We thank her work and hope you also enjoy the look of the notobooks with this style */ El estilo se ha aplicado =)