3 Autovalores, Autovectores y Descomposición Singular

Hasta ahora hemos resuelto $Ax=b$. Pero en modelización predictiva (ej. PageRank, PCA, estabilidad dinámica), nos interesa resolver: \[A \mathbf{v} = \lambda \mathbf{v}\] Donde $\lambda$ es el autovalor y $\mathbf{v}$ el autovector. Para matrices gigantes ($N > 10^5$), no podemos usar eigen(A) de R base. Necesitamos algoritmos iterativos.

¿Por qué “Eigen”?

La palabra eigen es alemana y significa “propio” o “característico”. El término fue popularizado por David Hilbert en 1904, aunque el concepto se remonta a los trabajos de Leonhard Euler sobre la rotación de cuerpos rígidos y de Augustin-Louis Cauchy sobre las formas cuadráticas. En español, a veces se les llama “valores propios” o “autovalores”.

3.1 El método de las potencias (the power method)

Este método busca encontrar el autovalor dominante (el de mayor valor absoluto, $|\lambda_1|$). Es fundamental porque es el motor matemático detrás del algoritmo PageRank de Google.

3.1.1 Algoritmo

Dado un vector inicial aleatorio $\mathbf{b}_0$: \[\mathbf{b}_{k+1} = \frac{A \mathbf{b}_k}{||A \mathbf{b}_k||}\] La sucesión converge al autovector dominante.

power_method <- function(A, epsilon = 1e-6, max_iter = 1000) {
  # 1. Vector inicial aleatorio (normalizado)
  n <- nrow(A)
  x <- rnorm(n)
  x <- x / sqrt(sum(x^2))
  
  lambda_old <- 0
  
  for(i in 1:max_iter) {
    # 2. Multiplicación (aprovecha la dispersión de A si es sparse)
    x_new <- as.vector(A %*% x)
    
    # 3. Estimación del autovalor (Cociente de Rayleigh simple)
    # Si x está normalizado, lambda es aproximadamente t(x) %*% A %*% x
    # Aquí simplificamos usando la norma del nuevo vector antes de normalizar
    lambda <- sqrt(sum(x_new^2)) 
    
    # Recuperar signo si es necesario (para autovalores reales)
    # if(sum(x_new * x) < 0) lambda <- -lambda
    
    # 4. Normalización
    x_new <- x_new / sqrt(sum(x_new^2))
    
    # 5. Convergencia
    if(abs(lambda - lambda_old) < epsilon) {
      return(list(
        eigenvalue = lambda, 
        eigenvector = x_new, 
        iterations = i
      ))
    }
    
    x <- x_new
    lambda_old <- lambda
  }
  warning("El método no convergió")
  return(list(eigenvalue = lambda, eigenvector = x, iterations = max_iter))
}

Probemos el método con una matriz dispersa grande pero sencilla.

# Creamos una matriz dispersa simétrica de 500x500
set.seed(2026)
n <- 500
M <- rsparsematrix(n, n, density = 0.01, symmetric = TRUE)

# Método manual
resultado <- power_method(M)

# Comparación con función base 
# (que usa LAPACK, mucho más lento en muy alta dimensión)
real <- eigen(as.matrix(M), only.values = TRUE)$values[1]

cat("Método Potencias: ", resultado$eigenvalue, "\n")

Método Potencias:  5.307119

cat("Eigen de R base:  ", real, "\n")

Eigen de R base:   5.291721

cat("Iteraciones:      ", resultado$iterations)

Iteraciones:       795

3.1.2 Determinación del autovalor (el cociente de Rayleigh)

En el algoritmo anterior, hemos visto cómo el vector converge a la dirección del autovector $\mathbf{v}$. Pero, ¿cómo obtenemos el valor escalar $\lambda$?

La forma más precisa de estimarlo es mediante el Cociente de Rayleigh. Si $A \mathbf{v} \approx \lambda \mathbf{v}$, multiplicando por la izquierda por $\mathbf{v}^T$:

\[\mathbf{v}^T A \mathbf{v} \approx \lambda \mathbf{v}^T \mathbf{v}\]

