

In this section we propose novel data driven machine learning (ML) methods for finding \emph{anomalous events} in simulated and real data sets. These methods will enable us to learn what these anomalies could be due to and we will be able to propose follow-up observations or model these spots in such a way that they can be removed from maps that are used in standard 2-point analyses.

We propose methods for two different settings: i) \emph{detection of known structures} and ii) \emph{anomalous event detection}. In the ``detection of known structures'' task, we propose algorithms for the case in which some interesting events (e.g. Lagrangian coherent structures, vortices) have already been identified, and their goal is to find similar ones in the provided datasets.
In the ``anomalous event detection'' problem the methods work in a completely unsupervised way; they require only real or simulated turbulence data and using this data they will try to identify the most anomalous parts automatically without requiring a labeled set of anomalous events. The methods proposed here will be applied to the 3D simulations. Once the anomalous regions have been detected in the 3D simulations, we can analyze the ensuing 2D maps to understand the characteristic signatures that are imprinted in the polarization data and then look for those signatures in the {\it Planck} data.

\emph{Detecting Known Structures in Turbulent Data.} In many scientific problems, scientists have a collected set of labeled examples and their goal is to find similar examples in the rest of the data. For example, when our goal is to find vortices or other known structures in a turbulent dataset, we might already have a labeled collection of vortices and another labeled set of examples not containing vortices. Our goal is to develop a detection algorithm that can learn the similarities and dissimilarities between these classes and to use this acquired knowledge to detect more vortices (or other structures) in the rest of the dataset.

This problem belongs to the family of classification problems. There are numerous algorithms proposed for classification problems and they are very well-studied in ML. One crucial challenging issue, however, is that most ML algorithms developed for similar classification problems can operate only on simple finite dimensional feature vectors, and they cannot be directly applied to the simulations, %or the{\it Planck} data,
where the inputs are not simple finite dimensional vectors. Advances in functional data analysis \cite{Ramsay05FDA}, nonparametric statistics, and ML, however, recently have started to study the question of how to develop learning systems that can process these complex data sets. This is a very rapidly developing field, with many open questions both in theory and practice.

\begin{figure}[h]
	\centering
	\subfigure[ Positive]{
		\includegraphics[height=2.8cm]{fig/turb_quiver_vortex}
		\label{fig:turb-train-vortex} % 25 28 121
	}
	\subfigure[ Negative]{
		\includegraphics[height=2.8cm]{fig/turb_quiver_simpleneg}
		\label{fig:turb-train-simpleneg} % 230 227 166
	}
	\subfigure[ Negative]{
		\includegraphics[height=2.8cm]{fig/turb_quiver_diamondneg}
		\label{fig:turb-train-diamondneg} % 115 146 128
	}
\hspace*{1cm}
 \subfigure[ Detected
 anomaly]{\includegraphics[height=2.85cm]{fig/turb_vort_hd}}
 \vspace*{-1mm}
	\caption{a,b,c: Training examples for the vortex classifier. d: Detected anomaly in the JHU\citep{perlman:07:turb} turbulence data. See Page 10 for further description.}
	\label{fig:turb-examples}
\end{figure}


As a specific example, suppose we are given a turbulence data set and we already identified a few interesting events in this data (e.g., vortices), and our goal is to find similar events. (See Fig.~\ref{fig:turb-examples}). In this case, we assume that these collected events are given as vector field patches, where each vector in the vector field patch has a start and an end position. We can represent any vector field patch as a set of vectors also known as a point cloud. Our goal can be rephrased as a supervised machine learning problem where given a set of positively and negatively labeled point clouds as training data (Fig.~\ref{fig:turb-examples} (a,b,c)), we want to find other events (i.e. other point clouds) in our dataset that are most similar to the positively labeled examples (vortex events) and dissimilar to the negatively labeled examples (non-vortex). Co-I Poczos has recently proposed Reproducing Kernel Hilbert Space  and Deep Neural Network based methods for this problem. We briefly summarize them below.


%\vni
\emph{Reproducing Kernel Hilbert Space Based Estimators.}
Poczos introduced Reproducing Kernel Hilbert Space based methods (Support Vector Machines and Gaussian Processes) to the domain of vector fields, point cloud data, and distributions \cite{poczos12CVPR}. The main idea is that kernel based methods can operate on any input domains as long as we can define and evaluate a positive semi-definite kernel function $k(p,q) $ on this domain. If the input instances are point clouds, then we can also think of each point cloud as a sample set sampled from an unknown distribution with unknown density $p$. Several kernel functions have been proposed that can operate on distributions and point clouds. They include the maximum mean discrepancy (MMD), polynomial kernel ($k(p,q)=(1+\int pq)^{d}$, Gaussian kernel $k(p,q)=c_1\exp(-\int (p-q)^2)/c_2$), and many others.
A crucial challenge is that in practice we do not know the density functions $p$ and $q$ of the point cloud data. All we have is some sample points sampled from $p$ and $q$ and we need to estimate $k(p,q)$. Recently our group has demonstrated that under some conditions $k(p,q)$ can be estimated efficiently from point cloud data. The estimators are completely non-parametric; we do not need to assume that $p$ or $q$ would have a Gaussian or other parametric distribution. One might also replace the $L_2$ distance in the Gaussian kernel to other divergences $D(p,q)$, such as Jensen-Shannon divergence, symmetric Renyi divergence or others. In this case the Gaussian kernel between two distributions becomes $k(p,q)=c_1\exp(-D(p,q)/c_2)$.

