Algebra Lineal con Python

Esta notebook fue creada originalmente como un blog post por Raúl E. López Briega en Mi blog sobre Python. El contenido esta bajo la licencia BSD.

Introducción#

Una de las herramientas matemáticas más utilizadas en machine learning y data mining es el Álgebra lineal; por tanto, si queremos incursionar en el fascinante mundo del aprendizaje automático y el análisis de datos es importante reforzar los conceptos que forman parte de sus cimientos.

El Álgebra lineal es una rama de las matemáticas que es sumamente utilizada en el estudio de una gran variedad de ciencias, como ser, ingeniería, finanzas, investigación operativa, entre otras. Es una extensión del álgebra que aprendemos en la escuela secundaria, hacia un mayor número de dimensiones; en lugar de trabajar con incógnitas a nivel de escalares comenzamos a trabajar con matrices y vectores.

El estudio del Álgebra lineal implica trabajar con varios objetos matemáticos, como ser:

• Los Escalares: Un escalar es un solo número, en contraste con la mayoría de los otros objetos estudiados en Álgebra lineal, que son generalmente una colección de múltiples números.

• Los Vectores:Un vector es una serie de números. Los números tienen una orden preestablecido, y podemos identificar cada número individual por su índice en ese orden. Podemos pensar en los vectores como la identificación de puntos en el espacio, con cada elemento que da la coordenada a lo largo de un eje diferente. Existen dos tipos de vectores, los vectores de fila y los vectores de columna. Podemos representarlos de la siguiente manera, dónde f es un vector de fila y c es un vector de columna:

  $$f=\begin{bmatrix}0&1&-1\end{bmatrix} ; c=\begin{bmatrix}0\\1\\-1\end{bmatrix}$$

• Las Matrices: Una matriz es un arreglo bidimensional de números (llamados entradas de la matriz) ordenados en filas (o renglones) y columnas, donde una fila es cada una de las líneas horizontales de la matriz y una columna es cada una de las líneas verticales. En una matriz cada elemento puede ser identificado utilizando dos índices, uno para la fila y otro para la columna en que se encuentra. Las podemos representar de la siguiente manera, A es una matriz de 3x2.

• Los Tensores:En algunos casos necesitaremos una matriz con más de dos ejes. En general, una serie de números dispuestos en una cuadrícula regular con un número variable de ejes es conocido como un tensor.

Sobre estos objetos podemos realizar las operaciones matemáticas básicas, como ser adición, multiplicación, sustracción y división, es decir que vamos a poder sumar vectores con matrices, multiplicar escalares a vectores y demás.

Librerías de Python para álgebra lineal#

Los principales módulos que Python nos ofrece para realizar operaciones de Álgebra lineal son los siguientes:

• Numpy: El popular paquete matemático de Python, nos va a permitir crear vectores, matrices y tensores con suma facilidad.

• numpy.linalg: Este es un submodulo dentro de Numpy con un gran número de funciones para resolver ecuaciones de Álgebra lineal.

• scipy.linalg: Este submodulo del paquete científico Scipy es muy similar al anterior, pero con algunas más funciones y optimaciones.

• Sympy: Esta librería nos permite trabajar con matemática simbólica, convierte a Python en un sistema algebraico computacional. Nos va a permitir trabajar con ecuaciones y fórmulas simbólicamente, en lugar de numéricamente.

• CVXOPT: Este módulo nos permite resolver problemas de optimizaciones de programación lineal.

• PuLP: Esta librería nos permite crear modelos de programación lineal en forma muy sencilla con Python.

Operaciones básicas#

Vectores#

Un vector de largo n es una secuencia (o array, o tupla) de n números. La solemos escribir como \(x=(x1,…,xn)\) o \(x=[x1,…,xn]\)

En Python, un vector puede ser representado con una simple lista, o con un array de Numpy; siendo preferible utilizar esta última opción.

