Descomposición Cholesky en matriz semidefinida positiva

Question

Al intentar aplicar el algoritmo descrito en el artículo:

Algoritmo de metrópolis adaptativo robusto con tasa de aceptación coercitiva (2011), Matti Vihola

He utilizado la descomposición de Cholesky para encontrar $S_n$ :

$S_n S_n^T = S_{n-1} D_n S_{n-1}^T$

El artículo asegura que la matriz (conocida) del lado derecho será siempre positiva definida. Sin embargo, después de muchas iteraciones, parece que la imprecisión de la coma flotante hace que sólo sea semidefinida positiva. Me gustaría saber si es posible (¿y cómo?) obtener $S_n$ utilizando un descomposición Cholesky robusta de una matriz con pivote que funciona cuando el Cholesky estándar falla.

El problema es que al usarlo, el resultado obviamente no es el $S_n$ matriz, sino una $X_n$ de la forma:

$P^T X_n Y X_n^T P = S_{n-1} D_n S_{n-1}^T$

Por lo que ha dicho @J.M. ici parece posible:

La descomposición de Cholesky puede fallar en coma flotante cuando se da una matriz simétrica semidefinida positiva. Sin embargo, se puede modificar Cholesky para hacer un pivote simétrico de modo que la matriz se factorice para "siempre que la matriz parezca definida positiva. Tendrá que modificar tu fórmula de Kalman si adoptas esto, sin embargo.

Por lo tanto, parece que tengo que hacer un poco de álgebra aquí para obtener $S_n$ No sé cómo hacerlo, ya que soy bastante novato en la materia.

Answer 1

Si $\mathbf{A}$ es semipositiva, tenemos que la secuencia $\{\mathbf{A}_k\} = \{\mathbf{A}+\frac{1}{k}\mathbf{I}\}$ consiste en matrices positivas 1 . Es más, $\mathbf{A}_k \to \mathbf{A}$ en la norma del operador, y con $\mathbf{A}_k = \mathbf{L}_k\mathbf{L}_k^T$ tenemos que $\mathbf{L}_k$ converge a un límite que denominamos $\mathbf{L}$ . Finalmente, $\mathbf{L}$ tiene las propiedades deseadas, es decir $\mathbf{A} = \mathbf{L}\mathbf{L}^T$ .

Por lo tanto, un algoritmo eficiente sería añadir una pequeña diagonal $\varepsilon \mathbf{I}$ à $\mathbf{A}$ para que $\mathbf{A} + \varepsilon \mathbf{I}$ es numéricamente definida positiva, y luego realizar la descomposición de Cholesky. Por último, observamos que el cálculo de los valores propios de $\mathbf{A}$ pide más operaciones en coma flotante.

---

1 esquema de prueba en wikipedia

Answer 2

Me encontré con un problema similar mientras implementaba un algoritmo del siguiente documento.

"Online and Batch Learning of Pseudo-Metrics" (Aprendizaje en línea y por lotes de pseudométricas) Shai Shalev-Shwartz, Yoram Singer y Andrew Y. Ng, ICML 2004.

En este caso también el algoritmo garantiza que la matriz que se calcula será positiva definida. Para mí, la cuestión parece ser computacional. Así que decidí encontrar la matriz más cercana que me permitiera continuar el cálculo. En el siguiente artículo se describe cómo se puede hacer esto.

"Computing a nearest symmetric positive semidefinite matrix", Nicholas J. Higham, Linear Algebra and its Applications, volumen 103, mayo de 1988, Páginas 103-118

1. Parte simétrica de $A$ dado por $B=(A+A^T)/2$
2. Eigendecomposición de $B$ da los vectores propios $Z$ y los valores propios $\lambda_i$ .
3. Aproximación positiva de $A$ se da como $Z~\text{diag}(\text{max}(0,\lambda_i))Z^T$

Código Matlab equivalente:

 [Z,L] = eig(0.5*(A+A'));
 A_out = Z*max(L,0)*Z';

El código disponible en el siguiente sitio web que resuelve un problema similar también utiliza la misma técnica.

OASIS: Aprendizaje en línea a gran escala de la similitud de imágenes a través de la clasificación