An Improved Incremental Singular Value Decomposition and New Error Bounds

Resumen Técnico: Una Descomposición en Valores Singulares Incremental Mejorada y Nuevas Cotas de Error

Enunciado del Problema
La Descomposición en Valores Singulares (SVD) Incremental, propuesta originalmente por Brand [2, 3] y extendida a productos internos ponderados por Fareed et al. [10], es un método estándar para la aproximación de rango bajo en escenarios de datos en flujo donde no es posible almacenar la matriz completa. El algoritmo actualiza una SVD truncada ($U \approx Q\Sigma R^\top$) a medida que llegan nuevas columnas. Sin embargo, en aritmética de punto flotante, cada actualización multiplica la base singular en ejecución $Q$ por un factor ortogonal pequeño. Si bien las multiplicaciones individuales preservan la ortogonalidad hasta el error de redondeo, estos errores se acumulan a lo largo de la longitud del flujo $n$. Eventualmente, $Q$ pierde la ortogonalidad numérica, lo que hace que los valores singulares sean inexactos y los criterios de truncamiento poco fiables.

El remedio estándar es la reortogonalización periódica (por ejemplo, mediante Gram–Schmidt Modificado) activada cuando la ortogonalidad se degrada. Sin embargo, la frecuencia necesaria para mantener un nivel de precisión prescrito ha permanecido como una pregunta abierta. Como señaló Brand [3, p. 28], no estaba claro con qué frecuencia es necesaria la reortogonalización para garantizar un cierto nivel general de precisión numérica, particularmente en entornos ponderados (por ejemplo, Descomposición en Modos Propios de Elementos Finitos) donde la reortogonalización es computacionalmente costosa.

Metodología
Los autores proponen una reestructuración del algoritmo de SVD incremental que pospone las actualizaciones de los factores singulares en ejecución. La metodología central se basa en dos observaciones:

1. Dominancia de No Incremento de Rango: En aplicaciones de rango bajo, la mayoría de las columnas entrantes se encuentran dentro del subespacio actual (residuo < $\text{tol}$). Estas columnas de "no incremento de rango" no expanden el subespacio.
2. Agrupación (Batching): El efecto acumulativo de $s$ columnas sucesivas de no incremento de rango puede reproducirse mediante una única SVD delgada de una matriz buffer corta y ancha que contiene los coeficientes de proyección, en lugar de realizar $s$ SVD de matrices bordeadas por separado.

El algoritmo propuesto (Algoritmo 5) mantiene un buffer para los coeficientes de proyección de las columnas de no incremento de rango. El factor en ejecución $Q$ permanece sin cambios durante estos pasos. Cuando llega una columna de "incremento de rango" (residuo $\ge \text{tol}$), el algoritmo:

• Vacía el buffer calculando una única SVD delgada de los coeficientes acumulados y actualizando $Q$ y $R$ en un lote.
• Procesa la nueva columna, aplicando potencialmente una segunda prueba de truncamiento sobre el valor singular más pequeño de la matriz bordeada local para descartar el rango despreciable.
• Aplica un Gram–Schmidt Clásico con dos pasadas (CGS-2) solo en estos eventos de incremento de rango para restaurar la $W$-ortogonalidad de la columna añadida.

Esta reestructuración reduce el número de grandes multiplicaciones ortogonales aplicadas a $Q$ de $O(n)$ (la longitud del flujo) a $O(r + J_{E2})$, donde $r$ es el rango numérico y $J_{E2}$ es el número de eventos de incremento de rango que provocan truncamiento.

Contribuciones Clave
El artículo proporciona dos avances teóricos primarios y una mejora algorítmica práctica:

1. Cota de Ortogonalidad Independiente de la Longitud del Flujo:
   Los autores demuestran (Teorema 4.1) que bajo un modelo probabilístico estándar de error de redondeo, la pérdida de $W$-ortogonalidad del factor izquierdo calculado $Q$ está acotada por:
   $$\|I - Q^\top W Q\|_2 \leq c_1 (r + J_{E2}) r \sqrt{m} \log(2/\delta) \varepsilon_{\text{mach}} + O(\varepsilon_{\text{mach}}^2)$$
   Crucialmente, esta cota es independiente de la longitud del flujo $n$. Depende de $m$ (dimensión de la columna) solo a través de un factor $\sqrt{m}$. Esto resuelve la pregunta abierta de Brand: la reortogonalización no es necesaria en cada paso; la estrategia de agrupación propuesta garantiza que la ortogonalidad se mantenga a precisión de máquina uniformemente en $n$.

2. Cota de Error Directo Aguda:
   El artículo establece una nueva cota de error directo para la factorización en la norma del operador ponderada por $W$ (Teorema 4.5):
   $$\|U - Q\Sigma R^\top\|_W \leq \sqrt{n} \, \text{tol} + (\text{términos de redondeo uniformes en } n)$$
   Esto afina la cota existente de Fareed y Singler [9], que era $n \, \text{tol}$. Los autores demuestran mediante un ejemplo constructivo que la tasa $\sqrt{n}$ es alcanzable, mostrando que la cota anterior era laxa por un factor de $\sqrt{n}$.

3. Eficiencia Algorítmica:
   Al reducir el número de costosos pasos de reortogonalización y agrupar las actualizaciones, el algoritmo propuesto reduce significativamente la sobrecarga computacional.

Resultados
Se realizaron experimentos numéricos en un problema de prueba parabólico 2D (tanto en configuraciones de elementos finitos sin ponderar como ponderadas) y en datos sintéticos.

• Ortogonalidad: Mientras que los algoritmos existentes (actualización directa de Brand y actualización reortogonalizada de Fareed) mostraron una degradación de la ortogonalidad que crecía con $n$ (desviándose hasta $10^{-3}$ en casos ponderados), el algoritmo propuesto mantuvo la ortogonalidad a precisión de máquina ($\sim 10^{-14}$) a lo largo de longitudes de flujo de hasta $n=10,000$ y dimensiones de columna de hasta $m=100,000$.
• Cotas de Error: Los experimentos confirmaron la tasa de truncamiento $\sqrt{n} \, \text{tol}$. En un ejemplo constructivo, el error medido coincidió con la predicción $\sqrt{n} \, \text{tol}$, mientras que la cota de Fareed–Singler ($n \, \text{tol}$) fue laxa por un factor que crece como $\sqrt{n}$.
• Rendimiento: El algoritmo propuesto se ejecutó 4.5$\times$ a 34$\times$ más rápido que sus competidores más cercanos. En el caso ponderado, el algoritmo existente pasó el 94% de su tiempo en reortogonalización, mientras que el algoritmo propuesto redujo este costo en tres órdenes de magnitud.

Significado
El artículo afirma resolver la pregunta abierta sobre la frecuencia de la reortogonalización en la SVD incremental. Al reestructurar el algoritmo para agrupar las actualizaciones de no incremento de rango, logra:

1. Ortogonalidad Uniforme: La base en ejecución permanece numéricamente ortogonal independientemente de la longitud del flujo, eliminando la necesidad de comprobaciones heurísticas y periódicas de reortogonalización que carecen de garantías a priori.
2. Análisis de Error Más Estricto: Proporciona una cota de error matemáticamente rigurosa y aguda ($\sqrt{n} \, \text{tol}$) que mejora al estado del arte en un factor de $\sqrt{n}$.
3. Aceleración Práctica: El método ofrece aceleraciones sustanciales en aplicaciones dominadas por los costos de reortogonalización, como la POD de elementos finitos, sin sacrificar la precisión.

Los autores señalan que el algoritmo ya se ha aplicado en versiones preliminares a EDPs de orden fraccionario en el tiempo, ecuaciones integro-diferenciales y optimización con restricciones de EDP, confiando en las garantías de ortogonalidad y error ahora establecidas. Se sugiere trabajo futuro para afinar el análisis de la deriva de ortogonalidad (actualmente conservador por $\sim 10^3$) y optimizar el paso de proyección, que actualmente domina el tiempo de ejecución.