# Vector como lista de Python
v1 = [2, 4, 6]
v1

[2, 4, 6]

# Vectores con numpy
import numpy as np

v2 = np.ones(3) # vector de solo unos.
v2

array([1., 1., 1.])

v3 = np.array([1, 3, 5]) # pasando una lista a las arrays de numpy
v3

array([1, 3, 5])

v4 = np.arange(1, 8) # utilizando la funcion arange de numpy
v4

array([1, 2, 3, 4, 5, 6, 7])

Representación gráfica#

Tradicionalmente, los vectores son representados visualmente como flechas que parten desde el origen hacia un punto.

Por ejemplo, si quisiéramos representar graficamente a los vectores \(v1=[2, 4]\), \(v2=[-3, 3]\) y \(v3=[-4, -3.5]\), podríamos hacerlo de la siguiente manera.

import matplotlib.pyplot as plt
from warnings import filterwarnings

%matplotlib inline
filterwarnings('ignore') # Ignorar warnings

def move_spines():
    """Crea la figura de pyplot y los ejes. Mueve las lineas de la izquierda y de abajo
    para que se intersecten con el origen. Elimina las lineas de la derecha y la de arriba.
    Devuelve los ejes."""
    fix, ax = plt.subplots()
    for spine in ["left", "bottom"]:
        ax.spines[spine].set_position("zero")
    
    for spine in ["right", "top"]:
        ax.spines[spine].set_color("none")
    
    return ax

def vect_fig(): 
    """Genera el grafico de los vectores en el plano"""
    ax = move_spines()
    
    ax.set_xlim(-5, 5)
    ax.set_ylim(-5, 5)
    ax.grid()
    vecs = [[2, 4], [-3, 3], [-4, -3.5]] # lista de vectores
    for v in vecs:
        ax.annotate(" ", xy=v, xytext=[0, 0],
                   arrowprops=dict(facecolor="blue",
                                  shrink=0,
                                  alpha=0.7,
                                  width=0.5))
        ax.text(1.1 * v[0], 1.1 * v[1], v)

vect_fig() # crea el gráfico

Operaciones con vectores#

Las operaciones más comunes que utilizamos cuando trabajamos con vectores son la suma, la resta y la multiplicación por escalares.

Cuando sumamos dos vectores, vamos sumando elemento por elemento de cada vector.

De forma similar funciona la operación de resta.

La multiplicación por escalares es una operación que toma a un número \(\gamma\), y a un vector \(x\) y produce un nuevo vector donde cada elemento del vector \(x\) es multiplicado por el número \(\gamma\).

En Python podríamos realizar estas operaciones en forma muy sencilla:

# Ejemplo en Python
x = np.arange(1, 5)
y = np.array([2, 4, 6, 8])
x, y

(array([1, 2, 3, 4]), array([2, 4, 6, 8]))

# sumando dos vectores numpy
x + y

array([ 3,  6,  9, 12])

# restando dos vectores
x - y

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

# multiplicando por un escalar
x * 2

array([2, 4, 6, 8])

y * 3

array([ 6, 12, 18, 24])

Producto escalar o interior#

El producto escalar de dos vectores se define como la suma de los productos de sus elementos, suele representarse matemáticamente como < x, y > o x’y, donde x e y son dos vectores.

Dos vectores son ortogonales o perpendiculares cuando forman ángulo recto entre sí. Si el producto escalar de dos vectores es cero, ambos vectores son ortogonales.

Adicionalmente, todo producto escalar induce una norma sobre el espacio en el que está definido, de la siguiente manera:

En Python lo podemos calcular de la siguiente forma:

# Calculando el producto escalar de los vectores x e y
x @ y

60

# o lo que es lo mismo, que:
sum(x * y), np.dot(x, y)

(60, 60)

# Calculando la norma del vector X
np.linalg.norm(x)

5.477225575051661

# otra forma de calcular la norma de x
np.sqrt(x @ x)

5.477225575051661

