Categoría: Matematicas

Uno de los mayores errores al aprender álgebra lineal es pensar que trata de números, fórmulas y cálculos. Esa es solo la superficie. En realidad, lo que estás aprendiendo es algo mucho más poderoso: cómo interpretar y combinar información.

Para entenderlo, hay que hacer un pequeño cambio de perspectiva. Una ecuación no es una operación. Es una afirmación sobre la realidad. Es una sentencia, \( x+y=10 \), está diciendo: “hay dos cantidades cuya suma es 10”.

En base a esta afirmación, piensa en una ecuación como una oración que describe al mundo:

• El perro es negro.
• El gato es naranja.

Ambas son sentencias que aportan información y al combinarlas tienes un conjunto de afirmaciones que combinadas te permiten deducir información al igual que un sistema de acuaciones.

Cuando aparece el sistema: combinar información

Imagina que estás resolviendo un problema como si fuera un pequeño caso de investigación. Tienes varias pistas y cada una aporta información parcial. Una pista puede decirte algo genérico, otra puede acotar un poco más, y juntas empiezan a dibujar una imagen más precisa. Eso es exactamente lo que hace un sistema de ecuaciones: combinar piezas de información para reducir la incertidumbre.

Pero no todas las combinaciones de información funcionan igual. Y aquí es donde el álgebra deja de ser mecánica y empieza a ser lógica. Cuando juntas varias sentencias, solo pueden pasar tres cosas:

1. Que cada nueva ecuación aporte algo distinto y coherente con lo anterior. En ese caso, el sistema converge hacia una única solución. Este es el escenario ideal. En términos matemáticos, el sistema es consistente y determinado. En términos prácticos, tienes suficiente información bien estructurada.
2. Las sentencias no se contradicen, pero tampoco aportan nada nuevo. Es como si alguien repitiera la misma pista con otras palabras. En matemáticas esto se traduce en ecuaciones dependientes, y el resultado es que no hay una única solución, sino infinitas. Sabes que algo cumple ciertas condiciones, pero no puedes identificarlo de forma única.
3. Las sentencias se contradicen. Una dice una cosa y otra dice lo contrario. En ese momento, el sistema deja de tener sentido. Has construido una realidad imposible.

🔹 3. Los tres comportamientos posibles (y por qué importan)

Aquí es donde quiero que prestes atención, porque esto es el núcleo conceptual.

Cuando combinas sentencias, solo pueden pasar tres cosas. No cuatro. No cinco. Tres.

Primer caso: las sentencias se complementan

Cada ecuación aporta información nueva. No se repiten, no se contradicen. Se refuerzan.

En ese caso, el sistema converge hacia una única solución.

Esto significa que has descrito la realidad con suficiente precisión como para identificar un único resultado.

Geométricamente, si quieres visualizarlo, es el punto donde dos líneas se cruzan.

Pero no te quedes con la imagen. Quédate con la idea:

Información suficiente y coherente produce una respuesta única.

Segundo caso: las sentencias se repiten

Aquí ocurre algo más sutil. Las ecuaciones no están mal. No hay contradicción. Pero una no aporta nada nuevo respecto a la otra.

Es como si alguien te dijera dos veces lo mismo con palabras distintas.

En ese caso, el sistema no falla, pero tampoco se define completamente. Tienes infinitas soluciones.

No porque el sistema sea incorrecto, sino porque no tienes suficiente información independiente.

Y esto es una idea muy importante: no importa cuántas ecuaciones tengas, sino cuánta información diferente contienen.

Tercer caso: las sentencias se contradicen

Aquí el sistema se rompe.

Una ecuación afirma algo, la otra afirma lo contrario. No existe ningún valor que pueda satisfacer ambas simultáneamente.

Y esto no es un problema de cálculo. Es un problema lógico.

Has construido un sistema imposible.

🔹 4. Qué ocurre cuando pasamos a matrices

Ahora, muchos de vosotros habéis visto matrices como algo técnico: filas, columnas, operaciones…

Pero quiero que lo veáis de otra forma.

Una matriz no es más que una forma de organizar sentencias.

• Cada fila → una ecuación → una afirmación
• Cada columna → una variable → una dimensión del problema

Y cuando operas con matrices, lo que estás haciendo en realidad es manipular información estructurada.

🔹 5. El concepto de rango (explicado sin fórmulas)

Aquí aparece una palabra que suele intimidar: el rango de una matriz.

Olvida la definición formal por un momento.

Quiero que lo entiendas así:

El rango es el número de sentencias realmente útiles que tienes.

No las que escribiste. Las que aportan información nueva.

Puedes tener 10 ecuaciones, pero si 8 son combinaciones de las otras 2, en realidad solo tienes 2 piezas de información.

Eso es el rango.

Y esto es lo que determina si un sistema tiene solución única, infinitas o ninguna.

🔹 6. Por qué esto es crucial fuera del aula

Ahora viene la parte que normalmente no se explica en matemáticas, pero es donde todo cobra sentido.

Esto mismo ocurre en ciencia de datos.

Un dataset no es más que un conjunto de sentencias:

• Cada fila → una observación
• Cada columna → una característica
• El modelo → una forma de combinar esa información

Si tus datos contienen información repetida, estás en el segundo caso: redundancia.

Si contienen contradicciones, estás en el tercero: inconsistencia.

Y si están bien construidos, estás en el primero: coherencia.

🔹 7. El error más común

Muchos estudiantes —y también muchos profesionales— creen que cuando algo falla, el problema está en el modelo.

No.

En muchísimos casos, el problema está en las sentencias.

Un modelo no puede arreglar un sistema de información mal construido.

Solo puede trabajar con lo que le das.

🔹 8. Cierre de la clase

Quiero que te quedes con esta idea, porque es la que marca la diferencia entre alguien que aplica fórmulas y alguien que entiende lo que hace:

Antes de resolver un sistema, pregúntate:

• ¿Estas ecuaciones aportan información nueva?
• ¿Se están repitiendo?
• ¿Se contradicen?

Si sabes responder eso, resolver el sistema es casi lo de menos.

Porque habrás entendido lo importante:

Las matemáticas no consisten en calcular resultados.
Consisten en interpretar información.

NumPy (Numerical Python) es un paquete de código abierto ampliamente utilizado en ciencia e ingeniería.

Realiza una amplia variedad de operaciones matemáticas avanzadas con alta eficiencia. En este artículo se verán varias funciones clave de NumPy como la creación de arrays, segmentación (slicing), indexación, cambio de forma (reshape) y apilamiento.

Arrays con Numpy

Los arrays son una de las estructuras de datos fundamentales de la librería NumPy, esenciales para organizar tus datos. Puedes pensar en ellos como una cuadrícula de valores, todos del mismo tipo.

Si has utilizado listas de Python anteriormente, recordarás que son convenientes, ya que puedes almacenar diferentes tipos de datos. Sin embargo, las listas de Python tienen funciones limitadas y ocupan más espacio y tiempo de procesamiento que los arreglos de NumPy.

NumPy proporciona un objeto array que es mucho más rápido y compacto que las listas de Python. A través de su amplia integración de API, la librería ofrece muchas funciones incorporadas que facilitan mucho el cómputo con solo unas pocas líneas de código. Esto puede ser una ventaja enorme al realizar operaciones matemáticas en grandes conjuntos de datos.

El objeto de array en NumPy se llama ndarray, que significa “array n-dimensional”. Un arreglo 1-D representa una lista estándar de valores en una dimensión. Recuerda que en NumPy, todos los elementos dentro de un array son del mismo tipo.

Crear arrays

Puedes crear un array de una dimensión utilizando simplemente la función array() la cual recibe una lista de valores como argumento.

a = np.array([1, 2, 3])
print(a)

a = np.array([1, 2, 3])
print(a)

[1 2 3]

[1 2 3]

Otra forma de implementar un array es utilizando np.arange(). Esta función devolverá un arreglo de valores espaciados uniformemente dentro de un intervalo dado.

b = np.arange(3)
print(b)

b = np.arange(3)
print(b)

[0 1 2]

[0 1 2]

Otro ejemplo usando las características de np.arange():

# Crear un array que comienza por 1, termina en 20, incrementando valores por 3.
c = np.arange(1, 20, 3)
print(c)

# Crear un array que comienza por 1, termina en 20, incrementando valores por 3.
c = np.arange(1, 20, 3)
print(c)

[ 1  4  7 10 13 16 19]

[ 1  4  7 10 13 16 19]

La función np.linspace() puede generar valores espaciados uniformemente en un intervalo. En este ejemplo el intervalo es entre 0 y 100 y el número de elementos en el array, 5.

lin_spaced_arr = np.linspace(0, 100, 5)
print(lin_spaced_arr)

lin_spaced_arr = np.linspace(0, 100, 5)
print(lin_spaced_arr)

[  0.  25.  50.  75. 100.]

[  0.  25.  50.  75. 100.]

Por defecto los valores devueltos por np.linspace son de punto flotante (np.float64). Puedes especificar fácilmente tu tipo de datos usando dtype.

c_int = np.arange(1, 20, 3, dtype=int)
print(c_int)

c_int = np.arange(1, 20, 3, dtype=int)
print(c_int)

[ 1  4  7 10 13 16 19]

[ 1  4  7 10 13 16 19]

Más arrays de NumPy

Puedes crear arreglos fácilmente con funciones integradas tales como: 

• np.ones(shape) – Devuelve un nuevo arreglo estableciendo los valores en uno.
• np.zeros(shape) – Devuelve un nuevo arreglo estableciendo los valores en cero.
• np.empty(shape) – Devuelve un nuevo arreglo sin inicializar.
• np.random.rand(shape) – Devuelve un nuevo arreglo con valores randoms. 

Arrays multidimensionales

Con NumPy también puedes crear arreglos con más de una dimensión. En los ejemplos anteriores, trabajaste con arreglos 1-D, donde puedes acceder a sus elementos utilizando un único índice. Un arreglo multidimensional tiene más de una columna. Piensa en un arreglo multidimensional como una hoja de Excel donde cada fila/columna representa una dimensión.

# Crear array de dos dimensiones (2-D)
two_dim_arr = np.array([[1,2,3], [4,5,6]])
print(two_dim_arr)

# Crear array de dos dimensiones (2-D)
two_dim_arr = np.array([[1,2,3], [4,5,6]])
print(two_dim_arr)

[[1 2 3]
 [4 5 6]]

[[1 2 3]
 [4 5 6]]

