public double empiricalHSICincompleteCholesky(TetradMatrix Gy,TetradMatrix Gx,TetradMatrix Gz,int m){
  int ky=Gy.columns();
  int kx=Gx.columns();
  int kz=Gz.columns();
  TetradMatrix H=KernelUtils.constructH(m);
  TetradMatrix Gcy=H.times(Gy);
  TetradMatrix Gcx=H.times(Gx);
  TetradMatrix Gcz=H.times(Gz);
  TetradMatrix A=new TetradMatrix(ky,kx);
  TetradMatrix Gcyt=Gcy.transpose();
  A=Gcyt.times(Gcx);
  TetradMatrix B=Gcy.times(A);
  TetradMatrix Kyx=new TetradMatrix(m,m);
  TetradMatrix Gcxt=new TetradMatrix(kx,m);
  Gcxt=Gcx.transpose();
  Kyx=B.times(Gcxt);
  double empHSIC=0.0;
  double xy=0.0;
  for (int i=0; i < m; i++) {
    empHSIC+=matrixProductEntry(B,Gcxt,i,i);
  }
  TetradMatrix Gytz=Gcyt.times(Gcz);
  TetradMatrix Gczt=Gcz.transpose();
  TetradMatrix Gztx=Gczt.times(Gcx);
  TetradMatrix Gztz=Gczt.times(Gcz);
  TetradMatrix Gztzr=Gztz.copy();
  for (int i=0; i < kz; i++) {
    Gztzr.set(i,i,Gztz.get(i,i) + this.regularizer);
  }
  TetradMatrix ZI=Gztzr.inverse();
  TetradMatrix ZIzt=ZI.times(Gczt);
  TetradMatrix Gzr=Gcz.copy();
  for (int i=0; i < m; i++) {
    for (int j=0; j < kz; j++) {
      Gzr.set(i,j,Gcz.get(i,j) * (-1.0 / this.regularizer));
    }
  }
  TetradMatrix Zinv=Gzr.times(ZIzt);
  for (int i=0; i < m; i++) {
    Zinv.set(i,i,Zinv.get(i,i) + (1.0 / this.regularizer));
  }
  TetradMatrix Gztzinv=Gczt.times(Zinv);
  TetradMatrix Gzinvz=Zinv.times(Gcz);
  TetradMatrix Gztinv2z=Gztzinv.times(Gzinvz);
  TetradMatrix Gytzztzinv2z=Gytz.times(Gztinv2z);
  TetradMatrix Gytzztzinv2zztx=Gytzztzinv2z.times(Gztx);
  TetradMatrix Gyytzztzinv2zztx=Gcy.times(Gytzztzinv2zztx);
  double second=0.0;
  for (int i=0; i < m; i++) {
    second+=matrixProductEntry(Gyytzztzinv2zztx,Gcxt,i,i);
  }
  empHSIC-=2 * second;
  TetradMatrix Gxtz=Gcxt.times(Gcz);
  TetradMatrix Gxtzztinv2z=Gxtz.times(Gztinv2z);
  TetradMatrix Gyytzztzinv2zztxxtzztinv2z=Gyytzztzinv2zztx.times(Gxtzztinv2z);
  for (int i=0; i < m; i++) {
    empHSIC+=matrixProductEntry(Gyytzztzinv2zztxxtzztinv2z,Gczt,i,i);
  }
  double betaz=0.0;
  for (int i=0; i < (m - 1); i++) {
    for (int j=(i + 1); j < m; j++) {
      betaz+=Math.pow(matrixProductEntry(Gcz,Gczt,i,j),2);
      betaz+=Math.pow(matrixProductEntry(Gcz,Gczt,j,i),2);
    }
  }
  empHSIC*=(m / (betaz * (m - 1)));
  return empHSIC;
}

