question: 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 ; }
target:  Empirical unconditional Hilbert-Schmidt Dependence Measure for X and Y given Z using incomplete Cholesky decomposition to approximate Gram matrices