Una forma alternativa de crear un arreglo multidimensional es cambiando la forma del arreglo 1-D inicial. Utilizando np.reshape() puedes reorganizar los elementos del arreglo anterior en una nueva forma.

# 1-D array 
one_dim_arr = np.array([1, 2, 3, 4, 5, 6])

# Multidimensional array using reshape()
multi_dim_arr = np.reshape(
                one_dim_arr, # the array to be reshaped
               (2,3) # dimensions of the new array
              )
# Print the new 2-D array with two rows and three columns
print(multi_dim_arr)

# 1-D array 
one_dim_arr = np.array([1, 2, 3, 4, 5, 6])

# Multidimensional array using reshape()
multi_dim_arr = np.reshape(
                one_dim_arr, # the array to be reshaped
               (2,3) # dimensions of the new array
              )
# Print the new 2-D array with two rows and three columns
print(multi_dim_arr)

[[1 2 3]
 [4 5 6]]

[[1 2 3]
 [4 5 6]]

Encontrar el tamaño, la forma y la dimensión.

Necesitarás saber cómo encontrar el tamaño, la dimensión y la forma de un array. Todos estos son atributos de un ndarray y se puede acceder a ellos de la siguiente manera:

• ndarray.ndim – Almacena el número de dimensiones del arreglo.
• ndarray.shape – Almacena la forma del arreglo. Cada número en la tupla denota la longitud de cada dimensión correspondiente.
• ndarray.size – Almacena el número de elementos en el arreglo.

Operaciones matemáticas con arrays

En esta sección, verás que NumPy te permite realizar rápidamente sumas, restas, multiplicaciones y divisiones elemento por elemento, tanto para arreglos 1-D como multidimensionales.

Las operaciones se realizan utilizando el símbolo matemático correspondiente: ‘+’, ‘-‘ y ‘*’. Recuerda que la suma de listas de Python funciona de manera completamente diferente, ya que concatena las listas creando una lista más larga; además, la resta y la multiplicación de listas de Python no funcionan.

arr_1 = np.array([2, 4, 6])
arr_2 = np.array([1, 3, 5])

# Adding two 1-D arrays
addition = arr_1 + arr_2
print(addition)

# Subtracting two 1-D arrays
subtraction = arr_1 - arr_2
print(subtraction)

# Multiplying two 1-D arrays elementwise
multiplication = arr_1 * arr_2
print(multiplication)

arr_1 = np.array([2, 4, 6])
arr_2 = np.array([1, 3, 5])

# Adding two 1-D arrays
addition = arr_1 + arr_2
print(addition)

# Subtracting two 1-D arrays
subtraction = arr_1 - arr_2
print(subtraction)

# Multiplying two 1-D arrays elementwise
multiplication = arr_1 * arr_2
print(multiplication)

[ 3  7 11]
[1 1 1]
[ 2 12 30]

[ 3  7 11]
[1 1 1]
[ 2 12 30]

Multiplicación de un vector por un escalar (broadcasting)

Supón que necesitas convertir millas a kilómetros. Para hacerlo, puedes usar las funciones de arreglos de NumPy que has aprendido hasta ahora. Puedes lograrlo realizando una operación entre un arreglo (millas) y un solo número (la tasa de conversión, que es un escalar). Dado que 1 milla = 1.6 km, NumPy calcula cada multiplicación dentro de cada celda. 

Este concepto se llama broadcasting (transmisión), el cual te permite realizar operaciones específicamente en arreglos de diferentes formas.

vector = np.array([1, 2])
vector * 1.6

vector = np.array([1, 2])
vector * 1.6

array([1.6, 3.2])

array([1.6, 3.2])

Indexación y segmentación (slicing)

La indexación es muy útil ya que te permite seleccionar elementos específicos de un array. También te permite seleccionar filas, columnas o planos completos, como verás en futuras tareas para arreglos multidimensionales. 

Indexación

Seleccionemos elementos específicos de los arreglos dados. 

# Seleccionar el tercer elemento del array.
a = ([1, 2, 3, 4, 5])
print(a[2])

# Seleccionar el primer elemento.
print(a[0])

# Seleccionar el tercer elemento del array.
a = ([1, 2, 3, 4, 5])
print(a[2])

# Seleccionar el primer elemento.
print(a[0])

3
1

3
1

Para arrays multidimensionales de forma \(n\), para indexar un elemento específico, debes ingresar \(n\) índices, uno por cada dimensión.

# Indexing on a 2-D array
two_dim = np.array(([1, 2, 3],
          [4, 5, 6], 
          [7, 8, 9]))

# Select element number 8 from the 2-D array using indices i, j.
print(two_dim[2][1])

# Indexing on a 2-D array
two_dim = np.array(([1, 2, 3],
          [4, 5, 6], 
          [7, 8, 9]))

# Select element number 8 from the 2-D array using indices i, j.
print(two_dim[2][1])

8

8

Segmentación (Slicing)

La segmentación te devuelve una sublista de los elementos que especifiques del array. La notación de segmentación define un valor de inicio y uno de fin, y copia la lista desde el inicio hasta el final, pero sin incluir este último (el final es exclusivo).

La sintaxis es:

arreglo[inicio:fin:paso]

Si no se pasa ningún valor al inicio, se asume que inicio = 0; si no se pasa ningún valor para el fin, se asume que fin = longitud del arreglo y si no se pasa ningún valor para el paso, se asume que paso = 1.

# Slice the array a to get the array [2,3,4]
sliced_arr = a[1:4]
print(sliced_arr)

# Slice the array a to get the array [2,3,4]
sliced_arr = a[1:4]
print(sliced_arr)

[2, 3, 4]

[2, 3, 4]

Apilamiento (Stacking)

Finalmente, el apilamiento es una característica de NumPy que permite una mayor personalización de los arrays. Consiste en unir dos o más arreglos, ya sea de forma horizontal o vertical, lo que significa que se realiza a lo largo de un nuevo eje.

• np.vstack() – apila verticalmente.
• np.hstack() – apila horizontalmente.
• np.hsplit() – divide un arreglo en varios arreglos más pequeños.

a1 = np.array([[1,1], 
               [2,2]])
a2 = np.array([[3,3],
              [4,4]])
print(f'a1:\n{a1}')
print(f'a2:\n{a2}')

a1 = np.array([[1,1], 
               [2,2]])
a2 = np.array([[3,3],
              [4,4]])
print(f'a1:\n{a1}')
print(f'a2:\n{a2}')

a1:
[[1 1]
 [2 2]]
a2:
[[3 3]
 [4 4]]

a1:
[[1 1]
 [2 2]]
a2:
[[3 3]
 [4 4]]

Apilar los arrays verticalmente

# Stack the arrays vertically
vert_stack = np.vstack((a1, a2))
print(vert_stack)

# Stack the arrays vertically
vert_stack = np.vstack((a1, a2))
print(vert_stack)

[[1 1]
 [2 2]
 [3 3]
 [4 4]]

[[1 1]
 [2 2]
 [3 3]
 [4 4]]

Apilar horizontalmente:

# Stack the arrays horizontally
horz_stack = np.hstack((a1, a2))
print(horz_stack)

# Stack the arrays horizontally
horz_stack = np.hstack((a1, a2))
print(horz_stack)

[[1 1 3 3]
 [2 2 4 4]]

[[1 1 3 3]
 [2 2 4 4]]

El ajuste por mínimos cuadrados es uno de los pilares del análisis de datos. Nos permite encontrar patrones y relaciones entre variables incluso cuando los datos no son perfectos. La idea esencial es siempre la misma:

Buscar los coeficientes que minimicen el error entre las observaciones reales y las predicciones del modelo.

A partir de aquí, la regresión lineal se convierte en la base de modelos más complejos de machine learning, donde la idea de “ajustar” parámetros para minimizar errores sigue siendo el núcleo de todo el proceso.

De los datos a la recta

Imaginemos que tenemos dos variables, x y y, y sospechamos que están relacionadas de forma lineal. Por ejemplo:

x = [1, 2, 3, 4, 5]
y = [7, 11, 15, 19, 23]

x = [1, 2, 3, 4, 5]
y = [7, 11, 15, 19, 23]

import matplotlib.pyplot as plt

# Datos
x = [1, 2, 3, 4, 5]
y = [7, 11, 15, 19, 23]

# Crear gráfico de dispersión
plt.scatter(x, y, color='blue', marker='o', label='Datos')

# Etiquetas y título
plt.title('Gráfico de dispersión: relación entre X y Y')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.legend()

# Mostrar gráfico
plt.show()

import matplotlib.pyplot as plt

# Datos
x = [1, 2, 3, 4, 5]
y = [7, 11, 15, 19, 23]

# Crear gráfico de dispersión
plt.scatter(x, y, color='blue', marker='o', label='Datos')

# Etiquetas y título
plt.title('Gráfico de dispersión: relación entre X y Y')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.legend()

# Mostrar gráfico
plt.show()

El modelo lineal

Si las variables son linealmente dependientes, podemos escribir:

$$y = \beta_0 + \beta_1 x$$

donde:

• β₀ (beta cero) es el intercepto, el valor de y cuando x = 0.
• β₁ (beta uno) es la pendiente, que indica cuánto cambia y por cada unidad de x.

Nuestro objetivo es determinar estos coeficientes (β₀ y β₁) de manera que la recta se ajuste lo mejor posible a los datos.

Representando el modelo como un sistema lineal

Podemos escribir la ecuación anterior en forma matricial.
Para cada observación de \(x\), formamos una fila en una matriz \(M\) que tiene dos columnas: una con los valores de \(x\) y otra con unos.

$$M =\begin{bmatrix}1 & 1 \\ 2 & 1 \\ 3 & 1 \\ 4 & 1 \\ 5 & 1 \end{bmatrix} , \beta = \begin{bmatrix} \beta_1 \\ \beta_0 \end{bmatrix}$$

Entonces el modelo se puede escribir como:

$$y = M \cdot \beta$$

donde \(y\) es el vector de observaciones reales.

Si tenemos más ecuaciones que incógnitas (más datos que parámetros), se trata de un sistema sobredeterminado, y no habrá una solución exacta. En esos casos, buscamos una solución aproximada, aquella que minimiza los errores o residuos.

Resolviendo con el método de los mínimos cuadrados

La solución óptima en el sentido de mínimos cuadrados se obtiene resolviendo la ecuación normal:

$$\hat{\beta} = (M^T M)^{-1} M^T y$$

Aquí:

• \( M^T \) es la traspuesta de \(M\).
• \( (M^T M)^{-1} \) es su inversa.
• \( \hat{\beta} \) son las estimaciones de los coeficientes.

