The \(d\)-variate \(t\) probability density function with \(\nu>0\) degrees of freedom (see and for an extensive review of this model) is given by \[ f_d({\bf x}\vert {\mu},{\Sigma},\nu) = \dfrac{\Gamma\left(\dfrac{\nu+d}{2}\right)}{\Gamma\left(\dfrac{\nu}{2}\right)\sqrt{(\pi\nu)^d\vert{\Sigma}\vert}} \left(1+\dfrac{({\bf x}-{\mu})^{T}{\Sigma}^{-1}({\bf x}-{\mu})}{\nu}\right)^{-\frac{\nu+d}{2}}, \]
The Kullback Leibler divergence between two \(d\)-variate \(t\) densities, \(f_{d,\nu}\) and \(f_{d,\nu^{\prime}}\), satisfies \[ D_{KL}(f_{d}(\cdot\vert {\mu},{\Sigma},\nu)\vert\vert f_{d,}(\cdot\vert {\mu},{\Sigma},\nu^{\prime})) = D_{KL}(f_{d}(\cdot\vert {\bf 0},{\bf I},\nu)\vert\vert f_{d}(\cdot\vert {\bf 0},{\bf I},\nu^{\prime})) =\int_{{\mathbb R}^n} f_{d}({\bf x}\vert {\bf 0},{\bf I},\nu) \log \dfrac{f_{d}({\bf x}\vert {\bf 0},{\bf I},\nu)}{f_{d}({\bf x}\vert {\bf 0},{\bf I},\nu^{\prime})} \,d{\bf x} \\ \] [1] found that this Kullback Leibler divergence can be written as \(D_{KL}(f_{d}(\cdot\vert {\mu},{\Sigma},\nu)\vert\vert f_{d,}(\cdot\vert {\mu},{\Sigma},\nu^{\prime})) =\) \[ \log \dfrac{K(d,\nu)}{K(d,\nu^{\prime})} -\dfrac{\nu+d}{2}\left[ \Psi\left(\dfrac{\nu+d}{2}\right) - \Psi\left(\dfrac{\nu}{2}\right) \right] + \dfrac{\nu^{\prime}+d}{2}K(d,\nu) \dfrac{\pi^{\frac{d}{2}}}{\Gamma\left(\frac{d}{2}\right)} \int_0^{\infty} \left(1+\dfrac{t}{\nu}\right)^{-\frac{\nu+d}{2}} t^{\frac{d}{2}-1}\log\left(1+\dfrac{t}{\nu^{\prime}}\right) dt, \] where \[K(d,\nu)=\dfrac{\Gamma\left(\dfrac{\nu+d}{2}\right)}{\Gamma\left(\dfrac{\nu}{2}\right)\sqrt{(\pi\nu)^d}}.\] Thus, this expression only involves a 1-dimensional integral, making the calculation of the Kullback Leibler divergence scalable and tractable for any dimension \(d=1,2,3,\dots\). The following R code shows the implementation of this divergence for several dimensions.
References
rm(list=ls())
# Package required for a nice format of the tables
library(knitr) ## Warning: package 'knitr' was built under R version 3.3.3# Normalising constant
K <- function(d,nu) (gamma(0.5*(nu+d))/( gamma(0.5*nu)*sqrt((pi*nu)^d) ))
# Kullback Liebler divergence
DKL <- function(nu,nup){
val1 <- log(K(d,nu)/K(d,nup)) -0.5*(nu+d)*(digamma(0.5*(nu+d))-digamma(0.5*nu))
tempf <- Vectorize(function(t) (1+t/nu)^(-0.5*(nu+d))*t^(0.5*d-1)*log(1+t/nup))
int<- integrate(tempf,0,Inf,rel.tol = 1e-9)$value
val2 <- (K(d,nu)*pi^(0.5*d)*0.5*(nup+d)/gamma(0.5*d))*int
return(val1+val2)
}
# Kullback Leibler in one dimension
d=1 # dimension
D <- matrix(0,ncol=3,nrow=29)
D[1,] <- c(1,0,DKL(1,2))
for(i in 2:29) D[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(D) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(D) <- 1:29
print(kable(D,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.000000e+00 1.130965e-01
## 2 6.209980e-02 1.916666e-02
## 3 1.363916e-02 5.897187e-03
## 4 4.700441e-03 2.411720e-03
## 5 2.047079e-03 1.169834e-03
## 6 1.032711e-03 6.364048e-04
## 7 5.768258e-04 3.760650e-04
## 8 3.472837e-04 2.365812e-04
## 9 2.214986e-04 1.563224e-04
## 10 1.478908e-04 1.074607e-04
## 11 1.024925e-04 7.631950e-05
## 12 7.326122e-05 5.570451e-05
## 13 5.375104e-05 4.161431e-05
## 14 4.032643e-05 3.171703e-05
## 15 3.084435e-05 2.459879e-05
## 16 2.399309e-05 1.937253e-05
## 17 1.894314e-05 1.546500e-05
## 18 1.515480e-05 1.249589e-05
## 19 1.226799e-05 1.020705e-05
## 20 1.003705e-05 8.419554e-06
## 21 8.290989e-06 7.007113e-06
## 22 6.908659e-06 5.879059e-06
## 23 5.802797e-06 4.969316e-06
## 24 4.909621e-06 4.229056e-06
## 25 4.181876e-06 3.621740e-06
## 26 3.584116e-06 3.119700e-06
## 27 3.089448e-06 2.701765e-06
## 28 2.677253e-06 2.351571e-06
## 29 2.331567e-06 2.056351e-06# Kullback Leibler in two dimensions
rm(d)
d=2 # dimension
D <- matrix(0,ncol=3,nrow=29)
D[1,] <- c(1,0,DKL(1,2))
for(i in 2:29) D[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(D) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(D) <- 1:29
print(kable(D,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.000000e+00 1.415927e-01
## 2 7.944154e-02 2.732554e-02
## 3 1.956127e-02 9.139055e-03
## 4 7.282898e-03 3.961126e-03
## 5 3.352862e-03 2.005447e-03
## 6 1.763603e-03 1.127459e-03
## 7 1.017649e-03 6.838292e-04
## 8 6.288552e-04 4.393997e-04
## 9 4.097277e-04 2.954756e-04
## 10 2.784711e-04 2.061445e-04
## 11 1.959103e-04 1.482667e-04
## 12 1.418505e-04 1.094059e-04
## 13 1.052417e-04 8.251520e-05
## 14 7.973102e-05 6.342075e-05
## 15 6.151054e-05 4.955533e-05
## 16 4.821469e-05 3.928773e-05
## 17 3.832778e-05 3.155168e-05
## 18 3.085192e-05 2.563259e-05
## 19 2.511425e-05 2.104082e-05
## 20 2.065128e-05 1.743421e-05
## 21 1.713759e-05 1.456933e-05
## 22 1.434074e-05 1.227021e-05
## 23 1.209212e-05 1.040777e-05
## 24 1.026762e-05 8.886050e-06
## 25 8.774727e-06 7.632841e-06
## 26 7.543654e-06 6.593187e-06
## 27 6.521166e-06 5.724832e-06
## 28 5.666240e-06 4.994967e-06
## 29 4.946970e-06 4.377892e-06# Kullback Leibler in three dimensions
rm(d)
d=3 # dimension
D <- matrix(0,ncol=3,nrow=29)
D[1,] <- c(1,0,DKL(1,2))
for(i in 2:29) D[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(D) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(D) <- 1:29
print(kable(D,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.000000e+00 1.551520e-01
## 2 8.850926e-02 3.208384e-02
## 3 2.313346e-02 1.129171e-02
## 4 9.021193e-03 5.087129e-03
## 5 4.306589e-03 2.654500e-03
## 6 2.331819e-03 1.528908e-03
## 7 1.377596e-03 9.458675e-04
## 8 8.680423e-04 6.178895e-04
## 9 5.749090e-04 4.213451e-04
## 10 3.962234e-04 2.974995e-04
## 11 2.821191e-04 2.162038e-04
## 12 2.064156e-04 1.609913e-04
## 13 1.545540e-04 1.223979e-04
## 14 1.180430e-04 9.474667e-05
## 15 9.172706e-05 7.450584e-05
## 16 7.236639e-05 5.940842e-05
## 17 5.786327e-05 4.795876e-05
## 18 4.682359e-05 3.914604e-05
## 19 3.829915e-05 3.227232e-05
## 20 3.163169e-05 2.684642e-05
## 21 2.635566e-05 2.251663e-05
## 22 2.213633e-05 1.902718e-05
## 23 1.872936e-05 1.618941e-05
## 24 1.595391e-05 1.386235e-05
## 25 1.367446e-05 1.193942e-05
## 26 1.178826e-05 1.033911e-05
## 27 1.021657e-05 8.998528e-06
## 28 8.898470e-06 7.868622e-06
## 29 7.786376e-06 6.910837e-06# Kullback Leibler in four dimensions
rm(d)
d=4 # dimension
D <- matrix(0,ncol=3,nrow=29)
D[1,] <- c(1,0,DKL(1,2))
for(i in 2:29) D[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(D) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(D) <- 1:29
print(kable(D,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.000000e+00 1.635764e-01
## 2 9.453489e-02 3.534458e-02
## 3 2.566131e-02 1.288059e-02
## 4 1.032586e-02 5.966897e-03
## 5 5.058505e-03 3.184816e-03
## 6 2.798317e-03 1.868943e-03
## 7 1.683192e-03 1.174484e-03
## 8 1.076901e-03 7.775137e-04
## 9 7.226231e-04 5.362981e-04
## 10 5.036958e-04 3.824496e-04
## 11 3.622049e-04 2.803748e-04
## 12 2.673282e-04 2.103906e-04
## 13 2.017137e-04 1.610568e-04
## 14 1.551273e-04 1.254411e-04
## 15 1.212923e-04 9.919120e-05
## 16 9.622739e-05 7.948973e-05
## 17 7.733310e-05 6.446368e-05
## 18 6.286835e-05 5.283823e-05
## 19 5.164053e-05 4.372748e-05
## 20 4.281619e-05 3.650418e-05
## 21 3.580232e-05 3.071672e-05
## 22 3.017012e-05 2.603498e-05
## 23 2.560495e-05 2.221429e-05
## 24 2.187277e-05 1.907101e-05
## 25 1.879743e-05 1.646569e-05
## 26 1.624477e-05 1.429133e-05
## 27 1.411160e-05 1.246499e-05
## 28 1.231775e-05 1.092180e-05
## 29 1.080038e-05 9.610588e-06# Kullback Leibler in fifty dimensions
rm(d)
d=50 # dimension
D <- matrix(0,ncol=3,nrow=29)
D[1,] <- c(1,0,DKL(1,2))
for(i in 2:29) D[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(D) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(D) <- 1:29
print(kable(D,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.0000000000 0.2074545331
## 2 0.1255888908 0.0552570704
## 3 0.0414442942 0.0243147651
## 4 0.0199175133 0.0133356785
## 5 0.0114607688 0.0082832698
## 6 0.0073359426 0.0055738121
## 7 0.0050388786 0.0039660586
## 8 0.0036393210 0.0029408104
## 9 0.0027292267 0.0022508204
## 10 0.0021074666 0.0017666052
## 11 0.0016659117 0.0014152297
## 12 0.0013423871 0.0011531719
## 13 0.0010991686 0.0009532154
## 14 0.0009123396 0.0007976672
## 15 0.0007661720 0.0006746422
## 16 0.0006499973 0.0005759313
## 17 0.0005563838 0.0004957252
## 18 0.0004800343 0.0004298258
## 19 0.0004170955 0.0003751422
## 20 0.0003647146 0.0003293608
## 21 0.0003207454 0.0002907246
## 22 0.0002835504 0.0002578810
## 23 0.0002518640 0.0002297764
## 24 0.0002246969 0.0002055813
## 25 0.0002012673 0.0001846359
## 26 0.0001809516 0.0001664106
## 27 0.0001632478 0.0001504769
## 28 0.0001477488 0.0001364853
## 29 0.0001341216 0.0001241487