Kernel Wasserstein Distance

TL;DR

This paper introduces a kernel-based method to compute the Wasserstein distance in a nonlinear feature space, enabling improved clustering of medical images with artifacts.

Contribution

It develops the first kernel method for Wasserstein distance, allowing nonlinear data analysis with applications in medical image artifact detection.

Findings

Kernel Wasserstein distance outperforms classical methods in artifact detection.

The approach effectively clusters CT slices with dental artifacts.

Improves analysis of corrupted data using kernel techniques.

Abstract

The Wasserstein distance is a powerful metric based on the theory of optimal transport. It gives a natural measure of the distance between two distributions with a wide range of applications. In contrast to a number of the common divergences on distributions such as Kullback-Leibler or Jensen-Shannon, it is (weakly) continuous, and thus ideal for analyzing corrupted data. To date, however, no kernel methods for dealing with nonlinear data have been proposed via the Wasserstein distance. In this work, we develop a novel method to compute the L2-Wasserstein distance in a kernel space implemented using the kernel trick. The latter is a general method in machine learning employed to handle data in a nonlinear manner. We evaluate the proposed approach in identifying computerized tomography (CT) slices with dental artifacts in head and neck cancer, performing unsupervised hierarchical clustering on the resulting Wasserstein distance matrix that is computed on imaging texture features extracted from each CT slice. Our experiments show that the kernel approach outperforms classical non-kernel approaches in identifying CT slices with artifacts.

Tables

Table 1. Table 1: GLCM-based 25 texture features used in this study.

No	Features	No	Features
1	Auto-correlation	14	Inverse Difference Moment
2	Joint Average	15	First Informal Correlation
3	Cluster Prominence	16	Second Informal Correlation
4	Cluster Shade	17	Inverse Difference Moment Normalized
5	Cluster Tendency	18	Inverse Difference Normalized
6	Contrast	19	Inverse Variance
7	Correlation	20	Sum Average
8	Difference Entropy	21	Sum Entropy
9	Dissimilarity	22	Sum Variance
10	Difference Variance	23	Haralick Correlation
11	Joint Energy	24	Joint Maximum
12	Joint Entropy	25	Joint Variance
13	Inverse Difference

Equations

$$W_p^p\left(\mu ,\nu \right)={\inf }_{\pi \in \Pi (\mu ,\nu )}{\int }_{{\mathrm{R}}^d\times {\mathrm{R}}^d}{\Vert x-y\Vert }^{p}d\pi (x,y),$$

$$\mathbf{K}={\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi },$$

$$k({{\bf x}_{i},{\bf x}_{j}})={\rm exp}\big{(}-\gamma||{{\bf x}_{i}-{\bf x}_{j}}||^{2}\big{)},$$

$$\mu =\frac{1}{n}\sum _{i=1}^n\varphi ({\mathbf{x}}_i)=\mathbf{\Phi s},\mathrm{and}\mathbf{\Sigma }=\frac{1}{n}\sum _{i=1}^n(\varphi ({\mathbf{x}}_i)-\mu )(\varphi ({\mathbf{x}}_i)-\mu {)}^{\mathrm{T}}=\mathbf{\Phi J}{\mathbf{J}}^{\mathrm{T}}{\mathbf{\Phi }}^{\mathrm{T}},$$

$$\mathbf{W}=\mathbf{\Phi J}=\frac{1}{\sqrt{n}}[(\varphi ({\mathbf{x}}_1)-\mu ),\cdots ,(\varphi ({\mathbf{x}}_n)-\mu )].$$

$$W_2({\nu }_1,{\nu }_2{)}^2=\Vert {\mathbf{m}}_1-{\mathbf{m}}_2{\Vert }^2+\mathrm{tr}({\mathbf{C}}_1+{\mathbf{C}}_2-2({\mathbf{C}}_1^{\frac{1}{2}}{\mathbf{C}}_2{\mathbf{C}}_1^{\frac{1}{2}}{)}^{\frac{1}{2}}),$$

$$\mathrm{tr}({\mathbf{C}}_1+{\mathbf{C}}_2-2({\mathbf{C}}_2{\mathbf{C}}_1{)}^{\frac{1}{2}}).$$

$$W_2(k{\nu }_1,k{\nu }_2{)}^2=\Vert {\mu }_1-{\mu }_2{\Vert }^2+\mathrm{tr}({\mathbf{\Sigma }}_1+{\mathbf{\Sigma }}_2-2({\mathbf{\Sigma }}_1^{\frac{1}{2}}{\mathbf{\Sigma }}_2{\mathbf{\Sigma }}_1^{\frac{1}{2}}{)}^{\frac{1}{2}}),$$