# vectores ortogonales
v1 = np.array([3, 4])
v2 = np.array([4, -3])

v1 @ v2

0

Matrices#

Las matrices son una forma clara y sencilla de organizar los datos para su uso en operaciones lineales.

Una matriz n × k es una agrupación rectangular de números con n filas y k columnas; se representa de la siguiente forma:

En la matriz A, el símbolo \(a_{nk}\) representa el elemento n-ésimo de la fila en la k-ésima columna. La matriz A también puede ser llamada un vector si cualquiera de n o k son iguales a 1. En el caso de n=1, A se llama un vector fila, mientras que en el caso de k=1 se denomina un vector columna.

Las matrices se utilizan para múltiples aplicaciones y sirven, en particular, para representar los coeficientes de los sistemas de ecuaciones lineales o para representar transformaciones lineales dada una base. Pueden sumarse, multiplicarse y descomponerse de varias formas.

Operaciones con matrices#

Al igual que con los vectores, que no son más que un caso particular, las matrices se pueden sumar, restar y la multiplicar por escalares.

Multiplicacion por escalares:

Suma de matrices:

Resta de matrices:

Para los casos de suma y resta, hay que tener en cuenta que solo se pueden sumar o restar matrices que tengan las mismas dimensiones, es decir que si tengo una matriz A de dimensión 3x2 (3 filas y 2 columnas) solo voy a poder sumar o restar la matriz B si esta también tiene 3 filas y 2 columnas.

# Ejemplo en Python
A = np.array([[1, 3, 2],
              [1, 0, 0],
              [1, 2, 2]])

B = np.array([[1, 0, 5],
              [7, 5, 0],
              [2, 1, 1]])

# suma de las matrices A y B
A + B

array([[2, 3, 7],
       [8, 5, 0],
       [3, 3, 3]])

# resta de matrices
A - B

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

# multiplicando matrices por escalares
A * 2

array([[2, 6, 4],
       [2, 0, 0],
       [2, 4, 4]])

B * 3

array([[ 3,  0, 15],
       [21, 15,  0],
       [ 6,  3,  3]])

# ver la dimension de una matriz
A.shape

(3, 3)

# ver cantidad de elementos de una matriz
A.size

9

Multiplicacion o Producto de matrices#

La regla para la multiplicación de matrices generaliza la idea del producto interior que vimos con los vectores; y esta diseñada para facilitar las operaciones lineales básicas. Cuando multiplicamos matrices, el número de columnas de la primera matriz debe ser igual al número de filas de la segunda matriz; y el resultado de esta multiplicación va a tener el mismo número de filas que la primer matriz y el número de la columnas de la segunda matriz. Es decir, que si yo tengo una matriz A de dimensión 3x4 y la multiplico por una matriz B de dimensión 4x2, el resultado va a ser una matriz C de dimensión 3x2.

Algo a tener en cuenta a la hora de multiplicar matrices es que la propiedad connmutativa no se cumple. AxB no es lo mismo que BxA.

Veamos los ejemplos en Python.

# Ejemplo multiplicación de matrices
A = np.arange(1, 13).reshape(3, 4) #matriz de dimension 3x4
A

array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])

B = np.arange(8).reshape(4,2) #matriz de dimension 4x2
B

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

# Multiplicando A x B
A @ B #resulta en una matriz de dimension 3x2

array([[ 40,  50],
       [ 88, 114],
       [136, 178]])

# Multiplicando B x A
B @ A

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-28-b55e34ad9c31> in <module>()
      1 # Multiplicando B x A
----> 2 B @ A


ValueError: shapes (4,2) and (3,4) not aligned: 2 (dim 1) != 3 (dim 0)

Este ultimo ejemplo vemos que la propiedad conmutativa no se cumple, es más, Python nos arroja un error, ya que el número de columnas de B no coincide con el número de filas de A, por lo que ni siquiera se puede realizar la multiplicación de B x A.

