Differential Similarity in Higher Dimensional Spaces: Theory and Applications

TL;DR

This paper extends the theory of differential similarity to higher dimensions, developing algorithms for clustering and coding that integrate geometric and probabilistic models, with applications to datasets like MNIST and CIFAR-10.

Contribution

It generalizes the differential similarity framework from 3D to n-dimensional spaces and introduces new solution and parameter estimation strategies.

Findings

Effective algorithms for n-dimensional differential similarity

Successful application to MNIST and CIFAR-10 datasets

Improved clustering and coding performance

Abstract

This paper presents an extension and an elaboration of the theory of differential similarity, which was originally proposed in arXiv:1401.2411 [cs.LG]. The goal is to develop an algorithm for clustering and coding that combines a geometric model with a probabilistic model in a principled way. For simplicity, the geometric model in the earlier paper was restricted to the three-dimensional case. The present paper removes this restriction, and considers the full $n$-dimensional case. Although the mathematical model is the same, the strategies for computing solutions in the $n$-dimensional case are different, and one of the main purposes of this paper is to develop and analyze these strategies. Another main purpose is to devise techniques for estimating the parameters of the model from sample data, again in $n$ dimensions. We evaluate the solution strategies and the estimation techniques by applying them to two familiar real-world examples: the classical MNIST dataset and the CIFAR-10 dataset.

Figures

Tables

Table 1. Table 1. Radius in the original Euclidean space of the 32,000 point Data Sphere and the 8,000 point Coordinate Sphere, for three prototypical clusters.

	(a)	(b)
	Radius of	Radius of
	Data Sphere	Coordinate Sphere
Prototype 09	2.95906	1.75322
Prototype 24	3.15226	2.30487
Prototype 27	2.97573	1.95364

Table 2. Table 2. Distances along the ρ 𝜌 \rho coordinate curves for the principal axes from a point on the Coordinate Sphere, both inwards and outwards. Compare Figures 13 and 14 .

	(c)	(d)	(e)
	Euclidean	Riemannian	Riemannian
	Distance	Distance	Distance
	Inwards	Inwards	Outwards
Prototype 09	2.40866	0.90783	4.01243
Prototype 24	2.90017	0.764126	6.93029
Prototype 27	2.0754	0.38668	3.41653

Table 3. Table 3. Riemannian and Euclidean distances along the Θ Θ \Theta coordinate curves for the maximal infinitesimal eigenvectors, from the principal axis to the Euclidean angles π / 2 𝜋 2 \pi/2 and − π / 2 𝜋 2 -\pi/2 . Compare Figure 15 .

	Positive	Negative
	Direction	Direction
Prototype 09	2.22261	4.16981
Prototype 24	2.71304	3.58147
Prototype 27	2.10106	2.46165

Table 4. Table 4. Prototype 09: Riemannian and Euclidean distances along the Θ Θ \Theta coordinate curves to π / 2 𝜋 2 \pi/2 and − π / 2 𝜋 2 -\pi/2 for three minimal infinitesimal eigenvectors. Compare Figure 16 .

	Positive	Negative	Total
	Direction	Direction	Distance
Curve 16	2.61453	3.16092	5.77545
Curve 17	2.68326	3.17738	5.86064
Curve 09	2.68949	3.11284	5.80233

Table 5. Table 5. Prototype 09: The ten Θ Θ \Theta coordinate curves with the least total Riemannian and Euclidean distances from the principal axis to π / 2 𝜋 2 \pi/2 and − π / 2 𝜋 2 -\pi/2 .

	Positive	Negative	Total
	Direction	Direction	Distance
Curve 42	2.80639	2.77794	5.58433
Curve 23	2.77045	2.81585	5.58631
Curve 30	2.79071	2.79673	5.58744
Curve 37	2.77144	2.81653	5.58799
Curve 29	2.8927	2.69908	5.59179
Curve 44	2.83735	2.75788	5.59523
Curve 11	2.70893	2.88847	5.59741
Curve 36	2.80777	2.79319	5.60096
Curve 26	2.71708	2.88505	5.60213
Curve 01	2.69536	2.90715	5.60251

Table 6. Table 6. Radius, in 3X49D and 49D Euclidean space, of the 8,000 point Data Sphere and the 2,000 point Coordinate Sphere, for seven prototypical clusters.

	Radius of		Radius of
	Data Sphere		Coordinate Sphere
	3X49D	49D	3X49D	49D
Prototype 14	2.11000	1.21821	1.42727	0.824038
Prototype 20	1.91678	1.10665	1.31318	0.758166
Prototype 23	1.97005	1.13741	1.34494	0.776500
Prototype 32	1.78403	1.03001	1.31377	0.758503
Prototype 34	2.06485	1.19214	1.47952	0.854199
Prototype 38	1.91457	1.10538	1.42930	0.825207
Prototype 40	1.76567	1.01941	1.19463	0.689722

Table 7. Table 7. Number of data points in the Θ Θ \Theta -Coordinate Clusters from Prototype 14 (see Figure 26 ) and Prototype 32 (see Figure 27 ).

	$\Theta$-Coordinate Cluster
	1	2	3	4	5	6	Total
Prototype 14	149	180	127	108	121	214	899
Prototype 32	217	99	135	172	223	219	1065

Table 8. Table 8. Riemannian distance ratios along the surface of the ellipsoidal approximations for the 6 Θ Θ \Theta -Coordinate Clusters from Prototype 14 and Prototype 32.

				$\Theta$-Coordinate Clusters
	N		Optimal	Within	Across
		R	0.794145	0.396452	0.312464
Prototype 14	899	G	0.855848	0.436969	0.333822
		B	0.794182	0.385526	0.301911
		R	1.171935	0.436817	0.380668
Prototype 32	1065	G	0.825405	0.345734	0.298927
		B	1.014526	0.393368	0.347713

Table 9. Table 9. Riemannian distance ratios for the optimal 10-dimensional coordinates in the Θ Θ \Theta -Coordinate Clusters from Prototype 14 and Prototype 32.

	$\Theta$-Coordinate Cluster
	1	2	3	4	5	6
Prototype 14	1.26189	1.14097	1.32156	1.31458	1.26402	1.23068
Prototype 32	1.33712	1.13318	1.04891	2.15558	1.17315	1.32071

Equations

$$\mathcal{L}$$

$$\frac{\partial }{\partial t}w(t,\mathbf{x})$$

$$=\mathcal{L}w(t,\mathbf{x}),\text{mbox}\mathrm{withinitialcondition}w(0,\mathbf{x})=f(\mathbf{x}),$$

$$\frac{\partial }{\partial t}p(t,\mathbf{x})$$

$$={\mathcal{L}}^{\ast }p(t,\mathbf{x}),$$

$$\frac{1}{2}\sum _{i,j}{a}^{ij}(\mathbf{x})\frac{{\partial }^2}{\partial x^i\partial x^j}p(\mathbf{x})-\sum _i\left(\frac{\partial }{\partial x^i}{b}^i(\mathbf{x})\right)p(\mathbf{x})-\sum _i{b}^i(\mathbf{x})\left(\frac{\partial }{\partial x^i}p(\mathbf{x})\right)$$

$$u=u(x,y,z),v=v(x,y,z),w=w(x,y,z),$$

$$\frac{\partial w}{\partial t}=\frac{1}{2}\Delta w+\nabla U(\mathbf{x})\cdot \nabla w,\mathrm{with}w(0,\cdot )=f.$$

$$\frac{\partial u}{\partial t}=\frac{1}{2}\Delta u-V(\mathbf{x})u,\mathrm{with}u(0,\cdot )=g,$$

$$u(t,\mathbf{x})={\int }_{\Omega }g(X_t)\exp \left[-{\int }_0^{t}V(X_s)ds\right]{\mathcal{W}}_{\mathbf{x}}(dX),$$

$$\frac{\hslash }{i}\frac{\partial \psi (t,\mathbf{x})}{\partial t}=\frac{{\hslash }^2}{2m}\Delta \psi (t,\mathbf{x})-V(\mathbf{x})\psi (t,\mathbf{x}),$$

$$X(t)=X(0)+{\int }_0^t\sigma (s,\omega )d\mathcal{B}(s,\omega )+{\int }_0^{t}b(s,\omega )ds,$$

$$dX(t)=\sigma (t,\omega )d\mathcal{B}(t,\omega )+b(t,\omega )dt.$$

$$dX(t)=\left(\begin{array}{ccc}{\sigma }_1^1& \dots & {\sigma }_d^1\\ ⋮& & ⋮\\ {\sigma }_1^n& \dots & {\sigma }_d^n\end{array}\right)\left(\begin{array}{c}d{\mathcal{B}}_1(t)\\ ⋮\\ d{\mathcal{B}}_d(t)\end{array}\right)+\left(\begin{array}{c}{b}^1\\ ⋮\\ b^n\end{array}\right)dt$$

$$dX(t)=\left(\begin{array}{ccc}{\sigma }_1^1& \dots & {\sigma }_d^1\\ ⋮& & ⋮\\ {\sigma }_1^n& \dots & {\sigma }_d^n\end{array}\right)\circ \left(\begin{array}{c}d{\mathcal{B}}_1(t)\\ ⋮\\ d{\mathcal{B}}_d(t)\end{array}\right)+\left(\begin{array}{c}{\tilde{b}}^1\\ ⋮\\ {\tilde{b}}^n\end{array}\right)dt$$

$${\tilde{b}}^i=b^i-\frac{1}{2}\sum _{k=1}^d\sum _{j=1}^n\frac{\partial {\sigma }_k^i}{\partial x^j}{\sigma }_k^j.$$

$${\rho }^{\prime }(t)$$

$$\rho (0)$$

$$\rho (t)={\mathbf{x}}_0+{\int }_0^t\frac{\nabla U(\rho (s))}{|\nabla U(\rho (s))|}ds,0\le t$$

$$\nabla U(\mathbf{x})=(P_0(\mathbf{x}),P_1(\mathbf{x}),\dots ,P_{n-1}(\mathbf{x})),$$

$$\nabla U(\mathbf{x})$$

$${\mathbf{V}}_1(\mathbf{x})$$

$${\mathbf{V}}_2(\mathbf{x})$$

$$\dots$$

$${\mathbf{V}}_{n-2}(\mathbf{x})$$

$${\mathbf{V}}_{n-1}(\mathbf{x})$$

$$\left[{\mathbf{V}}_i\partial ,{\mathbf{V}}_j\partial \right]$$

$$\displaystyle=\;\frac{1}{P_{0}^{2}}\left(\begin{array}[]{c}P_{0}\left[\frac{\partial P_{i}}{\partial x^{j+1}}-\frac{\partial P_{j}}{\partial x^{i+1}}\right]\\[8.61108pt] +\;P_{i}\left[\frac{\partial P_{j}}{\partial x^{1}}-\frac{\partial P_{0}}{\partial x^{j+1}}\right]\\[8.61108pt] +\;P_{j}\left[\frac{\partial P_{0}}{\partial x^{i+1}}-\frac{\partial P_{i}}{\partial x^{1}}\right]\end{array}\right)\frac{\partial}{\partial x^{1}}$$

$$\left[f\mathbf{V}\partial ,g\mathbf{W}\partial \right]$$

$$=f(\mathbf{V}\partial g)\mathbf{W}\partial +fg\mathbf{V}\partial \circ \mathbf{W}\partial$$

$$-g(\mathbf{W}\partial f)\mathbf{V}\partial -gf\mathbf{W}\partial \circ \mathbf{V}\partial$$

$$=fg\left[\mathbf{V}\partial ,\mathbf{W}\partial \right]+f(\mathbf{V}\partial g)\mathbf{W}\partial -g(\mathbf{W}\partial f)\mathbf{V}\partial$$

$$\mathbf{X}\partial$$

$$\mathbf{Y}\partial$$

$$\left(\begin{array}{c}\\ g_{i,j}(\mathbf{x})\\ \end{array}\right)=$$

$$\displaystyle\left(\begin{array}[]{cccccc}|\nabla U|^{2}&0&0&\ldots&0&0\\[4.30554pt] 0&P_{0}^{2}+P_{1}^{2}&P_{2}P_{1}&\ldots&P_{n-2}P_{1}&P_{n-1}P_{1}\\[4.30554pt] 0&P_{1}P_{2}&P_{0}^{2}+P_{2}^{2}&\ldots&P_{n-2}P_{2}&P_{n-1}P_{2}\\[4.30554pt] \ldots&&&&&\ldots\\[4.30554pt] 0&P_{1}P_{n-2}&P_{2}P_{n-2}&\ldots&P_{0}^{2}+P_{n-2}^{2}&P_{n-1}P_{n-2}\\[4.30554pt] 0&P_{1}P_{n-1}&P_{2}P_{n-1}&\ldots&P_{n-2}P_{n-1}&P_{0}^{2}+P_{n-1}^{2}\\[4.30554pt] \end{array}\right)$$

$$\sum _{i=1}^{n-1}{v}^i(t){\mathbf{V}}_i(\mathbf{x})$$

Full text

Differential Similarity in

Higher Dimensional Spaces:

Theory and Applications

L. Thorne McCarty

Department of Computer Science

Rutgers University

(Date: September, 2021.  ©  L. Thorne McCarty)

\markleft

L. Thorne McCarty

ABSTRACT: This paper presents an extension and an elaboration of the theory of differential similarity, which was originally proposed in [McC14, McC18]. The goal is to develop an algorithm for clustering and coding that combines a geometric model with a probabilistic model in a principled way. For simplicity, the geometric model in the earlier paper was restricted to the three-dimensional case. The present paper removes this restriction, and considers the full $n$-dimensional case. Although the mathematical model is the same, the strategies for computing solutions in the $n$-dimensional case are different, and one of the main purposes of this paper is to develop and analyze these strategies. Another main purpose is to devise techniques for estimating the parameters of the model from sample data, again in $n$ dimensions. We evaluate the solution strategies and the estimation techniques by applying them to two familiar real-world examples: the classical MNIST dataset and the CIFAR-10 dataset.

1. Introduction.

This paper presents an extension and an elaboration of the theory of differential similarity, which was originally proposed in [McC14, McC18]. The goal is to develop an algorithm for clustering and coding that combines a geometric model with a probabilistic model in a principled way. The geometric model is a Riemannian manifold with a Riemannian metric, ${g}_{ij}({\bf x})$, which is interpreted as a measure of dissimilarity. The probabilistic model consists of a stochastic process with an invariant probability measure that matches the density of the sample input data. The link between the two models is a potential function, $U({\bf x})$, and its gradient, $\nabla U({\bf x})$. Since the components of the gradient appear in the definition of the dissimilarity metric, the measure of dissimilarity will depend on the probability measure. Roughly speaking, the dissimilarity will be small in a region in which the probability density is high, and vice versa. Finally, the dissimilarity metric is used to define a coordinate system on the embedded Riemannian manifold, which leads to an “optimal” lower dimensional encoding of the original data.

For simplicity, the geometric model in [McC14, McC18] was restricted to the three-dimensional case, and this was the main deficiency of the theory. The present paper removes this restriction, and considers the full $n$-dimensional case. Although the mathematical model is the same, the strategies for computing solutions in the $n$-dimensional case are different, and one of the main purposes of this paper is to develop and analyze these strategies. Another main purpose is to devise techniques for estimating the parameters of the model from sample data, again in $n$ dimensions. We evaluate the solution strategies and the estimation techniques by applying them to two familiar real-world examples: the classical MNIST dataset [LBBH98] and the CIFAR-10 dataset [Kri09].

In the broadest terms, this work is an exploration of the manifold hypothesis in deep learning [RDV*+*12] [BCV13] [FMN16]. It is a common observation that real-world data in high-dimensional spaces tends to be concentrated on low-dimensional nonlinear subspaces, and this phenomenon seems to contribute to the success of deep neural networks in image recognition, speech recognition, and other similar tasks.

