Modified independent component analysis for extracting Eigen-Modes of a quantum system

Dynamical probes are popular experimental tools for studying quantum systems nowadays. Though a well-defined eigen-mode will manifest itself as a single frequency oscillation in the dynamical probe, the real experimental data can be quite messy because a probe usually excites multiple modes which are also mixed together with random white noise. Moreover, the oscillations usually quickly damp out due to finite decoherence time. These make it difficult to extract the frequencies of these oscillation measurement data, that is, the eigen-energies of these eigen-modes. Here we develop an unsupervised machine learning algorithm to solve this problem. Our method is inspired by the independent component analysis method and its application to the ‘cocktail party problem’, where the goal is to recover each voice from detectors that detect signals of many mixed voices. We demonstrate the advantage of our method by an example of analyzing the collective mode of a Bose–Einstein condensate of atomic gases. We believe that this method can find broad applications in analyzing data of dynamical experiments in quantum systems of different fields.


Introduction
In quantum systems, eigen-modes have discrete eigen-energies, and each eigen-mode with a fixed eigen-energy manifests itself as a single-frequency oscillation in a dynamical probe. The frequency of the oscillation, up to an ℏ, is the eigen-energy of the eigen-mode. Therefore it is quite a common task in many quantum physics experiments to extract a single-frequency oscillation from dynamical measurements.
If the measurement can be performed for a sufficiently long time, the frequency can be determined by the Fourier transformation. However, in many cases it is not possible to perform a long time measurement, either because the oscillation quickly damps out because of decoherence, or because the lifetime of the quantum system itself is limited. Hence, the frequency resolution of the Fourier transform method is always limited by 2π/T, where T is the total duration for taking data. An alternative way commonly used in practice is to fit the dynamics with a damped harmonic oscillator. Nevertheless, it is quite often that the probe excites multiple modes with different frequencies, and some times the oscillations are also embedded in noisy signals, therefore the fitting becomes not so reliable. Thus, when the measurement time is short, and the data contains multiple frequencies and is noisy, the task of extracting single-frequency oscillation becomes quite challenging. The purpose of this work is to design a machine learning algorithm to solve this problem.
To start with, we point out that this problem bears a lot of similarities with the famous 'cocktail party problem' that has been widely studied in neuroscience and computer science for decades, and the independent component analysis (ICA) is a quite efficient algorithm to solve the cocktail party problem. Let us first briefly review the ICA method [1] and its application to the classical 'cocktail party problem' [2].
The 'cocktail party problem' considers N-number of sources of voices, each of which is a sequence of data denoted by s i (t), (i = 1, . . . , N) and M-number of detectors with M > N, each of which obtains a signal x j (t), (j = 1, . . . , M) as a linear combination of all s i (t). That is to say, there is a M × N matrix A such that The matrix A is unknown and is assumed to be independent of t, and the goal is to find out A −1 such that one can recover each voice s i (t) from the signals obtained by the detectors {x j (t), (j = 1, . . . , M)}. The aforementioned problem in a quantum system is actually the same as the 'cocktail party problem' . Each single frequency oscillation from a given eigen-mode can be viewed as a source s i (t), and several different measurements x j (t) obtain several different superpositions of s i (t). To extract the frequency or eigen-energy of each mode, one first needs to recover s i (t) from the measurements {x j (t)}. What lies behind the ICA application of this problem is the central limit theorem [3]. For each signal, when we perform statistics over a certain duration of t, we can obtain the distribution for each signal, which always possesses certain features. The central limit theorem says that when one adds up many such signals together, the distribution of x j (t) will approach a Gaussian distribution and looks like a noise. Thus, if the inverting the transformation does not yield a single independent signal, in other words, if the inverting transformation still yields a summation of different signals, the distribution is still a Gaussian distribution. Therefore, in order for the inverted signal to be an independent signal, the distribution should be far away from Gaussian distribution as possible. The signal with Gaussian distribution displays maximum entropy for given mean and variance. Thus, the working principle of ICA is to find out A −1 such that the distribution of s i (t) = (A −1 ) ij x j (t) deviates from a Gaussian distribution as far as possible, or in a more quantitative description, that the entropy of s i (t) = ∑ j (A −1 ) ij x j (t) is as small as possible.
To be more precise, we first normalize the data x j (t) such that its mean valuex j (t) = 0 and its covariance matrix X † X = I, where the statistics is performed over a sufficiently long duration of t. We further require A † A = I, which ensures the mean value of s i (t) is also zero and its variance also equals one. By performing statistics over t, for a given sequence s(t), one can obtain its distribution P [s, u], where u is the range of s(t). Then, we can compute the entropy of a given sequence s(t) as If s Gau is a sequence obeying a Gaussian distribution with zero mean and unity variance, is the entropy maximum for sequences with same mean and variance. Hence, the ICA method is to find out is maximized. In practices, since it is hard to directly compute P[s, u], several formulae have been proposed to approximate H[s] [4].