$$\Vert {\mu }_1-{\mu }_2{\Vert }^2=\Vert {\mu }_1{\Vert }^2-2{\mu }_1^{\mathrm{T}}{\mu }_2+\Vert {\mu }_2{\Vert }^2.$$

$$\Vert {\mu }_1-{\mu }_2{\Vert }^2=\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^{n}k\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)-\frac{2}{nm}\sum _{i=1}^n\sum _{j=1}^{m}k\left({\mathbf{x}}_i,{\mathbf{y}}_j\right)+\frac{1}{m^2}\sum _{i=1}^m\sum _{j=1}^{m}k\left({\mathbf{y}}_i,{\mathbf{y}}_j\right).$$

$$=$$

$${\mathbf{\Phi }}_2\mathbf{R}{\mathbf{\Phi }}_1^{\mathrm{T}}\tilde{\mathbf{P}}=\tilde{\mathbf{P}}\tilde{\mathbf{\Lambda }}.$$

$$\mathrm{tr}{\left({\mathbf{\Phi }}_2\mathbf{R}{\mathbf{\Phi }}_1^{\mathrm{T}}\right)}^{\frac{1}{2}}=\mathrm{tr}{\left({\mathbf{\Phi }}_1^{\mathrm{T}}{\mathbf{\Phi }}_2\mathbf{R}\right)}^{\frac{1}{2}}=\mathrm{tr}{\left({\mathbf{K}}_{12}\mathbf{R}\right)}^{\frac{1}{2}}=\mathrm{tr}{\left({\mathbf{K}}_{12}{\mathbf{J}}_2{\mathbf{J}}_2^{\mathrm{T}}{\mathbf{K}}_{21}{\mathbf{J}}_1{\mathbf{J}}_1^{\mathrm{T}}\right)}^{\frac{1}{2}}.$$

$$=\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^{n}k\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)-\frac{2}{nm}\sum _{i=1}^n\sum _{j=1}^{m}k\left({\mathbf{x}}_i,{\mathbf{y}}_j\right)+\frac{1}{m^2}\sum _{i=1}^m\sum _{j=1}^{m}k\left({\mathbf{y}}_i,{\mathbf{y}}_j\right)+$$

$$\mathrm{tr}\left({\mathbf{J}}_1{\mathbf{J}}_1^{\mathrm{T}}{\mathbf{K}}_{11}\right)+\mathrm{tr}\left({\mathbf{J}}_2{\mathbf{J}}_2^{\mathrm{T}}{\mathbf{K}}_{22}\right)-2\mathrm{tr}{\left({\mathbf{K}}_{12}{\mathbf{J}}_2{\mathbf{J}}_2^{\mathrm{T}}{\mathbf{K}}_{21}{\mathbf{J}}_1{\mathbf{J}}_1^{\mathrm{T}}\right)}^{\frac{1}{2}}.$$

$$D_{\mathrm{KL}}(P\Vert Q)={\int }_{-\infty }^{\infty }p(x)\log \frac{p(x)}{q(x)}dx,$$

$$\displaystyle D_{\rm KL}({N}_{1}||{N}_{2})=\frac{1}{2}\big{\{}({\bf\mu}_{1}-{\bf\mu}_{2})^{\bf{\rm T}}{\bf\Sigma}_{2}^{-1}({\bf\mu}_{1}-{\bf\mu}_{2})+{\rm log}\frac{|{\bf\Sigma}_{2}|}{|{\bf\Sigma}_{1}|}+{\rm tr}[{\bf\Sigma}_{1}{\bf\Sigma}_{2}^{-1}]-l\big{\}},$$

$$\mathbf{H}=\mathbf{\Sigma }+\rho {\mathbf{I}}_l=\mathbf{\Phi J}{\mathbf{J}}^{\mathrm{T}}{\mathbf{\Phi }}^{\mathrm{T}}+\rho {\mathbf{I}}_l=\mathbf{W}{\mathbf{W}}^{\mathrm{T}}+\rho {\mathbf{I}}_l=\mathbf{\Phi S}{\mathbf{\Phi }}^{\mathrm{T}}+\rho {\mathbf{I}}_l,$$

$$(\mathbf{A}+\beta {\mathbf{U}{\mathbf{V}}^{\mathrm{T}})}^{-1}={\mathbf{A}}^{-1}-\beta {\mathbf{A}}^{-1}\mathbf{U}{\mathbf{\Sigma }}^{-1}{\mathbf{V}}^{\mathrm{T}}{\mathbf{A}}^{-1}.$$

$${\mathbf{H}}^{-1}$$