En el ejemplo de datos anteriores, este procedimiento nos da:

$$\beta_0 = 3, \quad \beta_1 = 4$$

por lo que la ecuación ajustada es:

$$y = 3 + 4x$$

Si probamos estos valores, obtenemos una coincidencia exacta con los datos: la relación es perfectamente lineal.

Cuando los datos no son perfectamente lineales

En la práctica, los datos reales rara vez se ajustan perfectamente a una línea.
Supongamos que ahora tenemos un conjunto diferente:

import matplotlib.pyplot as plt

# Datos
x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
y = [8, 10, 14, 18, 21, 23, 26, 27, 29, 33]

# Crear gráfico de dispersión
plt.scatter(x, y, color='blue', marker='o', label='Datos')

# Etiquetas y título
plt.title('Gráfico de dispersión: relación entre X y Y')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.legend()

# Mostrar gráfico
plt.show()

import matplotlib.pyplot as plt

# Datos
x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
y = [8, 10, 14, 18, 21, 23, 26, 27, 29, 33]

# Crear gráfico de dispersión
plt.scatter(x, y, color='blue', marker='o', label='Datos')

# Etiquetas y título
plt.title('Gráfico de dispersión: relación entre X y Y')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.legend()

# Mostrar gráfico
plt.show()

Vemos que no todos caen exactamente sobre una misma recta. Aun así, queremos una línea que represente la tendencia general: esa es la llamada línea de mejor ajuste (line of best fit).

La idea de los residuos

Cada punto tiene una pequeña distancia vertical hasta la recta estimada. Esa diferencia se llama residuo y se define como:

$$\text{residuo}_i = y_i – \hat{y}_i$$

donde \(\hat{y}_i\) es el valor predicho por el modelo. El método de mínimos cuadrados busca los parámetros β₀ y β₁ que minimicen la suma de los cuadrados de esos residuos:

$$\text{min} \sum_{i=1}^{n} (y_i – \hat{y}_i)^2$$

En forma desarrollada o explícita:

$$\text{min:} (y_1 – \hat{y}_1)^2 + \cdots $$

$$ + (y_n – \hat{y}_n)^2$$

Al hacerlo, la recta resultante será aquella que “pasa más cerca” de todos los puntos en promedio.

Interpretación visual

• β₀ desplaza la recta hacia arriba o abajo.
• β₁ modifica su inclinación.

Podemos imaginar que “movemos” y “rotamos” la línea hasta que la suma de los residuos sea lo más pequeña posible. En ese punto, hemos encontrado la línea de mejor ajuste.

Extensión a modelos polinomiales

La misma lógica se puede aplicar cuando la relación no es lineal. Por ejemplo, si creemos que \(y \) depende de \(x²\), podemos ajustar un modelo cuadrático:

$$y = \beta_0 + \beta_1 x + \beta_2 x^2$$

En este caso, la matriz \(M \) tendrá tres columnas: una para \(x^2\), una para \(x\) y una para los unos.

$$M = \begin{bmatrix} x_1^2 & x_1 & 1 \\ x_2^2 & x_2 & 1 \\ \vdots & \vdots & \vdots \\ x_n^2 & x_n & 1 \end{bmatrix}$$

Y nuevamente, los coeficientes se obtienen con la misma fórmula general:

$$\hat{\beta} = (M^T M)^{-1} M^T y$$

Si aplicamos este método a un conjunto de datos donde la relación es cuadrática, obtendremos una curva que se adapta mucho mejor a los puntos observados que una simple línea.

Generalización

Esta metodología puede ampliarse para ajustar polinomios de grado superior o otros tipos de funciones (logarítmicas, exponenciales, etc.), siempre que podamos expresar el modelo en forma lineal respecto a los coeficientes β.

Por ejemplo, para un polinomio cúbico \(y=\):

$$\beta_0 + \beta_1 x + \cdots + \beta_3 x^3$$

solo tendríamos que añadir una columna más a \(M\)con los valores de (x^3).

Consideremos un sistema sobredeterminado, donde la matriz \( M \) no es cuadrada, sino rectangular: tiene dimensiones \(n \times p \), siendo \( n > p \). Esto ocurre, por ejemplo, cuando tenemos más observaciones que variables.

Podemos imaginar un ejemplo sencillo:

• \( x \) –> representa los precios de diferentes productos del supermercado.
• \( M\) –> es una matriz que contiene la cantidad de cada producto que compraron distintas personas.
• \( y\)–> es el vector del coste total de la compra de cada individuo.

Hasta aquí, todo parece un sistema de ecuaciones lineales clásico: \( Mx = y \).
Sin embargo, en la vida real los precios no son constantes: cambian entre tiendas, promociones o regiones. Así que incorporamos esa variabilidad aleatoria en el modelo.

La idea probabilística

Si cada individuo experimenta ligeras variaciones en los precios, no hay un único valor “verdadero” para cada producto. Por eso, nuestro objetivo no será encontrar un precio exacto, sino una estimación promedio de los precios:

$$ \hat{x}$$

donde el símbolo del “sombrero” \(ˆ\) indica que se trata de una estimación. De este modo, cada elemento de \( \hat{x} \) representa el precio medio estimado de un producto.

Con esa estimación podemos obtener un coste total estimado multiplicando:

$$\hat{y} = M \hat{x}$$

Si nuestras estimaciones son buenas, entonces \(\hat{y} \) debería aproximarse lo más posible a los valores reales \( y \).

Qué significa “lo más posible”

Para medir qué tan cerca está nuestro resultado, podemos calcular la diferencia entre los valores reales y los estimados:

$$y – \hat{y}$$

y luego obtener su norma (la “distancia” entre ambos vectores). Nuestro objetivo es minimizar esa distancia. Este enfoque, de minimizar los errores cuadráticos, es justamente lo que hace la regresión OLS.

Cómo funciona OLS

En lugar de resolver directamente \( Mx = y \), multiplicamos ambos lados por la transpuesta de \( M \):

$$M^T M x = M^T y$$

Con esto obtenemos una matriz cuadrada \(p \times p\), que puede invertirse (si las columnas de \( M \) son linealmente independientes).

De ahí obtenemos la fórmula del vector estimado:

$$\hat{x} = (M^T M)^{-1} M^T y$$

Este vector \( \hat{x} \)contiene los coeficientes de regresión, también llamados parámetros estimados.

Comprobando la calidad del modelo

Una vez tenemos \( \hat{x} \), podemos comparar las predicciones \( \hat{y} \) con los valores reales \( y \).
Hay varias maneras de medir la calidad del ajuste:

1. La norma del error: mide la magnitud de las diferencias \( y – \hat{y}\).
2. El coeficiente de determinación \( R^2 \):
   Este valor, entre 0 y 1, indica qué proporción de la variación real de los datos está explicada por el modelo.
   • \( R^2 = 1 \): el modelo explica toda la variación (ajuste perfecto).
   • \(R^2 = 0\): el modelo no explica nada (mal ajuste).

En esencia, \(R^2 \) compara cuánto mejoran nuestras predicciones respecto a usar solo el promedio de los datos.

El concepto de residuo

El residuo es la diferencia entre el valor real y el valor estimado para cada observación individual:

$$\text{residuo}_i = y_i – \hat{y}_i$$

Gráficamente, si representamos los puntos reales y los estimados, cada residuo es la distancia vertical entre ellos. Un buen modelo tendrá residuos pequeños y distribuidos al azar.

Ejemplo: precios de supermercado

Supongamos que tenemos 9 productos (manzanas, naranjas, pan, etc.) y los registros de compras de 450 personas. Cada persona pudo haber pagado precios ligeramente diferentes por los mismos artículos, dependiendo de la tienda o el día.

• Si no hubiera variaciones (todos pagan el mismo precio), entonces \( \hat{y} = y \): el modelo reproduce exactamente los costes reales.
• Si hay variaciones, los puntos estimados y los reales ya no coincidirán perfectamente, y aparecerán residuos.

Aun así, la regresión OLS nos permite encontrar los promedios óptimos que minimizan los errores totales, proporcionando una estimación realista del comportamiento de compra.

Paso 1: importar librerías y definir matrices

import numpy as np

# Matriz M: cantidades compradas (filas = personas, columnas = productos)
M = np.array([
    [2, 3],   # persona 1: 2 manzanas, 3 naranjas
    [3, 5]    # persona 2: 3 manzanas, 5 naranjas
])

# Vector y: total gastado por cada persona
y = np.array([13, 21.6])

import numpy as np

# Matriz M: cantidades compradas (filas = personas, columnas = productos)
M = np.array([
    [2, 3],   # persona 1: 2 manzanas, 3 naranjas
    [3, 5]    # persona 2: 3 manzanas, 5 naranjas
])

# Vector y: total gastado por cada persona
y = np.array([13, 21.6])

Paso 2: aplicar la fórmula de OLS

$$\hat{x} = (M^T M)^{-1} M^T y$$

# Cálculo del estimador OLS
x_hat = np.linalg.inv(M.T @ M) @ M.T @ y
print("Estimaciones de precios:", x_hat)

# Cálculo del estimador OLS
x_hat = np.linalg.inv(M.T @ M) @ M.T @ y
print("Estimaciones de precios:", x_hat)

Salida esperada:

Estimaciones de precios: [0.2 4.2]

Estimaciones de precios: [0.2 4.2]

Esto indica que:

• Precio estimado de una manzana ≈ $0.20
• Precio estimado de una naranja ≈ $4.20

Paso 3: calcular las predicciones y los residuos

# Predicciones del modelo
y_hat = M @ x_hat

# Cálculo de los residuos
residuos = y - y_hat

print("Predicciones (ŷ):", y_hat)
print("Residuos:", residuos)

# Predicciones del modelo
y_hat = M @ x_hat

# Cálculo de los residuos
residuos = y - y_hat

print("Predicciones (ŷ):", y_hat)
print("Residuos:", residuos)

Salida esperada:

Predicciones (ŷ): [13.  21.6]
Residuos: [0. 0.]

Predicciones (ŷ): [13.  21.6]
Residuos: [0. 0.]

El modelo reproduce exactamente los valores observados, aunque —como vimos— los precios estimados no reflejan la realidad del sistema, porque el modelo no contempla variaciones individuales.

Paso 4: medir la calidad del ajuste (R²)

# Cálculo del coeficiente de determinación
R2 = 1 - np.sum((y - y_hat)**2) / np.sum((y - np.mean(y))**2)
print("Coeficiente R²:", R2)