There are now quite a few algorithms for manifold learning: [SSM98] [TB99] [RS00] [TSL00] [BN03] [Bra03] [DG03] [HR03] [ZZ04] [CL06] [WS06] [CSP*+*10] [CZF10] [YZG10]. These algorithms fall into several categories. Some use global methods: [TSL00] [WS06], while others are primarily local. Some use probabilistic models: [TB99] [HR03] [CSP*+*10], while others are primarily geometric. Among the local geometric algorithms, some are based on the Laplacian [BN03] or the Hessian [DG03], and some are based explicitly on a diffusion process [CL06]. For an analysis of this latter category, see [LW10]. Often, a discrete stochastic process is defined initially on a finite graph (e.g., as a random walk) and the limiting case is shown to be a diffusion on a manifold. See, e.g., [BN05] [HAvL07] [THJ10]. In such algorithms, nonlinear dimensionality reduction is usually achieved by a spectral decomposition of the Laplacian on the graph.

In a recent preprint on nonlinear dimensionality reduction [TJ18], Ting and Jordan develop a general theory for the class of local spectral methods: “These methods construct a matrix using only information in local neighborhoods and take a spectral decomposition to find a nonlinear embedding.” The class includes: [RS00] [BN03] [DG03] [ZZ04] [CL06]. Ting and Jordan’s general framework specifies a differential operator on a compact manifold, with a variety of boundary conditions, and they then analyze how each method converges from a set of conditions on a local neighborhood graph to an eigenproblem for the differential operator. Distinct methods correspond to distinct operators and distinct boundary conditions. The theory also leads to a ranking of the various algorithms, and the authors conclude that Local Tangent Space Alignment (LTSA) by Zhang and Zha [ZZ04] is the best.

However, if we are investigating the manifold hypothesis in the context of deep learning, Bengio and his colleagues argue in [BCV13] that we need a method with very different properties. One problem is the neighborhood graph, which has quadratic complexity. More significantly, any manifold learning algorithm based solely on local neighborhoods is not likely to generalize very well, in a deep network, beyond the intial training data. A better algorithm for deep learning would construct a parametric coordinate mapping that takes into account remote data, as well as local neighborhood data. Among the existing algorithms, [BCV13] singles out Local Coordinate Coding (LCC) by Yu, Zhang and Gong [YZG10], which has some similarities to sparse coding [OF96]. But most of the discussion in [BCV13] focuses on network architectures that can learn embedded manifolds directly from the structure of the data density itself, known as Regularized Auto-Encoders. Two types are considered: Denoising Auto-Encoders (DAEs) [Vin11] and a specialized form of Contractive Auto-Encoders (CAEs) [AB14], both of which can be shown to compute the gradient of the log of the input probability density.

The theory of differential similarity [McC14, McC18] matches the desiderata advocated by [BCV13] more closely than do the algorithms analyzed in [TJ18]. First, it is not based on neighborhood graphs: It is defined from the start on Euclidean ${\bf R}^{n}$. Second, the lower dimensional encodings in the theory are constrained globally as well as locally: In the geometric model, the relationship between the local tangent bundle and the global integral manifold is strictly determined by a classical theorem in differential geometry. In the probabilistic model, the relationship between the local diffusion process and the global probability density is strictly determined by a classical theorem on stochastic processes. The diffusion equation in the theory of differential similarity has a drift term as well as a Laplacian term, which means that the diffusion has an invariant probability measure, or a stationary probability density, unlike the diffusion in [LW10]. And the gradient of the log of the stationary probability density, which is computed by a DAE [Vin11] or a CAE [AB14] in the deep learning framework, is precisely the vector field, $\nabla U({\bf x})$, in the theory of differential similarity. We will return to this point in our discussion of “Future Work” in Section 8.

The balance of the paper is organized as follows: Section 2 is a review of “The Probabilistic Model” and Section 3 is a review of “The Geometric Model” from [McC14, McC18]. If the reader is willing to accept on authority (with citations) a few basic results on stochastic processes and differential geometry, this material should be accessible to anyone with a knowledge of linear algebra and advanced calculus. Section 2 is short, and it includes a simple proof of the basic theorem that we will need, without the extended discussion of stochastic processes in [McC14, McC18]. Section 3 is longer, because it is necessary to extend the geometric model to $n$ dimensions. We use a form of prototype coding for the coordinate system, measuring the distance from the origin (i.e., the “prototype”) in $n-1$ specified directions. Thus we define a radial coordinate, $\rho$, and the directional coordinates $\theta^{1},\theta^{2},\ldots,\theta^{n-1}$, collectively denoted by $\Theta$. We refer to this as the $\rho$, $\Theta$, coordinate system.

The paper then turns from theory to applications: “How to Estimate $\nabla U({\bf x})$ from Sample Data” in Section 4 and “Computing the Geodesic Coordinate Curves” in Section 5. To work with sample data, we borrow a technique from the literature on the mean shift algorithm [FH75] [Che95] [CM02]. Once we have an estimate of $\nabla U({\bf x})$, everything else in the model can be calculated from its components. In particular, since the coordinate curves are defined by geodesics on the embedded Riemannian manifold, they are determined by the Euler-Lagrange equations for the minimization of the energy functional over the Riemannian metric, ${g}_{ij}({\bf x})$. This is a large system of differential and algebraic equations, in a high-dimensional space, but it can be solved numerically in Mathematica.

Finally, Sections 6 and 7 demonstrate how the theory works when applied to real-world examples. Section 6, on the MNIST dataset [LBBH98], shows that our calculations lead to intuitively reasonable results, and Section 7, on the CIFAR-10 dataset [Kri09], adds another wrinkle: We show how to use quotient manifolds to build invariance into the geometric model, and we show how to use product manifolds to combine low-dimensional solutions into a higher dimensional problem space, so that our dimensionality reduction techniques can be applied recursively.

2. The Probabilistic Model.

The probabilistic model is known in the literature as Brownian motion with a drift term. More precisely, it is a diffusion process generated by the following differential operator:

[TABLE]

where $\Delta$ is the standard Laplacian expressed in Cartesian coordinates and $U(\mathbf{x})$ is a scalar potential function. Brownian motion, by itself, is generated by the differential operator $\frac{1}{2}\Delta$. But Brownian motion “dissipates,” that is, it has no invariant probability measure except zero. When we add a drift term, which is given here by $\nabla U(\mathbf{x})\cdot\nabla$, the invariant probability measure turns out to be finite and proportional to $e^{\,2\,U(\mathbf{x})}$. This means that $\nabla U(\mathbf{x})$ is proportional to the gradient of the log of the stationary probability density.

There are several ways to analyze this diffusion process, and establish this result. One classical approach is to use Kolmogorov’s backward and forward equations. See [Kol31]. Kolmogorov’s backward equation is:

[TABLE]

in which $\mathbf{a}(\mathbf{x})$ is a matrix of diffusion coefficients and $\mathbf{b}(\mathbf{x})$ is a vector of drift coefficients. Notice that the operator $\mathcal{L}$ in (1) is a special case of the operator $\mathcal{L}$ in (2). Kolmogorov’s forward equation is:

[TABLE]

in which $\mathcal{L}^{\mathbf{*}}$ is the formal adjoint of $\mathcal{L}$, and $p(t,\mathbf{x})$ is a probability density equal, in the limit, as $t\rightarrow 0$, to the unit probability mass at x. To find the stationary probability density of our diffusion process, we need to specialize the operator $\mathcal{L}^{\mathbf{*}}$ in (3) to the formal adjoint of the operator $\mathcal{L}$ in (1), and then set ${\partial p}/{\partial t}=\mathcal{L}^{\mathbf{*}}p=0$. But if $\mathbf{a}(\mathbf{x})$ is a constant matrix, then the right-hand side of (3) can be expanded and simplified to

[TABLE]

If we now set $\mathbf{a}(\mathbf{x})$ equal to the identity matrix and $\mathbf{b}(\mathbf{x})=\nabla U(\mathbf{x})$, it is a straightforward computation to verify that $\mathcal{L}^{\mathbf{*}}(e^{\,2\,U(\mathbf{x})})=0$.

Figure 1 is a three-dimensional example borrowed from [McC18]. The potential function $U({\bf x})$ is defined here as a quadratic polynomial over the variables $u$, $v$, and $w$, where

[TABLE]

is defined as a cubic polynomial coordinate transformation from $(x,y,z)$ to $(u,v,w)$. Thus $U({\bf x})$ is a sixth-degree polynomial in $x$, $y$, and $z$, and the gradient, $\nabla U({\bf x})$, is a fifth-degree polynomial. We call this the curvilinear Gaussian potential. One important property of this potential function is the fact that $\nabla U(x,y,z)=(0,0,0)$ at the origin, which means that $U(0,0,0)=0$ is an extremal point, a maximal point, in fact. Figure 1 is a contour plot for the surface at $U(x,y,z)=-10$. Figure 2(a) shows a StreamPlot of the gradient vector field generated by $\nabla U(x,y,z)$ at $z=-10$, and Figure 2(b) shows a stack of such stream plots, at the values $z=10$, $z=0$ and $z=-10$. Notice how the drift vector twists and turns to counteract the dissipative effects of the diffusion term, and maintain an invariant probability measure.

Looking at Figure 2, an interesting idea comes to mind: Could we use this gradient vector field to define a three-dimensional, $\rho$, $\Theta$, coordinate system? The radial coordinate, $\rho$, would follow the gradient vector, $\nabla U$, and the directional coordinates, $\theta^{1}$, $\theta^{2}$, would be orthogonal to $\rho$. We will see how to do this in Section 3.

The simple formula in (1) is all the reader needs to know about stochastic processes in order to understand the rest of this paper, including the examples in Sections 6 and 7, infra. However, Equation (1) is part of a much broader and deeper mathematical subject, discussed in Section 2 of [McC18]. Here is a summary, which could be skipped on a first reading:

Section 2.1 of [McC18] discusses the connection between Equation (1) and the Feynman-Kac formula [Fey48] [Kac49].

Specializing the operator $\mathcal{L}$ in (2) to the operator $\mathcal{L}$ in (1), we have:

[TABLE]

Now consider the following partial differential equation:

[TABLE]

where $V({\bf x})\;=\;\frac{1}{2}\left(\Delta U({\bf x})+|\nabla U({\bf x})|^{2}\right)$. We can show, by a straightforward calculation, using the definition of $V({\bf x})$ in terms of $U({\bf x})$, that $w(t,{\bf x})$ is a solution to (4) if and only if $e^{U(\bf x)}w(t,{\bf x})$ is a solution to (5) with initial value $u(0,\cdot)=e^{U}f$. See Lemma 1 in [McC18]. The Feynman-Kac formula associated with (5) is:

[TABLE]

where $X_{t}\equiv X(t,\omega)$ denotes a continuous path in ${\bf R}^{n}$, and $\mathcal{W}_{\bf x}$ denotes Wiener measure over all such paths beginning at $X_{0}={\bf x}$. Assuming mild regularity conditions, the theorem of Kac [Kac49], inspired by Feynman [Fey48], states that $u(t,{\bf x})$ as defined by (6) is a solution to (5). Furthermore, because of the relationship between (4) and (5), there is a similar integral, involving both $U({\bf x})$ and $V({\bf x})$, that provides a solution to (4). See Theorem 1 in [McC18].

These equations all have their origins in physics. For example, Equation (5) is a real-valued version of the Schrödinger equation:

[TABLE]

and Equation (6) with an $i$ in the exponent is Feynman’s famous “path integral” interpretation of quantum mechanics. (A more familiar version of Schrödinger’s equation can be obtained by multiplying both sides of Equation (7) by $i^{2}=-1$.) The Kolmogorov forward equation, Equation (3), is known to physicists as the Fokker-Planck equation. Also well known, to physicists, is Chapter 10 of [FH65], in which Feynman and Hibbs analyze a representation of the statistical density matrix in quantum statistical mechanics by means of a real-valued path integral in the form of Equation (6). The existence of these mathematical models in physics leads to a speculative conjecture: Could there be a physical device, at the molecular level, perhaps, that could compute analog solutions for various quantities associated with Equation (1)?

Section 2.2 of [McC18] discusses the interpretation of Equations (1) and (2) as stochastic differential equations, following the theories of both Itô [Ito51] and Stratonovich [Str66].

Itô’s theory starts with the definition of a stochastic integral in the following form:

[TABLE]

where the first integral is an *Itô integral *defined with respect to the Brownian motion $\mathcal{B}(t,\omega)$. In differential notation, this would be:

[TABLE]

Extending the notation to $n$ dimensions, let $\mathcal{B}_{1}(t,\omega),\dots,\mathcal{B}_{d}(t,\omega)$ be $d$ independent Brownian motion processes, and define the $n$-dimensional Itô process as follows:

[TABLE]

In this equation, ${\mathbf{\sigma}}\colon{\bf R}^{n}\to{\bf R}^{n\times d}$ is the “square root” of a, that is, ${\bf a}(\mathbf{x})={\mathbf{\sigma}}(\mathbf{x})\mathbf{\sigma}(\mathbf{x})^{T}$. One basic result of Itô’s theory is that Equation (8) defines the same stochastic process as the operator $\mathcal{L}$ in Equation (2). See Theorem 2 in [McC18].

For our purposes, however, the Itô process has a defect: It is not invariant under coordinate transformations. An alternative is to use the stochastic integral proposed by Stratonovich [Str66]. (Technically, in the discretization of $t$ that leads to the definition of the integral for $d\mathcal{B}(t,\omega)$, Itô’s theory evaluates the integrand at the initial point of the interval $[t_{j},t_{j+1}]$, while Stratonovich’s theory evaluates it at the mid point.) A common notation for this alternative is:

[TABLE]

It turns out that the Stratonovich integral satisfies a formula for the “chain rule” that is consistent with the Newton-Leibniz calculus, and thus Equation (9) behaves properly under coordinate transformations. Fortunately, we can use both theories, and translate back and forth between the two of them, because (8) and (9) define the same stochastic process whenever

[TABLE]

See Lemma 2 in [McC18]. This translation can therefore be used to rewrite in a nonlinear $\rho$, $\Theta$, coordinate system any stochastic process that was initially defined by Equation (1) in Euclidean ${\bf R}^{n}$.

In particular, Section 6 in [McC18] shows how to convert the example in Figure 1 from an Itô equation in Euclidean ${\bf R}^{3}$ into a Stratonovich equation in the coordinates ($\rho,\Theta$), and then back into an Itô equation with coefficients $\alpha^{ij}(\rho,\Theta)$ and $\beta^{i}(\rho,\Theta)$. One interesting consequence of these conversions is the calculation of the “drift correction vector field” illustrated in Figure 18 in Section 6 of [McC18].

Section 2.3 of [McC18] discusses integral curves and martingales on manifolds, and develops another interpretation of Equation (1) based on Stroock’s Theorem 7.3.10 in [Str93]. See Theorem 4 in [McC18]. This interpretation is not actually used in [McC18] to justify additional calculations. However, it is likely that Stroock’s work will be useful if we want to advance our theoretical understanding of how the stochastic process generated by Equation (1) interacts with the geometric model that we will construct in Section 3. Our coordinate system for the geometric model is based on integral curves, as we will see, and thus the papers of Stroock and Taniguchi [ST94][ST96] are highly relevant.

3. The Geometric Model.

To implement the idea of prototype coding in our geometric model, we need to define a radial coordinate, $\rho$, and the directional coordinates, $\theta^{1},\theta^{2},\ldots,\theta^{n-1}$, where $n$ is the dimensionality of the initial Euclidean space. But what we really want is a lower dimensional subspace, a $k$-dimensional subspace, say, where $k<n$. Somehow, we would like to choose $k-1$ out of the $n-1$ directional coordinates, and project our diffusion process onto the resulting $k-1$ dimensional space, which can then be combined with our one-dimensional radial coordinate to give us a $k$-dimensional subspace. How should these coordinate systems be defined?

First, we want the radial coordinate, $\rho$, to follow the drift vector, $\nabla U({\bf x})$. To do this, we define $\rho(t)$ to be the integral curve of the vector field $\nabla U({\bf x})$, starting at some initial point ${\bf x}_{0}$. More specifically, we define $\rho(t)$ to be the solution to the following differential equation:

[TABLE]

or, equivalently, the solution to the following integral equation:

[TABLE]