$$2D_{\mathrm{KL}}(N_1||N_2)=({\mu }_1-{\mu }_2{)}^{\mathrm{T}}{\mathbf{H}}_2^{-1}({\mu }_1-{\mu }_2)+\log \frac{|{\mathbf{H}}_2|}{|{\mathbf{H}}_1|}+\mathrm{tr}[{\mathbf{H}}_1{\mathbf{H}}_2^{-1}]-l.$$

$$({\mu }_1-{\mu }_2{)}^{\mathrm{T}}{\mathbf{H}}_2^{-1}({\mu }_1-{\mu }_2)={\mu }_1^{\mathrm{T}}{\mathbf{H}}_2^{-1}{\mu }_1+{\mu }_2^{\mathrm{T}}{\mathbf{H}}_2^{-1}{\mu }_2-{\mu }_1^{\mathrm{T}}{\mathbf{H}}_2^{-1}{\mu }_2-{\mu }_2^{\mathrm{T}}{\mathbf{H}}_2^{-1}{\mu }_1.$$

$${\mu }_i^{\mathrm{T}}{\mathbf{H}}_j^{-1}{\mu }_k$$

$$|{\mathbf{H}{}_1}^{-1}|$$

$$|{\mathbf{H}}_1|=\frac{1}{|H{{}_1}^{-1}|}=\frac{{\rho }^l}{|{\mathbf{I}}_n-{\mathbf{K}}_{11}{\mathbf{B}}_1|}.$$

$$log|{\mathbf{H}}_1|=\log \frac{{\rho }^l}{|{\mathbf{I}}_n-{\mathbf{K}}_{11}{\mathbf{B}}_1|}=l\log \rho -\log |{\mathbf{I}}_n-{\mathbf{K}}_{11}{\mathbf{B}}_1|.$$

$$\log \frac{|{\mathbf{H}}_2|}{|{\mathbf{H}}_1|}=log|{\mathbf{H}}_2|-log|{\mathbf{H}}_1|=\log |{\mathbf{I}}_n-{\mathbf{K}}_{11}{\mathbf{B}}_1|-\log |{\mathbf{I}}_m-{\mathbf{K}}_{22}{\mathbf{B}}_2|.$$

$$\mathrm{tr}[{\mathbf{H}}_1{\mathbf{H}}_2^{-1}]$$

$$={\rho }^{-1}({\theta }_{121}+{\theta }_{222}-{\theta }_{122}-{\theta }_{221})+\log |{\mathbf{I}}_n-{\mathbf{K}}_{11}{\mathbf{B}}_1|-\log |{\mathbf{I}}_m-{\mathbf{K}}_{22}{\mathbf{B}}_2|$$

$$+{\rho }^{-1}\mathrm{tr}[{\mathbf{S}}_1{\mathbf{K}}_{11}]-{\rho }^{-1}\mathrm{tr}[{\mathbf{S}}_1{\mathbf{K}}_{12}{\mathbf{B}}_2{\mathbf{K}}_{21}]-\mathrm{tr}[{\mathbf{B}}_2{\mathbf{K}}_{22}].$$

Full text

1 Introduction

Optimal mass transport (OMT) theory is an old research area with its roots in civil engineering (Monge 1781) and economics (Kantorovich 1942) [1]. Recently there has been an ever increasing growth in OMT research both theoretically and practically, with impact on numerous fields including medical imaging analysis, statistical physics, machine learning, and genomics [2, 3, 4, 5]. The classical OMT problem formulated by Monge in 1781 concerns finding the optimal way via the minimization of a transportation cost required to move a pile of soil from one site to another [6, 7, 8, 9]. Let $X$ and $Y$ denote two probability spaces with measures $\mu$ and $\nu$, respectively, and let $c(x,y)$ denote the transportation cost for moving one unit of mass from $x\in X$ to $y\in Y$. Then the OMT problem seeks to find a (measurable) transport map $T:X\rightarrow Y$ that minimizes the total transportation cost $\int_{X}c(x,T(x)){\rm\mu}(dx)$. In 1942, Kantorovich proposed a relaxed formulation that transforms the Monge’s nonlinear problem to a linear programming problem [8]. Based on the Kantorovich’s formulation, the ${\it L^{p}}$- Wasserstein distance between $\mu$ and $\nu$ on $\mathcal{R}^{d}$ is defined as:

[TABLE]

where $\Pi(\mu,\nu)$ is the set of all joint probability measures $\pi$ on $X\times Y$ whose marginals are $\mu$ and $\nu$. In particular, in this study, we focus on $\it L^{\rm 2}$- Wasserstein distance in which the squared Euclidean distance $c(x,y)={\|{x}-{y}\|}^{2}$ is the cost function [10]. Before introducing our proposed kernel Wasserstein distance, we first review some background on the kernel method.