Despejando $\lambda$:

\[\lambda = \frac{\mathbf{v}^T A \mathbf{v}}{\mathbf{v}^T \mathbf{v}}\]

Si nuestro autovector está normalizado ($||\mathbf{v}|| = 1$), la fórmula se simplifica elegantemente a: \[\lambda = \mathbf{v}^T A \mathbf{v}\]

3.1.3 Uso de la librería `matlib`

Aunque implementar el algoritmo manualmente es vital para entenderlo, R dispone de herramientas didácticas en el paquete matlib. La función powerMethod() no solo calcula el resultado, sino que muestra la evolución de las iteraciones, lo cual es fantástico para visualizar la velocidad de convergencia.

# Instalación si no está presente: install.packages("matlib")
library(matlib)

# Definimos una matriz sencilla 2x2
A <- matrix(c(4, 1, 
              2, 3), nrow = 2, byrow = TRUE)

# Ejecutamos el método de la potencia visualizando pasos
# plot = FALSE para no generar gráfico en este PDF (opcional)
resultado_matlib <- powerMethod(A, v = c(1, 1), plot = FALSE)

print(resultado_matlib)

$vector_iterations
     v1        v2
[1,]  1 0.7071068
[2,]  1 0.7071068

$iter
[1] 2

$vector
          [,1]
[1,] 0.7071068
[2,] 0.7071068

$value
[1] 5

Advertencia de Ingeniería: matlib vs Matrix

La librería matlib está diseñada con fines pedagógicos y trabaja con matrices densas estándar. Para los sistemas reales de ingeniería (matrices sparse de dimensión $10^5$), NO debe usar matlib. En su lugar, utilice implementaciones manuales optimizadas con el paquete Matrix (como la función power_method que definimos anteriormente) o librerías de alto rendimiento como RSpectra.

Ejercicio 3.1 (Autovalores de una matriz)

Crea una matriz de 3x3 con valores aleatorios.
Calcula los autovalores de la matriz usando la función eigen() de R base.
Calcula los autovalores de la matriz usando la función powerMethod() de la librería matlib.
Compara los resultados y discute las diferencias.

3.2 El Algoritmo QR (espectro completo)

Mientras que el método de las potencias encuentra solo el autovalor dominante, el Algoritmo QR es el “caballo de batalla” del álgebra lineal numérica moderna para encontrar todos los autovalores de una matriz.

3.2.1 La Lógica del Algoritmo

El método se basa en una iteración sorprendentemente simple llamada iteración QR. Comenzamos con $A_0 = A$ y en cada paso $k$:

Factorizamos la matriz actual en una ortogonal $Q_k$ y una triangular $R_k$ tal que $A_k = Q_k R_k$.
Multiplicamos en orden inverso para obtener la siguiente matriz: $A_{k+1} = R_k Q_k$.

Esta “mezcla” preserva los autovalores (las matrices son semejantes), pero empuja los valores hacia la diagonal principal. Si la matriz es simétrica, $A_k$ converge a una matriz diagonal donde los elementos $a_{ii}$ son los autovalores.

qr_simple_algo <- function(A, iter = 20) {
  A_k <- A
  # Guardamos la evolución para visualizar (opcional)
  historia <- numeric(iter)
  
  for(i in 1:iter) {
    # 1. Descomposición QR: A = Q * R
    descomp <- qr(A_k)
    Q <- qr.Q(descomp)
    R <- qr.R(descomp)
    
    # 2. Actualización: A_new = R * Q
    A_k <- R %*% Q
    
    # Guardamos el valor de la esquina superior (para ver convergencia)
    historia[i] <- A_k[1,1]
  }
  
  return(list(matriz_final = A_k, autovalores = diag(A_k)))
}

# --- PRUEBA DEL CONCEPTO ---
# Matriz simétrica 3x3
M <- matrix(c(5, 2, 0,
              2, 3, 1,
              0, 1, 1), nrow = 3)

resultado <- qr_simple_algo(M)

print("Matriz final (casi diagonal):")

[1] "Matriz final (casi diagonal):"