Since the vector field in (10) or (11) is normalized, the integral curve that solves these equations will be parameterized by Euclidean arc length. However, the parameterization that we choose is just a matter of convenience, and what we really want is a generalization of the concept of an integral curve, known as an integral manifold. A one-dimensional integral manifold is, roughly speaking, just the image of an integral curve without the parameterization, and it always exists, for any vector field.

For the directional coordinates, $\theta^{1},\theta^{2},\ldots,\theta^{n-1}$, the obvious generalization would be an integral manifold of dimension $n-1$, orthogonal to the integral manifold for $\rho$. But, for $k\geq 2$, a $k$-dimensional integral manifold exists if and only if certain conditions are satisfied, known as the Frobenius integrability conditions. Fortunately, as we will see, if we are looking for an integral manifold orthogonal to a vector field that is proportional to the gradient of a potential function, such as $\nabla U({\bf x})$, then the Theorem of Frobenius gives us the results that we want. Our analysis here is based on the standard literature in differential geometry. See, e.g., [Spi99], Chapter 6; [BG68], Chapter 3; [AM77], Chapter 8; [BC01], Chapter 1; [Lan95], Chapter VI.

Let’s consider an $n-1$ dimensional tangent subbundle, $E$, in ${\bf R}^{n}$ at some point ${\bf x}$ along the integral curve $\rho(t)$. We will initially use the Cartesian coordinates from the ambient space ${\bf R}^{n}$ to define a set of basis vectors for $E$, which suggests that one axis should be used to “center” the coordinate system and the other $n-1$ axes should be used to specify alternative directions in the vector space. To simplify both the exposition and our later calculations, we will always “center” our coordinate system on $x^{1}$ and simply permute the coordinate axes whenever we wish to make a different choice. It will be convenient to establish a special notation for the components of $\nabla U({\bf x})$ that reflects this convention. Thus we define:

[TABLE]

and observe that the term $P_{0}({\bf x})=\partial U({\bf x})/\partial x^{1}$ will play a special role because of our centering convention. We now define the basis vectors for $E$ as follows:

[TABLE]

It is straightforward to verify that $\nabla U({\bf x})$ is orthogonal to each ${\bf V}_{i}({\bf x})$, but we need to analyze the tangent subbundle more carefully to verify the Frobenius integrability conditions.

It is standard in differential geometry to think of a vector field as a differential operator, essentially the directional derivative with respect to a given vector ${\bf V}$. We will write this in shorthand notation as ${\bf V}\partial$. Let’s now consider the vector fields defined by ${\bf V}_{i}={\bf V}_{i}({\bf x})/P_{0}({\bf x})$ and ${\bf V}_{j}={\bf V}_{j}({\bf x})/P_{0}({\bf x})$, with $i\neq j$, and let’s compute the Lie bracket of ${\bf V}_{i}\partial$ and ${\bf V}_{j}\partial$. By a straightforward (but tedious) calculation, we have:

[TABLE]

We now substitute the definitions $P_{0}=\partial U({\bf x})/\partial x^{1}$, $P_{i}=\partial U({\bf x})/\partial x^{i+1}$ and $P_{j}=\partial U({\bf x})/\partial x^{j+1}$, and we note that the terms in the square brackets vanish identically by virtue of the equality of mixed partial derivatives. Thus $\left[{\bf V}_{i}\partial,{\bf V}_{j}\partial\right]=0$, which means that the vector fields ${\bf V}_{i}\partial$ and ${\bf V}_{j}\partial$ commute.

To formulate the Theorem of Frobenius, we need several definitions. We say that the vector field ${\bf V}\partial$ belongs to the tangent subbundle $E$ if ${\bf V}({\bf x})$ is an element of $E$ at each point ${\bf x}$ of the domain. Then, if $\left[{\bf X}\partial,{\bf Y}\partial\right]$ belongs to $E$ whenever ${\bf X}\partial$ belongs to $E$ and ${\bf Y}\partial$ belongs to $E$, for arbitrary ${\bf X}$ and ${\bf Y}$, we say that $E$ is involutive. If the tangent subbundle $E$ can be extended to a full integral manifold, we say that $E$ is integrable. These two concepts are related by the following:

Theorem 1** **(Frobenius).

A tangent subbundle, $E$, is integrable if and only if it is involutive.

Proof.

See [BC01], Sections 1.4 and 1.6, and Theorems 5, 6 and 7; [Lan95], Chapter VI, $\S$1 – $\S$4, and Theorems 1.1 and 1.2. ∎

We can now show that the tangent subbundle defined above by the basis vectors ${\bf V}_{1}({\bf x}),{\bf V}_{2}({\bf x}),\ldots,{\bf V}_{n-1}({\bf x})$, is involutive, hence integrable.

If ${\bf V}\partial$ and ${\bf W}\partial$ are vector fields and $f$ and $g$ are differentiable real-valued functions, we have the following identity for the expansion of Lie brackets:

[TABLE]

Now let ${\bf X}\partial$ and ${\bf Y}\partial$ be two arbitrary vector fields that belong to $E$, and write each of them in terms of their basis vectors:

[TABLE]

To compute $\left[{\bf X}\partial,{\bf Y}\partial\right]$, we apply (17) and use the fact that all terms in the form $\left[{\bf V}_{i}\partial,{\bf V}_{j}\partial\right]$ vanish because of (12), to show that the remaining terms form a linear combination of the basis vectors of E. Thus $\left[{\bf X}\partial,{\bf Y}\partial\right]$ belongs to $E$, which means that $E$ is involutive.

Figure 3 shows two views of the integral manifold for the three-dimensional example that was depicted in Figures 1 and 2 in Section 2. The view on the left is based on a coordinate system centered on the $x$ axis, and it shows several integral curves in the $xy$ plane (with $z$ constant) and the $xz$ plane (with $y$ constant). The view on the right is based on a coordinate system centered on the $y$ axis, which was computed by a permutation of the axes resulting in the definition $P_{0}({\bf x})=\partial U({\bf x})/\partial x^{2}$. It shows the same integral curves in the $xy$ plane (with $z$ constant) along with several new integral curves in the $yz$ plane (with $x$ constant). It should be clear that both sets of integral curves are tracing out the same two-dimensional integral manifold. Can this integral manifold be defined without reference to a specific coordinate system?

In [McC18], we reversed the procedure that we have been following here. Theorem 5 in [McC18] asserts the existence of a two-dimensional integral manifold orthogonal to any vector field that is proportional to the gradient of a scalar potential function, i.e., any vector field in the form $N({\bf x})\nabla U({\bf x})$. The theorem, as stated, makes use of the vector cross product and the “curl,” which is a three-dimensional concept, but it is actually a special case of a general result in ${\bf R}^{n}$ which follows from the dual version of the Theorem of Frobenius, expressed in terms of differential forms. Thus, if we wanted to, we could develop the theory of our integral manifold in a more “intrinsic” way, without reference to a special coordinate system. We would still have to introduce a coordinate system, of course, when we wanted to do computations, as we did in [McC18], but this would not be our starting point.

We could use either set of integral curves in Figure 3 to define a curvilinear coordinate system on the Frobenius integral manifold, but it would be a global coordinate system, since it follows the global Cartesian coordinates from the ambient space ${\bf R}^{n}$. This would not be particularly useful if we are looking for an “optimal” $k-1$ dimensional subspace. But another approach is to use these global vector fields to construct a local coordinate system. Any linear combination of ${\bf V}_{1}({\bf x}),{\bf V}_{2}({\bf x}),\ldots,{\bf V}_{n-1}({\bf x})$, could be taken as one of the basis vectors for the tangent subbundle, and we can vary this linear combination as we move around the integral manifold. To implement this idea, it is useful to define a Riemannian metric on the integral manifold, and the most natural way to do this is to define a metric tensor on all of ${\bf R}^{n}$, using the inner products of $\nabla U({\bf x})$, ${\bf V}_{1}({\bf x})$, ${\bf V}_{2}({\bf x})$, $\ldots$ , and ${\bf V}_{n-1}({\bf x})$, in that order. We thus define:

[TABLE]

If we want a uniform coordinate notation in place of $\rho$ and $\Theta$, we let $i$ and $j$ range over [math], $1$, $2$, $\ldots$, $n-1$, and we stipulate that $u^{0}=\rho$, and $u^{i}=\theta^{i}$, for $i=1,\ldots,n-1$. We will see later why this formula makes sense as a measure of dissimilarity, but for now we will simply adopt it as the definition of our Riemannian dissimilarity metric.

What is the purpose of this Riemannian dissimilarity metric? The main application of the metric is to compute geodesics on the surface of the integral manifold orthogonal to $\nabla U({\bf x})$. Since any linear combination of ${\bf V}_{1}({\bf x}),{\bf V}_{2}({\bf x}),\ldots,{\bf V}_{n-1}({\bf x})$, yields a vector in the tangent subbundle, $E$, we can construct vector fields in $E$ in the form

[TABLE]

for arbitrary functions

[TABLE]

For a geodesic, we are looking for a curve $\gamma(t)$ with values in ${\bf R}^{n}$ which minimizes the “energy” functional:

[TABLE]

subject to the constraint:

[TABLE]

This variational problem leads to a system of Euler-Lagrange equations for the curves $\gamma(t)$ and $\mathbf{v}(t)$, plus $n$ Lagrange multipliers, and by the existence and uniqueness theorems for ordinary differential equations, the resulting system will have a solution if we specify the initial conditions $\gamma(0)$ and $\mathbf{v}(0)$.

Figure 4 shows a local coordinate system for the three-dimensional example that was depicted in Figures 1 and 2 in Section 2, and for which we computed the global coordinate systems in Figure 3. We will use this illustration to explain the necessary calculations, but we are interested in generalizing from three dimensions to $n$ dimensions, and we will see that the higher dimensional case sometimes requires a slightly different treatment. In either case, there are three steps:

• •

Step One:

To find a principal axis for the $\rho$ coordinate, we minimize the Riemannian distance, ${g}_{i,j}({\bf x})$, along the drift vector.

• •

Step Two:

To choose the principal directions for the $\theta^{1},\theta^{2},\ldots,$ $\theta^{k-1}$ coordinates, we diagonalize the Riemannian matrix, $\left(\,{g}_{i,j}({\bf x})\,\right)$, at a fixed point, ${\bf x}_{0}$, along the principal axis.

• •

Step Three:

To compute the $\Theta$ coordinate curves, we follow the geodesics of the Riemannian metric, ${g}_{i,j}({\bf x}_{0})$, in each of the $k-1$ principal directions.

It turns out that Step Two requires the greatest modifications in the higher dimensional case, so we will defer that discussion and consider Step One and Step Three first.

In Step One, we compute the principal axes, which are the curves in Figure 4 that extend from the black dots to the red dot. One way to do this is to find a point at a fixed Euclidean distance from the origin for which the integral curve, $\rho(t)$, as defined in Equation (10) or (11), has minimal Riemannian length to the origin. However, if $|\nabla U({\bf x})|$ is monotonic, we obtain approximately the same results (i.e., within the accuracy of our numerical approximations) by simply minimizing $|\nabla U({\bf x})|^{2}$ on a sphere at a constant distance from the origin. For example, to find the black dot in the lower right corner of Figure 4, we first minimize $|\nabla U(x,y,z)|^{2}$ on the sphere $x^{2}+y^{2}+z^{2}=500$ to locate a point and a curve, $\rho(t)$, with Riemannian length $6.30873$. We then search within the neighborhood of this solution to find a point and a curve with a slightly smaller Riemannian length: $6.30863$. In Figure 4, we have chosen the latter solution, but in an $n$-dimensional space, the second solution is usually too complex, computationally, and we therefore stick with the simple minimization of $|\nabla U({\bf x})|^{2}$.

Let’s now jump ahead to Step Three. Setting $\gamma(0)={\bf x}_{0}$, which is the black dot on the principal axis in Figure 4, and setting $\mathbf{v}(0)$ equal to one of the principal directions from Step Two, we solve the Euler-Lagrange equations. The results are the coordinate curves $\theta^{1}$ and $\theta^{2}$ in Figure 4. These curves have interesting geometric properties. Since they are geodesics on the surface of the Frobenius integral manifold, they are always a constant Riemannian distance from the origin. Note that we have drawn several $\rho$ coordinate curves in Figure 4 from the $\theta^{2}$ curve to the origin. The Riemannian length of these curves is approximately $6.30863$, even though their Euclidean length varies considerably. On the other hand, the Riemannian distance along the $\theta^{1}$ and $\theta^{2}$ curves is the same as the Euclidean distance, and we can therefore use the Riemannian/Euclidean arc length to parameterize these curves, just as we have done with the $\rho$ coordinate curves.

In practice, when we try to solve the Euler-Lagrange equations in Step Three of our procedure, we often encounter a singularity, as the value of $P_{0}({\bf x})$ approaches zero. The solution is to shift the “center” of our coordinate system to a different coordinate axis, and then continue the computation with a new value for $P_{0}({\bf x})$. To validate this solution, though, we need to show that our definition of the Riemannian metric, ${g}_{i,j}({\bf x})$, is independent of our choice of a “centered” coordinate axis. For the three-dimensional case, the proof appears at the end of Section 4 in [McC18], and the proof for $n$ dimensions is a straightforward generalization. Let $\mathbf{u}(\mathbf{x})$ denote a $\rho$, $\Theta$, coordinate system centered on the $x^{1}$ axis, and let $\mathbf{\bar{u}}(\mathbf{x})$ denote a $\rho$, $\Theta$, coordinate system centered on the axis $x^{p}$, for some $p\neq 1$. (For example, look at the second image in Figure 3, in which $p=2$ and the centered coordinate axis is $x^{2}$.) The Jacobian matrix of the coordinate transformation from $\mathbf{\bar{u}}(\mathbf{x})$ to $\mathbf{u}(\mathbf{x})$ can be computed as follows:

[TABLE]

Now let ${g}_{i,j}({\bf x})$ and ${\bar{g}}_{k,l}({\bf x})$ denote the dissimilarity metric based on the $u^{i}$ and $\bar{u}^{k}$ coordinates, respectively. We can verify by a straightforward computation (compare the three-dimensional example in [McC18]) that

[TABLE]

But this is just an instantiation of the transformation law for a type (0,2) tensor. Thus, on an integral manifold of dimension $n-1$, for a fixed $\rho$, the dissimilarity metric, ${g}_{i,j}({\bf x})$, is independent of the global coordinate system used to define it.

Once we have computed the geodesic coordinate curves $\theta^{1},\theta^{2},\ldots,$ $\theta^{k-1}$, all of which emanate from ${\bf x}_{0}$, we also need to construct a system of transverse coordinate curves, which can emanate from any point, ${\bf x}$, on the integral manifold. For this purpose, we define the following flows:

[TABLE]

Here, $\hat{\theta}_{{\bf x}}(t)$ is an integral curve starting at ${\bf x}$, as in (11), and it has an equivalent definition by a differential equation, as in (10). We have one such flow equation for each geodesic, $\theta^{1},\theta^{2},\ldots,$ $\theta^{k-1}$, computed in Step Three, in which the coefficient vectors, $\mathbf{v}_{\theta^{1}}(s),\mathbf{v}_{\theta^{2}}(s),\ldots,\mathbf{v}_{\theta^{k-1}}(s)$, are the functions, $\mathbf{v}(t)$, that were computed at the same time. When ${\bf x}={\bf x}_{0}$, of course, the integral curve in Equation (25) coincides with the original geodesic coordinate curve. In the $n$-dimensional case, we usually apply Equation (25) in a fixed order to a sequence of coordinates to specify a point on the integral manifold. Thus, we start with the flow $\vec{\theta}^{\,1}_{t}({\bf x}_{0})$ and follow it for a distance $t_{1}$ to the point ${\bf x}_{1}$; we then follow the flow $\vec{\theta}^{\,2}_{t}({\bf x}_{1})$ for a distance $t_{2}$ to the point ${\bf x}_{2}$; and so on. In the three-dimensional case, as shown in Figure 4, it is convenient to define an additional geodesic coordinate curve, $\phi$, which is orthogonal to both $\theta^{1}$ and $\theta^{2}$. We then have three flows, $\vec{\theta}^{\,1}_{t}({\bf x})$, $\vec{\theta}^{\,2}_{t}({\bf x})$, and $\vec{\phi}_{t}({\bf x})$, and we can compose them either as $\vec{\phi}_{s}\circ\vec{\theta}^{\,1}_{t}({\bf x}_{0})$ or as $\vec{\phi}_{s}\circ\vec{\theta}^{\,2}_{t}({\bf x}_{0})$. (We will see later that this choice leads to an interesting comparative analysis.) Whatever choices we make, though, we want to make sure that our coordinate system covers the entire integral manifold.