Results
When we apply the ICA method to analyze the simulated experimental results of a dynamical probe, we will see that the ICA method can work but the outcome is not ideal. The reason that it does not work well can also be understood. Because the ICA only assumes that the signal has a certain feature but does not fully explore what exactly the feature is. This is good for the original 'cocktail party problem' because it does not require prior knowledge of each voice. However, as discussed above, in our quantum problem each eigen-mode has a fixed energy and manifests as a single frequency oscillation, but this feature of being single frequency oscillation is not utilized in the ICA method above. Therefore, the main result of this work is to present a modified ICA method that aims at finding out a signal of single frequency oscillation, short-noted as s-ICA.

The loss function of s-ICA
In short, let us consider a reference signal s ref = √ 2 cos(ωt). Here √ 2 is chosen to ensure that the variance equals one. The ICA method is to find out an A −1 such that each s i = ∑ j (A −1 ) ij x j is away from s Gau as much as possible; and our s-ICA method is to find out an A −1 such that each To quantify how s i is close to s ref , instead of using entropy we consider a quantity called the cumulants [5]. The cumulants is defined as where ⟨. . . ⟩ means performing average over a certain duration of t. For s ref , it is straightforward to compute where I 0 is the Bessel function of the first kind, and Therefore, our s-ICA method is to find out

Iteration in s-ICA
There is one subtlety in s-ICA. The analytical form equation (7) crucially relies on the distribution function equation (5) for s ref , but equation (5) is correct only when the statistics is performed over a duration that is an integer times of the period 2π/ω. Thus, when we calculate K[ ∑ j (A −1 ) ij x j , z], the statistics for all x j also needs to be carried out for the time interval being integer times of 2π/ω, otherwise K[s, z] always can not perfectly converge to K[s ref , z]. Hence, it enters a paradox. Since the goal is to separate out a single frequency oscillation to determine the frequency ω, ω is not known before analyzing. To solve this problem, our s-ICA method requires performing the minimization of L[s] iteratively. That is to say, we first choose an arbitrary period of t to perform statistics for x j , with which we minimize L[ ∑ j (A −1 ) ij x j ] to find out a s (1) i (t). Although s (1) i (t) is not a perfect single frequency oscillation because of this reason, we can still roughly determine a frequency ω (1) i by fitting the results with a single frequency sinusoidal oscillation. Then we perform statistics for x j with time duration 2π/ω i (t). Because this time, the statistics is performed in a time interval more close to integer times of the actual period, the cumulants is more close to the reference function equation (7). Therefore, s (2) i (t) will be more close to a single frequency oscillation than s (1) i (t), from which we can determine a frequency ω (2) i . We can continue the procedure and the accuracy of the frequency will be continuously improved, until we obtain a very good single frequency signal after k-steps and the frequency determined at each step also converges. This completes our s-ICA method. A comparison between ICA and our modified s-ICA method is shown in figure 1.

