unsupervised machine learning for physical concepts

Physics is built on many different concepts, such as force, entropy, and Hamiltonian. Scientists derive meaningful concepts from observations by their ingenuity and then use formulas to connect them, constituting the theory. As a principle, this theory is always as simple as possible. Since the beginning of the 21st century, machine learning has developed rapidly Liu et al. (2018) and has been widely used in various fields Teichert et al. (2019); Bourilkov (2019); Mishra et al. (2016), including machine-assisted scientific research Vamathevan et al. (2019); Valletta et al. (2017); Schmidt et al. (2019); Carleo et al. (2019). A natural question arises: can machines propose scientific theories by themselves? It is undoubtedly a fundamental and challenging goal. Many works have studied this from different aspects Ha and Jeong (2021); Zhang and Lin (2018); Boullé et al. (2021); Zobeiry and Humfeld (2019); Farina et al. (2020); Wu and Tegmark (2019); Zheng et al. (2018); D’Agnolo and Wulzer (2019). In general, the establishment of the theory can be divided into two steps: identify the critical variables from observations and connect them by formula. A technique known as symbolic regression has been developed for the second step Udrescu and Tegmark (2020); Udrescu et al. (2020). The authors propose a network named AI Feymann to automatically find the formula the system obeys. To increase the success rate of symbolic regression, the critical variables identified by the first step should be as few as possible.

In this work, we focus on the first step of establishing a theory. We suggest using both Topological Data Analysis(TDA)Wasserman (2018); chazal2017introduction; Murugan and Robertson (2019) and variational autoencoder(VAE) Kingma and Welling (2019, 2013); An and Cho (2015); Khobahi and Soltanalian (2019); Rezende et al. (2014); Higgins et al. (2016); Burgess et al. (2018) to extract meaningful variables from extensive observations. TDA is a recent and fast-growing tool to analyze and exploit the complex topological and geometric structures underlying data. This tool is necessary for our protocol if specific structures such as circles and sphere are in experimental data. VAE is a deep learning technique for learning latent representations. This technique has been widely used in many problemsLuchnikov et al. (2019); Cerri et al. (2019) as a generative model. It can also be seen as a manifold learning network for dimensionality reduction and unsupervised learning tasks.

Our protocol has two stages. Firstly, we use TDA to infer relevant topological features for experimental data. In the simplest case, where the manifold has a low dimension, an essential feature for us is the Betti numbers, topology invariants. Once we get the topological features, we can design the proper architecture and loss function of VAE. As we will show later, the latent variables of VAE need to form a manifold homeomorphic to the manifold composed of observations. After the training of the VAE network, the latent variables represent the key variables discovered by this machine. Thanks to the structure of the VAE network, the latent variables, and the observations are in one-to-one correspondence. Using the trained VAE network, one can calculate the latent variables and the observations from each other. That means the formula derived by symbolic regression connecting the latent variables can predict the experimental phenomena.

We test our protocol on three toy models. They have different topological features. The first is a classical coupled harmonic oscillator, where the observations constitute a circle embedded in three-dimensional Euclidean space. Another example is two balls rotating around the same center, with different radius and angular velocities. With the other ball as the reference system, the observations are the Cartesian coordinates of the ball. The observations constitute a lemniscate curve. The third is a two-level system, and the observations are the expected value of some physical quantity,i.e., some hermitian matrices. The observations constitute a sphere.

This paper is organized as follows. In section II, we describe the works been done before and the problem we want to solve in more detail. In section III, we introduce the architecture of neuron networks and argue why we need the manifold of latent variants should be homeomorphic to that of observations. In section IV, we show the performance of this protocol on three toy models. We compare the observations and latent variables to show the relation between them.

Data-driven scientific discoveries are not new ideas. It follows from the revolutionary work of Johannes Kepler and Sir Isaac Newton Roscher et al. (2020). Unlike the 17th century, we now have higher quality data and better calculation tools. People have done a lot of research on how to make machines help people make scientific discoveries in different contexts Carleo et al. (2019); Cichos et al. (2020); Karpatne et al. (2018); Wetzel (2017); Rodriguez-Nieva and Scheurer (2019); Wu and Tegmark (2019). In the early days, people paid more attention to the symbolic regression Kim et al. (2020); Udrescu and Tegmark (2020); Udrescu et al. (2020); Lample and Charton (2019). Another challenging direction is to let the machine design experiments. In  Melnikov et al. (2018); Krenn et al. (2016); Fösel et al. (2018); Bukov et al. (2018), authors designed automated search techniques and a reinforcement-learning-based network to generate new experimental setups. In condensed physics, machine learning has been used to characterize phase transitions Van Nieuwenburg et al. (2017); Uvarov et al. (2020); Carrasquilla and Melko (2017); Wang (2016); Ch’Ng et al. (2017); Yu and Deng (2021).