# Cálculo del coeficiente de determinación
R2 = 1 - np.sum((y - y_hat)**2) / np.sum((y - np.mean(y))**2)
print("Coeficiente R²:", R2)

Salida esperada:

Coeficiente R²: 1.0

Coeficiente R²: 1.0

El modelo tiene un ajuste perfecto con solo dos observaciones, pero los parámetros no son realistas: este es un ejemplo ideal para entender que un buen ajuste no siempre significa un modelo correcto.

Conclusión

La regresión de mínimos cuadrados ordinarios (OLS) nos permite resolver sistemas lineales cuando los datos no encajan exactamente, buscando el mejor ajuste promedio posible.
En este ejemplo, aunque el modelo predice perfectamente los costes totales, los parámetros estimados (los precios) son poco realistas porque los datos iniciales no cumplían las condiciones necesarias para una solución exacta.

En la práctica, el OLS es una herramienta central en ciencia de datos: es la base de la regresión lineal, de muchos métodos de modelado predictivo y de análisis estadístico que permiten capturar relaciones entre variables reales y ruidosas.

Imagina que tenemos una situación un poco diferente a la habitual. Supón que registramos el comportamiento de compra de varias personas en un supermercado. Tenemos los precios de los productos, los gastos totales de cada cliente, y una matriz que muestra cuántas unidades de cada producto compró cada uno.

Cada persona compra distintas cantidades de unos pocos artículos —por ejemplo, tres productos diferentes—, y tenemos los datos de 100 clientes. Esto significa que nuestro número de filas (n = 100) es mucho mayor que el número de columnas (p = 3). En otras palabras, la matriz que representa esas compras ya no es cuadrada, sino rectangular.

¿Y cuál es el problema?
Que las matrices rectangulares no tienen inversa, por lo que no podemos usar el método de “retroceso” (backsolving) que aplicábamos antes para resolver el sistema. En vez de rendirnos, vamos a usar un truco matemático muy útil: el transpuesto de la matriz.

Qué significa transponer una matriz

Transponer una matriz significa reflejar sus elementos respecto de la diagonal principal (la línea que va del vértice superior izquierdo al inferior derecho).
Por ejemplo:

• Si tienes un vector fila, su transpuesto será un vector columna.
• Si tienes una matriz cuadrada, transponerla intercambia sus filas por columnas.
• Si tienes una matriz rectangular de 3×2, su transpuesta será de 2×3.

Una propiedad interesante es que el transpuesto del transpuesto te devuelve la matriz original.

Qué ocurre al multiplicar una matriz por su transpuesta

Supón que tienes una matriz \( M \) de tamaño \( n \times p \) (por ejemplo, 100 × 3). Si calculas el producto de su transpuesta por ella misma, es decir \( M^T \times M\), el resultado sí es una matriz cuadrada de tamaño \( p \times p \) (en este ejemplo, 3 × 3).
Eso significa que podemos intentar invertirla.

En resumen:

• Si \( M \) es \(n \times p \),
• entonces \( M^T \) es \(p \times n\),
• y \( M^T M\) es \( p \times p \), una matriz cuadrada.

Esto nos abre una puerta para resolver el sistema que antes parecía imposible.

Aplicando el truco: resolver el sistema con el transpuesto

En lugar de intentar resolver directamente \( Mx = y \) (lo cual no podemos hacer porque \(M\) no tiene inversa), multiplicamos ambos lados por \( M^T \):

$$M^T Mx = M^T y$$

Ahora, el sistema tiene una matriz cuadrada \(M^T M \).
Si esta matriz tiene inversa, podemos calcular:

$$x = (M^T M)^{-1} M^T y$$

Y así obtenemos una solución para x.

Cuándo existe esa inversa

La clave está en la independencia lineal de las columnas de \(M\). Si las columnas son linealmente independientes (es decir, ninguna es combinación de las otras), entonces \( M^T M \) tendrá una inversa.

Si no lo son —por ejemplo, si una columna es el doble de otra—, entonces \( M^T M \) no será invertible y este método no funcionará.

Cómo hacerlo en Python

En NumPy puedes obtener el transpuesto de una matriz simplemente con .T.
Por ejemplo:

import numpy as np

M = np.array([[1, 2],
              [3, 4],
              [5, 6]])

M_trans = M.T
print(M_trans)

import numpy as np

M = np.array([[1, 2],
              [3, 4],
              [5, 6]])

M_trans = M.T
print(M_trans)

El resultado será:

[[1 3 5]
 [2 4 6]]

[[1 3 5]
 [2 4 6]]

Relación con la regresión lineal

Lo que acabamos de hacer manualmente es, en esencia, una regresión lineal por mínimos cuadrados. En lugar de usar fórmulas estadísticas, usamos álgebra matricial para encontrar los coeficientes (precios) que mejor explican los datos observados (costes totales).

En Python, este mismo procedimiento lo implementa la clase LinearRegression del paquete scikit-learn. Basta con:

from sklearn.linear_model import LinearRegression

modelo = LinearRegression(fit_intercept=False)
modelo.fit(M, y)
x_hat = modelo.coef_

from sklearn.linear_model import LinearRegression

modelo = LinearRegression(fit_intercept=False)
modelo.fit(M, y)
x_hat = modelo.coef_

En resumen

Cuando la matriz no es cuadrada, no podemos invertirla directamente. Pero si multiplicamos por su transpuesta, obtenemos una matriz cuadrada que sí puede tener inversa.
Ese es el truco del transpuesto, y es la base de la regresión lineal y de muchos métodos de optimización en ciencia de datos.

El backsolving es un método para resolver sistemas de ecuaciones lineales una vez que ya hemos reducido el sistema a forma triangular superior mediante eliminación gaussiana (o cualquier otra técnica de factorización como LU).

Trabajamos el mismo problema del articulo anterior:

$$ M = \begin{bmatrix} 3 & 6 & 2 \\ 10 & 3 & 8 \\ 1 & 7 & 5 \end{bmatrix} , p = \begin{bmatrix} p_1 \\ p_2 \\ p_3 \end{bmatrix}, c = \begin{bmatrix} 34 \\ 69 \\ 48 \end{bmatrix} $$

Queremos resolver los precios \( (p) \) de cada producto. En este caso utilizaremos eliminación gaussiana y luego sustitución hacia atrás — eso es, reducir el sistema a triangular superior y luego encontrar las incógnitas empezando por la última.

Eliminación gaussiana (reducción hacia triangular superior)

Escribimos el sistema en sus ecuaciones explícitas (producto fila por columna):

$$ \begin{cases} 3p_1 + 6p_2 + 2p_3 = 34 \\ 10p_1 + 3p_2 + 8p_3 = 69 \\ 1p_1 + 7p_2 + 5p_3 = 48 \end{cases}$$

Paso 1 — usar la primera ecuación para eliminar \(p_1\) de las filas 2 y 3.

Despejamos \(p_1\) (temporalmente) o eliminamos directamente:

Multiplicador para la fila 2: \(\frac{10}{3}=2\).
Hacemos \(R_2 \leftarrow R_2 – 2R_1\):

• Columna 1: \(10 – 2\cdot 3 = 0\).
• Columna 2: \(3 – 2\cdot 6 = 3 – 12 = -9\).

Multiplicador para la fila 3: \(\frac{1}{3}\).
Hacemos \(R_3 \leftarrow R_3 – \frac{1}{3}R_1\).

Después de eliminar \(p_1\) en filas 2 y 3 (operaciones sobre las filas):

$$\begin{pmatrix}3 & 6 & 2 & | & 34\\0 & -17 & 4 & | & 1\\0 & \frac{4}{3} & \frac{13}{3} & | & \frac{76}{3} \end{pmatrix}$$

Nota: la fila 2 quedó con coeficiente (-17) en la columna de \(p_2\). La fila 3 tiene fracciones — eso es normal en eliminación.

Paso 2 — usar la fila 2 para eliminar la componente \(p_2\) de la fila 3.

Multiplicador para \(R_3\): \(f=\frac{\frac{4}{3}}{-17}=-\frac{4}{51}\).
Hacemos \(R_3 \leftarrow R_3 – f R_2 = R_3 + \frac{4}{51}R_2\).

Al realizar esa operación la segunda columna de \(R_3\) queda 0 y obtenemos una ecuación sólo en \(p_3\):

$$\begin{pmatrix}3 & 6 & 2 & | & 34 \\ 0 & -17 & 4 & | & 1 \\ 0 & 0 & \frac{241}{51} & | & \frac{1296}{51} \end{pmatrix}$$

Paso 3 — despejar \(p_3\) y luego usar sustitución hacia atrás.

De la tercera fila:

$$\frac{241}{51},p_3 = \frac{1296}{51} \quad\Rightarrow\quad p_3 = \frac{1296}{241}$$

Sustituyendo \(p_3\) en la segunda fila:

$$-17,p_2 + 4,p_3 = 1 \quad\Rightarrow\quad p_2 = \frac{1 – 4p_3}{-17}$$

Sustituyendo el valor numérico de \(p_3\) (o directamente resolviendo con las fracciones) se obtiene \(p_2=3\).

Finalmente sustituimos \(p_2\) y \(p_3\) en la primera ecuación para obtener (p_1), que da (p_1=2).

Conclusión del método manual (eliminación + sustitución):

$$\mathbf p = \begin{pmatrix}2 \\ 3 \\ 5\end{pmatrix}$$

¿Qué significa calcular la inversa \(M^{-1}\) y cómo se obtiene con matriz aumentada?

Otra forma más “directa” (y que explica el porqué de la eliminación) es calcular la inversa de \(M\). Si \(M\) es invertible entonces

$$\mathbf p = M^{-1}\mathbf c$$

Para encontrar \(M^{-1}\) se construye la matriz aumentada \([,M \mid I,]\) y se aplican operaciones elementales por filas hasta transformar el lado izquierdo en la identidad; entonces el lado derecho habrá quedado como \(M^{-1}\).

Verificaciones numéricas y redondeo

En la práctica, cuando calculas \(M^{-1}M\) en ordenador verás que los elementos fuera de la diagonal no son exactamente cero sino números muy pequeños (por ejemplo \(-5\times 10^{-17})\). Eso es error numérico de punto flotante. Si redondeas (por ejemplo a 10 decimales) obtendrás la matriz identidad exacta en la presentación.

Resumen conceptual