Para una explicación más detallada del proceso de multiplicación de matrices, pueden consultar el siguiente tutorial.

La matriz identidad, la matriz inversa, la matriz transpuesta y el determinante#

La matriz identidad es el elemento neutro en la multiplicación de matrices, es el equivalente al número 1. Cualquier matriz multiplicada por la matriz identidad nos da como resultado la misma matriz. La matriz identidad es una matriz cuadrada (tiene siempre el mismo número de filas que de columnas); y su diagonal principal se compone de todos elementos 1 y el resto de los elementos se completan con 0. Suele representase con la letra I

Por ejemplo la matriz identidad de 3x3 sería la siguiente:

Ahora que conocemos el concepto de la matriz identidad, podemos llegar al concepto de la matriz inversa. Si tenemos una matriz A, la matriz inversa de A, que se representa como \(A^{-1}\) es aquella matriz cuadrada que hace que la multiplicación \(A\)x\(A^{-1}\) sea igual a la matriz identidad I. Es decir que es la matriz recíproca de A.

Tener en cuenta que esta matriz inversa en muchos casos puede no existir.En este caso se dice que la matriz es singular o degenerada. Una matriz es singular si y solo si su determinante es nulo.

El determinante es un número especial que puede calcularse sobre las matrices cuadradas. Se calcula como la suma de los productos de las diagonales de la matriz en una dirección menos la suma de los productos de las diagonales en la otra dirección. Se represente con el símbolo |A|.

Por último, la matriz transpuesta es aquella en que las filas se transforman en columnas y las columnas en filas. Se representa con el símbolo \(A^\intercal\)

Ejemplos en Python:

# Creando una matriz identidad de 2x2
I = np.eye(2)
I

array([[1., 0.],
       [0., 1.]])

# Multiplicar una matriz por la identidad nos da la misma matriz
A = np.array([[4, 7],
              [2, 6]])
A

array([[4, 7],
       [2, 6]])

A @ I # AxI = A

array([[4., 7.],
       [2., 6.]])

# Calculando el determinante de la matriz A
np.linalg.det(A)

10.000000000000002

# Calculando la inversa de A.
A_inv = np.linalg.inv(A)
A_inv

array([[ 0.6, -0.7],
       [-0.2,  0.4]])

# A x A_inv nos da como resultado I.
A @ A_inv

array([[1., 0.],
       [0., 1.]])

# Trasponiendo una matriz
A = np.arange(6).reshape(3, 2)
A

array([[0, 1],
       [2, 3],
       [4, 5]])

np.transpose(A)

array([[0, 2, 4],
       [1, 3, 5]])

Sistemas de ecuaciones lineales#

Una de las principales aplicaciones del Álgebra lineal consiste en resolver problemas de sistemas de ecuaciones lineales.

Una ecuación lineal es una ecuación que solo involucra sumas y restas de una variable o mas variables a la primera potencia. Es la ecuación de la línea recta.Cuando nuestro problema esta representado por más de una ecuación lineal, hablamos de un sistema de ecuaciones lineales. Por ejemplo, podríamos tener un sistema de dos ecuaciones con dos incógnitas como el siguiente:

La idea es encontrar el valor de \(x\) e \(y\) que resuelva ambas ecuaciones. Una forma en que podemos hacer esto, puede ser representando graficamente ambas rectas y buscar los puntos en que las rectas se cruzan.

En Python esto se puede hacer en forma muy sencilla con la ayuda de matplotlib.

# graficando el sistema de ecuaciones.
x_vals = np.linspace(0, 5, 50) # crea 50 valores entre 0 y 5
plt.plot(x_vals, (1 - x_vals)/-2) # grafica x - 2y = 1
plt.plot(x_vals, (11 - (3*x_vals))/2) # grafica 3x + 2y = 11
plt.axis(ymin = 0)

(-0.25, 5.25, 0, 5.875)