Let’s now return to Step Two, where we encounter a surprising mathematical fact. When we diagonalize the Riemannian matrix, $\left(\,{g}_{i,j}({\bf x})\,\right)$, in $n$ dimensions, we find only two distinct eigenvalues. The largest eigenvalue has multiplicity $2$: $\lambda_{0}=\lambda_{1}=|\nabla U|^{2}$, with corresponding eigenvectors:

[TABLE]

The smallest eigenvalue has multiplicity $n-2$: $\lambda_{2}=\lambda_{3}=\ldots=\lambda_{n-2}=\lambda_{n-1}=P_{0}^{2}$, with corresponding eigenvectors $\xi_{2},\;\xi_{3},\;\ldots,\;\xi_{n-2},\;\xi_{n-1}$, as follows:

[TABLE]

This means that the three-dimensional case is special, since it only has three eigenvectors:

[TABLE]

See Section 4 of [McC18]. We should think of these $\{\xi_{i}\}$ as infinitesimal eigenvectors. They tell us the maximal and minimal initial directions for the integrand of the energy functional:

[TABLE]

However, we are primarily interested in minimizing geodesic curves over finite distances, and there is no guarantee that minimizing the initial directions of the geodesics in the Euler-Lagrange equations will achieve this result.

Our solution to this problem in [McC18] is to rotate the infinitesimal eigenvectors around ${\bf x}_{0}$, and then compute the Riemannian/Euclidean distances along each rotated geodesic, up to some specified point, for example, up to the Euclidean angle $\pi/2$ from the origin, so that we can determine the global minimum or maximum. In Figure 4, the orange arrows depict the minimal infinitesimal eigenvector, $\xi_{2}$, in the positive $y$-direction and the negative $y$-direction, respectively. But the minimal geodesics over a finite distance, labelled as $\theta^{2}$, were obtained by a counter-clockwise rotation through the angle $\alpha=0.952169$ in the positive direction, and the angle $\alpha=1.12681$ in the negative direction. For the maximal infinitesimal eigenvector, $\xi_{1}$, the maximal geodesics over a finite distance, labelled as $\theta^{1}$, were obtained by a counter-clockwise rotation through the angle $\alpha=0.114166$. See Section 5.2 of [McC18] for the detailed calculations.

Now look at Figure 5. This figure includes the same view of the $\theta^{1}$ and $\theta^{2}$ coordinate curves as in Figure 4. It also includes $1000$ data points generated according to the curvilinear Gaussian probability distribution from Figure 1, and projected along the $\rho$ coordinate curve to the Frobenius integral manifold. Qualitatively, the density of the data is higher near the $\theta^{2}$ coordinate curve than it is near the $\theta^{1}$ coordinate curve, and we can quantify this observation by computing a value for the “reconstruction error” adapted to our curvilinear coordinate system. For this purpose, we use the geodesic coordinate curve, $\phi$, which is aligned with the green arrow in Figure 4, and its flow, $\vec{\phi}_{s}({\bf x})$. We want to compare two curvilinear coordinate systems for each data point: $(\rho,\theta^{1},\phi)$, which is defined by following the flow $\vec{\phi}_{s}\circ\vec{\theta}^{\,1}_{t}({\bf x}_{0})$ on the integral manifold, and $(\rho,\theta^{2},\phi)$, which is defined by following the flow $\vec{\phi}_{s}\circ\vec{\theta}^{\,2}_{t}({\bf x}_{0})$. What happens when we drop the $\phi$ coordinate, in each case? We can measure the “reconstruction error” by computing the Euclidean distance along the flow, $\vec{\phi}_{s}$, and scaling this value down, proportionately, by the position of the data point along the $\rho$ coordinate curve. We can then compute the root-mean-squared (RMS) reconstruction error for the 1000 data points in each coordinate system. It turns out that the RMS error for the truncation from $(\rho,\theta^{1},\phi)$ to $(\rho,\theta^{1})$ is 5.9431, and the RMS error for the truncation from $(\rho,\theta^{2},\phi)$ to $(\rho,\theta^{2})$ is 4.82787. Thus, according to this criterion, the “optimal” lower dimensional encoding is $(\rho,\theta^{2})$. See Section 5.2 of [McC18] again for the detailed calculations.

Recall that the $\theta^{2}$ coordinate curve was derived from the minimal infinitesimal eigenvector, $\xi_{2}$, by rotating it a small amount in order to minimize the length of the geodesic over a finite Euclidean angle. This is a purely geometric calculation. However, our analysis of the reconstruction error shows that projection of the data onto the $\theta^{2}$ subspace has better statistical properties than projection onto the $\theta^{1}$ subspace. Thus the link between the geometric model and the probabilistic model has computational implications. Here is an informal explanation:

$\dots$ The geodesic curves on the Frobenius integral manifold tend to follow the modes of the probability distribution. First, the origin of the coordinate system is a point at which $\nabla U({\bf x})=(0,0,0)$, which maximizes the probability density. Second, to compute the principal axis, we are looking for a point with a minimal Riemannian distance for a fixed Euclidean distance, or a maximal Euclidean distance for a fixed Riemannian distance. Under either formulation, this is an axis that maximizes probability. Third, for the directional coordinates, we are looking for a geodesic curve on the Frobenius integral manifold that covers a minimal Riemannian distance for a fixed angular Euclidean distance, or a maximal angular Euclidean distance for a fixed Riemannian distance. Under either formulation, again, this is a curve that maximizes probability. Thus, in general, we are minimizing dissimilarity and maximizing probability. This is the primary intuition behind the claim that we are constructing an “optimal” lower dimensional coordinate system.

See [McC18], Section 5.2.

For a more formal analysis of the link between the geometric model and the probabilistic model, and a possible mathematical explanation of why the clustering of data points in Figure 5 occurs, see Section 6 of [McC18].

What do we need to do to extend these ideas to $n$ dimensions? Rotating the eigenvectors to find the optimal initial directions is not likely to work in $n$ dimensions. In Figure 4, we rotated the minimal infinitesimal eigenvector in a two-dimensional plane, and we were able to compute and compare the geodesics under each rotation. These rotations and calculations would not be feasible in an $n-1$ dimensional hyperplane. However, the more important principle suggested by this example is that we should search among a set of candidates for a minimal geodesic over a finite distance. In $n$ dimensions, we can apply this principle to the set of minimal infinitesimal eigenvectors. First, note that any linear combination of the eigenvectors, $\xi_{2},\;\xi_{3},\;\ldots,\;\xi_{n-2},\;\xi_{n-1}$, is also an eigenvector associated with the eigenvalue $\lambda_{2}=P_{0}^{2}$. We can therefore define new eigenvectors in the form:

[TABLE]

where the coefficients, $c_{i}$, are either $+1$ or $-1$. (There are $2^{n-2}$ sequences of such coefficients, of course, but in practice we can generate a large subset randomly.) If we evaluate the integrand of the energy functional on these new eigenvectors, we have:

[TABLE]

Now, order the eigenvectors, $\{\zeta\}$, by their norms and apply the Gram-Schmidt orthogonalization procedure to a subsequence of length $n-2$ or more. We will then have an orthonormal basis for the $n-2$ dimensional tangent subbundle at ${\bf x}_{0}$, in which each basis vector represents a minimal initial direction for the solution of the Euler-Lagrange equations, by (81). Finally, from among these solutions, our initial strategy is simply to select the $k$ coordinate curves that have the minimal Riemannian length over a fixed angular Euclidean distance. We will see how well this works, using real data, in Sections 6 and 7, infra.

We will also analyze in Sections 6 and 7 a representation of clusters, using real data. The basic idea is illustrated for the three-dimensional case in Figure 6, which shows two copies of the curvilinear Gaussian potential from Figure 1. One copy has been translated from $(0,0,0)$ to $(20,20,-10)$. The other copy has been translated from $(0,0,0)$ to $(-20,-20,10)$ and rotated by $\pi/2$ around a line parallel to the $y$-axis. But the probability density is a mixture. If $U_{1}({\bf x})$ is the potential function for the first copy and $U_{2}({\bf x})$ is the potential function for the second copy, then the invariant probability density is given by:

[TABLE]

modulo an appropriate normalization factor. The advantage of this representation lies in the fact that our calculations for each copy will be almost independent of each other. Observe that the effective potential function for the mixture will be:

[TABLE]

Thus the gradient of $U({\bf x})$ in a neighborhood of $(20,20,-10)$ will be almost identical to the gradient of $U_{1}({\bf x})$ computed by itself, and the gradient of $U({\bf x})$ in a neighborhood of $(-20,-20,10)$ will be almost identical to the gradient of $U_{2}({\bf x})$ computed by itself. Or, in terms of our dissimilarity metric, the two clusters in Figure 6 will be exponentially far apart.

4. How to Estimate $\nabla U({\bf x})$ from Sample Data.

To apply the theory of differential similarity to real data, we need to estimate the quantities that appear in the equations for the $\rho$, $\Theta$, coordinate system. Fortunately, this is not hard to do: We will borrow a technique from the literature on the mean shift algorithm [FH75] [Che95] [CM02].

Consider a kernel density estimator with a Gaussian kernel:

[TABLE]

in which $\mathbf{s}_{k}$ is a sample data point and $\beta$ is a smoothing parameter. We can approximate a probability density by taking the average over these kernels:

[TABLE]

Now recall that $\nabla U({\bf x})$ is the gradient of the log of the stationary probability density in our theory. So we can differentiate explicitly:

[TABLE]

to obtain an estimate for $\nabla U({\bf x})$.

Figure 7 shows an example of Mathematica code that implements Equation (82) and calls it pmeanshift. The expression “mean shift” refers to the fact that the first term inside the square brackets computes the weighted mean of the sample data points, $\{\mathbf{s}_{k}\}$, around the point $\mathbf{x}$, and the second term subtracts $\mathbf{x}$ from this value to define a shift. This function is typically used in a gradient ascent algorithm, such as pgradascent. (In Figure 7, byteArray7X7 and zeroArray are tailored to the examples in Sections 6 and 7, infra.) In this case, the gradient ascent algorithm will find the mode of the probability density, $\hat{\mu}(\mathbf{x})$, for a given Sample. It is often applied iteratively: Choose a sample around $\mathbf{x}$, ascend to the mode of the probability density to find a new $\mathbf{x}$, choose another sample, and repeat.

We can certainly apply the code in Figure 7 to find the origin of an admissible coordinate system in our theory, and we will see examples of this kind of an application in Sections 6 and 7. However, Equation (82) has a much broader application than this, because every quantity that enters into the definition of the $\rho$ coordinates and the $\Theta$ coordinates depends on $\nabla U({\bf x})$. We will actually modify (82) slightly to simplify these calculations. The quantity $\sum_{k=1}^{n}K(\mathbf{s}_{k},\mathbf{x})$ in the denominator of the first term is a normalization factor, and we can multiply the formula by this factor (and divide out the factor $2\beta$) to write down an equation for the gradient without normalization:

[TABLE]

The second derivatives are now much easier to compute:

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta function. These modifications will not alter the specification of our coordinate system, for two reasons: (1) the definition of the $\rho$ coordinate in Equations (10) and (11) is already normalized, so the factor $\sum_{k=1}^{n}K(\mathbf{s}_{k},\mathbf{x})$ would cancel out; and (2) the basis vectors $\{{\bf V}_{i}({\bf x})\}$ can be divided by $P_{0}({\bf x})$ without affecting the Frobenius integral manifold, so again the factor $\sum_{k=1}^{n}K(\mathbf{s}_{k},\mathbf{x})$ would cancel out. The code for $\mathtt{DUjF}$ and $\mathtt{DiDUjF}$ is shown in Figure 8. The variables ${\beta}\mathtt{eta}$ and $\mathtt{SamplePoints}$ are intended to be defined externally, and the vector variable $\mathtt{VexpK}$ is expected to be instantiated by the kernel function $\mathtt{VkernelF}$.

5. Computing the Geodesic Coordinate Curves.

The main data structures in our theory are defined in Figure 9. But the only data structure needed to write down the Euler-Lagrange equations is the gradient: $\nabla U({\bf x})\,=\,(\,P_{0}({\bf x}),P_{1}({\bf x}),\ldots,P_{n-1}({\bf x})\,)$. Here are the Euler-Lagrange equations:

[TABLE]

This is a system of ordinary differential equations for $\mathbf{x}(t)$ and for the Lagrange multipliers, $\{\lambda_{i}(t)\}$, plus a system of algebraic equations for $\mathbf{v}(t)$, in which $d$ is the dimensionality of the space. In addition, we have to include the set of kernel equations, which enter into the definition of the gradient. So this is a very large system of equations, approximately $3*\mathtt{Ndims}+\mathtt{Nsamples}$.

Figure 10 shows the Mathematica code for the right-hand sides of the Euler-Lagrange equations, in three groups. Tracing this code back from Figure 10 to Figure 8, it is easy to verify that everything depends on the set of $\mathtt{SamplePoints}$ and the constant ${\beta}\mathtt{eta}$. Thus we should be able to apply this code to a real dataset. We will see how this works in Section 6 on the MNIST dataset, and in Section 7 on the CIFAR-10 dataset.

6. Example: 7$\times$7 Patches in the MNIST Dataset.

Figure 11 displays an architecture for deep learning on the MNIST dataset [LBBH98], based on several examples in the recent literature [Ran09] [Coa12]. The process starts in the lower-left corner and follows the arrows to the upper-right corner. The first step is to scan and randomly sample the 60,000 28$\times$28 images to extract a collection of 7$\times$7 “patches” from each one. Choosing a sampling rate of 10 scans per image, which is approximately 2%, we end up with 600,000 patches, each one represented as a point in a 49-dimensional space. The original image intensity at each pixel is represented by an integer in the range $[0,255]$, but we have scaled these values down to a real number in the range $[0,1]$. Thus the greatest distance between any two points in our 49-dimensional hypercube is $\sqrt{49}=7.0$. We will be using the theory of differential similarity, in this section of the paper, to reduce the dimensionality of the space to 12 dimensions, as shown in the lower-right corner of Figure 11. Note that the process continues upwards (following the blue arrow) by assembling four adjacent 7$\times$7 patches into a 2$\times$2 matrix and then resampling the image using the larger 14$\times$14 patch, but we will not pursue this analysis in the present paper. We are thus looking only at the first step in the process, which is an example of the classical problem of unsupervised feature learning.

There are several parameters that control the behavior of the algorithms defined in Sections 4 and 5. The smoothing parameter, $\beta$, can be set to different values in different circumstances. (In the traditional notation for a Gaussian, $\beta=1\,/\,2\sigma^{2}$.) We have experimented with a range of values, but we have found that the two values, $\beta=1/8$ (or $\sigma=2.0$) and $\beta=1$ (or $\sigma=0.707107$), are sufficient for most purposes. A typical strategy is to run a computation with $\beta=1/8$, for a coarse approximation, and then to refine the result by running the computation again with $\beta=1$. Another parameter is the variable $\mathtt{SamplePoints}$. Mathematica provides a function called Nearest which, when applied to a large dataset, returns a NearestFunction which maintains an efficient data structure that can find the $n$ points that are the nearest to any specified point in the dataset. Again, we have experimented with several values of $n$, but we have found that a Data Sphere of 32,000 points around each prototype yields good results. Within each Data Sphere, we draw two random 800 point samples and use these as inputs to the algorithms in Sections 4 and 5. A typical strategy is to run the computations separately on each 800 point sample, and if the results are qualitatively the same and quantitatively within the range that we would expect from random sampling, then we run a final computation on the union of the two samples. Note that a 1600 point sample is $5\%$ of a 32,000 point Data Sphere.