• Eliminación gaussiana: reduce el sistema a triangular superior (eliminación hacia adelante) y luego despeja las variables (sustitución hacia atrás).
• Inversa por matriz aumentada: hacer operaciones elementales sobre \([M\mid I]\) hasta convertir \(M\) en \(I\); lo que quede a la derecha será \(M^{-1}\).
• Comprobación: puedes verificar que \(M^{-1}\mathbf c = \mathbf p\) y \(M^{-1}M = I\) (hasta errores numéricos muy pequeños).

Código Python

import numpy as np

# Definición de la matriz y vectores (el sistema del ejemplo)
M = np.array([[3, 6, 2],
              [10, 3, 8],
              [1, 7, 5]], dtype=float)

c = np.array([34, 69, 48], dtype=float)

# 1) Resolver el sistema M p = c
p = np.linalg.solve(M, c)
print("Solución (np.linalg.solve):", p)

# 2) Calcular la inversa de M
M_inv = np.linalg.inv(M)
print("\nInversa M^{-1}:\n", M_inv)

# 3) Verificar que M^{-1} * c = p
p_from_inv = M_inv @ c
print("\nM^{-1} @ c =", p_from_inv)

# 4) Verificar que M_inv @ M ≈ I (mostramos la matriz y una versión redondeada)
I_approx = M_inv @ M
print("\nM^{-1} @ M (aprox):\n", I_approx)
print("\nM^{-1} @ M (redondeado a 10 decimales):\n", np.round(I_approx, 10))

import numpy as np

# Definición de la matriz y vectores (el sistema del ejemplo)
M = np.array([[3, 6, 2],
              [10, 3, 8],
              [1, 7, 5]], dtype=float)

c = np.array([34, 69, 48], dtype=float)

# 1) Resolver el sistema M p = c
p = np.linalg.solve(M, c)
print("Solución (np.linalg.solve):", p)

# 2) Calcular la inversa de M
M_inv = np.linalg.inv(M)
print("\nInversa M^{-1}:\n", M_inv)

# 3) Verificar que M^{-1} * c = p
p_from_inv = M_inv @ c
print("\nM^{-1} @ c =", p_from_inv)

# 4) Verificar que M_inv @ M ≈ I (mostramos la matriz y una versión redondeada)
I_approx = M_inv @ M
print("\nM^{-1} @ M (aprox):\n", I_approx)
print("\nM^{-1} @ M (redondeado a 10 decimales):\n", np.round(I_approx, 10))

Salida:

Solución (np.linalg.solve): [2. 3. 5.]

Inversa M^{-1}:
 [[ 0.17012448  0.06639004 -0.17427386]
 [ 0.17427386 -0.05394191  0.01659751]
 [-0.2780083   0.06224066  0.21161826]]

M^{-1} @ c = [2. 3. 5.]

M^{-1} @ M (aprox):
 [[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [-1.56125113e-16  1.00000000e+00 -5.55111512e-17]
 [ 5.55111512e-17  0.00000000e+00  1.00000000e+00]]

M^{-1} @ M (redondeado a 10 decimales):
 [[ 1.  0.  0.]
 [-0.  1. -0.]
 [ 0.  0.  1.]]

Solución (np.linalg.solve): [2. 3. 5.]

Inversa M^{-1}:
 [[ 0.17012448  0.06639004 -0.17427386]
 [ 0.17427386 -0.05394191  0.01659751]
 [-0.2780083   0.06224066  0.21161826]]

M^{-1} @ c = [2. 3. 5.]

M^{-1} @ M (aprox):
 [[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [-1.56125113e-16  1.00000000e+00 -5.55111512e-17]
 [ 5.55111512e-17  0.00000000e+00  1.00000000e+00]]

M^{-1} @ M (redondeado a 10 decimales):
 [[ 1.  0.  0.]
 [-0.  1. -0.]
 [ 0.  0.  1.]]

Qué observar al ejecutar:

• p debe salir como [2. 3. 5.].
• M_inv @ c debe coincidir con p (igual numéricamente, dentro de tolerancias de punto flotante).
• M_inv @ M será numéricamente la identidad; si ves valores muy pequeños fuera de la diagonal, es normal: redondea para comprobar que esos valores son efectivamente cero dentro de la precisión de la máquina.

Resumen

La eliminación y la inversa son dos caras de la misma moneda: la primera resuelve un sistema concreto, la segunda construye la “función inversa” que, aplicada a cualquier vector de costes, devolvería los precios.

En aplicaciones (ciencia de datos, econometría, ingeniería), normalmente no calculamos la inversa si solo queremos resolver un sistema: usamos solve (algoritmos de factorización) por razones de estabilidad numérica y eficiencia. Sin embargo, obtener la inversa es útil para entender conceptualmente la transformada inversa y para comprobaciones.

El Formato Matricial

La expresión general:

$$y=a_1​x_1​+a_2​x_2​+⋯+a_n​x_n​$$

Es el punto de partida que lleva directamente a la forma matricial de los sistemas lineales. El mismo sistema de ecuaciones de las frutas puede escribirse y resolverse en forma matricial. En la ciencia de datos los datasets son matrices, y resolver sistemas lineales ayuda a encontrar relaciones entre variables.

Partimos del sistema:

$$\begin{cases} 2x + 3y = 30 \\ x + y = 12 \end{cases}$$

En formato matricial:

$$\begin{bmatrix} 2 & 3 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 30 \\ 12 \end{bmatrix}$$

Significado de las partes:

La forma vectorial compacta del sistema matricial se expresa como:

$$A\vec{x} = \vec{b}$$

Donde:

• \(A\) –> es la matriz de coeficientes
• \(\vec{x}\) –> es el vector de incógnitas
• \(\vec{b}\) –> es el vector de términos independientes

Por tanto:

$$A = \begin {bmatrix} 2 & 3 \\ 1 & 1 \end{bmatrix}, \vec{x} = \begin{bmatrix} x \\ y \end{bmatrix}, \vec{b} = \begin{bmatrix} 30 \\ 12 \end{bmatrix}$$

En ciencia de datos y machine learning esta forma es fundamental porque:

• Permite expresar relaciones entre miles de variables en una sola ecuación.
• Es la base de la regresión lineal, la reducción de dimensionalidad, y el entrenamiento de modelos.
• Toda la manipulación se hace con álgebra matricial, que es más eficiente computacionalmente.

Resolviendo problemas reales.

En artículos anteriores vimos operaciones sobre vectores y matrices desde la perspectiva matemática. Vamos a ver cómo esas mismas operaciones pueden resolver problemas reales, como descubrir cuánto cuesta cada producto de un supermercado.

Escenario

Volvemos a nuestro supermercado de frutas. Imaginemos que tenemos tres personas:
Bob, Alice y Tim, que han ido al supermercado y compraron manzanas, naranjas y peras. Podemos representar sus compras en una matriz \(M\):

$$M = \begin{bmatrix} 3 & 6 & 2 \\ 10 & 3 & 8 \\ 1 & 7 & 5 \end{bmatrix}$$

Cada fila representa una persona, y cada columna representa un producto:

El coste total

Sabemos cuánto pagó cada persona:

$$C =\begin{bmatrix} 34 \\ 69 \\ 48 \end{bmatrix}$$

Y queremos averiguar cuánto cuesta cada producto. Presentemos los importes que gasto cada uno en un vector que llamamos \( P \) , vector de precios:

$$P = \begin{bmatrix} p_1 \\ p_2 \\ p_3 \end{bmatrix} $$

Donde:

• \(p_1\) –> precio de la manzana
• \(p_2\) –> precio de la naranja
• \(p_3\) –> precio de la pera

Sabemos que:

$$M \cdot P = C$$

o, dicho en palabras: la cantidad comprada multiplicada por el precio unitario da el coste total.

Escribiendo las ecuaciones

Multiplicar la matriz \(M\) por el vector \(P\) nos da tres ecuaciones:

$$ \begin{cases} 3p_1 + 6p_2 + 2p_3 = 34 \\ 10p_1 + 3p_2 + 8p_3 = 69 \\ 1p_1 + 7p_2 + 5p_3 = 48 \end{cases}$$

Esto es un sistema de ecuaciones lineales con tres incógnitas: \( p_1, p_2, p_3 \).

Representar el sistema de ecuaciones en formato matricial:

$$ \begin{bmatrix} 3 & 6 & 2 \\ 10 & 3 & 8 \\ 1 & 7 & 5 \end{bmatrix} \cdot \begin{bmatrix} p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} 34 \\ 69 \\ 48 \end{bmatrix} $$

Recordando, tenemos aqui:

• Matriz de coeficientes: contiene los números que multiplican a las variables (las cantidades de las frutas).
• Vector de incógnitas: las variables que queremos encontrar (precio de las frutas).
• Vector de resultados: los valores en el lado derecho (total de dinero gastado).

Cómo transformar el sistema para poder resolverlo

Hasta ahora has visto qué es un sistema de ecuaciones y cómo representarlo en forma matricial. Pero hay un problema práctico: tal como están escritas, las ecuaciones no siempre resultan fáciles de resolver.

En otros articulos hablamos de las caracteristicas del un sistema, como la singularidad. Aqui veremos como transformar el sistema sin cambiar su significado, haciendo mas simple su solucion.

Transformar sin perder información

Existe una serie de operaciones que podemos aplicar a las ecuaciones sin alterar sus soluciones:

• Intercambiar ecuaciones
• Multiplicar o dividir una ecuación por un número distinto de cero
• Sumar o restar ecuaciones entre sí

Estas operaciones son la base del Algoritmo de eliminación de Gauss. Su objetivo es reorganizar el sistema en una forma mucho más sencilla de interpretar.

Ejemplo conceptual del proceso

Partimos de un sistema cualquiera y aplicamos operaciones de fila para:

1. Elegir un elemento clave (pivote)
2. Convertir en cero todo lo que hay debajo
3. Repetir el proceso en las siguientes filas

El resultado es una matriz en la que la estructura del sistema se vuelve transparente.

Es como resolver un sistema de ecuaciones común por eliminación, pero aplicado directamente a los números de una matriz (llamada matriz de coeficientes). El objetivo es limpiar la matriz para que las respuestas sean evidentes.

Forma escalonada por filas (row echelon form).

Se consigue eliminando la información redundante y aislando la información útil. Cada paso reduce el sistema a algo más simple, manteniendo exactamente las mismas soluciones.

Una matriz está en esta forma cuando:

• Cada fila empieza más a la derecha que la anterior
• Debajo de cada elemento clave (pivote) hay ceros
• Las filas sin información (todos ceros) quedan al final

$$A = \begin{bmatrix}1 & * & * \\ 0 & 1 & * \\ 0 & 0 & 1 \end{bmatrix}$$

• Es útil porque te deja el sistema listo para “sustitución hacia atrás”.

Forma escalonada reducida: la solución directa

Si continúas el proceso un paso más, llegas a la forma escalonada reducida (reduced row echelon form) Aquí ocurre algo muy potente, en su columna solo hay ceros (arriba y abajo)

Cada pivote es 1

$$A = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

• Es el nivel máximo de limpieza.
• Aquí la solución es directa, cada fila te dice el valor de cada variable.

¿Qué pasa si el sistema es “Singular”?

A veces, al hacer la eliminación gaussiana, una fila completa se convierte en ceros \( 0 = 0 \)). Esto significa que:

• El sistema es singular (las ecuaciones eran dependientes o redundantes).
• En la matriz, verás un 0 en la diagonal principal en lugar de un 1.
• Esto indica que el sistema puede tener infinitas soluciones o ninguna.

Método de eliminación gaussiana

La eliminación gaussiana (o reducción de filas) es el proceso de manipular una matriz mediante operaciones matemáticas simples para simplificarla y resolver sistemas de ecuaciones.

Pasos en la solución.

1 – Crear la matriz aumentada:

$$
\left(
\begin{array}{ccc|c}
3 & 6 & 2 & 34 \\
10 & 3 & 8 & 69 \\
1 & 7 & 5 & 48
\end{array}
\right)
$$

2 – Elegir un buen pivote y hacer intercambio de filas.

En este ejemplo la fila 3 tiene coeficiente 1 y la cambiamos por la 1.

$$
F_1 \leftrightarrow F_3
$$

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
10 & 3 & 8 & 69 \\
3 & 6 & 2 & 34
\end{array}
\right)
$$

3- Eliminación columna 1

Queremos eliminar las variables \( p_1 \) de las filas 2 y 3.

Buscamos eliminar las \(p_1\), lo primero es restar la fila 1 con la fila 2. Para esto hacemos una normalización de la fila dos, buscando que se elimine la variable. En notación de matrices se expone como:

$$
F_2 \leftarrow F_2 – 10F_1
$$

Lo que quiere decir es que la nueva fila \( F_2 \) será el resultado de la resta de la fila \( F_2 \) y la fila \( F_1 \) multiplicada por 10. Esto hace que se elimine \(p_1\) en la operación \(10p_1 – 10p_1\).

La segunda eliminación de la variable es creando una nueva fila 3 producto de la resta entre las filas \( F_3 \) y la \( F_1 \) multiplicando esta por 3.

$$
F_3 \leftarrow F_3 – 3F_1
$$

La matriz resultante es:

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
0 & -67 & -42 & -411 \\
0 & -15 & -13 & -110
\end{array}
\right)
$$

4 – Normalización pivote 2

De la matriz actual buscamos reducir \(p_2\). Normalizamos la fila 2, dividiéndola entre 67, llevando \(p_2\) a 1.

$$
F_2 \leftarrow \frac{1}{-67}F_2
$$

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
0 & 1 & \frac{42}{67} & \frac{411}{67} \\
0 & -15 & -13 & -110
\end{array}
\right)
$$

5 – Eliminación columna 2

Normalizamos la \(F_3\), dividiendo entre 15, y la restamos con \(F_2\) eliminando la variable.

$$
F_3 \leftarrow F_3 + 15F_2
$$

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
0 & 1 & \frac{42}{67} & \frac{411}{67} \\
0 & 0 & \frac{-241}{67} & \frac{-1205}{67}
\end{array}
\right)
$$

6 – Normalización pivote 3

Multiplicamos \(F_3\) por \( \frac{-67}{241} \) para llevar \(p_3\) a 1.

$$
F_3 \leftarrow \frac{-67}{241}F_3
$$

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
0 & 1 & \frac{42}{67} & \frac{411}{67} \\
0 & 0 & 1 & \frac{1205}{241}
\end{array}
\right)
$$

7 – Sustitución hacia atrás

Calcular el valor de \(p_3\).

$$
p_3 = \frac{1205}{241} = 5
$$

Sustituir \(p_3\) en \(F_2\)

$$
p_2 + \frac{42}{67} \cdot 5 = \frac{411}{67}
$$

$$
p_2 = 3
$$

Sustituir \(p_3\) y \(p_2\) en \(F_1\)

$$
p_1 + 7 \cdot 3 + 5 \cdot 5 = 48
$$

$$
p_1 = 2
$$

Resultado final

$$
\begin{cases}
p_1 = 2 \\
p_2 = 3 \\
p_3 = 5
\end{cases}
$$

Paso 4 — verificación (comprobación directa)

Multiplicamos la matriz por el vector de precios para ver si obtenemos los costes:

$$M = \begin{bmatrix} 3 & 6 & 2 \\ 10 & 3 & 8 \\ 1 & 7 & 5 \end{bmatrix}, p = \begin{bmatrix} 2 \\ 3 \\ 5 \end{bmatrix}$$

Producto \(Mp\):

• Fila 1: \( 3 \cdot 2 + 6 \cdot 3 + 2 \cdot 5 = 34\)
• Fila 2: \( 10 \cdot 2 + 3 \cdot 3 + 8 \cdot 5 = 69\)
• Fila 3: \( 1 \cdot 2 + 7 \cdot 3 + 5 \cdot 5 = 48\)

Método de Gauss-Jordan: llevar el sistema hasta la matriz identidad

Una vez hemos aplicado la sustitución hacia atrás y ya conocemos los valores de las variables, puede surgir una pregunta lógica: ¿Y si en lugar de calcular las variables “a mano”, dejamos la matriz ya resuelta automáticamente? Aquí es donde entra el método de Gauss-Jordan.

¿Qué aporta Gauss-Jordan?

El método de Gauss (el que hemos usado) termina en una matriz escalonada, y a partir de ahí necesitamos hacer sustitución hacia atrás. El método de Gauss-Jordan va un paso más allá. Sigue aplicando operaciones hasta convertir la matriz en la matriz identidad

💡 En ese momento:

• Ya no necesitas sustituir
• Los valores de las variables aparecen directamente en la última columna

¿Por qué es útil llegar a la identidad?

Porque transforma el sistema en algo trivial:

• Cada fila representa directamente una variable
• El sistema queda completamente desacoplado
• Es el método que se usa en:
  • cálculo matricial
  • inversa de matrices
  • implementaciones en programación

Continuamos el ejemplo

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 5 & 48 \\
0 & 1 & \frac{42}{67} & \frac{411}{67} \\
0 & 0 & 1 & \frac{1205}{241}
\end{array}
\right)
$$

Paso 1: Eliminar la columna 3 de la fila 2 usando la fila 3:

Operación \( F_2 \rightarrow F_2 – \frac{42}{67}F_3\)

• Columna 3:

$$ \frac{42}{67} – \frac{42}{67} \cdot 1 = 0$$

• Termino independiente:

$$ \frac{411}{67} – \frac{42}{67} \cdot 1 \frac{1205}{241}$$

$$ \frac{42 \cdot 1205}{67 \cdot 241} = \frac{50610}{16147}$$

Ahora llevamos \( \frac{411}{67}\) al mismo denominador:

$$ \frac{411}{67} = \frac{411 \cdot 241}{67 \cdot 241} = \frac{99051}{16147}$$

Restamos:

$$ \frac{99051}{16147} – \frac{50610}{16147} = \frac{48441}{16471} = 3$$

Nueva \(F_2\):

$$\left(\begin{array}{ccc|c}0 & 1 & 0 & 3 \end{array}\right)$$

Paso 2: eliminar la columna 3 de la fila 1 usando la fila 3

Operación \( F_1 \rightarrow F_1 – 5F_3\)

• Columna 3:

$$5 – 5 = 0$$

• Termino independiente:

$$ 48 – 5 \cdot \frac{1205}{241}$$

$$ 5 \cdot \frac{1205}{241} = \frac{6025}{241}$$

Pasamos 48 a fracción:

$$ 48 = \frac{11568}{241}$$

Restamos:

$$ \frac{11568}{241} – \frac{6025}{241} = \frac{5543}{241} = 23$$

Calculamos el término independiente de la fila 3 que no se usará nuevamente:

$$ \frac{1205}{241} = 5$$

Matriz tras paso 1:

$$
\left(
\begin{array}{ccc|c}
1 & 7 & 0 & 23 \\
0 & 1 & 0 & 3 \\
0 & 0 & 1 & 5
\end{array}
\right)
$$

PASO 2: eliminar la columna 2 (usando la fila 2)

Operación \( F_1 \rightarrow F_1 – 7F_2\)

Columna 2:

$$7 – 7 = 0$$

Termino independiente:

$$23 – 7 \cdot 3 = 2$$

Matriz de identidad:

$$
\left(
\begin{array}{ccc|c}
1 & 0 & 0 & 2 \\
0 & 1 & 0 & 3 \\
0 & 0 & 1 & 5
\end{array}
\right)
$$

Solucion final

$$
\begin{cases}
p_1 = 2 \\
p_2 = 3 \\
p_3 = 5
\end{cases}
$$

Solución con Python

Para solucionar con python usamos numpy.linalg.solve, que aplica internamente eliminación de Gauss con pivotado, y es la forma más eficiente y segura de resolver sistemas lineales \(A \cdot x = b \) en Python.

import numpy as np

# Matriz de coeficientes
A = np.array([
    [3, 6, 2],
    [10, 3, 8],
    [1, 7, 5]
], dtype=float)

# Vector de resultados
b = np.array([34, 69, 48], dtype=float)

# Resolver el sistema A·x = b
sol = np.linalg.solve(A, b)

# Mostrar resultados
p1, p2, p3 = sol
print(f"p1 = {p1:.2f}, p2 = {p2:.2f}, p3 = {p3:.2f}")

import numpy as np

# Matriz de coeficientes
A = np.array([
    [3, 6, 2],
    [10, 3, 8],
    [1, 7, 5]
], dtype=float)

# Vector de resultados
b = np.array([34, 69, 48], dtype=float)

# Resolver el sistema A·x = b
sol = np.linalg.solve(A, b)

# Mostrar resultados
p1, p2, p3 = sol
print(f"p1 = {p1:.2f}, p2 = {p2:.2f}, p3 = {p3:.2f}")

Salida:

p1 = 2.00, p2 = 3.00, p3 = 5.00

p1 = 2.00, p2 = 3.00, p3 = 5.00