The kernel space method is based on the following idea. Suppose that we are given a data set of $n$ samples in a native space, denoted by ${\bf X}=[{\bf x}_{1},{\bf x}_{2},\cdots,{\bf x}_{\it n}]\in\mathcal{R}^{d}$. The input data can be mapped (transformed) into a higher dimensional feature space (called the kernel space) via a nonlinear mapping function $\phi$ [11, 12]. Let ${\bf\Phi}$ be ${\bf\Phi}_{l\times n}=[\phi({\bf x}_{1}),\phi({\bf x}_{2}),\cdots,\phi(\bf x_{\it n})]$, i.e., the transformed data, where $l$ is the number of features in the feature (kernel) space with $l>d$. To avoid complex data handling in the feature space, and to avoid the explicit computation of the mapping funciton $\phi$, typically one applies the kernel trick. More precisely, given any positive definite kernel function $k$, one can find an associated mapping function $\phi$ such that $k({\bf x},{\bf y})=<\phi{({\bf x})},\phi{({\bf y})}>$ with ${\bf x},{\bf y}\in\mathcal{R}^{d}$ [13, 14]. The resulting kernel Gram matrix ${\bf K}$ is defined as:

[TABLE]

where the $ij$th element, $\phi({\bf x}_{i})^{\rm T}\phi({\bf x}_{j})$, is computed using a kernel function $k({{\bf x}_{i},{\bf x}_{j}})=<\phi{({\bf x}_{i})},\phi{({\bf x}_{j})}>$. Common choices of kernel functions are the polynomial and radial basis function (RBF) kernels. In this study, the following RBF kernel will be employed:

[TABLE]

where $\gamma>0$ controls the kernel width. We will fix $\gamma=1$ in what follows. The mean and the covariance matrix in the feature space are estimated as:

[TABLE]

where ${\bf s}_{n\times 1}=\frac{1}{n}{\bf\vec{1}}^{\rm T}$, $\bf{J}=\frac{\rm 1}{\sqrt{\it{n}}}(\bf I_{\it n}-s{\bf\vec{1}})$, and ${\bf\vec{1}}=[1,1,\cdots,1]$. Then in Eq. (4), denoting $\bf{\Phi J}$ by ${\bf W}$ we have

[TABLE]

Of note, these equations are used to compute the kernel Wasserstein distance as described in the following section.

2 Methods

In this section, we introduce the classical L2-Wasserstein distance between Gaussian measures and then propose a novel approach to compute L2-Wasserstein distance in the kernel space (denoted as kernel L2-Wasserstein distance). For comparison, we also provide a brief review of the Kullback-Leibler distance in the kernel space that we proposed in [15]. The code was implemented using MATLAB R2018b.

2.1 L2-Wasserstein Distance

For two Gaussian measures, $\nu_{1}=N_{1}({\bf m}_{\rm 1},{\bf C}_{\rm 1})$ and $\nu_{2}=N_{2}({\bf m}_{\rm 2},{\bf C}_{\rm 2})$ on $\mathcal{R}^{d},$ the L2-Wasserstein distance between the two distributions may be computed as follows [16]:

[TABLE]

where ${\rm tr}$ is the trace. The term ${\rm tr}({\bf C}_{\rm 1}+{\bf C}_{\rm 2}-2({\bf C}_{\rm 1}^{\frac{1}{2}}{\bf C}_{\rm 2}{\bf C}_{\rm 1}^{\frac{1}{2}})^{\frac{1}{2}})$ can be expressed as follows [17, 18, 19]:

[TABLE]

For convenience, we sketch the proof of this fact.

Proof.