print(round(resultado$matriz_final, 3))

     [,1]  [,2]  [,3]
[1,] 6.29 0.000 0.000
[2,] 0.00 2.294 0.000
[3,] 0.00 0.000 0.416

cat("\nAutovalores estimados (Diagonal):", 
    round(resultado$autovalores, 4), "\n")


Autovalores estimados (Diagonal): 6.2899 2.2943 0.4158

cat("Autovalores reales (función eigen):", round(eigen(M)$values, 4))

Autovalores reales (función eigen): 6.2899 2.2943 0.4158

Nota de Ingeniería

En R, para matrices dispersas, raramente calculamos todos los autovalores (sería demasiado costoso). Usamos métodos como Arnoldi / Lanczos (disponibles en paquetes como RSpectra) que generalizan el método de las potencias para encontrar los primeros autovalores.

Ejercicio 3.2 (Algoritmo QR)

Crea una matriz de 3x3 con valores aleatorios.
Calcula los autovalores de la matriz usando la función eigen() de R base.
Calcula los autovalores de la matriz usando la función qr_simple_algo() que definimos anteriormente.
Compara los resultados y discute las diferencias.

3.2.2 Ejemplo Real: Dinámica de Poblaciones (Matrices de Leslie)

P.H. Leslie y su Matriz

Patrick Holt Leslie fue un ecólogo británico que, en 1945, introdujo este modelo matricial para describir el crecimiento de poblaciones. Aunque parece una aplicación puramente biológica, su estructura matemática es la misma que encontramos en los sistemas de recomendación y en la propagación de información en redes sociales.

Una de las aplicaciones más elegantes de los autovalores en biología es el estudio de poblaciones estructuradas por edad. La matriz de Leslie permite predecir cómo evolucionará una población a lo largo del tiempo.

Consideremos una especie con tres etapas de vida: jóvenes, pre-adultos y adultos. Su dinámica se rige por:

Fecundidad: Cuántos descendientes produce cada individuo de cada edad.
Supervivencia: Qué fracción de individuos pasa a la siguiente etapa.

Intuición Matemática

Desde un punto de vista teórico, el éxito de este análisis depende del Teorema de Perron-Frobenius. Este teorema garantiza que si una matriz tiene todas sus entradas positivas (o es irreducible y primitiva, como la mayoría de matrices de Leslie), posee un autovalor real único $\lambda$ que es estrictamente mayor que el valor absoluto de cualquier otro autovalor. Este es el “ritmo” dominante que dicta el futuro de la población.

Interpretación:

Si $\lambda > 1$, la población crece indefinidamente.
Si $\lambda < 1$, la población tiende a la extinción.
El autovector asociado al autovalor dominante nos indica la distribución estable por edades: a largo plazo, la proporción de individuos en cada etapa se mantendrá constante según las componentes de este vector.

# Matriz de Leslie L:
# Juveniles -> Pre-adultos -> Adultos
L <- matrix(c(0,   0.5, 2.0,  # Fecundidad (solo adultos producen mucho)
              0.1, 0,   0,    # 10% de juveniles sobrevive a pre-adultos
              0,   0.5, 0),   # 50% de pre-adultos sobrevive a adultos
            nrow = 3, byrow = TRUE)

# Autovalores
eigen_L <- eigen(L)
lambda_dominante <- Re(eigen_L$values[1])

cat("Tasa de crecimiento a largo plazo (lambda):", round(lambda_dominante, 3), "\n")

Tasa de crecimiento a largo plazo (lambda): 0.5

# Distribución estable (normalizada)
v_estable <- Re(eigen_L$vectors[,1])
v_estable <- v_estable / sum(v_estable)
names(v_estable) <- c("Jóvenes", "Pre-adultos", "Adultos")
print(v_estable)

    Jóvenes Pre-adultos     Adultos 
  0.7142857   0.1428571   0.1428571

3.2.3 Ejemplo IT: Centralidad de Vector Propio (Eigenvector Centrality)

En redes sociales, seguridad informática o análisis de infraestructuras, no todos los nodos son iguales. La Centralidad de Vector Propio propone que un nodo es importante si está conectado a otros nodos importantes.

