Numpy: Álgebra lineal

NumPy ofrece una biblioteca robusta para realizar operaciones de álgebra lineal, comparable a sistemas como MATLAB o Octave, pero con la flexibilidad de Python. El módulo numpy.linalg proporciona todas las funciones fundamentales para trabajar con vectores y matrices, desde operaciones básicas hasta descomposiciones matriciales avanzadas. Esta capacidad hace de NumPy una herramienta esencial para campos como el aprendizaje automático, procesamiento de señales, física computacional y muchas otras áreas de la ciencia e ingeniería.

Vectores y matrices en NumPy

En NumPy, tanto los vectores como las matrices se representan mediante arrays. Un vector se implementa como un array unidimensional, mientras que una matriz es un array bidimensional. Esto permite una sintaxis clara y consistente para todas las operaciones de álgebra lineal.

Para crear un vector en NumPy, podemos utilizar la función array():

import numpy as np

# Crear un vector
v = np.array([1, 2, 3, 4])
print(v)  # [1 2 3 4]
print(v.shape)  # (4,)

Para crear una matriz, simplemente pasamos una lista de listas a la función array():

# Crear una matriz
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(A)
# [[1 2 3]
#  [4 5 6]
#  [7 8 9]]
print(A.shape)  # (3, 3)

NumPy ofrece varias funciones especializadas para crear matrices comunes en álgebra lineal:

# Matriz de ceros
zeros = np.zeros((3, 3))
print(zeros)

# Matriz de unos
ones = np.ones((2, 4))
print(ones)

# Matriz identidad
eye = np.eye(3)
print(eye)

# Matriz diagonal
diag = np.diag([1, 2, 3, 4])
print(diag)

# Matriz aleatoria
random = np.random.rand(3, 3)
print(random)

Operaciones matriciales básicas

Suma y resta de matrices

La suma y resta de matrices en NumPy se realiza elemento a elemento utilizando los operadores + y -:

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Suma de matrices
C = A + B
print(C)
# [[ 6  8]
#  [10 12]]

# Resta de matrices
D = A - B
print(D)
# [[-4 -4]
#  [-4 -4]]

Estas operaciones requieren que las matrices tengan dimensiones compatibles. En el caso de la suma y la resta, las matrices deben tener exactamente la misma forma.

Multiplicación escalar

Para multiplicar una matriz por un escalar, simplemente utilizamos el operador *:

# Multiplicación escalar
scalar = 2
E = scalar * A
print(E)
# [[2 4]
#  [6 8]]

Multiplicación de matrices

La multiplicación matricial es una operación fundamental en álgebra lineal y se puede realizar de varias formas en NumPy:

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Multiplicación matricial usando el operador @
C = A @ B
print(C)
# [[19 22]
#  [43 50]]

# Forma alternativa usando la función matmul
C = np.matmul(A, B)
print(C)
# [[19 22]
#  [43 50]]

# Forma alternativa más antigua usando la función dot
C = np.dot(A, B)
print(C)
# [[19 22]
#  [43 50]]

El operador @ fue introducido en Python 3.5 específicamente para la multiplicación de matrices, y es la forma más clara y recomendada para realizar esta operación.

Es importante recordar que para que la multiplicación de matrices A×B sea posible, el número de columnas de A debe ser igual al número de filas de B. Si A es una matriz de forma (m, n) y B es una matriz de forma (n, p), el resultado será una matriz de forma (m, p).

Multiplicación elemento a elemento

A veces necesitamos multiplicar matrices elemento por elemento en lugar de utilizar la multiplicación matricial. En NumPy, esto se puede hacer con el operador *:

# Multiplicación elemento a elemento
D = A * B
print(D)
# [[ 5 12]
#  [21 32]]

Potencia de matrices

La potencia de una matriz se puede calcular con la función matrix_power() del módulo linalg:

import numpy.linalg as la

# Potencia de una matriz
A = np.array([[1, 2], [3, 4]])
A_squared = la.matrix_power(A, 2)
print(A_squared)
# [[ 7 10]
#  [15 22]]