Luego de haber graficado las funciones, podemos ver que ambas rectas se cruzan en el punto (3, 1), es decir que la solución de nuestro sistema sería \(x=3\) e \(y=1\). En este caso, al tratarse de un sistema simple y con solo dos incógnitas, la solución gráfica puede ser de utilidad, pero para sistemas más complicados se necesita una solución numérica, es aquí donde entran a jugar las matrices.

Ese mismo sistema se podría representar como una ecuación matricial de la siguiente forma:

Lo que es lo mismo que decir que la matriz A por la matriz \(x\) nos da como resultado el vector b.

En este caso, ya sabemos el resultado de \(x\), por lo que podemos comprobar que nuestra solución es correcta realizando la multiplicación de matrices.

# Comprobando la solucion con la multiplicación de matrices.
A = np.array([[1., -2.],
              [3., 2.]])
x = np.array([[3.],[1.]])

A @ x

array([[ 1.],
       [11.]])

Para resolver en forma numérica los sistema de ecuaciones, existen varios métodos:

• El método de sustitución: El cual consiste en despejar en una de las ecuaciones cualquier incógnita, preferiblemente la que tenga menor coeficiente y a continuación sustituirla en otra ecuación por su valor.

• El método de igualacion: El cual se puede entender como un caso particular del método de sustitución en el que se despeja la misma incógnita en dos ecuaciones y a continuación se igualan entre sí la parte derecha de ambas ecuaciones.

• El método de reduccion: El procedimiento de este método consiste en transformar una de las ecuaciones (generalmente, mediante productos), de manera que obtengamos dos ecuaciones en la que una misma incógnita aparezca con el mismo coeficiente y distinto signo. A continuación, se suman ambas ecuaciones produciéndose así la reducción o cancelación de dicha incógnita, obteniendo una ecuación con una sola incógnita, donde el método de resolución es simple.

• El método gráfico: Que consiste en construir el gráfica de cada una de las ecuaciones del sistema. Este método (manualmente aplicado) solo resulta eficiente en el plano cartesiano (solo dos incógnitas).

• El método de Gauss: El método de eliminación de Gauss o simplemente método de Gauss consiste en convertir un sistema lineal de n ecuaciones con n incógnitas, en uno escalonado, en el que la primera ecuación tiene n incógnitas, la segunda ecuación tiene n - 1 incógnitas, …, hasta la última ecuación, que tiene 1 incógnita. De esta forma, será fácil partir de la última ecuación e ir subiendo para calcular el valor de las demás incógnitas.

• El método de Eliminación de Gauss-Jordan: El cual es una variante del método anterior, y consistente en triangular la matriz aumentada del sistema mediante transformaciones elementales, hasta obtener ecuaciones de una sola incógnita.

• El método de Cramer: El cual consiste en aplicar la regla de Cramer para resolver el sistema. Este método solo se puede aplicar cuando la matriz de coeficientes del sistema es cuadrada y de determinante no nulo.

La idea no es explicar cada uno de estos métodos, sino saber que existen y que Python nos hacer la vida mucho más fácil, ya que para resolver un sistema de ecuaciones simplemente debemos llamar a la función solve().

Por ejemplo, para resolver este sistema de 3 ecuaciones y 3 incógnitas.

Primero armamos la matriz A de coeficientes y la matriz b de resultados y luego utilizamos solve() para resolverla.

# Creando matriz de coeficientes
A = np.array([[1, 2, 3],
              [2, 5, 2],
              [6, -3, 1]])
A

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

# Creando matriz de resultados
b = np.array([6, 4, 2])
b

array([6, 4, 2])

# Resolviendo sistema de ecuaciones
x = np.linalg.solve(A, b)
x

array([0., 0., 2.])

# Comprobando la solucion
A @ x == b

array([ True,  True,  True])

Programación lineal#

La programación lineal estudia las situaciones en las que se exige maximizar o minimizar funciones que se encuentran sujetas a determinadas restricciones.