In addition to the definition of the Data Sphere, it is also useful to define a Coordinate Sphere for each prototype. Recall that the principal axis is defined by a point at a fixed Euclidean distance from the origin that has a minimal Riemannian distance to the origin, as measured along the $\rho$ coordinate curve. The fixed Euclidean distance gives us a sphere, of course, and it is reasonable again to think of the size of this sphere in terms of the number of data points it contains. Combined with a 32,000 point Data Sphere, we have found that an 8,000 point Coordinate Sphere yields good results.

Figure 12 shows an initial selection of 35 prototypes for the 600,000 7$\times$7 patches. The prototypes outlined in blue were selected (subjectively) for a more detailed investigation. (Prototype 02, in which most of the pixels are black, was only used to provide a lower dimensional test case for the development of our algorithms) We will therefore focus our attention on the three examples: Prototype 09, Prototype 24, Prototype 27. The full set of 35 prototypes was constructed by a combination of: (i) pgradascent in Figure 7, with xstart generated randomly, with $\beta=1/8$, and with a 2000 point Sample retrieved by the NearestFunction at each iteration; plus (ii) a traditional clustering algorithm (FindClusters in Mathematica) applied to the output of pgradascent; and finally (iii) some manual pruning at the end of the process to eliminate redundancies. We then constructed a 32,000 point Data Sphere around each prototype, and used these for all the subsequent data analyses. Note that $35\times 32,000=1,120,000$, so we can potentially partition the entire dataset almost twice over. In fact, an analysis ex post shows that the union of the 32,000 point Data Spheres around these 35 prototypes covers $92.5\%$ of the distinct 7$\times$7 patches in our sample. So there is room for some improvement here, but not much.

Although each prototype is located at the center of its Data Sphere, by definition, it will not necessarily be located at a mode of the probability distribution, since we have altered the sample in pgradascent to include the nearest 32,000 data points. Thus we now compute the integral curve of $\nabla U$, starting at the original prototype, to arrive at a modified prototype, where $\nabla U({\bf x})=0$. The results are shown in the center column of Figure 13. We then draw a Coordinate Sphere around each modified prototype, and apply Step One of the procedure defined in Section 3 to compute the points at which the $\rho$ coordinate curves drawn inwards to the origin have minimal Riemannian length. These points are shown in the right column of Figure 13. The $\rho$ coordinate curves themselves are shown in Figure 14.

Table 1 shows the size of the Data Sphere and the Coordinate Sphere for the three prototypical clusters that we are analyzing. These values should be compared with the maximal distance in the hypercube ($\sqrt{49}=7.0$), and with the standard deviation of the kernel density estimator when the smoothing parameter $\beta=1$ ($\sigma=0.707107$).

The distances in Table 2 reveal some of the main properties of our geometric model. The Euclidean distance along each $\rho$ coordinate curve, Column (c) in Table 2, is greater than the Euclidean distance between its two end points, Column (b) in Table 1. This is an indicator of the fact that the principal axis is a curve, and not a straight line. The Riemannian distance inwards along the $\rho$ coordinate curves, Column (d) in Table 2, is much less than the Euclidean distance in Column (c). This is an indicator of the fact that the principal axis is traversing a region of high probability density. But Column (e) in Table 2, which displays the Riemannian distance computed outwards along the $\rho$ coordinate curve for approximately the same Euclidean distance, is much larger because the probability density in this region is much smaller.

It is also interesting to analyze Figures 13 and 14 qualitatively. These images were constructed by solving a system of differential equations on pixels in a high-dimensional image space, a very local process, but they seem to possess certain global coherence properties. For Prototype 09, we see a simple geometric shape modified by a shift transformation. The mapping from the left column to the center column in Figure 13 is basically a two-dimensional shift, one pixel up and one pixel to the right. The mapping from the right column to the center column, which is displayed in Figure 14 as a discretely sampled continuous curve, shifts the same geometric shape up by two pixels. For Prototype 27, we see another simple geometric shape modified by a sequence of global transformations. The mapping from the left column to the center column in Figure 13 is basically a counter clockwise rotation. The lower left corner is fixed, approximately, while the upper end of the geometric shape is shifted one pixel to the left. The mapping from the right column to the center column in Figure 13, which is displayed in Figure 14 as a discretely sampled continuous curve, shifts the same geometric shape one pixel to the left. Prototype 24 is somewhat more complex. Here the discretely sampled continuous mapping in Figure 14 shifts the left half of the geometric shape down by one pixel, while it bends the right half of the geometric shape up to construct a thick horizontal line.

Keeping these coherence properties in mind, let’s now look at Step Two and Step Three of the procedure defined in Section 3. In Step Two, we first diagonalize the Riemannian matrix, $\left(\,{g}_{i,j}({\bf x})\,\right)$, at the points on the principal axes shown in the right column of Figure 13. Taking Prototype 09 as an example, the two eigenvalues are 35.7443 and 13.572, and the infinitesimal eigenvectors are constructed to have unit norm. Thus the integrand of the energy functional for the maximal infinitesimal eigenvector, $\xi_{1}$, is 35.7443. (See Equation (78) in Section 3.) For the minimal infinitesimal eigenvectors, we adopt the following strategy: We generate 10,000 random linear combinations of the eigenvectors, $\xi_{2},\;\xi_{3},\;\ldots,\;\xi_{48}$, using Equation (80), we sort these linear combinations by their norms, and we select every 200th entry in the sorted list. We then apply the Gram-Schmidt orthogonalization procedure (using Othogonalize in Mathematica) to the selected eigenvectors to construct a set of 47 orthonormal basis vectors. For each of these basis vectors, the integrand of the energy functional is 13.572. (See Equation (81) in Section 3.) Furthermore, since $\xi_{1}$ was orthogonal to all of the eigenvectors, $\xi_{2},\;\xi_{3},\;\ldots,\;\xi_{48}$, we now have an orthonormal basis for the 48-dimensional tangent subbundle at the point on the principal axis that will serve as the origin of our $\Theta$ coordinate curves.

Figure 15 displays the $\Theta$ coordinate curves that were computed in Step Three by solving the Euler-Lagrange equations with the initial directions, $\mathbf{v}(0)$, set to the values of the maximal infinitesimal eigenvectors, $\xi_{1}$. (The first row for each prototype follows in the positive direction; the second row follows in the negative direction; but this is an arbitrary convention.) The curves are continued until they reach a point at a Euclidean angle of $\pi/2$ from the origin. Sometimes the computation encounters a singularity in the value of $P_{0}({\bf x})$, and it is necessary to shift the “center” of the coordinate system to a different coordinate axis in order to continue with the algorithm, but this does not happen more than once for each of these curves. Table 3 shows the distances along the coordinate curves to the Euclidean angles $\pi/2$ in the positive direction and $-\pi/2$ in the negative direction. Recall that the Riemannian and Euclidean distances computed along a geodesic curve on the Frobenius integral manifold are the same, since each point on the manifold is a constant Riemannian distance from the origin. Thus there is only a single entry for the distances in Table 3, for each direction.

For the minimal infinitesimal eigenvectors, we will work exclusively with Prototype 09, as an example, but we will compute all 47 $\Theta$ coordinate curves. (Note: All calculations up to this point have used a 1600 point sample, following a comparison of the results for two 800 point samples. However, when we double the size of the sample from 800 to 1600, we increase the size of the data structures in Figure 10 by a factor of four, and we also increase the CPU time to find a solution to the Euler-Lagrange equations, empirically, by a factor of four. Thus, since we are now computing 47 geodesics, in two directions, we will reduce the computational burden by restricting our analysis here to a single 800 point sample.) Figure 16 shows the three $\Theta$ coordinate curves with the smallest Riemannian and Euclidean distances in the positive direction, which also happen to be the curves that have the largest Riemannian and Euclidean distances in the negative direction. Table 4 shows these distances, up to the Euclidean angles $\pi/2$ and $-\pi/2$, respectively. The total distance is shown in the third column. Note that the total distance along the $\Theta$ coordinate curve for the maximal infinitesimal eigenvector in Prototype 09, according to Table 3, is considerably larger: 6.39242.

Figure 17 shows another way to understand the relationship between the curves in Figures 14, 15 and 16: a three-dimensional visualization. How to map 49 dimensions down to three? The black dot at the origin is the modified Prototype 09 in the center column of Figure 13. The 49-dimensional image space is rotated around the origin to align the $x$ axis with the principal axis, so that the coordinates at the tip of the red arrow are $(1.75322,0,0)$, where $1.75322$ is the radius of the Coordinate Sphere, as shown in Table 1. For the $y$ coordinate, we rotate the image space again around the tip of the red arrow to align the $y$ axis with the infinitesimal eigenvectors that are used to generate the $\Theta$ coordinate curves, but we do this separately for each curve. Thus, although the initial directions of the four curves depicted in Figure 17 are mutually orthogonal in the 49-dimensional image space, in the visualization they are all drawn in the same plane. Finally, we need to collapse the remaining 47 dimensions down to the $z$ axis, somehow. Our solution is simply to sum all the remaining coordinates after applying the $x$ and $y$ rotations to the first two, which means that the initial value of $z$ at the tip of the red arrow will be zero, and the length of the curves in the image space will correspond roughly to the length of the curves in the visualization. In the image space, the coordinate curves are extended to the Euclidean angles $\pi/2$ and $-\pi/2$, and this is true in the visualization as well. Note that the positive direction towards $\pi/2$ corresponds to a clockwise rotation in the visualization, as can be seen by comparing the length of the two segments of the black curve in Figure 17.

The green curves in Figure 17 also illustrate some of the main properties of our geometric model. These are the $\rho$ coordinate curves drawn inwards from the points on the $\Theta$ coordinate curve for the maximal infinitesimal eigenvector at the Euclidean angles $\pi/2$ and $-\pi/2$. The Euclidean distance along the curve from $\pi/2$ (on the left) is 1.50939 and the Riemannian distance is 0.947386. The Euclidean distance along the curve from $-\pi/2$ (on the right) is 2.1958 and the Riemannian distance is 0.922917. These values should be compared to the values in Table 2 for the corresponding distances along the principal axis for Prototype 09: The Euclidean distances in Figure 17 are less, but the Riemannian distances are approximately the same (within the tolerance of our numerical calculations). This is a reflection of the fact that the geodesics on the Frobenius integral manifold are a constant Riemannian distance from the origin.

We saw in our analysis of the three-dimensional curvilinear Gaussian example in Section 3 that the “optimal” $\Theta$ coordinate curve was the one with the minimal Riemannian or Euclidean length, when measured at a fixed Euclidean angle from the origin. What choices does this principle suggest for Prototype 09? Looking at Table 3 for the maximal infinitesimal eigenvector, we see that the distance in the positive direction for Prototype 09 is less than all of the distances in Table 4, but the distance in the negative direction is greater than all of the distances in Table 4. Thus the principle seems to be inconclusive in this case. In constructing Table 4 for the minimal infinitesimal eigenvectors, however, we chose the three curves with the smallest distances in the positive direction and the largest distances in the negative direction, and it is instructive to look at these choices as they are displayed in Figure 16. Qualitatively, all three curves, 16, 17, and 9, seem to be more compact and coherent in the positive direction than in the negative direction, which is consistent with the principle developed in Section 3. But in this case, we have additional choices. Table 5 lists the ten $\Theta$ coordinate curves for which the total distance is minimal. Although we do not have the space to display these examples, they are in both directions qualitatively similar to the curves in the positive direction in Figure 16.

Thus, for an “optimal” 12-dimensional coordinate system for Prototype 09, one possible choice would be: (i) the $\rho$ coordinate; (ii) the $\Theta$ coordinate in Table 3, which is depicted in Figure 15; and (iii) the ten $\Theta$ coordinates in Table 5. However, since these numbers are very close together, it is not clear if there is enough separation among them to make a difference.

Let us, therefore, proceed with a quantitative evaluation of the coordinates in Table 5. There are two issues: First, we need an efficient method to compute the $\rho$, $\Theta$, coordinates for any arbitrary point in the original 49-dimensional Euclidean space. Second, we need a higher-dimensional generalization of the “reconstruction error” that was defined for the three-dimensional case in Section 5.2 of [McC18], and applied to the projected data points in Figure 5, supra.

It is straightforward to compute the $\rho$ coordinate, and project a data point onto the Frobenius integral manifold. Then, to locate this projected point on the surface of the manifold, we need to follow the coordinate flows, given by Equation (25), in a fixed order. We can formulate this task as an optimization problem: Find a sequence of parameters, $t_{1},t_{2},t_{3},\ldots$, for the flows $\vec{\theta}^{\,1}_{t_{1}}({\bf x}_{0}),\vec{\theta}^{\,2}_{t_{2}}({\bf x}_{1}),\vec{\theta}^{\,3}_{t_{3}}({\bf x}_{2}),\ldots$, that takes us from the principal axis at ${\bf x}_{0}$ to a location as close as possible to the projected point on the manifold. (For an illustration of how this works in the three-dimensional case, see Figure 15 in Section 5.2 of [McC18].) However, if we tried to use the optimization code in Mathematica, say, to solve this problem directly, we would never get an answer. The complexity is too high: The parameters, $t_{1},t_{2},t_{3},\ldots$, are mutually dependent, and each step in the calculation requires a numerical solution of the integral equations in (25). Our approach, instead, is to first compute an ellipsoidal approximation to the manifold, based on the simple Gaussian example in Section 5.1 of [McC18]. (See, e.g., the paragraph in [McC18] that includes Equation (31).) The analytical expressions for the flows $\vec{\theta}_{s}({\bf x})$ and $\vec{\phi}_{t}({\bf x})$ in Section 5.1 of [McC18] can easily be extended to 49 dimensions, and it turns out that our optimization problem can be solved, on the 48-dimensional ellipsoidal surface, by FindMinimum in Mathematica. Once we have a solution for the ellipsoidal approximation, we can use it to guide the search for a solution to the integral equations in (25) on the Frobenius integral manifold itself.111For a detailed explanation and justification of this ellipsoidal approximation, and a discussion of how an analytical solution on the ellipsoidal surface can be converted into a sequence of coordinate flows on the manifold, see Section 7, infra, on the CIFAR-10 dataset. Even so, this last step is slow: In a 49-dimensional example, it takes approximately 20 times more CPU time than the initial calculation on the surface of the ellipsoid.

The definition of the reconstruction error in the three-dimensional case was illustrated in Figure 15 in [McC18]: We computed and compared the Riemannian/Euclidean distances along the two $\phi$ coordinate curves, one emanating from $\theta^{1}$, the other emanating from $\theta^{2}$, and both constrained to lie on the Frobenius integral manifold. Similarly, in the 49-dimensional case, we can compute the total Riemannian/Euclidean distance along the 10 coordinate curves in Table 5, and compare this to the total Riemannian/Euclidean distance along the remaining 37 coordinates. For this calculation, we will use 669 sample points. (Recall that we used a single 800 point sample to solve the Euler-Lagrange equations. For the reconstruction error, we will use a different 800 point sample, while excluding those points that lie in the opposite hemisphere from the principal axis. Also, for computational reasons, we will do the calculations on the ellipsoidal approximation rather than on the manifold itself.) The results are not good. The RMS distance for the 10 coordinates in Table 5 is 1.49091, and for the remaining 37 coordinates, the RMS distance is 5.48309. The ratio of these two numbers is 0.27191, which is almost exactly the ratio 10/37 (= 0.27027). Thus, choosing the 10 coordinate curves in Table 5 provides no advantage over a random shuffling of the coordinates.

However, when we look more closely at the calculation of the reconstruction error, we observe an interesting phenomenon: The total Riemannian/Euclidean distance along the 47 $\Theta$ coordinate curves, for each sample point, tends to be concentrated on a small number of coordinates. Thus, for each sample point, if we sort the $\Theta$ coordinate curves, and choose the 10 coordinates for which the optimized Riemannian/Euclidean distance computed along the flows towards the sample point is the greatest, we get a very different result. The RMS distance for the top 10 coordinates is now 3.48581, and for the remaining 37 coordinates it is 3.38028. The ratio of these two numbers is 1.03122. Furthermore, the top 10 coordinates tend to fall into clusters. We will return to this point at the end of Section 7, but for now we will simply illustrate the phenomenon using a conventional clustering algorithm (FindClusters in Mathematica). It turns out that we can partition the top 10 coordinates associated with the 669 sample points into 18 clusters, most of which are dominated by a single element, which is therefore the mode of the cluster distribution.