Si $A$ es la matriz de adyacencia de una red, queremos asignar una puntuación $x_i$ a cada nodo tal que: \[x_i = \frac{1}{\lambda} \sum_{j \in G} a_{ij} x_j\] En forma matricial: $A \mathbf{x} = \lambda \mathbf{x}$. ¡Es exactamente un problema de autovectores!

library(igraph)
# Creamos una red pequeña "star-like"
g <- make_star(5) 

plot(g)

A <- as_adjacency_matrix(g, sparse = FALSE)
print(A)

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    1    0    0    0    0
[3,]    1    0    0    0    0
[4,]    1    0    0    0    0
[5,]    1    0    0    0    0

# Calculamos centralidad por autovector
# Para ello, necesitamos la matriz de adyacencia transpuesta:
# De la columna j se va a la fila i: del nodo j se va al nodo i.
ev <- eigen(t(A))
centralidad <- Re(ev$vectors[,1])

# El nodo 1 (centro de la estrella) tendrá la mayor importancia
print(round(abs(centralidad), 3))

[1] 1 0 0 0 0

Esta es la base conceptual de lo que Google hace con PageRank, aunque con modificaciones sutiles para asegurar convergencia en la “web salvaje”.

3.3 Descomposición en Valores Singulares (SVD)

Aunque los autovalores son fundamentales para matrices cuadradas, la Descomposición en Valores Singulares (SVD) es la herramienta universal del álgebra lineal, aplicable a cualquier matriz rectangular.

Definición 3.1 (Teorema SVD) Toda matriz real $A$ de dimensión $m \times n$ puede factorizarse como: \[A = U \Sigma V^T\] Donde $U$ ($m \times m$) y $V$ ($n \times n$) son matrices ortogonales, mientras que $\Sigma$ ($m \times n$) es una matriz diagonal que contiene los valores singulares no negativos, ordenados decrecientemente ($\sigma_1 \ge \sigma_2 \ge \dots \ge 0$).

La SVD es el motor matemático detrás de numerosas aplicaciones modernas. En la compresión de imágenes, nos permite aproximar una matriz visual guardando solo los primeros $k$ valores singulares (Teorema de Eckart-Young). En los sistemas de recomendación, es la base del filtrado colaborativo matricial. Además, el Análisis de Componentes Principales (PCA) no es más que una SVD aplicada a la matriz de datos centrada en la media.

3.3.1 Implementación en R

Para realizar la SVD en R, utilizamos la función svd(), que descompone la matriz en sus tres componentes: $U$, $D$ (diagonal con los valores singulares) y $V$.

# Creamos una matriz sencilla (no necesariamente cuadrada)
A <- matrix(c(1, 2, 3, 
              4, 5, 6), nrow = 2, byrow = TRUE)

# Aplicamos SVD
res_svd <- svd(A)

# Mostramos los componentes
# Valores singulares (diagonal de Sigma)
print(res_svd$d)

[1] 9.5080320 0.7728696

# Matrices U y V
print(res_svd$u)

           [,1]       [,2]
[1,] -0.3863177 -0.9223658
[2,] -0.9223658  0.3863177

print(res_svd$v)

           [,1]       [,2]
[1,] -0.4286671  0.8059639
[2,] -0.5663069  0.1123824
[3,] -0.7039467 -0.5811991

3.3.1.1 Reconstrucción y Aproximación

Una propiedad fundamental es que podemos reconstruir la matriz original como $A = U D V^T$. Si decidimos usar solo los $k$ valores singulares más grandes, obtenemos la mejor aproximación de rango $k$ (Teorema de Eckart-Young).

# Reconstrucción completa: U %*% diag(d) %*% t(V)
A_rec <- res_svd$u %*% diag(res_svd$d) %*% t(res_svd$v)
print(A_rec)

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

# Aproximación de Rango 1 (usando solo el primer valor singular)
A_k1 <- res_svd$u[, 1, drop = FALSE] %*% (res_svd$d[1]) %*% t(res_svd$v[, 1, drop = FALSE])
print(A_k1)

         [,1]     [,2]     [,3]