By the property of trace, ${\rm tr}({\bf C}_{1}+{\bf C}_{2}-2({\bf C}_{2}{\bf C}_{1})^{\frac{1}{2}})={\rm tr}({\bf C}_{1})+{\rm tr}({\bf C}_{2})-2{\rm tr}({\bf C}_{2}{\bf C}_{1})^{\frac{1}{2}}.$ Therefore, we need to prove that ${\rm tr}({\bf C}_{2}{\bf C}_{1})^{\frac{1}{2}}={\rm tr}({\bf C}_{\rm 1}^{\frac{1}{2}}{\bf C}_{\rm 2}{\bf C}_{\rm 1}^{\frac{1}{2}})^{\frac{1}{2}}$. Note that ${\bf C}_{\rm 1}$ and ${\bf C}_{\rm 2}$ are symmetric positive semidefinite, and ${\bf C}_{\rm 2}{\bf C}_{\rm 1}$ is diagonalizable and has nonnegative eigenvalues [20]. The eigenvalue decomposition of ${\bf C}_{2}{\bf C}_{1}$ can be computed as ${\bf C}_{2}{\bf C}_{1}{\bf P}={\bf P\Lambda}$ where ${\bf P}$ and $\Lambda$ are the estimated eigenvector and eigenvalue matrices, respectively. Multiplying both sides by ${\bf C}_{1}^{\frac{1}{2}}$, we have ${\bf C}_{1}^{\frac{1}{2}}{\bf C}_{2}{\bf C}_{1}^{\frac{1}{2}}{\bf C}_{1}^{\frac{1}{2}}{\bf P}={\bf C}_{1}^{\frac{1}{2}}{\bf P\Lambda}.$ That is, ${\bf C}_{1}^{\frac{1}{2}}{\bf C}_{2}{\bf C}_{1}^{\frac{1}{2}}$ has an eigenvector matrix ${\bf C}_{1}^{\frac{1}{2}}{\bf P}$ and an eigenvalue matrix ${\bf\Lambda}$ which is the same as that of ${\bf C}_{2}{\bf C}_{1}$ [21]. Let $\lambda_{1},\lambda_{2},...,\lambda_{k}$ be the distinct eigenvalues of ${\bf C}_{2}{\bf C}_{1}$. Then the eigenvalues of $({\bf C}_{2}{\bf C}_{1})^{\frac{1}{2}}$ are $\sqrt{\lambda_{1}},\sqrt{\lambda_{2}},...,\sqrt{\lambda_{k}}$. Therefore, ${\rm tr}({\bf C}_{2}{\bf C}_{1})^{\frac{1}{2}}={\rm tr}({\bf C}_{\rm 1}^{\frac{1}{2}}{\bf C}_{\rm 2}{\bf C}_{\rm 1}^{\frac{1}{2}})^{\frac{1}{2}}$. ∎

In particular, when ${\bf C}_{1}={\bf C}_{2}$, we have $W_{2}(\nu_{1},\nu_{2})^{2}=\|{\bf m}_{\rm 1}-{\bf m}_{\rm 2}\|^{2}.$

2.2 L2-Wasserstein Distance in Kernel Space

Suppose that we are given two Gaussian measures, $k\nu_{1}$ and $k\nu_{2}\in\mathcal{R}^{\it l}$, in the kernel space with mean ${\bf\mu}_{\it i}$ and covariance matrix ${\bf\Sigma}_{\it i}$, for ${\it i}$=1 and 2, and two sets of sample data in the native space, ${\bf X}=[{\bf x}_{1},{\bf x}_{2},\cdots,{\bf x}_{\it n}],\;{\bf Y}=[{\bf y}_{1},{\bf y}_{2},\cdots,{\bf y}_{\it m}]\in\mathcal{R}^{d}$ associated with $k\nu_{1}$ and $k\nu_{2}$, respectively. Then, as in Eq. (6), the L2-Wasserstein distance between the two distributions is given by:

[TABLE]

where the definitions of ${\bf\mu}_{\it i}$ and ${\bf\Sigma}_{\it i}$ are shown in Eq. (4). Note that ${\rm tr}({\bf\Sigma}_{\rm 1}^{\frac{1}{2}}{\bf\Sigma}_{\rm 2}{\bf\Sigma}_{\rm 1}^{\frac{1}{2}})^{\frac{1}{2}}={\rm tr}\left({\bf\Sigma}_{2}{\bf\Sigma}_{1}\right)^{\frac{1}{2}}$. The first term, $\|{\bf\mu}_{\rm 1}-{\bf\mu}_{\rm 2}\|^{2}$, in Eq. (8) can be obtained as follows [22]:

[TABLE]

Via a simple computation, Eq. (9) can be expressed as:

[TABLE]

Now we simplify the last term in Eq. (8) by using Eq. (4) :

$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm tr}({\bf\Sigma}_{1}+{\bf\Sigma}_{2}-2({\bf\Sigma}_{2}{\bf\Sigma}_{1})^{\frac{1}{2}})={\rm tr}\left({\bf\Sigma}_{1}\right)+{\rm tr}\left({\bf\Sigma}_{2}\right)-2{\rm tr}\left({\bf\Sigma}_{2}{\bf\Sigma}_{1}\right)^{\frac{1}{2}}$

[TABLE]

where ${\bf R}={\bf J}_{2}{\bf J}_{2}^{\rm T}{\bf K}_{21}{\bf J}_{1}{\bf J}_{1}^{\rm T}$ and ${\bf K}_{\it ij}={\bf\Phi}_{i}^{\rm T}{\bf\Phi_{\it j}}$. Note that ${\bf\Sigma_{\rm 1}}$ and ${\bf\Sigma_{\rm 2}}$ are symmetric positive semidefinite. Therefore, ${\bf\Sigma_{\rm 2}}{\bf\Sigma_{\rm 1}}={\bf\Phi}_{2}{\bf R}{\bf\Phi}_{1}^{\rm T}$ is diagonalizable and has nonnegative eigenvalues. Suppose that the eigenvector and eigenvalue matrices of ${\bf\Phi}_{2}{\bf R}{\bf\Phi}_{1}^{\rm T}$ are $\bf\tilde{P}$ and ${\bf\tilde{\Lambda}}$:

[TABLE]

By multiplying both sides by ${\bf\Phi}_{1}^{\rm T}$, we have ${\bf\Phi}_{1}^{\rm T}{\bf\Phi}_{2}{\bf R}{\bf\Phi}_{1}^{\rm T}{\bf\tilde{P}}={\bf\Phi}_{1}^{\rm T}{\bf\tilde{P}\tilde{\Lambda}}$. That is, the eigenvalue matrix of ${\bf\Phi}_{1}^{\rm T}{\bf\Phi}_{2}{\bf R}$ is the same as that of ${\bf\Phi}_{2}{\bf R}{\bf\Phi}_{1}^{\rm T}$ with ${\bf\tilde{\Lambda}}$. Therefore, the following equations hold:

[TABLE]

Finally, using Eq. (10) and Eq. (13), the kernel L2-Wasserstein distance can be expressed as:

$W_{2}(k\nu_{1},k\nu_{2})^{2}$=$\|{\bf\mu}_{1}-{\bf\mu}_{2}\|^{2}+{\rm tr}({\bf\Sigma}_{1}+{\bf\Sigma}_{2}-2({\bf\Sigma}_{1}^{\frac{1}{2}}{\bf\Sigma}_{2}{\bf\Sigma}_{1}^{\frac{1}{2}})^{\frac{1}{2}})$

[TABLE]

In a special case, when ${\bf\Sigma}_{1}={\bf\Sigma}_{2}$, we have $W_{2}(k\nu_{1},k\nu_{2})^{2}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}{\it k}\left({\bf x}_{i},{\bf x}_{j}\right)-\frac{2}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}{\it k}\left({\bf x}_{i},{\bf y}_{j}\right)+\frac{1}{m^{2}}\sum_{i=1}^{m}\sum_{j=1}^{m}{\it k}\left({\bf y}_{i},{\bf y}_{j}\right)$.

2.3 Kullback-Leibler Divergence in Kernel Space

Kullback-Leibler (KL) divergence is another type of method to compare two probability distributions [23]. It is not a distance measure due to its asymmetric nature. Let $P$ and $Q$ be two continuous probability distributions with the corresponding probability densities $p(x)$ and $q(x)$, respectively. Then the KL divergence or relative entropy of $P$ and $Q$ over the same variable $x$ is defined as:

[TABLE]

where $D_{\mathrm{KL}}(P\|Q)$ equals zero if and only if $P=Q$. Given two Gaussian measures in the kernel space, ${N}_{1}({\bf\mu}_{1},{\bf\Sigma}_{1})$ and ${N}_{2}({\bf\mu}_{2},{\bf\Sigma}_{2})$ with both $\in\mathcal{R}^{\it l}$, the KL divergence may be computed to be:

[TABLE]

where $|{\bf\Sigma_{\it i}}|$ is the determinant of covariance matrix ${\bf\Sigma}_{\it i}$. Note that $l$ is the number of features in the kernel space and indeed this is unknown. Importantly, this variable is canceled out when Eq. (16) is completely solved, which will be explained later. In a special case, when {\bf\Sigma}={\bf\Sigma}_{1}={\bf\Sigma}_{2},D_{\rm KL}({N}_{1}||{N}_{2})=\frac{1}{2}\big{\{}({\bf\mu}_{1}-{\bf\mu}_{2})^{\bf{\rm T}}{\bf\Sigma}^{-1}({\bf\mu}_{1}-{\bf\mu}_{2})\}.

Singularity Problem

In many real problems, the number of samples is considerably smaller than the number of features, leading to the covariance matrix being singular, and therefore non-invertible. Typically, one deals with data in a higher dimensional space as in this study. To avoid the singularity problem, several methods have been proposed [24, 25]. In this work, we employ a simple method by adding some positive values to the diagonal elements of the covariance matrix [26]. Therefore, the modified covariance matrix is of full rank and invertible. Set

[TABLE]

where $\bf W=\Phi J$ as in Eq. (5), $\bf S=JJ^{\rm T}$, $l$ is the number of features, and ${\bf I}_{\it l}$ is an $l\times l$ identity matrix. In this study, $\rho=0.1$ is used.