In recent years, some works have contributed to letting the machine search for key physical concepts. They used the different networks to extract key physical parameters from the experimental data Iten et al. (2020); Zheng et al. (2018); Lu et al. (2020). In Zheng et al. (2018), authors propose an unsupervised method to extract physical parameters by interaction networks Battaglia et al. (2016). Another helpful tool is the VAE network, which has been widely used for similar goal Lu et al. (2020); Iten et al. (2020).

VAE network is a powerful tool, by which one can minimize the numbers of extracted variables by making the variables independent of each other  Higgins et al. (2016); Burgess et al. (2018). To do this, one can choose the prior distribution $P(z)$ as latent variables $z$. However, this method fails when the manifold of observations have some topological features. This paper aims to solve this problem.

Here we describe this goal in more detail. Suppose we have an experimental system $\mathcal{S}$. In this system $\mathcal{S}$, some physical variables change as the experiment progresses. We use $\mathcal{P}$ to denote the set consists of the value of these physical variables. We use $\mathcal{P}^{(k)}$ to represent the value of these variables in time $k$. From the system $\mathcal{S}$ we can derive an experimental data set $\mathcal{E}$. Every data point is a vector $\mathcal{E}^{(k)}\in\mathcal{E}$ where $\mathcal{E}^{(k)}$ denotes the data point belongs to time $k$. From the perspective of physics, the change of experimental data must be attributed to the change of key physical variables. Therefore there is a function $f:\mathcal{P}\rightarrow\mathcal{E}$ such that $f(\mathcal{E}^{(k)})=\mathcal{P}^{(k)}$ for any $k$. In this paper we aim to find the function $\tilde{f}$ and $\tilde{f}^{-1}$ such that $\tilde{\mathcal{P}}^{(k)}=\tilde{f}^{-1}(\mathcal{S}^{(k)})$ where $\tilde{\mathcal{P}}^{(k)}$ constitute a set $\tilde{\mathcal{P}}$. We call $\tilde{\mathcal{P}}^{(k)}$ the effective physical variables. We remark that $\tilde{\mathcal{P}}^{(k)}$ is not necessary to equal to $\mathcal{P}^{(k)}$, but their dimensions should be the same. Effective physical variables are enough to describe the experimental system. In fact, one can redefine the physical quantity can get a theory that looks very different but the predicted results are completely consistent with existing theories.

The problem arise if we try to use neuron network to find proper function $\tilde{f}$ and $\tilde{f}^{-1}$. All the functions neuron network can simulate is continuousLABEL:, so $\tilde{f}$ and $\tilde{f}^{-1}$ must be continuous. In a real physical system, $f$ and $f^{-1}$ don’t have to be like this. Therefore, in some cases, we can never find a $\tilde{\mathcal{P}}^{(k)}$ which has the same dimension with $\mathcal{P}^{(k)}$. For example, suppose we have a ball rotating around a center. The observable data is the location of the ball, denoted by Cartesian coordinates $\mathcal{E}^{(k)}=(x^{(k)},y^{(k)})$, which form a circle, and the simplest physical variable is the angle $\mathcal{P}^{(k)}=\theta^{(k)}\in[0,2\pi)$. In this case, $f:\theta\rightarrow(x,y)$ is continuous while $f^{-1}$ is not, which cannot be approximated by neural network. This is because $f$ is a periodic function of $\theta$. More generally, as long as the manifold composed of $\mathcal{E}$ has holes with any dimension, this problem arises. We call these physical variables topological physical variables(TPVs). Back to the last example, we can find the Betti numbers of the ball’s location is $[1,1,0]$, where the second $1$ means that it has one TPVs. Due to this reason, we suggest using TDA to identify the topological features firstly. After knowing the numbers of TPVs, we can design proper latent variables and the corresponding loss function $\mathcal{L}$. For the case we have two PPVs, we need two latent variables, named $x$ and $y$, and we add the topological term $|x^{2}+y^{2}-1|$ to $\mathcal{L}$ to restrict them.

In at least two cases, the manifold will have holes. One is the $\mathcal{E}$ of a conversed system,e.g., the classical harmonic oscillator. Another situation stems from the limits of the physical quantity itself, such as the single-qubit state, which forms a sphere. In addition to these two categories, sometimes the choice of observations and references also affects their topological properties.

As shown in figure  2, This network consists of two parts. The left is called encoder, and the right is called decoder. The encoder part is used to encode high-dimensional data into low-dimensional space, while the decoder part is used to recover the data, i.e., the input from low-dimensional data. The latent variables are low-dimensional data, and we hope they contain all the information needed to recover the input data. So that we require the output $\hat{x}$ is as close as possible to input $x$. Therefore, we can train this neuron network without supervision.

In general, the encoder and decoder can simulate arbitrary continuous function $f$ and $f^{-1}$. That imposes restrictions on latent variables:

This limitation means that the Betti numbers should keep the same. Suppose the inputs $x$ form a circle, as shown in 1, the easiest quantity to describe these data is the angle. According to 1, this neuron network can never find such a quantity because the angle will form a line segment that is not homeomorphic to a circle. On the other hand, VAE network can’t even find the Cartesian coordinate $\{x_{1},x_{2}\}$ because $x_{1}$ and $x_{2}$ are not irrelevant. To handle this case, we suggest analyzing the Betti numbers as a priori knowledge of input $x$ and tell the neuron network.

The general loss function is wirtten as Higgins et al. (2016):

Here $\theta,\phi$ are the parameters of decoder and encoder networks, respectively. $p(z)$ is the prior distribution about the latent variables $z$. $q_{\phi}(z|x)$ and $p_{\theta}(x|z)$ denote the conditional probability. $D_{KL}$ means the KL divergence of two distributions.

To calculate the loss function for a given output $\hat{x}$, one needs to set up a prior $p(z)$, and parameterize the distribution $p_{\theta}$ and $q_{\phi}$. In traditional VAE, the distributions $p(z)$, $p_{\theta}$ and $q_{\phi}$ are often supposed to be multidimensional Gaussian distribution. The mean of $p(z)$ is zero the and variance is $1$. The mean and variance of $p_{\theta}$ and $q_{\phi}$ is the output of encoder network.

In our cases, we choose the prior $p(\vec{z})$ as $p(\vec{z})=p(\vec{z_{g}})\times p(\vec{z_{t}})$. Where $\vec{z_{g}}$ denotes the GPVs, and $\vec{z_{t}}$ denotes the TPVs. With the same as traditional VAE, in order to minimize the number of GPVs, a convenient way is to choose $p(\vec{z_{g}})$ as Gaussion distribution $p(\vec{z_{g}})=N(0,I)$. According to the topological features, we choose a proper topological term $\mathcal{T}$. We choose the $p(\vec{z_{t}})$ as $p(\vec{z_{t}})=Ae^{-T}$. Here $A$ is a constant. The variance can be asorpted into the hyperparameters in loss function.

For the conditional distribution we choose as s $q_{\phi}(\mathbf{z}\mid\mathbf{x})=N(\vec{z},I)$, with constantly variance and means determined by the encoder network. $p_{\theta}(\mathbf{x}\mid\mathbf{z})$ has the same expression $p_{\theta}(\mathbf{x}\mid\mathbf{z})=N(x,I)$. In this assumption, the KL divergence can be written as

Thus the loss function can be written as a simplier expression:

where $\alpha,\beta,\gamma$ are hyperparameters and $\mathcal{T}$ is the topological term depending on the Betti numbers. For example, when the Betti numbers are $\{1,0,1\}$, $\mathcal{T}$ can be written as $\mathcal{T}=|z_{1}^{2}+z_{2}^{2}+z_{3}^{2}-1|$. When the Betti numbers are $\{1,2,0\}$, $\mathcal{T}$ can be written as $\mathcal{T}=|(z_{1}^{2}+z_{2}^{2})-2(z_{1}^{2}-z_{2}^{2})|$, which is known as a Lemniscate. Here we use $\vec{z}$ to denote the latent variables. $P(\vec{z_{g}})$ is the distribution of $\vec{z_{g}}$. We need to choose some latent variables as TPVs, and others as general physical variables(GPVs). In general, we know how many TPVs we need from the Betti numbers, but we usually don’t know how many GPVs we need. One solution is to set up as many GPVs as possible, and the redundant GPVs will be zero.

We test this neuron network by numerical simulation. Three toy models whose manifolds are different are considered here.

We use pytorch to implement neural networks. Our code can be found here. We use the same structure in the encoding and decoding network. They have two hidden layers, each with 20 neurons. In encoder network, we choose Tanh as the active function, while in decoder network, we choose ReLU as the active function. In numerical simulation, we first generate a database of $1000\times m$, where m is the dimension of the observations. We then use the TDA tools to analyze the Betti numbers of the manifold constituted by observations and set up the latent neurons. We used the Adam optimizer and set the learning rate to $0.0001$. At each training session, we randomly select 100 sets of data from the database as a batch, then calculate the average of its loss and reverse propagate it. All training can be done on the desktop in less than 6 hours.

In this work, we discuss the defect of the traditional VAE network and propose a simple solution. We can extract the minimum effective physical variables from the experimental data with the improved method by classifying the latent variables. We test our approach on three models. They represent three different situations that may arise that traditional VAE cannot handle. Some of them come from experimental restrictions, and some come from physical laws. We think the latter situation is more essential. However, in the more complex case, the Betti numbers are not the only useful topological features. Two manifolds may have the same Betti numbers but not Homeomorphic. In such a case, more topological features are needed to design reasonable restrictions on latent variables. When the manifold dimension is higher, it may be hard for TDA to calculate the Betti numbers. One efficient way is to firstly reduce the dimensions by traditional autoencoder, and then calculate the Betti numbers of latent variables. Another important question is how to ensure that the relation between meaningful variables and observations is simple. We leave it for the future work.