Analyzing collective modes of a BEC
As a demonstration of our method, we apply our method to analyze the collective oscillations of a Bose-Einstein condensate (BEC) of cold atoms in a harmonic trap, which is a very common experimental measurement in cold atom experiments [6][7][8][9]. We use this as an example simply because it is easy to generate simulated experimental data in this case. We first consider an equilibrium Thomas-Fermi density distribution where Θ(x) = max(0, x/|x|) is the Heaviside function, and for simplicity, we consider a two-dimensional geometry with a harmonic trap V(r) = mω 2 ⊥ (x 2 + y 2 )/2 and ω ⊥ being the trapping frequency. µ is the chemical potential that later will be adjusted to satisfy the total number of atoms conservation condition at any given time, and the total number of atoms is chosen as N = 10 3 in our simulation. g is the interaction strength, and we here set ℏ = 1, m = 1, ω ⊥ = 1 and g = 10. We consider that three different modes are excited, and they are which are the dipole mode, the quadrupole mode and the breathing mode, with their frequencies being ω 1 = 1, ω 2 = √ 2 and ω 3 = 2, respectively [10,11]. We also add noise that varies at different spatial locations and different time, denoted by ϵ(r, t) and rate τ to mimic the damping situation in the experiments. Then the simulated data of the real time density dynamics is given by (r, t)). (13) where f i (i = 1, 2, 3) denotes the amplitudes of the three modes.
Here we shall also note that, although x, x 2 − y 2 and x 2 + y 2 initially are three functions orthogonal to each other, due to the boundary condition imposed by the Heaviside function and the atom number conservation condition, they are no longer orthogonal. Thus, the method like the principle component analysis [12] also does not work very well in this case (See supplementary material for details (stacks.iop.org/MLST/1/025010/mmedia)). Now we will use the ICA method and our s-ICA method to separate out each δn i from data n(r, t) given by equation (13). As an example, we choose f 1 = 1, f 2 = 0.2, f 3 = 0.2 and τ = 20, and ε uniformly distributed in the range [−2, 2]. Here we only require short-time information of couple oscillation periods, and in fact, for this application another advantage is that we only need information from a few spatial points. Let us consider three locations denoted by r 1 , r 2 and r 3 , and the density dynamics at these points play the role as detectors. Their density dynamics all contain these three frequencies but the coefficients are different because it depends on the spatial locations. We plot n(r i , t) (i = 1, 2, 3) in figure 1(a)-(c), which show irregular temporal behaviors. The goal of ICA or s-ICA is to find out a proper combination of them that exhibits the single frequency oscillatory behavior.
where A 1 , A 2 , φ, ω are all fitting parameters and ω is what we want. Nevertheless, by fitting these signals, one can obtain ω Here we use the optimization processes to minimize the loss function equation (8) with Newton's method.
Here we should note that so far we only obtain the frequency but can determine neither the amplitude nor the damping rate. In fact, both ICA and our s-ICA have a problem that the amplitude can not be uniquely determined. In fact, signals figures 1(j)-(l) are all distributed between [− √ 2, √ 2] because the data n(r i , t) (i = 1, 2, 3) have been preprocessed to be zero mean and unity variance. For instance, to obtain the amplitude, we will perform another fitting. Here we consider density dynamics from all spatial points without pretreatment, that is Now all ω i are already known from the s-ICA analysis. For each r point, there are only four fitting parameters C i (r) (i = 0, . . . , 3) to fit data points of a sequence of t. The results of the fitting are shown in figure 3 which reveal the spatial information of each mode. In this case, figure 3(a) is the background. Figure 3(b) shows positive amplitude in one side and negative amplitude in another side, which changes the BEC's center-of-mass and is the dipole mode. Figure 3(c) has two nodal lines and concentrates at the surface, and it is the quadruple mode. Figure 3(d) shows negative amplitudes at the center and positive amplitudes in outside, which changes the BEC's side and is the breathing mode [11]. With similar fitting, we can also determine the damping rate of each mode.

Outlook
In summary, in this work we have developed a modified ICA method to extract the eigen-energy of a quantum system from a dynamical probe. We remark that our method and the Fourier transformation respectively work well in two different regimes. The Fourier transformation requires taking data over a sufficiently long time duration, but it does not require many data within each period. If the total time duration is short, the frequency resolution obtained from the Fourier transformation will be quite limited.
On the contrast, our method only needs data for a short time duration such as couple periods, but we require the accumulation of sufficient data during the time interval of a few oscillation periods that allows performing accurate statistics. Hence, our method and the Fourier transformation are very complementary, and therefore we believe our method can find wide application in future data analysis. Although the system considered here is a quantum one, the data is the expectation value of certain observable and is a class one. Hence, our method may even find its application beyond physics problems. Another remark is that it is of great interest to consider the quantum data and the quantum analogy of the 'cocktail party problem' , and to see whether the similar ICA method can work there [13].