Consiste en optimizar (minimizar o maximizar) una función lineal, denominada función objetivo, de tal forma que las variables de dicha función estén sujetas a una serie de restricciones que expresamos mediante un sistema de inecuaciones lineales.

Para resolver un problema de programación lineal, debemos seguir los siguientes pasos:

1. Elegir las incógnitas.

2. Escribir la función objetivo en función de los datos del problema.

3. Escribir las restricciones en forma de sistema de inecuaciones.

4. Averiguar el conjunto de soluciones factibles representando gráficamente las restricciones.

5. Calcular las coordenadas de los vértices del recinto de soluciones factibles (si son pocos).

6. Calcular el valor de la función objetivo en cada uno de los vértices para ver en cuál de ellos presenta el valor máximo o mínimo según nos pida el problema (hay que tener en cuenta aquí la posible no existencia de solución).

Veamos un ejemplo y como Python nos ayuda a resolverlo en forma sencilla.

Supongamos que tenemos la siguiente función objetivo:

y las siguientes restricciones:

Podemos resolverlo utilizando PuLP, CVXOPT o graficamente (con matplotlib) de la siguiente forma.

# Resolviendo la optimizacion con pulp
from pulp import *

# declarando las variables
x1 = LpVariable("x1", 0, 800)   # 0<= x1 <= 40
x2 = LpVariable("x2", 0, 1000) # 0<= x2 <= 1000

# definiendo el problema
prob = LpProblem("problem", LpMaximize)

# definiendo las restricciones
prob += x1+1.5*x2 <= 750 
prob += 2*x1+x2 <= 1000
prob += x1>=0
prob += x2>=0

# definiendo la funcion objetivo a maximizar
prob += 50*x1+40*x2

# resolviendo el problema
status = prob.solve(GLPK(msg=0))
LpStatus[status]

# imprimiendo los resultados
(value(x1), value(x2))

(375.0, 250.0)

# Resolviendo el problema con cvxopt
from cvxopt import matrix, solvers

A = matrix([[-1., -2., 1., 0.], # columna de x1
            [-1.5, -1., 0., 1.]]) # columna de x2
b = matrix([750., 1000., 0., 0.]) # resultados
c = matrix([50., 40.]) # funcion objetivo

# resolviendo el problema
sol=solvers.lp(c,A,b)

     pcost       dcost       gap    pres   dres   k/t
 0: -2.5472e+04 -3.6797e+04  5e+03  0e+00  3e-01  1e+00
 1: -2.8720e+04 -2.9111e+04  1e+02  2e-16  9e-03  2e+01
 2: -2.8750e+04 -2.8754e+04  1e+00  8e-17  9e-05  2e-01
 3: -2.8750e+04 -2.8750e+04  1e-02  4e-16  9e-07  2e-03
 4: -2.8750e+04 -2.8750e+04  1e-04  9e-17  9e-09  2e-05
Optimal solution found.

# imprimiendo la solucion.
print('{0:.2f}, {1:.2f}'.format(sol['x'][0]*-1, sol['x'][1]*-1))

375.00, 250.00

# Resolviendo la optimizacion graficamente.
x_vals = np.linspace(0, 800, 10) # 10 valores entre 0 y 800
plt.plot(x_vals, ((750 - x_vals)/1.5)) # grafica x1 + 1.5x2 = 750
plt.plot(x_vals, (1000 - 2*x_vals)) # grafica 2x1 + x2 = 1000
plt.axis(ymin = 0)
plt.show()

Como podemos ver en el gráfico, ambas rectas se cruzan en la solución óptima, x1=375 y x2=250.

Con esto termino esta introducción al Álgebra lineal con Python.

Espero que hayan disfurtado de este tutorial tanto como yo disfrute escribirlo!

Saludos!

Este post fue escrito utilizando IPython notebook. Pueden descargar este notebook o ver su version estática en nbviewer.