Décomposition de Cholesky générant des NaN en Java

Question

Je ne sais pas si c'est une maths.se ou une question SO, mais je vais avec SO car je pense que c'est lié à Java.

Je suis un livre de texte sur les processus gaussiens (R&W) et l'implémentation de quelques exemples en Java. Une étape commune pour plusieurs exemples est de générer une décomposition de Cholesky d'une matrice de covariance. Dans ma tentative je peux obtenir des résultats réussis pour des matrices jusqu'à une taille limitée (33x33). Cependant, pour tout ce qui est plus grand, un NaN apparaît dans la diagonale (en 32,32) et ainsi toutes les valeurs suivantes dans la matrice sont également NaNs.

Le code est indiqué ci-dessous et la source du NaN est indiquée dans la méthode cholesky. Essentiellement l'élément de covariance a[32][32] est 1.0, mais la valeur de sum est un peu au dessus (1.0000001423291431), donc la racine carrée est imaginaire. Mes questions sont les suivantes:

1. Est-ce un résultat attendu de l'algèbre linéaire, ou, par exemple, un artefact de ma mise en œuvre?
2. Comment éviter ce problème dans la pratique?

Notez que je ne recherche pas les recommandations de bibliothèques à utiliser. C'est simplement pour ma propre compréhension.

Toutes mes excuses pour la longueur, mais je l'ai essayé de fournir un MWE complet:

import static org.junit.Assert.assertFalse; 

import org.junit.Test; 

public class CholeskyTest { 

    @Test 
    public void testCovCholesky() { 
     final int n = 34; // Test passes for n<34 
     final double[] xData = getSpread(-5, 5, n); 
     double[][] cov = covarianceSE(xData); 
     double[][] lower = cholesky(cov); 
     for(int i=0; i<n; ++i) { 
      for(int j=0; j<n; ++j) { 
       assertFalse("NaN at " + i + "," + j, Double.isNaN(lower[i][j])); 
      } 
     } 
    } 

    /** 
    * Generate n evenly space values from min to max inclusive 
    */ 
    private static double[] getSpread(final double min, final double max, final int n) { 
     final double[] values = new double[n]; 
     final double delta = (max - min)/(n - 1); 
     for(int i=0; i<n; ++i) { 
      values[i] = min + i*delta; 
     } 
     return values; 
    } 

    /** 
    * Calculate the covariance matrix for the given observations using 
    * the squared exponential (SE) covariance function. 
    */ 
    private static double[][] covarianceSE (double[] v) { 
     final int m = v.length; 
     double[][] k = new double[m][]; 
     for(int i=0; i<m; ++i) { 
      double vi = v[i]; 
      double row[] = new double[m]; 
      for(int j=0; j<m; ++j) { 
       double dist = vi - v[j]; 
       row[j] = Math.exp(-0.5*dist*dist); 
      } 
      k[i] = row; 
     } 
     return k; 
    } 

    /** 
    * Calculate lower triangular matrix L such that LL^T = A 
    * Using Cholesky decomposition from 
    * https://rosettacode.org/wiki/Cholesky_decomposition#Java 
    */ 
    private static double[][] cholesky(double[][] a) { 
     final int m = a.length; 
     double[][] l = new double[m][m]; 
     for(int i = 0; i< m;i++){ 
      for(int k = 0; k < (i+1); k++){ 
       double sum = 0; 
       for(int j = 0; j < k; j++){ 
        sum += l[i][j] * l[k][j]; 
       } 
       l[i][k] = (i == k) ? Math.sqrt(a[i][i] - sum) : // Source of NaN at 32,32 
        (1.0/l[k][k] * (a[i][k] - sum)); 
      } 
     } 
     return l; 
    } 
} 

Answer 1

Hmm, je pense que je l'ai trouvé une réponse à ma propre question, du même manuel que je suivais. De R&W p.201:

Dans la pratique, il peut être nécessaire d'ajouter un petit multiple de la matrice $ \ epsilon identité I $ à la matrice de covariance pour numériques raisons. En effet, les valeurs propres de la matrice K peuvent décroître très rapidement [...] et sans cette stabilisation la décomposition de Cholesky échoue. L'effet sur les échantillons générés est d'ajouter bruit de la variance supplémentaire indépendant epsilon $.

Ainsi, le changement suivant semble être suffisante:

private static double[][] cholesky(double[][] a) { 
    final int m = a.length; 
    double epsilon = 0.000001; // Small extra noise value 
    double[][] l = new double[m][m]; 
    for(int i = 0; i< m;i++){ 
     for(int k = 0; k < (i+1); k++){ 
      double sum = 0; 
      for(int j = 0; j < k; j++){ 
       sum += l[i][j] * l[k][j]; 
      } 
      l[i][k] = (i == k) ? Math.sqrt(a[i][i]+epsilon - sum) : // Add noise to diagonal values 
       (1.0/l[k][k] * (a[i][k] - sum)); 
     } 
    } 
    return l; 
} 

Answer 2

Je viens de finir d'écrire ma propre version d'une routine de Cholesky décomposition en C++ et JavaScript. Au lieu de calculer L, il calcule U, mais je serais curieux de le tester avec la matrice qui provoque l'erreur NaN. Seriez-vous en mesure de poster la matrice ici, ou contactez-moi (info dans le profil.)