public double empiricalHSICincompleteCholesky ( TetradMatrix Gy , TetradMatrix Gx , TetradMatrix Gz , int m ) { int ky = Gy . columns ( ) ; int kx = Gx . columns ( ) ; int kz = Gz . columns ( ) ; TetradMatrix H = KernelUtils . constructH ( m ) ; TetradMatrix Gcy = H . times ( Gy ) ; TetradMatrix Gcx = H . times ( Gx ) ; TetradMatrix Gcz = H . times ( Gz ) ; TetradMatrix A = new TetradMatrix ( ky , kx ) ; TetradMatrix Gcyt = Gcy . transpose ( ) ; A = Gcyt . times ( Gcx ) ; TetradMatrix B = Gcy . times ( A ) ; TetradMatrix Kyx = new TetradMatrix ( m , m ) ; TetradMatrix Gcxt = new TetradMatrix ( kx , m ) ; Gcxt = Gcx . transpose ( ) ; Kyx = B . times ( Gcxt ) ; double empHSIC = 0.0 ; double xy = 0.0 ; for ( int i = 0 ; i < m ; i ++ ) { empHSIC += matrixProductEntry ( B , Gcxt , i , i ) ; } TetradMatrix Gytz = Gcyt . times ( Gcz ) ; TetradMatrix Gczt = Gcz . transpose ( ) ; TetradMatrix Gztx = Gczt . times ( Gcx ) ; TetradMatrix Gztz = Gczt . times ( Gcz ) ; TetradMatrix Gztzr = Gztz . copy ( ) ; for ( int i = 0 ; i < kz ; i ++ ) { Gztzr . set ( i , i , Gztz . get ( i , i ) + this . regularizer ) ; } TetradMatrix ZI = Gztzr . inverse ( ) ; TetradMatrix ZIzt = ZI . times ( Gczt ) ; TetradMatrix Gzr = Gcz . copy ( ) ; for ( int i = 0 ; i < m ; i ++ ) { for ( int j = 0 ; j < kz ; j ++ ) { Gzr . set ( i , j , Gcz . get ( i , j ) * ( - 1.0 / this . regularizer ) ) ; } } TetradMatrix Zinv = Gzr . times ( ZIzt ) ; for ( int i = 0 ; i < m ; i ++ ) { Zinv . set ( i , i , Zinv . get ( i , i ) + ( 1.0 / this . regularizer ) ) ; } TetradMatrix Gztzinv = Gczt . times ( Zinv ) ; TetradMatrix Gzinvz = Zinv . times ( Gcz ) ; TetradMatrix Gztinv2z = Gztzinv . times ( Gzinvz ) ; TetradMatrix Gytzztzinv2z = Gytz . times ( Gztinv2z ) ; TetradMatrix Gytzztzinv2zztx = Gytzztzinv2z . times ( Gztx ) ; TetradMatrix Gyytzztzinv2zztx = Gcy . times ( Gytzztzinv2zztx ) ; double second = 0.0 ; for ( int i = 0 ; i < m ; i ++ ) { second += matrixProductEntry ( Gyytzztzinv2zztx , Gcxt , i , i ) ; } empHSIC -= 2 * second ; TetradMatrix Gxtz = Gcxt . times ( Gcz ) ; TetradMatrix Gxtzztinv2z = Gxtz . times ( Gztinv2z ) ; TetradMatrix Gyytzztzinv2zztxxtzztinv2z = Gyytzztzinv2zztx . times ( Gxtzztinv2z ) ; for ( int i = 0 ; i < m ; i ++ ) { empHSIC += matrixProductEntry ( Gyytzztzinv2zztxxtzztinv2z , Gczt , i , i ) ; } double betaz = 0.0 ; for ( int i = 0 ; i < ( m - 1 ) ; i ++ ) { for ( int j = ( i + 1 ) ; j < m ; j ++ ) { betaz += Math . pow ( matrixProductEntry ( Gcz , Gczt , i , j ) , 2 ) ; betaz += Math . pow ( matrixProductEntry ( Gcz , Gczt , j , i ) , 2 ) ; } } empHSIC *= ( m / ( betaz * ( m - 1 ) ) ) ; return empHSIC ; }
empirical unconditional hilbert - schmidt dependence measure for x and y given z using incomplete cholesky decomposition to approximate gram matrices