[1,] 1.574546 2.080114 2.585681
[2,] 3.759361 4.966446 6.173530

En este ejemplo 2x3, la aproximación de rango 1 captura la “estructura principal” de la matriz con un solo escalar y dos vectores. En matrices de millones de píxeles, esto permite reducciones de tamaño masivas con pérdida visual mínima.

3.3.1.2 Ejemplo Práctico: Compresión de Imagen (Matriz)

Podemos visualizar la potencia de la SVD usando una matriz de datos que represente una superficie o imagen. Utilizaremos el dataset volcano de R base (una matriz de 87x61 que representa el volcán Maunga Whau).

# 1. Matriz original
V <- volcano
svd_volcano <- svd(V)

# Función para aproximar
aproximar <- function(res_svd, k) {
  res_svd$u[, 1:k] %*% diag(res_svd$d[1:k]) %*% t(res_svd$v[, 1:k])
}

# 2. Visualización
par(mfrow = c(1, 4), mar = c(1, 1, 2, 1))
image(V, main = "Original (Rango 61)", axes = FALSE, col = terrain.colors(100))
image(aproximar(svd_volcano, 2), main = "k = 2", axes = FALSE, col = terrain.colors(100))
image(aproximar(svd_volcano, 5), main = "k = 5", axes = FALSE, col = terrain.colors(100))
image(aproximar(svd_volcano, 15), main = "k = 15", axes = FALSE, col = terrain.colors(100))

Observa cómo con solo 5 valores singulares (de los 61 posibles), ya capturamos la silueta general del volcán. Con 15, la diferencia visual es mínima. Esto reduce el almacenamiento necesario de $87 \times 61 = 5307$ valores a solo $15 \times (87 + 1 + 61) = 2235$ valores (una reducción del 58%).

Ejercicio 3.3 (Compresión por SVD)

Crea una matriz de 10x10 usando outer(1:10, 1:10, "+") y súmale un poco de ruido aleatorio con rnorm(100, sd = 0.1).
Calcula su SVD.
Aproxima la matriz usando solo los 2 primeros valores singulares.
Calcula el “error de reconstrucción” (la norma de la diferencia entre la matriz original y la aproximada).

3.3.2 Resultados Técnicos y Estabilidad

Más allá de la aproximación, la SVD ofrece una visión profunda sobre la arquitectura de la matriz:

Relación con Eigendecomposition: Los valores singulares de $A$ son las raíces cuadradas de los autovalores de $A^T A$ (o $AA^T$). Los vectores en $V$ son los autovectores de $A^T A$.
Número de Condición ($\kappa$): Se define como la razón entre el valor singular más grande y el más pequeño: $\kappa(A) = \sigma_1 / \sigma_n$. Un número de condición muy alto indica que la matriz es “casi singular”, lo que significa que pequeños errores en los datos (o redondeo numérico) pueden causar errores gigantes en la resolución de $Ax=b$.
Pseudoinversa de Moore-Penrose: Permite “invertir” matrices que no tienen inversa (rectangulares o de rango deficiente). Se calcula fácilmente con SVD como $A^+ = V \Sigma^+ U^T$, donde $\Sigma^+$ contiene $1/\sigma_i$ en la diagonal.

# Matriz casi singular
A_enferma <- matrix(c(1, 1, 
                      1, 1.000001), nrow = 2)

s <- svd(A_enferma)$d
cat("Número de condición:", s[1] / s[2])

Número de condición: 4000002

3.3.2.1 Ejemplo: Cálculo de la Pseudoinversa

Cuando una matriz es rectangular ($2 \times 3$), no tiene inversa tradicional. Sin embargo, la SVD nos permite calcular su pseudoinversa de Moore-Penrose ($A^+$), esencial para resolver problemas de mínimos cuadrados.

A <- matrix(c(1, 2, 3, 
              4, 5, 6), nrow = 2, byrow = TRUE)

res <- svd(A)
# A+ = V * (1/D) * Ut
D_inv <- diag(1/res$d)