En el artículo anterior se trabajo los autovalores y autovectores, es común quedarse con la sensación de que son conceptos abstractos. Pero en realidad, están detrás de muchas de las herramientas más potentes del análisis de datos moderno.

Una de ellas es el Análisis de Componentes Principales (PCA), una técnica que permite comprimir información sin perder lo esencial.

En este artículo veremos cómo los eigenvectores sirven para reducir la dimensionalidad de los datos, y cómo este principio se aplica a algo tan visual como la compresión de imágenes de rostros humanos.

Cuando trabajamos con datos reales —imágenes, audio, sensores o registros de clientes— solemos tener miles de variables. Eso significa que el procesamiento se vuelve lento, costoso y muchas veces innecesario, porque buena parte de los datos son redundantes.

El PCA (Principal Component Analysis) se basa precisamente en encontrar las direcciones de mayor variabilidad dentro de los datos. Y esas direcciones no son otras que los autovectores de la matriz de covarianza del conjunto.

Visualizando el concepto con una imagen

Una imagen digital no es más que una matriz de números, donde cada valor representa la intensidad del color (por ejemplo, del blanco al negro si es en escala de grises).

Para este ejemplo usaremos el famoso dataset de rostros de celebridades incluido en scikit-learn.

from sklearn.datasets import fetch_lfw_people
import matplotlib.pyplot as plt

# Cargamos el dataset de caras
faces = fetch_lfw_people(min_faces_per_person=60)
image = faces.images[0]

plt.imshow(image, cmap='bone')
plt.title(f"Ejemplo de imagen: {faces.target_names[0]}")
plt.axis('off')
plt.show()

from sklearn.datasets import fetch_lfw_people
import matplotlib.pyplot as plt

# Cargamos el dataset de caras
faces = fetch_lfw_people(min_faces_per_person=60)
image = faces.images[0]

plt.imshow(image, cmap='bone')
plt.title(f"Ejemplo de imagen: {faces.target_names[0]}")
plt.axis('off')
plt.show()

Cada imagen tiene 62 filas × 47 columnas, es decir, 2.914 píxeles o variables.
PCA nos permitirá representar esta misma información con muchas menos variables.

Aplicando PCA paso a paso

PCA consiste básicamente en tres fases:

1. Calcular la matriz de covarianza del conjunto de datos.
2. Obtener sus autovectores y autovalores.
3. Proyectar los datos sobre las direcciones asociadas a los autovalores más grandes, que contienen la mayor parte de la varianza.

En código, podemos hacerlo así:

from sklearn.decomposition import PCA
import numpy as np

# Aplanamos las imágenes 2D a vectores 1D
X = faces.data

# Aplicamos PCA para conservar solo 15 componentes principales
pca = PCA(n_components=15)
X_pca = pca.fit_transform(X)

# Reconstruimos las imágenes desde los componentes principales
X_reconstructed = pca.inverse_transform(X_pca)

# Mostramos la imagen original y la reconstruida
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
axes[0].imshow(X[0].reshape(62, 47), cmap='bone')
axes[0].set_title("Imagen original")
axes[0].axis('off')

axes[1].imshow(X_reconstructed[0].reshape(62, 47), cmap='bone')
axes[1].set_title("Reconstruida (15 componentes)")
axes[1].axis('off')

plt.show()

from sklearn.decomposition import PCA
import numpy as np

# Aplanamos las imágenes 2D a vectores 1D
X = faces.data

# Aplicamos PCA para conservar solo 15 componentes principales
pca = PCA(n_components=15)
X_pca = pca.fit_transform(X)

# Reconstruimos las imágenes desde los componentes principales
X_reconstructed = pca.inverse_transform(X_pca)

# Mostramos la imagen original y la reconstruida
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
axes[0].imshow(X[0].reshape(62, 47), cmap='bone')
axes[0].set_title("Imagen original")
axes[0].axis('off')

axes[1].imshow(X_reconstructed[0].reshape(62, 47), cmap='bone')
axes[1].set_title("Reconstruida (15 componentes)")
axes[1].axis('off')

plt.show()

¿Qué estamos haciendo en realidad?

Cada componente principal es un autovector de la matriz de covarianza.
Cada uno representa una dirección en el espacio de datos donde la información cambia más.
Los autovalores indican cuánta “varianza” explica cada componente.

En este caso, al reconstruir la imagen con solo 15 componentes, vemos una versión algo borrosa, pero perfectamente reconocible. Hemos pasado de 2.914 variables a solo 15, manteniendo casi toda la información útil.

Si reducimos a 5 componentes, la imagen pierde detalle; si aumentamos a 50, se ve casi igual que la original. Todo depende del equilibrio que queramos entre tamaño y fidelidad.

Las “eigenfaces”

Una de las formas más visuales de entender PCA en imágenes es observar los autovectores directamente. En este caso se les conoce como eigenfaces, o “caras propias”.

Cada eigenface representa un patrón característico: sombras de ojos, contorno del rostro, luz lateral, expresión, etc.

Podemos verlos así:

fig, axes = plt.subplots(3, 5, figsize=(9, 6))
for i, ax in enumerate(axes.flat):
    ax.imshow(pca.components_[i].reshape(62, 47), cmap='bone')
    ax.set_title(f"Comp {i+1}")
    ax.axis('off')
plt.suptitle("Los 15 eigenfaces principales")
plt.show()

fig, axes = plt.subplots(3, 5, figsize=(9, 6))
for i, ax in enumerate(axes.flat):
    ax.imshow(pca.components_[i].reshape(62, 47), cmap='bone')
    ax.set_title(f"Comp {i+1}")
    ax.axis('off')
plt.suptitle("Los 15 eigenfaces principales")
plt.show()

Los eigenvectores son, en el fondo, una forma de entender la estructura interna de los datos.

Verás que, aunque una transformación lineal cambia la dirección de la mayoría de los vectores, existen algunos vectores especiales que no rotan, solo se estiran o encogen.
Esos vectores son los autovectores (eigenvector), y la cantidad por la que se estiran o encogen se llama autovalor (eigenvalues).

Intuición: vectores que no cambian de dirección

Imagina una matriz que representa una transformación lineal —por ejemplo, una matriz de “shear” (cizallamiento), que deforma un cuadrado en un paralelogramo.
En la mayoría de los casos, todos los vectores cambian de dirección al aplicar la transformación… excepto algunos muy particulares. Esos vectores especiales mantienen su dirección (aunque pueden cambiar su longitud). Por eso decimos que son autovectores.

Ejemplo de una transformación lineal:

import numpy as np

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

square = np.array([[0, 0, 1, 1, 0],
                   [0, 1, 1, 0, 0]])
transformed = M @ square

import numpy as np

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

square = np.array([[0, 0, 1, 1, 0],
                   [0, 1, 1, 0, 0]])
transformed = M @ square

fig, ax = plt.subplots(figsize=(8,8))
ax.plot(square[0], square[1], 'b--', label='Original')
ax.plot(transformed[0], transformed[1], 'r-', label=f'Shear con b={b}')

ax.set_xlim(-1, 3)
ax.set_ylim(-1, 3)
ax.axhline(0, color='black')
ax.axvline(0, color='black')
ax.set_aspect('equal')
ax.legend()
plt.show()

fig, ax = plt.subplots(figsize=(8,8))
ax.plot(square[0], square[1], 'b--', label='Original')
ax.plot(transformed[0], transformed[1], 'r-', label=f'Shear con b={b}')

ax.set_xlim(-1, 3)
ax.set_ylim(-1, 3)
ax.axhline(0, color='black')
ax.axvline(0, color='black')
ax.set_aspect('equal')
ax.legend()
plt.show()

Definición

Un vector \( \vec{v} \) es un autovector de una matriz \( M \) si, al aplicar la matriz transformadora sobre él, el resultado es el mismo vector escalado por un número \( \lambda \):

$$M \vec{v} = \lambda \vec{v}$$

Donde:

• \( M \) → es una matriz cuadrada \(n \times n\)
• \( \vec{v} \) → es un autovector
• \( \lambda \) → es el autovalor asociado

Este número \( \lambda \) indica cuánto se estira o encoge el vector cuando la matriz actúa sobre él.

Interpretación geométrica

• Si \( \lambda > 1 \): el vector se alarga.
• Si \( 0 < \lambda < 1 \): el vector se encoge.
• Si \( \lambda < 0 \): el vector cambia de sentido (se invierte).
• Si \( \lambda = 1 \): el vector queda igual.

Calculando los eigenvectors y eigenvalues

Continuamos con el ejemplo inicial:

import numpy as np

# Calcular eigenvalores y eigenvectores
w, V = np.linalg.eig(M)

print("\nEigenvalores:\n", w)
print("\nEigenvectores (columnas):\n", V)

import numpy as np

# Calcular eigenvalores y eigenvectores
w, V = np.linalg.eig(M)

print("\nEigenvalores:\n", w)
print("\nEigenvectores (columnas):\n", V)

Salida:

Eigenvalores:
 [1.38196601 3.61803399]

Eigenvectores (columnas):
 [[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

Eigenvalores:
 [1.38196601 3.61803399]

Eigenvectores (columnas):
 [[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8,8))

# Cuadrado original (azul) y transformado (rojo)
ax.plot(square[0], square[1], 'b--', label='Original (cuadrado unidad)')
ax.plot(transformed[0], transformed[1], 'r-', label='Transformado M·S')

# Ejes
ax.axhline(0, color='black', linewidth=1)
ax.axvline(0, color='black', linewidth=1)

# Eigenvectores
for i in range(len(w)):
    v = V[:, i]            # vector propio
    Mv = M @ v             # su transformación
    λ = w[i]

    # Dibuja el eigenvector (azul)
    ax.quiver(0, 0, v[0], v[1], angles='xy', scale_units='xy', scale=1,
              color='blue', width=0.007, label=fr'$v_{i+1}$')

    # Dibuja el vector transformado (rojo)
    ax.quiver(0, 0, Mv[0], Mv[1], angles='xy', scale_units='xy', scale=1,
              color='red', width=0.007, label=fr'$M v_{i+1} = \lambda_{i+1} v_{i+1}$')

    # Muestra texto con el eigenvalor
    ax.text(Mv[0]*1.1, Mv[1]*1.1, fr'$\lambda_{i+1}={λ:.2f}$', color='red', fontsize=10)

ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_aspect('equal', adjustable='box')
ax.legend(loc='upper left', facecolor='none', edgecolor='none')
ax.set_title("Transformación lineal y Eigenvectores", fontsize=14)
ax.set_xlabel("x")
ax.set_ylabel("y")
plt.grid(True)
plt.show()

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8,8))