# Matriz inversa (potencia -1)
A_inv = la.matrix_power(A, -1)
print(A_inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

Transposición de matrices

La transposición de una matriz consiste en reflejar sus elementos respecto a la diagonal principal, convirtiendo filas en columnas y viceversa. NumPy ofrece varias formas de transponer una matriz:

A = np.array([[1, 2, 3], [4, 5, 6]])
print(A)
# [[1 2 3]
#  [4 5 6]]

# Usando el atributo T
A_T = A.T
print(A_T)
# [[1 4]
#  [2 5]
#  [3 6]]

# Usando la función transpose
A_T = np.transpose(A)
print(A_T)
# [[1 4]
#  [2 5]
#  [3 6]]

Para matrices complejas, NumPy también proporciona la transposición conjugada (o adjunta), que además de transponer la matriz, calcula el conjugado complejo de cada elemento:

# Matriz compleja
C = np.array([[1+2j, 3-4j], [5+6j, 7-8j]])

# Transposición conjugada (adjunta)
C_H = C.conj().T  # o simplemente C.H en versiones recientes
print(C_H)

Determinante de una matriz

El determinante es un valor escalar derivado de los elementos de una matriz cuadrada. Es fundamental en muchas aplicaciones de álgebra lineal y proporciona información sobre las propiedades de la matriz, como su invertibilidad.

NumPy proporciona la función det() para calcular el determinante:

import numpy.linalg as la

A = np.array([[1, 2], [3, 4]])
det_A = la.det(A)
print(det_A)  # -2.0

B = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
det_B = la.det(B)
print(det_B)  # 0.0

Un determinante igual a cero indica que la matriz es singular (no invertible), lo que significa que sus columnas o filas son linealmente dependientes.

Rango de una matriz

El rango de una matriz es la dimensión del espacio vectorial generado por sus columnas (o filas). Es igual al número máximo de columnas (o filas) linealmente independientes.

NumPy proporciona la función matrix_rank() para calcular el rango:

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
rank_A = la.matrix_rank(A)
print(rank_A)  # 2

B = np.eye(3)  # Matriz identidad 3x3
rank_B = la.matrix_rank(B)
print(rank_B)  # 3

El rango de la matriz A es 2, lo que confirma que no todas sus columnas son linealmente independientes (de hecho, su determinante es cero). La matriz identidad B tiene rango completo (igual a su dimensión), lo que indica que todas sus columnas son linealmente independientes.

Matriz inversa

La matriz inversa es uno de los conceptos fundamentales en álgebra lineal. Para una matriz cuadrada A, su inversa A⁻¹ satisface la ecuación A · A⁻¹ = A⁻¹ · A = I, donde I es la matriz identidad.

NumPy proporciona la función inv() para calcular la matriz inversa:

A = np.array([[1, 2], [3, 4]])
A_inv = la.inv(A)
print(A_inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

# Verificación: A · A⁻¹ = I
I = np.dot(A, A_inv)
print(np.round(I, 10))  # Redondeo para manejar errores de punto flotante
# [[1. 0.]
#  [0. 1.]]

Es importante destacar que no todas las matrices tienen inversa. Una matriz es invertible si y solo si su determinante es distinto de cero. Si intentamos calcular la inversa de una matriz singular, NumPy lanzará un error LinAlgError:

B = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
try:
    B_inv = la.inv(B)
except la.LinAlgError as e:
    print(f"Error: {e}")
    # Error: Singular matrix

Para sistemas de ecuaciones en los que la matriz de coeficientes podría ser singular, es más apropiado utilizar métodos como la pseudoinversa o soluciones de mínimos cuadrados.

Pseudoinversa de Moore-Penrose

La pseudoinversa de Moore-Penrose es una generalización de la matriz inversa que existe para cualquier matriz, incluso si no es cuadrada o es singular. Es particularmente útil para resolver sistemas de ecuaciones sobredeterminados o subdeterminados.

NumPy proporciona la función pinv() para calcular la pseudoinversa:

A = np.array([[1, 2, 3], [4, 5, 6]])
A_pinv = la.pinv(A)
print(A_pinv)
# [[-0.94444444  0.44444444]
#  [-0.11111111  0.11111111]
#  [ 0.72222222 -0.22222222]]

# Verificación para una matriz rectangular
I_approx = np.dot(A, A_pinv)
print(np.round(I_approx, 10))
# [[1. 0.]
#  [0. 1.]]

La pseudoinversa es especialmente útil en aplicaciones como regresión lineal, procesamiento de señales y problemas de optimización.

Sistemas de ecuaciones lineales

Resolver sistemas de ecuaciones lineales es una de las aplicaciones más comunes del álgebra lineal. Un sistema de ecuaciones lineales se puede expresar en forma matricial como Ax = b, donde A es la matriz de coeficientes, x es el vector de incógnitas y b es el vector de términos independientes.

NumPy proporciona la función solve() para resolver sistemas de ecuaciones lineales:

# Sistema:
# 2x + y = 1
# 3x + 2y = 2

A = np.array([[2, 1], [3, 2]])
b = np.array([1, 2])

x = la.solve(A, b)
print(x)  # [0. 1.]

# Verificación
result = np.dot(A, x)
print(result)  # [1. 2.]

Para sistemas que no tienen solución exacta (sobredeterminados) o tienen infinitas soluciones (subdeterminados), podemos utilizar la función lstsq() para encontrar la solución de mínimos cuadrados:

# Sistema sobredeterminado:
# x + y = 2
# 2x + y = 3
# 3x - y = 2

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

x, residuals, rank, s = la.lstsq(A, b, rcond=None)
print(x)  # [1. 1.]
print(residuals)  # [1.]  # Suma de los cuadrados de los residuos

La función lstsq() busca el vector x que minimiza la norma euclidiana ‖Ax - b‖.

Valores y vectores propios

Los valores propios (o eigenvalores) y vectores propios (o eigenvectores) juegan un papel crucial en muchas aplicaciones de álgebra lineal. Para una matriz cuadrada A, un vector no nulo v es un vector propio con valor propio λ si Av = λv.

NumPy proporciona la función eig() para calcular valores y vectores propios:

A = np.array([[1, 2], [3, 4]])

# Cálculo de valores y vectores propios
eigenvalues, eigenvectors = la.eig(A)

print("Valores propios:")
pr