Recently we developed efficient nonparametric estimators for these divergences. They can be used to create a positive semi-definite kernel function and plug that kernel function into a support vector machine (or other kernel machines) that can be used for classification and regression.
We call this method support distribution machine (SDM), and have demonstrated in many practical problems (cosmology \cite{Ntampaka15Machine,Ntampaka16Dynamical}, computer vision \cite{poczos12CVPR}, and civil engineering \cite{Jung15Exploration}) that these algorithms can often outperform the state-of-the art. We have also generalized other ML tools (such as Lasso \cite{Oliva14a}, Ridge regression \cite{szabo16JMLR}, and Deep Convolutional Neural networks \cite{Ravanbakhsh16ICML} to distribution, set, function, and other complex inputs, and these algorithms can also be used to process turbulence data.

\emph{Nonparametric Divergence Estimation}. Below we briefly summarize how we can use point cloud data with unknown densities $p$ and $q$ to estimate the kernel value $k(p,q)$. A nice property of these estimators is that in many cases they are simple and can estimate these functionals of the density \emph{without} estimating the density itself.

Let $p$, $q$ be $\R^d \to \R$ density functions,
and let $\alpha \in \R\setminus\{1\}$.  The R\'enyi-$\alpha$ divergence ($D^{R}_\alpha$) is defined as the following functional of the densities $p$ and $q$:
$D^{R}_\alpha(p\|q)\doteq\frac{1}{\alpha-1} \log \int p^{\alpha}(x)q^{1-\alpha}(x) d x$.

Let $X_{1:n}\doteq(X_1,\ldots,X_n)$ be an
i.i.d.\ sample from a distribution with density $p$, and similarly let
$Y_{1:m}\doteq(Y_1,\ldots,Y_m)$ be an i.i.d.\ sample from a
distribution having density $q$.
%Let $X$, $Y$ $\in \M \subseteq \mathbb{R}^d$, and let
Let $\rho_k(i)$ denote the Euclidean distance of the $k^{th}$ nearest
neighbor of $X_i$ in the sample $X_{1:n}$, and similarly let
$\nu_k(i)$ denote the distance of the $k^{th}$ nearest neighbor of
$X_i$ in the sample $Y_{1:m}$. Let $\Gamma$ denote the gamma function.
Recently, Wang et al. \cite{Wang-Kulkarni-Verdu2009} derived an estimator for the
 \acro{KL}-divergence ($D^{KL}(p,q)\doteq\int p\log \frac{p}{q}$), and we proposed a consistent estimator for the R\'enyi-$\alpha$ divergence (P\'oczos et al. \cite{poczos11aistats}):
 \begin{align*}
 \widehat{D}^{KL}(X_{1:n}\|Y_{1:m})&\doteq\frac{d}{n} \sum_{i=1}^n \log \frac{\nu_k(i)}{\rho_k(i)} +\log\frac{m}{n-1}.\\
\widehat D^R_{\alpha}(X_{1:n}\|Y_{1:m})&\doteq \frac{1}{\alpha-1}\log
\frac{1}{n}\sum_{i=1}^n\left(\frac{(n-1)\rho_k^d(i)}{m\nu_k^d(i)}\right)^{1-\alpha}\frac{\Gamma(k)^2}{\Gamma(k-\alpha+1)\Gamma(k+\alpha-1)}.
\end{align*}

In Figure~\ref{fig:rho-nu}, we illustrate how to calculate $\rho_k(i)$ and $\nu_k(i)$ quantities that are needed for these estimators. They can be used to estimate the kernel function $k(p,q)$.
%It is worth noting that divergence estimators can also be used for mutual information estimation, since the mutual information by definition is the divergence between the joint density and the product of the marginal densities. The above estimators are simple; they only use certain $k$th-nearest neighbor distance based statistics that can be efficiently calculated using $k$-d trees \cite{Bentley:1975:MBS}. Recently,
%Sz\'ekely has proposed an estimator that uses all distances between the sample points (\ie the complete Euclidean graph of the sample points) to define a new quantity for measuring the dependence between random variables \cite{szekely07measuring}.

%\vspace*{-.3cm}
\begin{figure}[htb]
  \centering
    \includegraphics[width=6cm]{fig/divergence_demonstration3}
    \vspace*{-1cm}
    \caption{Calculating $\rho_k(i)$ and $\nu_k(i)$. Blue dots with X signs and red dots show samples $X_{1:n}$ and $Y_{1:m}$, respectively. We fixed $k=3$ in this example, i.e. we need to calculate the 3rd nearest neighbors of each $X_i$. }
  \label{fig:rho-nu}
\end{figure}


\emph{Methods for anomalous event detection}. In addition to finding instances of known patterns, our methods can also be generalized to find anomalous, rare events in a completely unsupervised way. One possible approach is to flag those point cloud regions anomalous whose distributions---using the above listed MMD, KL, Renyi-$\alpha$, or another divergence---is far from most of the distributions of the other regions of the point cloud data. Another possible approach is to use the so-called One-class SVM method \cite{NIPS1999_1723} with an appropriate kernel that can operate on point cloud data. The One-class SVM method is an unsupervised method for detecting anomalous points in datasets.\\
% part of the process of exploring the results of a large-scale simulation is looking for \emph{unexpected rare phenomena} that scientists can analyze further.\\


%{\emph Deep Neural Network Based Methods} ****\\

\emph{Preliminary demonstrations on detection of known structures.}
To show a preliminary result on the applicability of these methods to detect vortices in turbulence datasets, we performed the following experiment.
We trained an SDM classifier using turbulence data from the JHU Turbulence Data Cluster\footnote{\url{http://turbulence.pha.jhu.edu}} \citep{perlman:07:turb},
(TDC). TDC simulates fluid flow through time on a 3-dimensional grid, calculating 3-dimensional velocities and pressures of the fluid at each step.
We trained an SDM on a manually labeled training set of 11 positive (containing vortex) and 20 negative examples (not containing vortices).
Training examples are shown in Fig.~\ref{fig:turb-examples}(a--c).
This classifier was then used to evaluate groups along $z$-slices of the data. One slice of the resulting probability
estimates is shown in Fig.~\ref{fig:turb-scores}; the arrows represent the
mean velocity at each classification point. The high-probability region on the
left is a canonical vortex, while the slightly-lower probability region in the
upper-right deviates a little from the canonical form. Note that some other
areas show somewhat complex velocity patterns but mostly have low
probabilities. As we can see in the figure, our method can accurately detect vortices. Leave-one-out cross-validation yielded an
accuracy of 97\% (one mistake).

In this proposal we will extend these ideas and show that similar methods can also be used to detect other Lagrangian Coherent Structures in turbulent datasets. Using recent advances in ML on big data sets, we will also study how to scale these novel classification algorithms up to large datasets. In related function-to-function regression and distribution regression problems we were able to scale up regression methods to big data problems \cite{oliva2015fastf2f,Oliva2014fast}.\\



\emph{Preliminary results on unsupervised anomalous event detection}.
In the above described TDC simulation we consider the vertices in a
local cubic region as a group, and the goal is to find groups of vertices
whose velocity distributions (i.e.~moving patterns) are unusual and
potentially interesting. To demonstrate anomaly detection with SDMs, we trained a ``one-class SDM'' on 100
typical distributions from the turbulence data, with centers chosen uniformly at
random. This one-class SDM method then was asked to process the data and find regions that are most untypical and look the most different from the examples in the training set. As we can see in Fig.~\ref{fig:turb-anomaly}, the two vortices in the data set are picked out, but
the area with the highest score is a diamond-like velocity pattern; these may well be less
common and therefore more anomalous in the dataset than vortices. In Fig.~\ref{fig:turb-examples}(d) we show another detected anomaly in the same TDC dataset. As we can see, this is indeed a rare and interesting interaction among three vortices positioned in a straight line.

In this proposal we will extend these ideas to be able to process large datasets and find different types of anomalies in the turbulence datasets. We believe that SDM classifiers and anomaly detectors can serve as simulation exploration tools that will allow scientists to
iteratively look for anomalous phenomena and label some of them. Classifiers
could then find more instances of those phenomena and compute statistics about
their occurrence, while anomaly detection would be iteratively refined to
highlight only what is truly new.
%We will use these methods to find anomalies that can lead to the non-Kolmogorov contribution in atmospheric turbulence


\begin{figure*}
	\centering
	\subfigure[Classification probabilities]{
%\includegraphics[angle=0,width=3.1in]{fig/turb_scores_quiver}
\includegraphics[angle=90,width=3in]{fig/turb_scores_quiver}
		\label{fig:turb-scores}
	}
	\subfigure[Anomaly scores]{
		% \includegraphics[angle=0,width=3.1in]{fig/turb_anomaly_quiver} 
		 \includegraphics[angle=90,width=3in]{fig/turb_anomaly_quiver}
		\label{fig:turb-anomaly}
	}
 \vspace*{-1mm}
	\caption{Classification and anomaly scores along with velocities for one $54\times48$ slice of the turbulence data.}
 	\label{fig:turb-score-quivers}
\end{figure*}