Figure 18 is an illustration. The top row displays 6 sample points whose top 10 coordinates are included in the modes of 6 distinct clusters. The second row shows the projection of each sample point onto the Frobenius integral manifold. The first 5 columns represent the clusters with the largest modes, i.e., the clusters with the most number of repeated elements, in descending order. Notice that the projections of these 5 sample points onto the Frobenius integral manifold tend to form coherent and distinct geometric shapes. The last column represents the cluster to which the 10 coordinates in Table 5 would be assigned. This is a distinct cluster in the output of FindClusters, but it overlaps with the cluster in the second column, and the projections of the second and the last sample points onto the Frobenius integral manifold show a clear resemblance. The third row in Figure 18 displays the results of following the $\Theta$ coordinate curves from the principal axis for Prototype 09 towards the projected sample points, first, following the coordinate flow for the maximal infinitesimal eigenvector (see Table 3 and Figure 15), and then following the top 10 coordinate flows for the minimal infinitesimal eigenvectors. (Actually, we are displaying in the third row of Figure 18 the corresponding points on the ellipsoidal approximation.) We are thus looking at an 11-dimensional approximation of the 48-dimensional points from the Frobenius integral manifold.

For additional insight into how this coordinate system works, Figure 19 shows an expansion of the fifth column of Figure 18. The top row traces the path of the $\rho$ coordinate curve from the original 49-dimensional sample point to a point on the 48-dimensional Frobenius integral manifold. The bottom rows follow the $\Theta$ coordinate curves from the principal axis to a point on an 11-dimensional submanifold. In this example, the initial step in the direction of the maximal infinitesimal eigenvector is the largest, but this is not always the case. Note that the last step in Figure 19 is slightly different from the last step in Figure 18. This is because the coordinate curves in Figure 19 were computed on the Frobenius integral manifold itself, whereas in Figure 18 they were computed on the ellipsoidal approximation. In either case, though, if we continued to follow the $\Theta$ coordinate curves for an additional 37 steps, we would arrive at a point that is almost identical to the target point in the second row of Figure 18 or the first row of Figure 19. In this sense, we are justified in saying that we have found an optimal encoding for the sample points in Figure 18 on a low-dimensional nonlinear subspace of the original high-dimensional Euclidean space.

It may be instructive to compare Figures 14–16 and 18–19 to the results that would be obtained from Principal Component Analysis [Pea01]. Figure 20 shows the first six principal eigenvectors computed on the covariance matrix for the 32,000 point Data Sphere around Prototype 09. The first two eigenvectors show some similarities to the curves for Prototype 09 in Figures 14 and 15. But the subsequent eigenvectors are very different, and they do not seem to be related to any coherent global geometric shapes in the MNIST dataset.

7. Example: 7$\times$7 RGB Patches in the CIFAR-10 Dataset.

The second example that we will study is the CIFAR-10 dataset [Kri09], a classification task which is generally considered to be much harder than MNIST. The dataset was constructed from small, low-resolution, full-color photographs, in ten categories: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, truck. Each image is represented as a 32$\times$32 array of pixels in three color channels: RED, GREEN, BLUE. There are 50,000 images in the training set and 10,000 images in the test set. We performed global contrast normalization on the images in the training set, and then scaled the normalized image intensity down to a real number in the range [0, 1], but we did no other preprocessing, such as whitening [Kri09]. To facilitate a straightforward comparison to our results on the MNIST dataset, we defined a two-pixel border around each image in the training set and sampled the set of 7$\times$7 patches within the 28$\times$28 interior, drawing 12 samples per image. We thus started out with a sample that is comparable to the sample in Figure 11: 600,000 7$\times$7 patches. Our plan is to proceed along the same path as in Figure 11, but only as far as the lower-right corner, as we did with MNIST, to produce a 12-dimensional encoding of the patches. The main difference is that we are now working with three color channels instead of one monochrome channel.

We will also introduce two new constructs in this Section that were not part of our analysis of the MNIST dataset in Section 6: quotient manifolds and product manifolds.

Quotient manifolds are used to build invariance into the geometric model, consistent with known properties of the probabilistic model. For example, in our discussion of Figures 13 and 14 in Section 6, we noted that the mappings of Prototype 09 could be described by a shift transformation applied to a simple geometric shape, one or two pixels up or down, or left, or right. This makes sense, probabilistically. If a simple geometric shape is located at the mode of a probability distribution, then the horizontal and vertical translations of this shape will have approximately the same probability, and will also be located at a mode. But we know this, a priori, from the properties of images. Thus we do not need to compute the transformation from the data, as we did in Figures 13 and 14. We can build it into the geometric model from the start.

The classical technique to impose translation invariance on an image classification task is to use a Convolutional Neural Network [LBD*+*89] [LBBH98] [ZF14], since the convolution operation is equivariant with respect to translation. However, there are other possible techniques, based on more general transformation groups [ALR*+*16] [CW16]. Let $G$ be a group, and let $\{T_{g}\;|\;g\in G\}$ be a set of transformations on the space of images, parameterized by $g$. Suppose we hypothesize, for a particular choice of $G$, a particular choice of $T_{g}$, and a particular representation of the space of images, that the probability of $T_{g}(\textbf{x})$ is identical for all $g\in G$. We can then form a quotient space over the space of images, and work with a representative of this quotient in our probabilistic calculations. For the space of 7$\times$7 RGB patches in the CIFAR-10 dataset, the most useful transformations are: vertical translations, horizontal translations, and horizontal flips.

Figure 21 illustrates how this works. These are patches from a photograph of a white and gray cat sitting in front of a green vegetative background. Assume that our sampling procedure has selected the patch in the center of this grid, and we are looking at all the patches shifted zero-to-three pixels up and down, and zero-to-three pixels left and right. (At the edges of the 28$\times$28 images, we are using the discarded two pixel border and an additional one pixel border set to the mean value of the image intensities in each color channel.) We arrange the grid on a torus: The vertical shifts lie along the “vertical” axis of the torus, and the horizontal shifts are extended by a sequence of horizontal flips and lie along the “horizontal” axis of the torus. By our probabilistic hypothesis, every patch on this torus has equal probability, and thus the “toroidal grid” is the quotient space that we are looking for. How do we choose a representative of this quotient? What we want is an unambiguous procedure based on some salient property of the image. For the translations, we choose the patch that has the maximum variance in the three color channels, which is the patch outlined in blue in Figure 21. For the horizontal flips, we choose the orientation that maximizes the mean intensity on the left-hand side, which means that the patch outlined in blue in Figure 21 needs to be flipped horizontally. Note that our sampling procedure could initially select two distinct patches from a 28$\times$28 image and transform both of them into the same patch, but that is precisely what is required by our probabilistic model.

Product manifolds are used to combine low-dimensional solutions into a higher dimensional problem space, so that our dimensionality reduction techniques can be applied recursively. One example is shown in Figure 11 for the MNIST dataset with a single monochrome channel. Looking at the lower-right corner, assume that we have mapped the 7$\times$7 patches into a set of 12-dimensional submanifolds. Following the blue arrow, we assemble four adjacent 12-dimensional submanifolds into a product manifold and sample the images again using the 14$\times$14 patch. This means that each 14$\times$14 sample is mapped into four nonlinear submanifolds to produce a point in a 48-dimensional probability space. And, as indicated by the arrows pointing to the left in the middle of Figure 11, the task now is to apply the theory of differential similarity again to reduce the dimensionality of this space down to 12. Product manifolds — in particular, products of Grassmann manifolds — have been used successfully in recent research on the recognition of human actions. See, e.g., [Lui12]. They do not appear to have played a prominent role so far in work on the classification of static images.

We will not study the 14$\times$14 sample space in the present paper (but see Section 8 on “Future Work”). Instead, we will study a simpler example based on the RED/GREEN/BLUE color channels in the CIFAR-10 dataset. The dimensionality of a 7$\times$7 RGB patch is $3\times 49=147$, and we could either take this space as our starting point, or we could start with three separate 49-dimensional spaces. Let’s take the latter approach. We will apply the theory of differential similarity, separately, in each color channel, to produce three mappings into a set of 12-dimensional submanifolds, and then combine these into a set of 36-dimensional product manifolds. Several questions arise as we take the next step: Assume that we want to construct a lower-dimensional submanifold embedded in one of these product manifolds. How should we do this? If there are several distinct 12-dimensional submanifolds in each color channel, is it necessary to consider all possible combinations when we define the product space? How should we sample this space? We will suggest answers to these questions in the present paper, with the expectation that this will give us some insights into the analogous questions for the 14$\times$14 sample space.

For our initial selection of prototypes, however, we will start out with a set of data points in the full 147-dimensional space. Recall that our selection of prototypes for the MNIST dataset included a manual post-processing step. After some experimentation, we have developed a fully automated procedure for the CIFAR-10 dataset, with several adjustable parameters. There are two stages:

1. (1)

The first stage applies pgradascent from Figure 7 to the full dataset. We process the set of CIFAR-10 patches in five batches of 120,000 points each, generating 160 random choices of xstart in each batch, with $\beta=1/128$ (or $\sigma=8.0$) and with a 200 point Sample retrieved by the NearestFunction at each iteration. With these parameters, pgradascent converges very quickly to a mode very close to xstart. 2. (2)

The second stage computes and follows the integral curve of $\nabla U$, starting at one of the outputs from the first stage, to arrive at a new mode, ${\bf x}$, where $\nabla U({\bf x})=0$. In this computation, we set $\beta=1$ and we set $\mathtt{SamplePoints}$ to an 800 point random sample drawn from a 4,000 point data sphere around the starting point.

We also filter the output of each stage. The outputs from the first stage are sorted by the average of the variance in each color channel, then processed from max variance to min variance, and a mode is filtered out if it lies within the radius of a 4,000 point data sphere around one of its sorted predecessors. This step reduces the size of the output in each batch to: 58, 64, 59, 52, 59. The five batches are then merged, sorted, and filtered again, to produce a list of 94 modes as the input to the second stage. Similarly, the output from the second stage is filtered to remove any mode that lies within the radius of a 1,000 point data sphere around one of its sorted predecessors. The final result is a sorted list of exactly 50 prototypes.

Figure 22 shows a selection of these 50 prototypes, in three groups: high variance, medium variance, low variance. Seven prototypes of medium variance are shown in Figure 23, along with the two stages of their computation. Note that the value of xstart that led to the derivation of Prototype 14 is the patch outlined in blue in Figure 21, flipped horizontally. We will trace the analysis of these seven prototypes in the subsequent discussion, with particular attention to Prototypes 14 and 32.

Table 6 shows the size of an 8,000 point Data Sphere and a 2,000 point Coordinate Sphere for our seven prototypical clusters, and explicitly distinguishes between the 147-dimensional RGB space and the three 49-dimensional R/G/B spaces. In each case, we have computed the radius for 3X49D and then divided by $\sqrt{3}$ to get the radius for 49D, on the assumption that the spheres in each color channel should be the same size. Recall that the sizes of the Data Sphere and the Coordinate Sphere in the MNIST example were 32,000 and 8,000 points, respectively, and here we have chosen much smaller numbers, 8,000 and 2,000 points, respectively. We will also generally use sample sizes that are one-half the size of the samples in the MNIST example. This means that we will probably see clusters that are more coherent, and closer to a Gaussian, perhaps, but with more variability in the estimates. This is not to imply that these are better choices, just different choices. We should think of this analysis as one data point in an exploration of the parameter space for applications of the theory of differential similarity to various image processing tasks.

The next step is to compute the $\rho$ coordinate curves along the principal axis of each prototypical cluster. From this point on, as explained previously, we will do all computations separately in each color channel, although we will then combine the three R/G/B channels to generate full-color images. Figure 24 is analogous to Figure 14 in Section 6, showing the principal $\rho$ coordinate curves for Prototypes 14, 20, 23, 32, 34, 38, and 40. (Note: To locate the principal axes in each color channel, we did two global searches using NMinimize in Mathematica, with $\beta=1/128$ and with a 400 point sample defining $\nabla U({\bf x})$ for each search, followed by a local search using FindMinimum in Mathematica, which combined the two samples and set $\beta=1$.) The curves here are actually extended to the left, beyond the 2,000 point Coordinate Sphere where the principal axis is computed, to the 8,000 point Data Sphere. The color shifts along these curves are obvious, but there are also shifts in the shape of the images, as there were in the MNIST examples. For instance, for Prototype 14, reading from right to left, we see a counter clockwise rotation of the global geometric shape.

Figure 25 is analogous to Figure 15 in Section 6, for Prototypes 14 and 32. In each case, the first row shows the coordinate curve in the direction of the maximal infinitesimal eigenvector, to the Euclidean angle $\pi/2$, and the second row shows the coordinate curve in the opposite direction, to the Euclidean angle $-\pi/2$. (Note: To find a solution to the Euler-Lagrange equations in these examples, we set $\beta=1$ and used a single 400 point sample defining $\nabla U({\bf x})$.) The starting points, the images on the left, are points on the $\rho$ coordinate curves from Figure 24 at their intersections with the 2,000 point Coordinate Spheres. We could draw similar curves, of course, for the minimal infinitesimal eigenvectors,222For reasons that are not entirely clear, when we compute the eigenvectors of the matrix $\left(\,{g}_{i,j}({\bf x})\,\right)$ in Mathematica, the function Eigensystem sometimes returns a conjugate pair of complex eigenvalues. We ignore the complex eigenvectors when we compute the coordinate curves. Thus, for Prototype 14, GREEN, we have only 45 coordinate curves for the minimal infinitesimal eigenvectors, and for Prototype 32, RED and BLUE, we have only 43. This creates some minor problems when we try to measure the reconstruction error. See footnote 5, infra. analogous to Figure 16 in Section 6, but we first need to discuss our strategies for the computation of optimal lower-dimensional subspaces, a task to which we now turn.

In our work on optimal lower-dimensional subspaces for the MNIST example, in Section 6, we developed an ellipsoidal approximation to the $\Theta$ coordinates on the Frobenius integral manifold, and we defined a measure of the reconstruction error for this coordinate system in a higher dimensional space. We also generated a new sample (800 points, reduced to the 669 points in the positive hemisphere) to test these concepts. The ten coordinate curves in Table 5 did not fare very well, but we discovered that we could find an optimal set of ten coordinate curves for each individual sample point. Furthermore, these coordinate curves tended to fall into distinct clusters, which were illustrated in Figure 18. The main challenge in extending these ideas to the CIFAR-10 dataset is to generalize from a single monochrome channel to three R/G/B color channels, i.e., to apply our techniques to the product manifold.

Figure 26 for Prototype 14 and Figure 27 for Prototype 32 are analogous to Figure 18 in Section 6, showing the end results of our calculations. The six columns for each prototype represent six distinct $\Theta$ coordinate clusters. The first row in each figure shows six 147-dimensional sample points that are the nearest to the modes of each cluster, and the second row shows the projection of these points onto three 48-dimensional Frobenius integral manifolds. The third row shows three 11-dimensional R/G/B approximations, printed as full-color images. We will explain, step-by-step, how the sample points were selected and how the clusters were constructed, using as a running example a point from the upper left corner of Figure 26: Prototype 14, Cluster 01, Channel RED.

First, as promised in footnote 1, let’s look more closely at the ellipsoidal approximation. We start with an ellipsoid that roughly matches the Frobenius integral manifold along its semi-axes. For the principal axis, this is just the radius of the Coordinate Sphere, as shown in Table 6, which is 0.824038 for our running example. Let $a$ equal the reciprocal of the square of the length of this axis. For the coordinate curves aligned with the maximal infinitesimal eigenvectors, we saw in Table 3 and Figure 17 for MNIST that the distances along the axes in the positive and the negative directions were substantially different, and the same is true for CIFAR-10. In our running example, the matching semi-axis is 0.71319 in the positive direction and 0.589725 in the negative direction, so we split the space and use a different ellipsoid in each half. Let $b$ equal the reciprocal of the square of the length of the axis in each case. The equation for the coordinate flow, $\vec{\theta}^{\,1}_{t}({\bf x})$, along the surface of the ellipsoid can easily be generalized from the equations for the three-dimensional case in Section 5.1 of [McC18]. Let ${\bf x}=(x,y,z_{1},\ldots,z_{47})$. Then