# Cuadrado original (azul) y transformado (rojo)
ax.plot(square[0], square[1], 'b--', label='Original (cuadrado unidad)')
ax.plot(transformed[0], transformed[1], 'r-', label='Transformado M·S')

# Ejes
ax.axhline(0, color='black', linewidth=1)
ax.axvline(0, color='black', linewidth=1)

# Eigenvectores
for i in range(len(w)):
    v = V[:, i]            # vector propio
    Mv = M @ v             # su transformación
    λ = w[i]

    # Dibuja el eigenvector (azul)
    ax.quiver(0, 0, v[0], v[1], angles='xy', scale_units='xy', scale=1,
              color='blue', width=0.007, label=fr'$v_{i+1}$')

    # Dibuja el vector transformado (rojo)
    ax.quiver(0, 0, Mv[0], Mv[1], angles='xy', scale_units='xy', scale=1,
              color='red', width=0.007, label=fr'$M v_{i+1} = \lambda_{i+1} v_{i+1}$')

    # Muestra texto con el eigenvalor
    ax.text(Mv[0]*1.1, Mv[1]*1.1, fr'$\lambda_{i+1}={λ:.2f}$', color='red', fontsize=10)

ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_aspect('equal', adjustable='box')
ax.legend(loc='upper left', facecolor='none', edgecolor='none')
ax.set_title("Transformación lineal y Eigenvectores", fontsize=14)
ax.set_xlabel("x")
ax.set_ylabel("y")
plt.grid(True)
plt.show()

Matrices simétricas

Cuando una matriz es simétrica (es decir, \( M = M^T \)), se cumple algo muy útil:

• Todos sus autovectores son perpendiculares entre sí.
• Sus autovalores son reales.

Por ejemplo:

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

autovalores, autovectores = np.linalg.eig(A)
print("Autovalores:", autovalores)
print("Autovectores:\n", autovectores)

print("Producto punto entre los autovectores:", np.dot(autovectores[:,0], autovectores[:,1]))

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

autovalores, autovectores = np.linalg.eig(A)
print("Autovalores:", autovalores)
print("Autovectores:\n", autovectores)

print("Producto punto entre los autovectores:", np.dot(autovectores[:,0], autovectores[:,1]))

La salida mostrará que el producto punto es 0, confirmando que los autovectores son ortogonales.

Autovalores: [5. 3.]
Autovectores:
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
Producto punto entre los autovectores: 0.0

Autovalores: [5. 3.]
Autovectores:
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
Producto punto entre los autovectores: 0.0

Aplicaciones en ciencia de datos

Los eigenvectores y eigenvalores son fundamentales en data science:

Análisis de Componentes Principales (PCA)

• Qué hace: Reduce dimensiones de un dataset encontrando los eigenvectores de la matriz de covarianza.
• Ejemplo: Dataset Iris (150×4) → eigenvectores indican las direcciones de mayor variación (ej: longitud vs ancho de pétalo).

Redes Neuronales

• Eigenvectores optimizan transformaciones en capas de redes profundas.
• Ejemplo: Ajustar pesos para maximizar patrones.

Procesamiento de Imágenes

• Eigenvectores en eigenfaces identifican características faciales clave.
• Ejemplo: Reconocimiento facial en apps de seguridad.

Sistemas Dinámicos

• Modelan cómo cambian datos con el tiempo (ej: predicciones financieras).
• Eigenvalores determinan la estabilidad de los sistemas.

Historia real: En 2024, Spotify usó eigenvectores en su sistema de recomendaciones para identificar patrones en gustos musicales, mejorando la precisión en un 10%.

Cómo encontrar los Eigenvalues y Eigenvectors

Para cualquier matriz \(M\):

1. Primero encontramos todos los eigenvalores λ.
2. Luego, usando esos eigenvalores, encontramos los eigenvectores correspondientes.

La ecuación base \(M \vec{v} = \lambda \cdot \vec{v}\) será la clave de todo.

Pero para trabajar con ella algebraicamente, vamos a transformarla un poco.

La matriz identidad

Necesitamos la matriz especial llamada la matriz identidad, que se denota por \(I\). La matriz identidad es una matriz cuadrada con unos en la diagonal principal (de arriba a la izquierda hasta abajo a la derecha) y ceros en todas las demás posiciones. Ya la vimos como ejemplo anteriormente en el articulo de las matrices.

Por ejemplo, la matriz identidad de tamaño 3×3 es:

$$I = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix}$$

Y se llama “identidad” porque, al multiplicarla por cualquier matriz o vector, el resultado es exactamente el mismo.

Podemos decir que: \(I \cdot M = M\) y tambien que: \(I \cdot v = v\)

Manipulamos la ecuación para llegar al sistema característico

Volvamos a nuestra ecuación:

$$M \vec{v} = \lambda \cdot \vec{v}$$

Podemos escribir el lado derecho como \(\lambda I \cdot v\), ya que \(I \cdot v.\).

Entonces la ecuación queda así:

$$M \vec{v} = \lambda I \cdot \vec{v}$$

Restamos el término del lado derecho en ambos lados:

$$M \cdot \vec{v} \ – \lambda I \cdot \vec{v} = 0 $$

Sacamos vvv como factor común (como si fuera un número):

$$(M \ – \lambda I) \cdot \vec{v} = 0 $$

Y esta ecuación es fundamental: se conoce como la ecuación característica.

La condición de existencia de eigenvectores

Esta ecuación tiene dos tipos de soluciones posibles:

1. La trivial, cuando \(v=0\).
   Pero ese vector no nos dice nada (no tiene dirección ni longitud), así que lo descartamos.
2. La no trivial, cuando \(v≠0\).
   Para que exista esta solución, debe cumplirse una condición especial.

Recuerda una propiedad del álgebra lineal:

Una ecuación de la forma \(A \cdot v=0\) solo tiene una solución no trivial \(v≠0\) si la matriz A es singular, es decir, su determinante es igual a cero.

Aplicando esto a nuestra ecuación \((M \ – \lambda I) \cdot \vec{v} = 0 \), la condición será: \(det(M \ – \lambda I) \cdot \vec{v} = 0 \)

Esta ecuación nos permitirá encontrar los valores de \( \lambda \), los eigenvalores.

Interpretación geométrica: el área que se hace cero

Detrás de todo esto hay una idea visual muy poderosa.

Imagina que tienes una matriz \(M\) que transforma puntos del plano \(XY\). Puedes representarla como una transformación del cuadrado unitario (un cuadrado de lado 1 anclado en el origen).

Cuando multiplicas \(M\) por las coordenadas de los vértices de ese cuadrado, obtienes un nuevo cuadrilátero. En general, este nuevo cuadrilátero es un paralelogramo.

• Si \(M\) solo escala el cuadrado, obtienes un rectángulo.
• Si además deforma el cuadrado (por ejemplo, en un cizallamiento), obtienes un paralelogramo inclinado.

El área de ese paralelogramo es una medida de cuánto “estira” o “encoge” la transformación de \(M\).
Y el determinante de la matriz está directamente relacionado con esa área.

Cuando el determinante es cero, significa que la transformación ha colapsado el plano en una línea o un punto, es decir, ha perdido una dimensión. Eso ocurre exactamente cuando estamos en un eigenvalor, porque en ese caso el espacio queda “aplastado” en una dirección particular: la del eigenvector.

Ejemplo visual con una matriz 2×2

Imaginemos una matriz genérica \(M\) de 2×2:

$$M = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$

Y tomemos el cuadrado unitario cuyos vértices son (0,0), (1,0), (1,1) y (0,1).

Cuando aplicamos \(M\) sobre esos puntos, el cuadrado se transforma en un paralelogramo cuyos vértices vienen dados por los productos matriciales:

$$M \cdot S$$

donde \(S\) contiene las coordenadas de los vértices del cuadrado.
El área de ese nuevo paralelogramo se puede calcular directamente como:

$$\text{Área} = ad – bc$$

que es precisamente el determinante de \(M\).

Así que cuando el determinante \(ad−bc=0\), el paralelogramo “se aplasta” y pierde área. Geométricamente, eso es lo que ocurre cuando hay un eigenvalor: la transformación ya no conserva el área y “colapsa” el espacio en una dirección.

Cálculo con Python: encontrar los eigenvalores

Ejemplo con matriz 2×2.

$$M = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$$

M = np.array([[2, 1],
              [1, 2]])

M = np.array([[2, 1],
              [1, 2]])

Calculamos eigenvalores y eigenvectores

eigenvalores, eigenvectores = np.linalg.eig(M)

print("Matriz M:")
print(M)

print("\nEigenvalores:")
print(eigenvalores)

print("\nEigenvectores (por columnas):")
print(eigenvectores)

eigenvalores, eigenvectores = np.linalg.eig(M)

print("Matriz M:")
print(M)

print("\nEigenvalores:")
print(eigenvalores)

print("\nEigenvectores (por columnas):")
print(eigenvectores)

Eigenvalores:
[1.38196601 3.61803399]

Eigenvectores (por columnas):
[[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

Eigenvalores:
[1.38196601 3.61803399]

Eigenvectores (por columnas):
[[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

Ejemplo con matriz 3×3.

M3 = np.array([
    [4, -2, 1],
    [0, 3, -1],
    [0, 0, 2]
])

eigenvalores3, eigenvectores3 = np.linalg.eig(M3)

print("Matriz M3:")
print(M3)

print("\nEigenvalores:")
print(eigenvalores3)

print("\nEigenvectores (por columnas):")
print(eigenvectores3)

M3 = np.array([
    [4, -2, 1],
    [0, 3, -1],
    [0, 0, 2]
])

eigenvalores3, eigenvectores3 = np.linalg.eig(M3)

print("Matriz M3:")
print(M3)

print("\nEigenvalores:")
print(eigenvalores3)

print("\nEigenvectores (por columnas):")
print(eigenvectores3)

Matriz M3:
[[ 4 -2  1]
 [ 0  3 -1]
 [ 0  0  2]]

Eigenvalores:
[4. 3. 2.]

Eigenvectores (por columnas):
[[1.         0.89442719 0.33333333]
 [0.         0.4472136  0.66666667]
 [0.         0.         0.66666667]]

Matriz M3:
[[ 4 -2  1]
 [ 0  3 -1]
 [ 0  0  2]]

Eigenvalores:
[4. 3. 2.]

Eigenvectores (por columnas):
[[1.         0.89442719 0.33333333]
 [0.         0.4472136  0.66666667]
 [0.         0.         0.66666667]]