public static TetradMatrix incompleteCholeskyGramMatrix(List<Kernel> kernels,DataSet dataset,List<Node> nodes,double precision){
  if (precision <= 0) {
    throw new IllegalArgumentException("Precision must be > 0");
  }
  int m=dataset.getNumRows();
  TetradMatrix G=new TetradMatrix(m,m);
  double[] Dadv=new double[m];
  int[] p=new int[m];
  for (int i=0; i < m; i++) {
    Dadv[i]=evaluate(kernels,dataset,nodes,i,i);
    p[i]=i;
  }
  int cols=m;
  for (int k=0; k < m; k++) {
    double best=Dadv[k];
    int bestInd=k;
    for (int j=(k + 1); j < m; j++) {
      if (Dadv[j] > best / .99) {
        best=Dadv[j];
        bestInd=j;
      }
    }
    if (best < precision) {
      cols=k - 1;
      break;
    }
    int pk=p[k];
    p[k]=p[bestInd];
    p[bestInd]=pk;
    double dk=Dadv[k];
    Dadv[k]=Dadv[bestInd];
    Dadv[bestInd]=dk;
    for (int j=0; j < k; j++) {
      double gk=G.get(k,j);
      G.set(k,j,G.get(bestInd,j));
      G.set(bestInd,j,gk);
    }
    double diag=Math.sqrt(Dadv[k]);
    G.set(k,k,diag);
    for (int j=(k + 1); j < m; j++) {
      double s=0.0;
      for (int i=0; i < k; i++) {
        s+=G.get(j,i) * G.get(k,i);
      }
      G.set(j,k,(evaluate(kernels,dataset,nodes,p[j],p[k]) - s) / diag);
    }
    for (int j=(k + 1); j < m; j++) {
      Dadv[j]-=Math.pow(G.get(j,k),2);
    }
    Dadv[k]=0;
  }
  TetradMatrix Gm=new TetradMatrix(m,cols);
  for (int i=0; i < m; i++) {
    for (int j=0; j < cols; j++) {
      Gm.set(i,j,G.get(i,j));
    }
  }
  return Gm;
}

public static TetradMatrix incompleteCholeskyGramMatrix ( List < Kernel > kernels , DataSet dataset , List < Node > nodes , double precision ) { if ( precision <= 0 ) { throw new IllegalArgumentException ( "Precision must be > 0" ) ; } int m = dataset . getNumRows ( ) ; TetradMatrix G = new TetradMatrix ( m , m ) ; double [ ] Dadv = new double [ m ] ; int [ ] p = new int [ m ] ; for ( int i = 0 ; i < m ; i ++ ) { Dadv [ i ] = evaluate ( kernels , dataset , nodes , i , i ) ; p [ i ] = i ; } int cols = m ; for ( int k = 0 ; k < m ; k ++ ) { double best = Dadv [ k ] ; int bestInd = k ; for ( int j = ( k + 1 ) ; j < m ; j ++ ) { if ( Dadv [ j ] > best / .99 ) { best = Dadv [ j ] ; bestInd = j ; } } if ( best < precision ) { cols = k - 1 ; break ; } int pk = p [ k ] ; p [ k ] = p [ bestInd ] ; p [ bestInd ] = pk ; double dk = Dadv [ k ] ; Dadv [ k ] = Dadv [ bestInd ] ; Dadv [ bestInd ] = dk ; for ( int j = 0 ; j < k ; j ++ ) { double gk = G . get ( k , j ) ; G . set ( k , j , G . get ( bestInd , j ) ) ; G . set ( bestInd , j , gk ) ; } double diag = Math . sqrt ( Dadv [ k ] ) ; G . set ( k , k , diag ) ; for ( int j = ( k + 1 ) ; j < m ; j ++ ) { double s = 0.0 ; for ( int i = 0 ; i < k ; i ++ ) { s += G . get ( j , i ) * G . get ( k , i ) ; } G . set ( j , k , ( evaluate ( kernels , dataset , nodes , p [ j ] , p [ k ] ) - s ) / diag ) ; } for ( int j = ( k + 1 ) ; j < m ; j ++ ) { Dadv [ j ] -= Math . pow ( G . get ( j , k ) , 2 ) ; } Dadv [ k ] = 0 ; } TetradMatrix Gm = new TetradMatrix ( m , cols ) ; for ( int i = 0 ; i < m ; i ++ ) { for ( int j = 0 ; j < cols ; j ++ ) { Gm . set ( i , j , G . get ( i , j ) ) ; } } return Gm ; }
approximates gram matrix using incomplete cholesky factorization