[TABLE]

For the coordinate curves aligned with the minimal infinitesimal eigenvectors, the matching semi-axes are very close together, in the range 0.610 to 0.625 for our running example, so we simply use the mean of the positive and negative values. Let $c_{i}$ equal the reciprocal of the square of the mean length of axis $i$. Then

[TABLE]

The most important property of these equations is the fact that they are closed-form analytical solutions for the coordinate flows defined in (25), and thus the optimization algorithms in Mathematica can handle them efficiently.

Figure 28 shows the optimal coordinate flows for our running example: Prototype 14, Cluster 01, Channel RED. The goal is to approximate the coordinates of the target point on a 48-dimensional Frobenius integral manifold, which is the point in the RED channel of the image in the second row in Figure 26. The first three coordinates are: $\rho$, $\theta^{1}$, $\theta^{2}$. It is difficult to see this in Figure 28, but the $\theta^{2}$ coordinate curve is tilted slightly inwards, and there is a small gap between the point at the end of this curve and the manifold itself. As we proceed along the remaining 46 coordinate curves, the gap grows larger, and the coordinates get closer to the target point. Let’s analyze this approximation procedure, quantitatively. When we first compute the coordinate flows, the best we can do is to use a random ordering of the $\theta$ coordinates for the minimal infinitesimal eigenvectors: $z_{1},\ldots,z_{47}$. Once we have a solution, though, we can measure the Riemannian/Euclidean distance along each coordinate curve and identify the 10 coordinates for which this distance is the greatest. For our running example, these are:

[TABLE]

We can then compute the coordinate flows for the minimal infinitesimal eigenvectors again, starting with these 10 coordinates. The difference is apparent in the following table, where we are displaying the distance in the ambient Euclidean space between specific points on the coordinate curves and the target point on the manifold:

[TABLE]

The distance at the end of coordinate curve $\theta^{1}$ does not change, and the distance at the end of coordinate curve $\theta^{48}$ is almost the same, but there is a big difference at the end of coordinate curve $\theta^{11}$. Another way to present this result is to compute the total Riemannian/Euclidean distance along the various segments of the coordinate curves, as we did in the MNIST example:

[TABLE]

We can then say that the ratio between the length of the first 10 and the last 37 coordinate segments is 0.212804 with a random ordering and 0.830623 with an optimal ordering of the coordinates.

To continue, we need a new sample of RGB patches and the optimal coordinates for each one. In the MNIST example, we restricted this new sample to the positive hemisphere, that is, to the half-space that extends from the principal axis to the Euclidean angles $\pi/2$ and $-\pi/2$. Although we can compute $\Theta$ coordinate curves beyond these points, they become less accurate when the angles become much larger. If we were developing a production model for an image processing task, the best approach would be to find the antipodal point in the negative hemisphere and construct additional $\Theta$ coordinate curves from there, as we did for the curvilinear Gaussian example in [McC18], Section 5.2. For our present purposes, however, we can glean enough information from one hemisphere, and to do more would be duplicative. We therefore generate 1200 new RGB samples apiece for Prototypes 14 and 32, and restrict these to their positive hemispheres, in two ways: (i) we can either take the intersection of the positive hemispheres for the three R/G/B channels, or (ii) we can take the union, and it turns out to be advantageous to compute both and use each one for slightly different purposes. For Prototype 14, the intersection contains 481 points and the union contains 899 points. For Prototype 32, the intersection contains 185 points and the union contains 1065 points. (Note: This comparison shows, incidentally, how the principal axes are splayed out much more in Prototype 32 than in Prototype 14.)

Let’s now compute the $\Theta$ coordinate clusters, as we did in the MNIST example. The first fact we discover is that FindClusters in Mathematica does not work for the CIFAR-10 data. The clusters are completely “flat” with modes of at most two points. But it turns out, fortuitously, that we can use the theory of differential similarity to find clusters in the space of optimal coordinates. Consider a point with a list of 10 optimal coordinates, such as the example above from Prototype 14, Cluster 01, Channel RED. Write this as a 47-dimensional binary vector encoding the positions of the coordinates:

[TABLE]

Do the same with the GREEN and the BLUE channels, and construct the 141-dimensional product space. We can now define $\nabla U({\bf x})$ on this space, treating the RGB product as a real-valued vector space with the sample points situated at the corners of a 141-dimensional unit cube. To find the cluster modes, we follow a variant of the strategy that we developed above to construct the CIFAR-10 prototypes. To allocate RGB patches to clusters, we compute the Riemannian distance along the integral curve of $\nabla U$ towards each cluster mode, and select the mode that is the shortest distance away.

Within this framework, there are several alternatives. Let’s assume that we want the cluster modes to lie in the positive hemisphere in all three color channels. This means that we need to work with the intersection of our sample points, 481 points for Prototype 14 and 185 points for Prototype 32. Let’s assume that we want six cluster modes for each prototype. After some experimentation, it is easy to find a set of parameters that produces these results. For Prototype 14, we run pgradascent with 200 random choices out of 481 data points for xstart, with $\beta=1/8$, and with a 40 point Sample retrieved by the NearestFunction at each iteration. For Prototype 32, we run pgradascent with all 185 data points for xstart, with $\beta=1/8$, and with a 30 point Sample retrieved by the NearestFunction at each iteration. The filtering strategy is the same for both prototypes: The outputs from pgradascent are sorted by their estimated probability values, then processed from max probability to min probability, and a mode is filtered out if it lies within one-half of the radius of a 20 point data sphere around one of its sorted predecessors. Note that the final six modes will be sorted by probability, but they will be continuous-valued, not discrete. If we want binary-valued outputs, we can set the top ten coordinates in each color channel to 1 and the rest to 0. The highest probability mode for Prototype 14 is then:

[TABLE]

which corresponds to Cluster 01 in Figure 26. Notice, though, that this is not quite the same as our running example, since the coordinates in the RED channel are slightly different. Instead, our running example is the sample point that is the nearest to the mode for Cluster 01, and it has the following optimal coordinates:

[TABLE]

This is also the criterion for the selection of the other patches in the first row of Figures 26 and 27.

When we use the coordinates for Mode 01 to compute distances on the Frobenius integral manifold and its ellipsoidal approximation, they are slightly less than optimal. For example, in this table:

[TABLE]

the distance at the end of coordinate curve $\theta^{11}$ is 0.639119 instead of 0.608948, and in this table:

[TABLE]

the ratio between the length of the first 10 and the last 37 coordinate segments is 0.771451 instead of 0.830623. But this is exactly what we want from cluster modes. We can compute optimal coordinates, as we have seen, but that would give us a different subspace for each sample point, which is not very useful. What we want is a small number of clusters that provide reasonable approximations to these optimal calculations.

For the next step, allocating RGB patches to clusters, let’s work with the union of our sample points: 899 points for Prototype 14 and 1065 points for Prototype 32. Given any RGB patch, we need to compute the Riemannian distance of its optimal coordinates along the integral curve of $\nabla U$ towards each of the cluster modes. For Prototype 14, we set $\beta=1/16$ and we set $\mathtt{SamplePoints}$ to the 40 points that are the nearest to the cluster mode in the space of optimal coordinates. For Prototype 32, we set $\beta=1/16$ and we set $\mathtt{SamplePoints}$ to the 30 points that are the nearest to the cluster mode in the space of optimal coordinates. The patch is then assigned to the mode from which the Riemannian distance is the least. Note that these calculations are using distances and probabilities in the space of optimal coordinates, not in the original image space. The sizes of the resulting clusters are shown in Table 7.

Do these clusters, in fact, provide reasonable approximations to the optimal calculations in the original image space? Using standard statistical tests, we can show that, within each cluster, the Riemannian/Euclidean distance along the first 10 coordinates is greater (statistically significant in most cases333For Prototype 32, there are three exceptions with $0.1>p>0.001$ and one exception with $p>0.1$. Total cases: 36 $=$ 2 $\times$ 6 $\times$ R/G/B. at $p=0.001$) and the Riemannian/Euclidean distance along the remaining coordinates is less (statistically significant in most cases444For Prototype 14, there is one exception with $p>0.1$. For Prototype 32, there are six exceptions with $0.1>p>0.01$ and four exceptions with $p>0.1$. Total cases: 36 $=$ 2 $\times$ 6 $\times$ R/G/B. at $p=0.01$) than the corresponding distances for points outside the cluster. For an indication of the magnitude of these differences, see Table 8. The numbers in this table are the ratios between the RMS of the Riemannian/Euclidean distances along the first 10 coordinates and the RMS of the Riemannian/Euclidean distances along the remaining coordinates, as developed in the MNIST example. For the column labelled Optimal, the distances are computed using the optimal coordinates for every data point in every cluster, but separately for each color channel. This column thus sets the scale for the highest possible ratios in each case.555As explained in footnote 2, Prototype 14, GREEN, has only 45 coordinate curves for the minimal infinitesimal eigenvectors, and Prototype 32, RED and BLUE, has only 43. Thus the ratios in these cases are artificially high. In the next two columns, the distances are computed using the coordinates for each cluster, in two different ways. For the column labelled Within, the RMS values are computed for the data points within each cluster, using just the coordinates for that cluster mode. For the column labelled Across, the RMS values are computed for all the data points using the coordinates for all the cluster modes. Comparing the ratios in these two columns to the ratios for the optimal coordinates, there is a consistent pattern for the two prototypes and the three color channels. The evidence therefore supports the claim that these are reasonable approximations.

The ellipsoidal approximation played a central role in the preceding analysis, but it actually has two somewhat different roles to play in the theory of differential similarity. First, it could be a true approximation, because we are primarily interested in finding points on the Frobenius integral manifold that satisfy various conditions. Second, it could be a target for the mapping to a lower-dimensional subspace, regardless of its relationship to the manifold.

For an example of the first role, consider again our running example: Prototype 14, Cluster 01, Channel RED. Using the coordinates for Mode 01, we saw that the distance to the target point on the Frobenius integral manifold from the point on the ellipsoidal approximation at the end of coordinate curve $\theta^{11}$ was 0.639119. We would like to transfer this result to the manifold itself. The first step is to solve Equations (25), the integral equations for the coordinate flows on the manifold, one segment at a time, matching the flows on the ellipsoid. For example, if ${\bf x}_{0}$ is the initial point on the principal axis, we start by solving for $\vec{\theta}^{\,1}_{t_{1}}({\bf x}_{0})$ and then computing the value of the parameter, $t_{1}$, that minimizes the distance from $\vec{\theta}^{\,1}_{t_{1}}({\bf x}_{0})$ to the end of coordinate curve $\theta^{1}$, calling this the baseparameter. Continuing in the same manner, we compute the parameters for $\vec{\theta}^{\,2}_{t_{2}}({\bf x}_{1}),\vec{\theta}^{\,3}_{t_{3}}({\bf x}_{2}),\ldots,\vec{\theta}^{\,11}_{t_{11}}({\bf x}_{10})$, calling these the flowparameters. These first steps are packaged into a function called InitializeFlowParameters. We now have a sequence of coordinate points on the manifold, but the end result is actually further away from the target point, by a small amount: 0.673622. The next step is to apply a function called UpdateFlowParameters, and to do so iteratively. This function starts with solutions $\vec{\theta}^{\,1}_{t_{1}}({\bf x}_{0}),\vec{\theta}^{\,2}_{t_{2}}({\bf x}_{1}),\ldots,\vec{\theta}^{\,11}_{t_{11}}({\bf x}_{10})$, for a fixed sequence of starting points, ${\bf x}_{0},{\bf x}_{1},\ldots,{\bf x}_{10}$, and computes the baseparameter and the flowparameters for a modified formula that minimizes the distance to the target point, but takes the coordinate curves off the manifold. This is not a problem, however, because we can now use the new baseparameter and flowparameters to solve the integral equations again for the coordinate flows. After several iterations of UpdateFlowParameters, our algorithm converges to a point on the manifold at a distance of 0.442212 from the target point.666Conditions for the convergence of this algorithm are still open, and need to be investigated further. This is the point that produces the RED channel in Cluster 01, in the third row of Figure 26. Similar calculations produce the GREEN and the BLUE channels in Cluster 01, as well as the other images in the third row of Figures 26 and 27.

The reader has probably noticed that there are high frequency variations in the color patterns in the third row of Figures 26 and 27. The explanation is simple: We are looking at an 11-dimensional submanifold in each color channel, with the remaining coordinates truncated and their values set to 0.0. Since the point of truncation is essentially random, and different for each color, we see random color variations. It is natural to refer to this phenomenon as “truncation noise.” However, this is just a problem with the presentation of the submanifold when it is embedded in the original 147-dimensional space, and if we could smooth out the noise we would see that the global shape and the overall color balance of the RGB patches has been preserved. There are several ways to do this. For example, one cheap cosmetic trick is to replace the truncated coordinate values with the corresponding coordinate values from the sample point nearest to the cluster mode, that is, from the first row of Figures 26 and 27.

In the second role for the ellipsoidal approximation, however, we don’t care about the mapping back to the Frobenius integral manifold because we are only interested in the mapping forward. At this point, we have three 12-dimensional manifolds defined by a $\rho$, $\Theta$, coordinate system, and we want to apply the theory of differential similarity again to construct a lower-dimensional submanifold of the 36-dimensional product manifold. But to do this, we need to work with sample data points in a Euclidean space, for two reasons: (i) the fundamental equation for Brownian motion with a drift term, Equation (1), is only valid in a Euclidean space with a Cartesian coordinate system; and (ii) the mean shift algorithm, Equation (82) and Figure 7, only works in a Euclidean space with a Cartesian coordinate system. The ellipsoidal approximation gives us exactly what we want, in the simplest possible form. Although the coordinate flows are defined on the surface of an ellipsoid, they are expressed in Cartesian coordinates: $(x,y,z_{1},\ldots,z_{47})$, initially, but then truncated to $(x,y,z_{1},\ldots,z_{10})$. In addition, the product manifold is just the ordinary Cartesian product.

Let’s see how this works for the data points in the $\Theta$ coordinate clusters in Table 7. We start with a 36-dimensional Cartesian coordinate system and our goal is to reduce it to 12 dimensions. For each of the 2 $\times$ 6 clusters, we proceed through the same steps we have seen several times before: find a prototype, choose a coordinate sphere, locate the principal axis and define the $\rho$ coordinate curves, compute the maximal and minimal infinitesimal eigenvectors, compute the geodesics for the $\theta$ coordinate curves, define the $\theta$ coordinate flows, and compute the optimal ellipsoidal approximation for the coordinates on the Frobenius integral manifold that are determined by the minimal infinitesimal eigenvectors. The results shown in Table 9 are analogous to the Optimal results in Table 8: These are the ratios between the RMS of the Riemannian/Euclidean distances along the 10 optimal coordinates and the RMS of the Riemannian/Euclidean distances along the remaining coordinates.777These numbers are higher than the numbers in Table 8 because there are 36 coordinates in this case, in total, and we are computing the ratio between the optimal 10 coordinates for the minimal infinitesimal eigenvectors and the remaining 24 coordinates. See also footnote 5. We can rescale the denominator in these ratios by $37/24$ to get an approximate comparison. For example, in Prototype 14, Cluster 01, the ratio 1.26189 would be rescaled to 0.818523.

We can now process any arbitrary 7$\times$7 RGB patches through these coordinate systems. As an illustration, let’s look at two more random points from Prototype 14 and Prototype 32. See Figure 29, first column. (Note: For each prototype, we generated 100 additional random patches from the 8,000 point Data Sphere, filtered out all patches with a color channel in the negative hemisphere, and selected the one located at the greatest distance from the prototype.) We have seen columns two and three before, in other examples. For the fourth column, for each prototype, we have computed the Riemannian distance in the space of optimal coordinates to the modes of all six clusters, and selected the nearest one, which happens to be Mode 04 in each case. Finally, for the fifth column, we computed the coordinates of the data point in the optimal 12-dimensional ellipsoidal approximation for Cluster 04, as described above. Notice that there is a slight smoothing of the truncation noise in these final 12D images.