In the next computation, we will employ the Woodbury formula [27], which we now review. Let $\bf A$ be a square $r\times r$ invertible matrix, and let $\bf U$ and $\bf V$ be matrices of size $r\times k$ with $k\leq r$. Assume that the $k\times k$ matrix ${\bf\Sigma}={\bf I}_{k}+\beta{\bf V^{\rm T}A^{\rm-1}U}$ is invertible, where $\beta$ is an arbitrary scalar. Then the Woodbury formula states that

[TABLE]

Accordingly, utilizing this formula, we can compute the inverse of $\bf H$:

[TABLE]

where $\bf B=JM^{\rm-1}J^{\rm T}$, ${\bf M}=\rho{\bf I}_{n}+{\bf W^{\rm T}W}=\rho{\bf I}_{n}+{\bf J^{\rm T}}{\bf\Phi}^{\rm T}{\bf\Phi}{\bf J}=\rho{\bf I}_{n}+{\bf J^{\rm T}KJ}$, and $n$ is the number of samples. In ${\bf H^{\rm-1}}$, some mapping functions are still left. These will be replaced with kernel functions when kernel KL divergence is completely solved.

Calculation of Kernel KL Divergence

In the Experiments section below, we will compare Wasserstein and KL divergence based kernel methods. Accordingly, we sketch the necessary theory for the kernel KL divergence approach.

Suppose that we are given two Gaussian measures in the kernel space, ${N}_{1}({\bf\mu}_{1},{\bf\Sigma}_{1}),\;{N}_{2}({\bf\mu}_{2},{\bf\Sigma}_{2})\in\mathcal{R}^{\it l}$, consisting of $n$ and $m$ samples, respectively. Assume that ${\bf\Sigma}_{1}$ and ${\bf\Sigma}_{2}$ are singular in the higher dimensional space. Let ${\bf H}_{1}$ and ${\bf H}_{2}$ denote the approximate covariance matrices for the two distributions as in Eq. (17). Then, the kernel KL divergence is expressed as follows:

[TABLE]

We now solve each term separately: (1) $(\mu_{1}-\mu_{2})^{\bf{\rm T}}{\bf H}_{2}^{-1}(\mu_{1}-\mu_{2})$, (2) ${\rm log}\frac{|{\bf H}_{2}|}{|{\bf H}_{1}|}$, and (3) ${\rm tr}[{\bf H}_{1}{\bf H}_{2}^{-1}]$. The first term consists of four sub-terms:

[TABLE]

Substituting Eq. (4) and Eq. (Singularity Problem) into each sub-term $\mu_{i}^{{\rm T}}{\bf H}_{j}^{-1}\mu_{k}$, we have

[TABLE]

where $\theta_{ijk}={\bf s}_{i}^{\rm T}{\bf K}_{ik}{\bf s}_{k}-{\bf s}_{i}^{\rm T}{\bf K}_{ij}{\bf B}_{j}{\bf K}_{jk}{\bf s}_{k}$. As a result, all the mapping functions in the first term can be replaced with kernel functions. For the second term, we should compute the determinant of ${\bf H}$. To accomplish this, we use a simple trick by computing the determinant of ${\bf H}^{-1}$ instead of ${\bf H}$.

[TABLE]

where ${\bf Q_{\rm 1}}={\bf\Phi}_{1}{\bf B}_{1}$. Now we compute $|{\bf H}_{1}|$ as follows:

[TABLE]

By taking $\mathrm{logarithm}$ of $|{\bf H}_{1}|$, we have

[TABLE]

Therefore, we have the second term composed of kernel functions:

[TABLE]

The third term can be replaced with kernel functions using properties of trace:

[TABLE]

Consequently, we have solved kernel KL divergence by replacing all mapping functions with kernel functions and it is expressed as:

$2D_{\rm{KL}}({N}_{1}||{N}_{2})=(\mu_{1}-\mu_{2})^{\bf{\rm T}}{\bf H}_{2}^{-1}(\mu_{1}-\mu_{2})+{\rm log}\frac{|{\bf H}_{2}|}{|{\bf H}_{1}|}+{\rm tr}[{\bf H}_{1}{\bf H}_{2}^{-1}]-{\it l}$

[TABLE]

Note that the $l$ is canceled out. Moreover, since $D_{\mathrm{KL}}(P\|Q)\neq D_{\mathrm{KL}}(Q\|P)$, in this study an average value of two KL measures is used: J_{\mathrm{KL}}(P\|Q)=\frac{1}{2}\big{\{}D_{\mathrm{KL}}(P\|Q)+D_{\mathrm{KL}}(Q\|P)\big{\}} [28].

3 Experiments

3.1 Data

