1. Fundamentos Matemáticos

Una matriz cuadrada A de orden n es diagonalizable si existe una matriz invertible P y una matriz diagonal D tal que:

A = P·D·P⁻¹

Donde D contiene los valores propios (eigenvalues) λ₁,…,λₙ en su diagonal y las columnas de P son los vectores propios (eigenvectors) correspondientes.

• La diagonalización es posible si y solo si A tiene n vectores propios linealmente independientes.
• Para matrices simétricas reales (A = Aᵀ) siempre se cumple la condición (teorema espectral).
• En caso contrario, se recurre a la forma de Jordan o a descomposiciones como Schur.

2. Algoritmo Paso a Paso

1. Calcular los valores propios de A resolviendo det(A-λI)=0.
2. Obtener los vectores propios resolviendo (A-λI)v=0 para cada λ.
3. Verificar que los vectores propios forman una base linealmente independiente (determinante de P ≠ 0).
4. Construir P con los vectores propios como columnas y D con los λ en la diagonal.
5. Si P es invertible, A está diagonalizada; de lo contrario, la matriz no es diagonalizable.

3. Implementación en Python

Utilizaremos numpy.linalg.eig y scipy.linalg.eig para obtener valores y vectores propios.

import numpy as np
from scipy import linalg
# Matriz de ejemplo (simétrica)
A = np.array([[4, 1, 2],
              [1, 3, 0],
              [2, 0, 5]], dtype=float)
# 1️⃣ Numpy: valores y vectores propios
vals, vecs = np.linalg.eig(A)
print('Valores propios:', vals)
print('Vectores propios (columnas):\n', vecs)
# 2️⃣ Verificación de diagonalización
P = vecs
D = np.diag(vals)
A_rec = P @ D @ np.linalg.inv(P)
print('Reconstrucción (error máximo):', np.max(np.abs(A - A_rec)))
# 3️⃣ Scipy (más robusto con matrices complejas)
vals_s, vecs_s = linalg.eig(A)
print('Scipy - Valores propios:', vals_s)

En la práctica, numpy.linalg.eig es suficiente para matrices pequeñas y bien condicionadas. Para tamaños mayores (< 10⁴) o matrices con números complejos, scipy.linalg.eig ofrece mayor estabilidad y opciones de balanceo.

4. Comparativa con Otras Descomposiciones

5. Rendimiento y Escalabilidad

Para matrices densas de gran tamaño, el costo O(n³) se vuelve prohibitivo. Algunas estrategias:

• Uso de bibliotecas optimizadas: Intel MKL, OpenBLAS y cuSOLVER (GPU).
• Reducción a forma tridiagonal: Algoritmos de Lanczos para matrices simétricas.
• Descomposición parcial: Calcular sólo los k mayores valores propios con scipy.sparse.linalg.eigsh.

Ejemplo de cálculo parcial:

from scipy.sparse import random as sparse_random
from scipy.sparse.linalg import eigsh
# Matriz dispersa 10000x10000 con densidad 0.001
A_sparse = sparse_random(10000, 10000, density=0.001, format='csr')
# Obtener los 5 valores propios más grandes
vals_k, vecs_k = eigsh(A_sparse, k=5, which='LM')
print('Top 5 valores propios:', vals_k)

6. Buenas Prácticas y Solución de Problemas

Numerical Stability

• Escala la matriz antes de la diagonalización (normalización).
• Utiliza scipy.linalg.eig con balanced=True para reducir el error de redondeo.

Cuando la matriz no es diagonalizable

• Comprueba el rango de P. Si np.linalg.matrix_rank(P) < n, la matriz es defectiva.
• Recurre a la descomposición de Jordan (no disponible directamente en NumPy, pero sí en SymPy).

Seguridad en entornos productivos

• Valida la entrada: la matriz debe ser cuadrada y numérica.
• Limita el tamaño máximo aceptado para evitar ataques de denegación de servicio (DoS) por cálculos O(n³).

7. Casos de Uso Reales

• Vibraciones mecánicas: Cálculo de modos propios de estructuras.
• Google PageRank: Diagonalización del operador de transición para encontrar el vector propio dominante.
• Compresión de imágenes (PCA): Diagonalización de la matriz de covarianza para obtener componentes principales.