Figure 29 also illustrates how the theory of differential similarity can be implemented recursively. The base step starts with a random sample of data points and constructs an optimal 12-dimensional submanifold in each R/G/B color channel. The first recursive step for product manifolds starts with another random sample of data points, and sends it through the pipeline to construct an optimal 12-dimensional submanifold in each $\Theta$ coordinate cluster. The next step would repeat this process for the 14$\times$14 sample space, and so on. We will not pursue these subsequent steps in the present paper, but the reader should be able to anticipate roughly how they would operate.

8. Future Work.

The motivating principle in the theory of differential similarity is simple: We want to construct a geometric model that can be tuned to match the probabilistic model, in any given domain.

This is easy to do for the three-dimensional examples in [McC18]. The origin of the coordinate system is the mode of the probability density; the principal axis is a $\rho$ coordinate curve with a minimal Riemannian length, measured inwards from a fixed Euclidean distance; the Frobenius integral manifold is a surface of constant probability, at a fixed Riemannian distance from the origin, parameterized by a pair of $\Theta$ coordinate curves. We saw in the present paper that the geometric models for the higher-dimensional MNIST and CIFAR-10 examples have the same properties, but with more $\Theta$ coordinate curves. Most significantly, as shown in Table 2 and Figure 17, the data points on the Frobenius integral manifold for MNIST are all (approximately) the same Riemannian distance from the origin, even though their Euclidean distances vary substantially.

When we proceed to the construction of a lower-dimensional subspace, however, the three-dimensional cases and the higher-dimensional cases begin to diverge. It makes sense to measure the reconstruction error along the surface of the Frobenius integral manifold, following the coordinate flows given by (25). These are curves of constant probability, after all, and their Riemannian length is the same as their Euclidean length. Figure 5 shows that the data points from the curvilinear Gaussian potential are clustered around the geodesic coordinate curve that was generated by the minimal infinitesimal eigenvector, and the calculation of the reconstruction error confirms this phenomenon in quantitative terms. (See also Figures 14 and 15 in [McC18].) We observed a similar phenomenon in the analysis of the rotated Gaussian example in Section 5.1 of [McC18] (Figures 7, 8, 9). These examples are based on synthetic data generated from a known probability density, and thus there must be an explanation of these facts about the reconstruction error in purely mathematical terms. There is a conjecture at the end of Section 6 of [McC18] (see Figures 17 and 18), but the question needs further study.

In the MNIST and CIFAR-10 examples, by contrast, we are working with empirical distributions, and when we project these data points onto the Frobenius integral manifold, we see empirical clusters. This is shown in Figure 18 for MNIST, and in Figures 26 and 27 for CIFAR-10. We have seen that we can encode these clusters using small subsets of the $\Theta$ coordinates, which we have fixed arbitrarily at ten coordinates. Thus each cluster is represented by a distinct direction from the prototype in the high-dimensional image space. With this perspective, it is not surprising that the ten $\Theta$ coordinates in Table 5 did not yield a satisfactory reconstruction error. They are essentially random, since they are not pointing coherently towards any of the empirical clusters.

To compute a representation of the $\Theta$ coordinate clusters, it was necessary to develop an extension of the basic theory. First, we showed that it is possible to compute the optimal coordinates for each data point. All we need to do is to measure the Riemannian/Euclidean distance along each segment of the coordinate flows to the target point on the Frobenius integral manifold, and then choose the ten segments with the maximal length. Second, we showed how to search for clusters in the space of optimal coordinates. In the MNIST case, we could do this using FindClusters in Mathematica, but in the CIFAR-10 case, FindClusters did not work, and we computed clusters using a novel adaptation of the theory of differential similarity. The basic idea is to encode the optimal coordinates in a binary vector, and then define $\nabla U({\bf x})$ on a real-valued vector space, with the binary sample points situated at the corners of a unit cube in that space. It is remarkable that this idea works, since we are mapping a discrete search problem into a continuous space, using $\nabla U({\bf x})$ to compute the modes of the probability density in that space, which is a continuous computation, and then truncating the solutions back to a set of binary vectors. It would be interesting to investigate this idea further, in different but related contexts.

There were two more additions to the basic theory in the discussion of CIFAR-10 in Section 7: quotient manifolds and product manifolds. Quotient manifolds were used to build invariance into the geometric model, consistent with known properties of the probabilistic model. The basic idea is to form a quotient space over the space of images, and then work with a representative of this quotient in our probabilistic calculations. This approach should be studied further, and compared to similar approaches in papers such as [ALR*+*16] [CW16]. Product manifolds were used to combine low-dimensional solutions into a higher dimensional problem space, so that our dimensionality reduction techniques could be applied recursively. The main example in Section 7 was the 147-dimensional RGB manifold constructed out of three 49-dimensional R/G/B manifolds. It is important to understand how our dimensionality reduction techniques, such as the search for $\Theta$ coordinate clusters, as discussed in the previous paragraph, are configured to work with product manifolds. Recall that we computed the $\Theta$ coordinate clusters and the Riemannian distances within them in a 141-dimensional product space constructed from the optimal coordinates for the RED, GREEN, BLUE, channels in combination, not separately. If we had instead computed $6$ clusters in each channel, we would have to consider $6^{3}$ clusters in the product manifold, and perhaps estimate the probabilities of each combination in order to decide which ones to retain and which ones to discard. The point of any dimensionality reduction technique is that we should expect to lose information, but we want (i) the simplest coding scheme that (ii) preserves as much information as possible. For these reasons, in a coding scheme that replaces optimal coordinates with cluster modes, it is better to compute the cluster modes in the product manifold than in the separate factor manifolds.

The same strategy can be applied to the 14$\times$14 sample space in Figure 11. If we computed, say, 40 prototypes in the 7$\times$7 sample space, we would have to consider $40^{4}=2,560,000$ prototypical clusters in the product manifold. Instead, we can compute 40, 50, 60, $\ldots$, prototypes directly in the 2$\times$2 product space. We can then apply the same strategies that we used for the R/G/B product manifolds in the CIFAR-10 example to construct a small number of $\Theta$ coordinate clusters for each prototype. Moreover, any prototype from the 7$\times$7 sample space that does not appear in the 2$\times$2 product space can be pruned from our construction, and this process of proliferating and pruning product manifolds can be continued, recursively, up through the hierarchy in Figure 11 and similar architectures. Our initial results on the MNIST dataset show that this idea works, but we cannot yet say how well it works. There are several hyperparameters that have to be set: the constant $\beta$, the size of $\mathtt{SamplePoints}$, the size of the Data Sphere and the Coordinate Sphere, etc., and these values can vary as we move through the architecture. The architecture itself is a variable: the size of the patches, the configuration of the product manifolds, the number of prototypes, the number of optimal $\Theta$ coordinates, etc., all have an impact on the outcomes. These are empirical questions, which are currently under investigation. When we have fully analyzed the experimental results, we will publish them in a paper with the working title: “Deep Learning with a Riemannian Dissimilarity Metric.”

The MNIST dataset has often been used as a testbed for experiments on unsupervised learning. Much of the recent work on this problem has been based on deep generative models implemented with neural networks, in particular, either Variational Auto-Encoders (VAEs) [KW14] [RMW14], Generative Adversarial Networks (GANs) [GPAM*+*14], or a combination of the two [MSJG16] [MNG17]. In the VAE framework, for example, if $\mathbf{z}$ is a latent variable (also called a “code”) and $\mathbf{x}$ is a data point, then the generator (also called the “decoder”) is the probability $p(\mathbf{x}|\mathbf{z};\theta)$ with prior $p(\mathbf{z};\theta)$, which is usually parameterized by a neural network. The true posterior, $p(\mathbf{z}|\mathbf{x};\theta)$, is assumed to be intractable, and so it is approximated by the recognizer $q(\mathbf{z}|\mathbf{x};\phi)$ (also called the “encoder”), which is itself parameterized by a neural network. The parameters $\phi$ and $\theta$ are estimated by maximizing a lower bound on the log-likelihood of the data, a process called variational inference. If we can learn these conditional probabilities, then there are several ways to build a semi-supervised classifier. For an example of how this works on the MNIST dataset with only 100 labelled examples, see [KRMW14]. The pure GAN framework does not have a recognizer, but the generator, $p(\mathbf{x}|\mathbf{z};\theta)$, is trained with an auxiliary discriminative network acting as an adversary. One way to use a GAN for unsupervised learning is to add a term to the loss function for the auxiliary discriminative network that maximizes the mutual information between the data distribution and the predicted class distribution, in an attempt to disentangle the code for the classes. See, e.g., [Spr16] [CDH*+*16]. Although deep generative models have achieved good results for unsupervised learning on MNIST, they do not seem to work as well on natural images. On CIFAR-10, for example, the best current results for unsupervised learning are instead based on a variety of ad hoc methods: [HPP*+*20] [VGVG*+*20] [PHK*+*21]. For a general critique of almost all recent work on the use of deep generative models to disentangle the factors of variation in a dataset of images, without supervision, see [LBL*+*19].

In Section 1, supra, we mentioned an earlier class of auto-encoders, known generically as Regularized Auto-Encoders. In particular, we noted that Denoising Auto-Encoders (DAEs) [Vin11] and a specialized form of Contractive Auto-Encoders (CAEs) [AB14] can be shown to compute the gradient of the log of the input probability density [BCV13]. In other words, these are network architectures that can be trained to compute the vector field $\nabla U({\bf x})$. This correspondence suggests several questions for further investigation:

1. (1)

The earliest papers on Denoising Auto-Encoders, [VLBM08] and [VLL*+*10], offer an intuitive justification of the DAE algorithm in terms of the manifold hypothesis. Consider Figure 2 in these papers, and the following informal explanation:

During denoising training, we learn a stochastic operator $p(X|\tilde{X})$ that maps a corrupted $\tilde{X}$ back to its uncorrupted $X$ …. Corrupted examples are much more likely to be outside and farther from the manifold than the uncorrupted ones. Thus stochastic operator $p(X|\tilde{X})$ learns a map that tends to go from lower probability points $\tilde{X}$ to nearby high probability points $X$, on or near the manifold. [VLL*+*10], at 3380.

Similar justifications are offered in [RBDV12] to explain how Contractive Auto-Encoders tend to generate samples along a data manifold. Can these explanations be formalized within the theory of differential similarity, using Brownian motion with the drift term $\nabla U(\mathbf{x})\cdot\nabla$ as a model? See Section 2.3 in [McC18], and see the derivations in Section 6 of [McC18]. 2. (2)

In another theoretical perspective on Denoising Auto-Encoders, [VLBM08] shows how to derive the DAE training criterion by maximizing a variational bound on a generative model that includes the corruption of the data: $X\rightarrow\tilde{X}$. This paper thus anticipates some of the later work on Variational Auto-Encoders [KW14] [RMW14]. The theory of differential similarity does not use a neural network, of course, but we are computing in our geometric model a low-dimensional code $\mathbf{z}=r(\mathbf{x})$ for the data point $\mathbf{x}$, with some probability density $q(\mathbf{z}|\mathbf{x})$. It would be interesting to compare these architectures, layer by layer. 3. (3)

Vincent reported in [Vin11] that there was an “unsuspected link” between the training of a Denoising Auto-Encoder and the score matching technique [Hyv05] [Hyv07] for learning the parameters of an unnormalized density model over continuous-valued data. This result was later generalized in [AB14] to a specialized form of Contractive Auto-Encoder. What Hyvärinen refers to here as the “score” is just the vector field $\nabla U({\bf x})$. Using his score matching technique, Hyvärinen avoids estimating $\nabla U({\bf x})$ directly, but we have taken a different approach within the theory of differential similarity. We use the Frobenius integral manifold as a model of the data distribution, and this gives us an interesting compromise between traditional parametric and traditional nonparametric statistical estimation. Note that the Frobenius integral manifold can represent almost any probability density, subject to mild regularity conditions, but its shape is strictly determined once we specify an orthonormal set of infinitesimal eigenvectors at a point along the principal axis. This appears to be a novel statistical estimation technique, but it needs to be studied more carefully. 4. (4)

We have used very small sample sizes in our work on the MNIST and CIFAR-10 datasets, as compared to the sample sizes needed to train neural networks. See, for example, the size of the $\Theta$ coordinate clusters in Table 7. However, empirically, these small samples seem to yield consistent and stable estimates for the $\rho$ and $\Theta$ coordinate curves, and it is important to understand, theoretically, why this is the case. A significant step in this direction is the work of Arias-Castro, Mason and Pelletier [ACMP16] on the consistency and stability of the mean shift algorithm, which answers some of our questions about the $\rho$ coordinate curves. Similar results for the $\Theta$ coordinate curves would be very useful.

Note that the $\rho$, $\Theta$, coordinate system imposes a strong inductive bias on both the model and the data, as advocated in [LBL*+*19], and this bias can be manipulated systematically by choosing different quotient manifolds and product manifolds in different domains.

There are two popular dimensionality reduction algorithms that compute and optimize probability distributions on low-dimensional manifolds: (i) SNE or t-SNE, and (ii) UMAP. Since these algorithms are used primarily for data visualization, the embedding spaces are usually only two- or three-dimensional. Stochastic Neighbor Embedding (SNE) [HR03] centers a Gaussian at each data point $\mathbf{x}_{i}$ in the original high dimensional space, and also at each data point $\mathbf{y}_{i}$ in the low dimensional space, and uses the conditional probabilities derived from these Gaussians to measure the similarities of $\mathbf{x}_{i}$ and $\mathbf{y}_{i}$ to each of their $k$ neighbors. The goal is to adjust the mapping $\{\mathbf{x}_{i}\}\mapsto\{\mathbf{y}_{i}\}$ so that the mismatch between the probability distributions in the high and low dimensional spaces is as small as possible, and this is achieved by minimizing the sums of the Kullback-Leibler divergences over all data points. An improvement of this algorithm, called t-SNE [VdMH08], replaces the Gaussian in the low dimensional space by a Student t-distribution with one degree of freedom, because it has much heavier tails. This is all somewhat ad hoc, and the Uniform Manifold Approximation and Projection (UMAP) algorithm [MHM20] achieves similar and arguably better results with a rigorous mathematical foundation. The basic idea of UMAP is to rescale the metric at each point $\mathbf{x}_{i}$ in the high dimensional space, in order to approximate a uniform probability density in that space, and then to glue these incompatible metric spaces together, using a clever construction from algebraic topology and category theory. The optimization algorithm can then minimize the cross-entropy between the topological representation of the data in the high dimensional space and the corresponding topological representation in its low dimensional projection. Here is a comment by the authors on their methodology:

…[W]e feel that strong theory and mathematically justified algorithmic decisions are of particular importance in the field of unsupervised learning. This is, at least partially, due to a plethora of proposed objective functions within the area. …UMAP’s design decisions were all grounded in a solid theoretic foundation and not derived through experimentation with any particular task focused objective function. [MHM20] at 2–3.

Although the goals are quite different, the theory of differential similarity also has a solid mathematical foundation, drawn from the fields of stochastic processes and differential geometry. See [McC14, McC18].

There is one operation that appears repeatedly in the diagram in Figure 11: We construct a product manifold consisting of four prototypical clusters, and then construct a submanifold which is itself a prototypical cluster. We can now take a big step: We can use this construction to define the semantics of an atomic formula in a logical language, that is, to define a predicate with four arguments. The general idea is to replace the standard semantics of classical logic, based on sets and their elements, with a semantics based on manifolds and their points. For a detailed discussion of how this works, see [McC15]. (These ideas emerge naturally from earlier work on the applications of AI to Law.) The natural setting for these developments is a logical language based on category theory, or what is known as a categorical logic. Thus, in [McC15], we analyze a categorical logic based on the category of differential manifolds (Man), which is weaker than a logic based on the category of sets (Set) or the category of topological spaces (Top). Let’s call this a manifold logic. In a manifold logic, as explained above, the prototypical clusters provide an interpretation of the atomic formulas, and the proof theory extends the differential manifold structure throughout the entire language. A technical paper addressed to the computational logic community on these ideas is currently in preparation.