A_plus <- res$v %*% D_inv %*% t(res$u)

# Verificación: A %*% A+ %*% A debe ser igual a A
print(round(A %*% A_plus %*% A, 3))

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

3.4 Ejercicios adicionales

Esta práctica está diseñada para realizarse en aproximadamente una hora y consolida los conceptos de este capítulo.

Ejercicio 3.4 (Dinámica de Población) Se estudia una colonia de aves con tres grupos: 0-1 año (juveniles), 1-2 años (sub-adultos) y >2 años (adultos).

Los juveniles no promedian descendientes. Los sub-adultos promedian 1 y los adultos promedian 2.
La probabilidad de que un juvenil pase a sub-adulto es 0.3.
La probabilidad de que un sub-adulto pase a adulto es 0.5.

Define la matriz de Leslie $L$ para este sistema.
Si empezamos con una población de 100 individuos en cada categoría, ¿cuántos habrá después de 10 años? (Usa el producto matricial repetido o potencias de la matriz).
Calcula el autovalor dominante. ¿Crece o decrece la población?
Encuentra la distribución estable por edades.

Ejercicio 3.5 (SVD y Ruido en Datos) La SVD se usa para “limpiar” datos. Vamos a simular una señal con mucho ruido.

Crea una matriz $A$ de 50x50 donde cada elemento $a_{ij} = i + j$. Esta es nuestra “señal pura” (tiene rango 2).
Súmale ruido: A_ruido <- A + matrix(rnorm(2500, sd = 5), 50, 50).
Calcula la SVD de A_ruido y grafica los valores singulares (plot(res$d)). ¿Cuántos valores singulares destacan sobre el “suelo” de ruido?
Reconstruye la matriz usando solo los 2 primeros valores singulares. Calcula la diferencia (error) respecto a la señal pura $A$ y compáralo con el error de la matriz ruidosa original.

Ejercicio 3.6 (Análisis de Semántica Latente) En procesamiento de lenguaje natural, la SVD se usa para encontrar “temas” en documentos. Tenemos 3 documentos y 3 palabras clave (R, Python, Datos).

El Análisis de Semántica Latente (LSA) es una técnica de procesamiento de lenguaje natural que permite descubrir la estructura “oculta” (latente) en una colección de textos. El problema del lenguaje es que es ambiguo (sinonimia y polisemia): diferentes palabras pueden significar lo mismo, y una palabra puede tener varios significados según el contexto.

La SVD resuelve esto tratando el lenguaje como un problema geométrico:

Reducción de Ruido: Al quedarnos solo con los valores singulares más grandes, eliminamos las variaciones “accidentales” del lenguaje (ruido) y nos quedamos con los patrones de co-ocurrencia más fuertes.
Espacio Semántico: Los documentos se proyectan en un espacio de baja dimensión donde la cercanía no depende de compartir las mismas palabras exactas, sino de compartir el mismo “contexto” o “tema”.
Descubrimiento de Conceptos: Las matrices $U$ y $V$ de la SVD nos dan, respectivamente, la relación de las palabras con los temas y de los documentos con los temas.

Para organizar la información, usaremos la denominada matriz de términos-documentos. Se define como una matriz donde las filas representan los términos y las columnas representan los documentos. En nuestro caso, los términos son “R”, “Python” y “Datos”, y los documentos son Doc1, Doc2 y Doc3.

Crea una matriz de términos-documentos (3x3) donde:
- Doc1: 5 veces “R”, 0 “Python”, 2 “Datos”.
- Doc2: 1 vez “R”, 4 “Python”, 1 “Datos”.
- Doc3: 0 vez “R”, 1 “Python”, 5 “Datos”.
Calcula la SVD de esta matriz.
Analiza la matriz $V$ (espacio de documentos). ¿Qué documentos están más “cerca” entre sí según sus componentes en los dos primeros valores singulares?
Calcula la matriz aproximada de rango 2. ¿Qué palabra ha ganado importancia en el Doc1 que antes no tenía (debido a la relación latente entre conceptos)?