We investigated the utility of kernel L2-Wasserstein distance to identify slices with dental artifacts in computerized tomography (CT) scans in head and neck cancer. Serious image degradation caused by metallic fillings or crowns in teeth is a common problem in CT images. We analyzed 1164 axial slices from 44 CT scans that were collected from 44 patients with head and cancer who were treated in our institution. This retrospective study was approved by the institutional review board and informed consent was obtained from all patients. Before the analysis, each CT slice was labeled as noisy or clean based on the presence of dental artifacts by a medical imaging expert, resulting in 276 noisy and 888 clean slices. Figure 1 shows representative noisy and clean CT slices from two different scans.

3.2 Texture Features

Intensity thresholding (at the 5th percentile) was performed to exclude air voxels on each CT slice. The Computational Environment for Radiological Research (CERR) radiomics toolbox was then used to calculate the gray-level co-occurrence matrix (GLCM) from the remaining voxels using 64 gray levels and with a neighborhood of 8 voxels across all 4 directions in 2D [29, 30]. A total of 25 scalar features were then extracted from the GLCM (listed in Table 1). For further information on the GLCM features used, see https://github.com/cerr/CERR/wiki/Radiomics. Each feature was normalized between 0 and 1 for further analysis.

3.3

Experimental Results

For 1164 CT slices, we computed kernel Wasserstein distance between each pair of slices on 25 GLCM-based texture features. After that, we conducted unsupervised hierarchical clustering using the resulting distance matrix. Figure 2(a) shows a heatmap of the symmetric distance matrix for each pair of slices and Figure 2(b) presents the hierarchical clustering result. We identified two clusters: Cluster 1 with blue lines and Cluster 2 with red lines, consisting of 666 and 498 slices, respectively. In Figure (b), the bar under the hierarchical graph indicates the actual labels with blue for clean slices and red for noisy slices. As a result, Cluster 1 has 658 clean and 8 noisy slices whereas Cluster 2 has 230 clean and 268 noisy slices (table in Figure 2(b)). That is, Cluster 1 and 2 were significantly enriched for clean and noisy slices, respectively, with a Chi-square test p-value < 0.0001. Prediction rates were 97.10% and 74.10% for noisy and clean slices, respectively, and overall prediction rate was 79.6%. In Figure 2(a), the order of slices is the same as that shown in Figure 2(b) and the two clusters were divided by the black dot lines; the left bottom block represents Cluster 1 and the right top block for Cluster 2. The areas with blue color indicate close distance between slices whereas the areas with red color indicate far distance. Not surprisingly, the blue areas are mostly shown in the two blocks that represent the distances within each cluster. On the other hand, other two blocks (in the left top and right bottom) mostly have red areas, implying far distance between two clusters (between noisy and clean slices).

We compared performance of kernel Wasserstein distance with other methods, including Wasserstein, KL, and kernel KL distances, using Chi-square statistic. Note that for KL and kernel KL we computed the average value of two KL measures, i.e., $J_{\mathrm{KL}}(P\|Q)$. We repeated the analysis process noted above for alternative methods. Of note, kernel methods (kernel Wasserstein distance and kernel KL distance) had the same accuracy with a Chi-square statistic of 436 (Figure 3). By contrast, non-kernel methods (Wasserstein distance and KL distance) had substantially lower Chi-square statistics with 21 and 91, respectively, showing the superiority of kernel methods in this application.

Figure 4 shows scatter plots for the correlation between kernel Wasserstein distance and Wasserstein distance in Cluster 1 (Figure 4(a)), Cluster 2 (Figure 4(b)), and between two clusters (Figure 4(c)). Compared to the correlation in Cluster 1 and between two clusters, the correlation in Cluster 2 that was enriched for noisy slices showed relatively more linear pattern. This may suggest that the impact of kernel Wasserstein distance in Cluster 1 and between two clusters is larger than classical Wasserstein distance, thereby leading to improved performance.

4 Conclusion

The Wasserstein distance is a powerful tool with a wide range of applications. Although extensively used, the method of computing this in the kernel space is lacking. In this paper, we proposed a computational method to solve L2-Wasserstein distance in the kernel space. We applied this method to a medical imaging problem in which CT scans are often degraded by artifacts arising from high-density materials. Our unsupervised method consisting of kernel L2-Wasserstein distance and hierarchical clustering showed a good level of performance in identifying noisy CT slices, outperforming conventional Wasserstein distance. Notably, kernel Kullback-Leibler distance also obtained comparable performance. This implies the nonlinearity of data and thus nonlinear analysis using kernel methods would be more likely to be essential. Future work will focus on further applications of kernel Wasserstein distance in imaging and biological data analysis.

Acknowledgments

This research was funded in part through National Institutes of Health/National Cancer Institute Cancer Center Support grant P30 CA008748.