We consider the general dimensionality reduction problem of locating in a high-dimensional data cloud, a $k$-dimensional non-Gaussian subspace of interesting features. We use a projection pursuit approach—we search for mutually orthogonal unit directions which maximise the 2-Wasserstein distance of the empirical distribution of data-projections along these directions from a standard Gaussian. Under a generative model, where there is a underlying (unknown) low-dimensional non-Gaussian subspace, we prove rigorous statistical guarantees on the accuracy of approximating this unknown subspace by the directions found by our projection pursuit approach. Our results operate in the regime where the data dimensionality is comparable to the sample size, and thus supplement the recent literature on the non-feasibility of locating interesting directions via projection pursuit in the complementary regime where the data dimensionality is much larger than the sample size.

§ INTRODUCTION

*Equal contribution.footnote
A central question in statistics and machine learning concerns the recovery of useful or interesting features from data. A huge body of literature exists that focuses on such feature extraction tasks. Often the statistician encounters high-dimensional data of which only a relatively low-dimensional subspace is of interest. A family of algorithms, often described by the umbrella term projection pursuit [5, 7], are particularly well-suited for such tasks. By restricting attention to low-dimensional subspaces, projection pursuit allows the statistician to evade the so-called “curse-of-dimensionality”, which plagues most classical learning algorithms in high-dimensional settings. Furthermore, projection pursuit helps the statistician to discard noisy and information-poor features. Some prominent members of this family of techniques include Principal Component Analysis (PCA), Independent Component Analysis (ICA), matching pursuit, etc.

Perhaps the simplest projection pursuit algorithm is PCA (see, e.g., [9],[8]), which  considers the subspace generated by the top $k$ eigenvectors of the sample covariance matrix. In effect, PCA tries to find linear combinations of the original features which explain the most variability. While very useful in its own regard, PCA is limited by the fact that it only considers variances. Thus it works very well when the superfluous features have considerably lower variance than the signal, e.g., in noise reduction problems.

On the other hand, consider a setup where the interesting components are non-Gaussian, while the rest are Gaussian of comparable variability. The Gaussian components cannot be treated as noise (in the sense of having smaller variance than the signal component) but are simply superfluous or “uninteresting”. PCA has difficulty separating the interesting non-Gaussian components in such scenarios.

To overcome this limitation, various methods conceptually similar to PCA have been proposed. A broad class of such methods goes by the name of ICA([10]). Broadly speaking, there are two families of ICA algorithms. One class of algorithms focuses more on ensuring that the signal directions are statistically independent, thus they minimize mutual information. The other focuses on finding directions in which the data is the “least Gaussian” (i.e. most interesting). In this paper, we are also interested in the latter objective.

Many approaches are possible for finding non-Gaussian directions, depending upon our definition of “non-Gaussian”. A natural way to quantify non-Gaussianity would be to measure the deviation of some aspect of a probability measure of interest from that of a standard Gaussian. For instance, one could use measures such as kurtosis([6]) or negentropy ([3],[12]). Alternatively, one could measure the deviation from a Gaussian using suitable probability metrics such as the Kolmogorov-Smirnov (KS) distance, the Wasserstein distance, etc.

We now state the general projection pursuit approach in the context of the problem of identifying non-Gaussian components with the following simple example. Consider a $p$-dimensional random vector $X$ which satisfies the following: there is a unknown direction $u_*$ such that $u_*^\top X$ is non-Gaussian, while $(I - u_*u_*^\top) X$, the distribution of $X$ in the orthogonal complement of $u_*$, is $(p - 1)$-dimensional standard Gaussian, and, further, the non-Gaussian component $u_*^\top X$ and the Gaussian part $(I - u_*u_*^\top)X$ are statistically independent.  Suppose $\dist(\nu, \nu')$ is some measure of quantifying the distance between two probability measures $\nu$ and $\nu'$. $\dist$ could be a divergence between probability measures (e.g., the Kullback-Liebler divergence) or a proper metric (e.g., KS distance).
Suppose we observe a sample $X_1, \ldots, X_n$ from $\nu$, the distribution of $X$. Our goal is to recover the unknown direction $u_*$. The main idea of projection pursuit is then to find a unit direction $\hat{u}$ such that the empirical distribution of the data projected on $\hat{u}$ (i.e. $\frac{1}{n}\sum_{i = 1}^n \delta_{\hat{u}^\top X_i}$) is the farthest from the standard Gaussian distribution with respect to $\dist$, i.e.

\[
\hat{u} = \argmax_{u \,:\, \norm{u} = 1} \dist\bigg(\frac{1}{n}\sum_{i = 1}^n \delta_{u^\top X_i}, \nu_g\bigg),
\]

where $\nu_g$ is the standard Gaussian measure. We will later formalise a version of this for general $k$. Our paper is interested in the case when $\dist$ is the 2-Wasserstein distance between probability measures with finite second moment. Specifically, we analyse the question of whether the recovered directions can be guaranteed (with high probability) to be from the signal space or not.

Some work in this regard has been done in [1] and [11]. [1] use the KS distance for $\dist$. They show that if the data is purely Gaussian (i.e. in a null model with no interesting directions), then two completely different phenomena occur according as whether the data-dimensionality-to-sample-size ratio $\frac{p}{n}$ goes to zero or infinity. In the former regime, all projections are Gaussians (in fact, this is known from the earlier work of [4]). On the other hand, in the latter regime, given any arbitrary probability distribution $\tilde{\nu}$, with high probability, one can find a (data-dependent) direction along which the data set is distributed as $\tilde{\nu}$. In other words, one can find directions along which the data is as far from Gaussianity (in the KS metric) as one desires. This means that projection-pursuit can spot fake signal amidst complete noise. When applied to actual data in this regime, there is no way of  knowing if the found direction came from some underlying signal space, or if it is a mirage of signal in a Gaussian desert!

[11] prove a similar result for the 2-Wasserstein distance. Both the above papers argue that under the null model of $N(0,\bbI_p)$, when $p/n$ converges to a sufficiently small constant, the empirical distributions of projections of the data points in every direction are close to the standard Gaussian distribution. This obviously begs the question if, under a spiked alternative, one could find directions along which the data projections are non-Gaussian. <cit.> study the question of obtaining such a signal direction under a specific model of supervised learning.

Our contributions. Our analysis is done in the context of an alternative model of unsupervised learning. We suppose that the sample comes from a spiked Gaussian model, i.e. there is a $k$-dimensional subspace in which the distribution is decidedly not Gaussian. We first show that under subgaussian tail assumptions, in every direction, the empirical distribution formed by the data-projections and the true marginal distribution in the same direction are uniformly close. This result is a substantial extension of similar results in [1] and [11] to a more general setting. Further, we also show, using a peculiar property of the $2$-Wasserstein distance, that one can recover an orthonormal set of vectors which form an approximate basis of the signal space. (This can be thought of as an instance of the general strategy of matching pursuit.) In particular, each recovered vector's component in the independent Gaussian space is inversely proportional to the signal-to-noise ratio. Finally, if the signal-to-noise ratio is sufficiently large, then we give a methodology to accurately estimate $k$, i.e. the dimension of the signal space. This allows one to use our sequential procedure even in cases where very little is known about the signal space.

§ SET-UP

Given a probability measure $\mu$ in $\bbR^p$, and a vector $v \in \bbR^p$, we define the action of $v$ on $\mu$, $v\sharp \mu$ to be the marginal density of $\mu$ in the direction $v$. In particular,

If $X$ is a random variable in $\bbR^p$ from the measure $\mu$ and $v^\top X$ is the dot product of $v$ and $X$, then $v\sharp \mu$ is defined to be the density function of the real-valued random variable $v^\top X$.

We also need the notion of the $2$-Wasserstein distance $d_{W_2}(\mu, \nu)$ between two probability measures $\mu$ and $\nu$, defined via

\[
d_{W_2}^2(\mu, \nu) \coloneqq \inf_{\substack{\text{all couplings $\pi$ of } X, Y \\ X \sim \mu, Y \sim \nu}} \bbE_\pi [X - Y]^2.
\]

We now introduce a formal set-up for the non-Gaussian component recovery problem which we will analyse.

Suppose that we have data $X_1, \ldots, X_n$ i.i.d. from a $\sigma$-subgaussian distribution $\nu$ on $\R^p$,[A random variable $X \in \bbR^p$ with mean $\mu$ is subgaussian with parameter $\sigma$ or the distribution is in $SG_p(\sigma)$ iff $P\left( \norm{X -\mu} \geq t \right)  \leq Ce^{-\frac{t^2}{2\sigma^2}}.$]  where $\Psi$ has the following structure: Suppose $X \sim \Psi$.

* $\bbE X = 0; \var(X) = I_p$[While for the purposes of our proof, we have assumed that the covariance matrix is identity throughout, this assumption is heuristically not much different from working with whitened data.].

* There is a (unknown) $k$-dimensional subspace $U$ such that $\Pi_U X$ has a sufficiently non-Gaussian distribution, and $\Pi_{U^{\perp}}X$ has a close-to-Gaussian distribution, in the sense that there exists constants $\kappa_1, \kappa_2$ with,

\[
\inf_{u \in U, \|u\| = 1} d_{W_2}(u \sharp \Psi , \Phi) > \kappa_1 > \kappa_2 >
\sup_{u \in U^{\perp}, \|u\| = 1}d_{W_2}(u \sharp \Psi, \Phi).
\]

* $\Pi_U X$ and $\Pi_{U^\perp} X$ are independent.

We will denote the gaussian subspace, $U^{\perp}$ by $W$.
Given the sample $X_1, \cdots, X_n$, our goal is to recover the space $U$. The quantities $\kappa_1$ and $\kappa_2$, as we will see below, dictate a separation condition necessary to distinguish the non-Gaussian signal components from the Gaussian part.

To further motivate our setup let us quickly look at a simple distribution satisfying the assumptions above :

Let the data $X$ be generated from a mixture of gaussians i.e. $\sum_{i=1}^k \frac{N(u_i, \bbI_p)}{k}$, where $u_i$'s are some vectors in $\bbR^p$. Note quickly that if $v$ is any norm $1$ vector orthogonal to all $u_i$, then $v^\top X$ follows $N(0,1)$ and is independent of $u_i^\top X$. Clearly then this is a specific example of our model, with $U = \textit{Span}\{u_1, ..., u_k\}$ and $W = U^{\perp}$. In fact the projection of $\Psi$ to the subspace $W$ is the distribution $N(0,I_{p-k})$.

Before moving onto the technical results we also quickly define the following notations we will be using throughout.

Given a $p$-dimensional distribution $\nu$ and a $1$-dimensional distribution $\mu$

\[
d(\nu, \mu) \coloneqq \sup_{u : \norm{u} = 1} d_{w_2}(u\sharp\nu, \mu).
\]

Note that when $p=1$, this is simply the $2$-Wasserstein distance between the $\nu$ and $\mu$. For larger $p$, when $\mu =\Phi$, our distance $d$ captures how non-gaussian the distribution $\nu$ can become in a particular direction.

Given a $p$-dimensional distribution $\nu$ and a $1$-dimensional distribution $\mu$

\[
d_{min}(\nu, \mu) \coloneqq \inf_{u : \norm{u} = 1} d_{w_2}(u\sharp\nu, \mu).
\]

In essence, when $\mu = \Phi$, $d_{min}$ gives a measure of separation from the $1$-dimensional gaussian.

Let $W$ be some subspace of $\bbR^p$, and let $q_1, ..., q_t$ be some orthonormal basis of $W$. If $Q$ is the matrix whose columns are given by $q_i$ then given a random variable $X \in \bbR^p$, where $X \sim \mu$, we define $\mu_{|Q}$ to be the distribution of $Q^TX$.

We remark that as the distance $d(\mu,\Phi)$ is rotationally invariant, given a fixed subspace $W$, $d(\mu_{|Q},\Phi)$ is the same regardless of what orthonormal basis one chooses. Thus one can consider the quantity $d(\mu_{|W},\Phi)$ unambiguously.

Finally for the sake of clarity of our conclusion we define a signal to noise ratio for the distribution $\Psi$, $SNR$ as

\[
SNR = \sqrt{\frac{d(\Psi,\Phi)^2 - d(\Psi_{|W},\Phi)^2}{d(\Psi, \Phi)^2 - d_{min}(\Psi_{|U},\Phi)^2}}.
\]

Note that based on our definition of $\kappa_1$ and $\kappa_2$, our SNR is always larger than $\sqrt{\frac{d(\Psi,\Phi)^2 - \kappa_2^2}{d(\Psi, \Phi)^2 - \kappa_1^2}}$.
We will show that we can construct with high probability $k$ orthonormal vectors  $v_1, ..., v_k$ such that

\[
\norm{\textit{Proj}_{W}(\hat{v}_j)} \leq \frac{2}{SNR}.
\]

We also show that we can estimate $k$ if $k \leq \frac{SNR^2}{4}$.

§ MAIN RESULTS

The following proposition is the central pivot granting us leverage to most of our results.

Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $n$,$p$ go to infinity in a way such that $p/n \rightarrow \gamma$. Then given a positive constant $\epsilon$, there exists a positive constant $\gamma_{\sigma,\epsilon}$ depending on $\sigma$ and $\epsilon$ such that when $\gamma \leq \gamma_{\sigma,\epsilon}$, we have

$$P\left(sup_{u \in S_{p-1}} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right) - \bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right]}  > \epsilon \right) < De^{-nc_{\sigma, \gamma,\epsilon}},$$

where $c_{\sigma,\gamma,\epsilon}$ is some positive constant dependent on $\sigma$, $\gamma$, and $\epsilon$.

Proposition <ref> uniformly bounds the difference between the data dependent (and thus random) quantity, $d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right),$ and the deterministic quantity,$ \bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right],$ dependent only on $u$.

Then given our assumptions <ref> on the distribution $\Psi$ above, we will state the following theorem (proved in Section <ref>):

Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $n$,$p$ go to infinity in a way such that $p/n \to \gamma$. Given an $\epsilon > 0$, there exists a constants $\gamma_{\sigma,\epsilon}$ dependent on $\sigma$ and $\epsilon$  and  $C_{\sigma}$  depending on $\sigma$ such that if $\gamma \leq \gamma_{\sigma,\epsilon}$, the following statement is true with high probability for all  unit vectors $u$ in $\bbR^p$ simultaneously :

\[\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right) - d_{W_2}\left(u\sharp\Psi, \Phi\right)}  \leq \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}.\]

We note that Theorem <ref> needs very little assumptions on the distribution $\Psi$. We only need $\Psi$ to be $\sigma$-subgaussian. The main upshot of the theorem is that it implies with uniform high probability that in every direction the empirical distribution of the projection is as far away from Gaussian, as the true marginal distribution in that direction. Thus heuristically if we want to find directions in which $\Psi$ is not gaussian it makes sense to maximise the quantity $d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right).$ We can now proceed to state conditions under which the recovered directions have a very small component in the Gaussian subspace, $W$. A proof of Theorem <ref>  is written in Section <ref>.

Let $W$ be the gaussian subspace of $\Psi$. Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $n$,$p$ go to infinity in a way such that $p/n \to \gamma$. Given $\epsilon > 0$, there exists a constant $\gamma_{\sigma,\epsilon}$ dependent on $\sigma$ and $\epsilon$ such that if $\gamma \leq \gamma_{\sigma,\epsilon}$, then with asymptotic high probability for any $u \in S_{p-1}$ such that $d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right) \geq \sqrt{1-\delta^2}d(\Psi,\Phi) + \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}, $

we have that
\[\norm{\textit{Proj}_{W}(u)} \leq \delta\frac{d(\Psi,\Phi)}{\sqrt{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}}.\]

Now that we have introduced most of our bulky technology, we can use it to prove the following simple Corollary. In the interest of space we have moved a detailed proof of the Corollary to the Appendix in Section <ref>. This in turn allows us to argue the validity of procedure in the vein of the general idea of matching pursuit.

Let $U$ be a $k$-dimensional sub-space of $\bbR^p$, where $k$ is a constant. Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $l < k$ be some integer. Let $v_1 , ..., v_l$ be some vectors in $\bbR^p$. Then there exists some constant $C_{\sigma}$, depending on $\sigma$ such that given $\epsilon > 0$ there exists with high probability a unit vector, $u$ which is orthonormal to all $v_i$ such that

\[
d_{W_2}\bigg(\frac{1}{n}\sum_{i=1}^n \delta_{u^\top X_i}, \Phi\bigg) \geq d_{min}(\Psi_{|U},\Phi) - \epsilon -  \frac{C_{\sigma}}{\sqrt[4]{n}}.
\]

Suppose now that the distribution $\Psi$ is such that there is a $k$ dimensional subspace $U$ such that $d_{min}(\Psi_{|U},\Phi) = \sqrt{1-\delta^2} d(\Psi,\Phi)$. That is the “top $k$ directions” are a constant factor far from gaussian as the maximum possible. Then using Corollary <ref> we can with high probability sequentially construct vectors $v_1, ..., v_k$ such that for every $1 \leq j \leq k$,

\[
d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v_j^TX_i}}{n}, \Phi\right) \geq \sqrt{1-\delta^2} d(\Psi,\Phi) - \epsilon - \frac{C_{\sigma}}{\sqrt[4]{n}}.
\]

Then setting $\epsilon = \frac{\delta^2 d(\Psi,\Phi)}{2}$, and $n$ large enough such that $\frac{4C_{\sigma}}{\sqrt[4]{n}} \leq  \delta^2 d(\Psi,\Phi)$, we have that

\[
d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v_j^TX_i}}{n}, \Phi\right) \geq \sqrt{1-\delta^2} d(\Psi,\Phi) - \epsilon - \frac{C_{\sigma}}{\sqrt[4]{n}} \geq \sqrt{1-4\delta^2} d(\Psi,\Phi) + \epsilon  +\frac{C_{\sigma}}{\sqrt[4]{n}}..
\]

Now we can use Theorem <ref> with $\epsilon$ as set above. Thus if $p/n$ converges to a sufficiently small constant $\gamma$, then for large enough $n$ with high probability we have that for every $j$,

\[
\norm{\textit{Proj}_{W}(v_j)} \leq 2\delta\frac{d(\Psi,\Phi)}{\sqrt{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}}.
\]

In other words the $k$-space that we found (i.e. the one spanned by $v_1, ..., v_k$) is mostly orthogonal to $W$, the subspace where the distribution is close to Gaussian. The above discussion then gives the following natural method to estimate $k$ vectors which are almost orthogonal to $W$. For $1 \leq j \leq k$, let

\[
\hat{v}_j =  \argmax_{\norm{\hat{v}_j} = 1; \forall t < j : \hat{v}_j^T\hat{v}_t = 0} d_{W_2}\bigg(\frac{1}{n}\sum_{i=1}^n \delta_{\hat{v}_j^TX_i},\Phi\bigg).
\]

The above discussion then implies that, with high probability for large enough $n$ we have that (by invoking Theorem <ref> with $\delta$ such that $d_{min}(\Psi_{|U},\Phi) = \sqrt{1-\delta^2} d(\Psi,\Phi)$):

\[
\norm{\textit{Proj}_{W}(\hat{v}_j)} \leq 2\sqrt{\frac{d(\Psi,\Phi)^2 - d_{min}(\Psi_{|U},\Phi)^2}{d(\Psi, \Phi)^2 - d(\Psi_{|W},\Phi)^2}} = \frac{2}{SNR}.
\]

A common problem that often occurs in such problems is that $k$ is unknown. To give some answer to this question we first consider the following corollary which is proved in detail in the Appendix in Section <ref>.

Given  integers $m > k+1$, let $\delta$ be a positive real number such that $4\delta^2 < \frac{1}{m} \left(1 - \frac{d(\Psi_{|W}, \Phi)^2}{d(\Psi, \Phi)^2}\right).$ Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $n$,$p$ go to infinity in a way such that $p/n \rightarrow \gamma$ Given $\epsilon > 0$ there is a $\gamma_{\sigma,\epsilon}$, where $\gamma_{\sigma,\epsilon}$ is a constant depending on $\sigma$, $\epsilon$, such that if $\gamma \leq \gamma_{\sigma,\epsilon}$ then with high probability there does not exist a set of $k+1$ orthonormal unit vectors $v_1, ..., v_{k+1}$ such that

\[
d_{W_2}\bigg(\frac{1}{n}\sum_{i=1}^n \delta_{v_j^TX_i}, \Phi\bigg) \geq \sqrt{1-4\delta^2} d(\Psi,\Phi) + \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}. \]

Continuing the discussion prior to the corollary, we consider $\delta$ such that $d_{min}(\Psi_{|U},\Phi) = \sqrt{1-\delta^2} d(\Psi,\Phi)$. Note that when $k+1 < \frac{1}{4}\frac{d(\Psi,\Phi)^2 - d(\Psi_{|W},\Phi)^2}{d(\Psi, \Phi)^2 - d_{min}(\Psi_{|U},\Phi)^2} = \frac{SNR^2}{4}$ the hypothesis of Corollary <ref> is true. This gives us a natural cutoff point for our sequential algorithm. We can stop at $\hat{k}$, if for $\epsilon = \frac{\delta^2 d(\Psi,\Phi)}{2}$, and $n$ large enough such that $\frac{4C_{\sigma}}{\sqrt[4]{n}} \leq  \delta^2 d(\Psi,\Phi)$ we have that

\[
d_{W_2}\bigg(\frac{1}{n}\sum_{i=1}^n \delta_{\hat{v}_{\hat{k}+1}^TX_i}, \Phi\bigg) < \sqrt{1-4\delta^2} d(\Psi,\Phi) + \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}. \]

Corollary <ref> then implies that this stopping rule ensures with high probability that $\hat{k} \leq k$. On the other hand, the discussion following Corollary <ref> means that the same stopping rule ensures that $\hat{k} \geq k$. In effect we have that if $k+1 < \frac{SNR^2}{4}$, then with high probability $\hat{k} = k$.

§ CONCLUSION

In this article, we have considered the problem of isolating a non-Gaussian independent component from a Gaussian counterpart under certain separability assumptions. We have theoretically analysed the approximation accuracy of a projection pursuit procedure. In contrast to more traditional procedures like PCA, we do not need the variances of the superfluous feature directions to be small. We only need a distributional gap between directions which are Gaussian and those which are not.

Since the proposed method involves optimisation of the objective function $d_{W_2}\big(\frac{1}{n}\sum_{i=1}^n \delta_{v^\top X_i}, \Phi\big)$ as $v$ varies over the unit sphere, two natural questions immediately come to mind. First of all, since our objective function is markedly non-convex, designing an efficient algorithm that can find a global minimum (or even good local minima) would be a significant addition to present work.

Secondly, it needs to be investigated if similar results are true for distances other than the $2$-Wasserstein distance. It is plausible that some distances would be more suitable both from a theoretical perspective and also the practical optimisation aspect. We leave the investigation of these questions for future work.

§ PROOFS

§.§ Proof of Proposition <ref>

Quickly noting that for any vector $u \in \bbR^p$, with $\norm{u}=1$ we have,
\[ P\left( \abs{u.(X-\mu)} \geq t \right) \leq P\left( \norm{X -\mu} \geq t \right),\]

we get the following simple proposition :

If $X \in \bbR^p$ is in $SG_p(\sigma)$ and $u \in \bbR^p$ be any unit norm vector (i.e. $\norm{u} = 1$), then $u^\top X \in SG_1(\sigma)$.

The following is a simple proposition which bounds the norms of the sample covariance matrix of subgaussian random variables. It is a slightly reworded version of Theorem 6.5 of [13]

Let $X_1, ..., X_n$ be iid sample from a $\sigma$ subgaussian distribution in $\bbR^p$ with covariance matrix $\bbI$. Then there exists universal constants $c_1,c_2,c_3$ such that we have for all $\delta > 0$ that

\[P\left( \norm{\frac{\sum X_iX_i^T}{n}}_2 \geq 1 + \sigma^2\left(c_1 \left(\sqrt{\frac{p}{n}} + \frac{p}{n} \right)  + \delta\right) \right) \leq c_2e^{-nc_3\min(\delta, \delta^2)}
\]

We will also use the following result (for $p=2$) on Wasserstein distances and sample convergences found in [2] as Corollary 7.17.

Let $p$ be some positive integer. Let $\mu$ be some distribution such that for some $s > p$ its $s$'th moment exists and is bounded. Then if $X_1, ..., X_n$ are iid random variables sampled from $\mu$, we have

\[
\bbE_{X_1, ..., X_n}\left[d_{W_p}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \mu \right)^p\right] \leq \frac{C}{\sqrt{n}},
\]

where $C$ is some absolute constant dependent upon the upper bound of the $s$'th moment.

To prove Proposition <ref>, we will finally be needing the following lemma (proved in the Appendix in section <ref>,) on the concentration of the 2-Wasserstein distance between a subgaussian measure $\mu$ and the empirical measure $\mu_n = \frac{1}{n}\sum_{i = 1}^n \delta_{X_i}$ of an i.i.d. sample $X_1, \ldots, X_n$ from $\mu$ (a similar result with a log-Sobolev assumption on $\mu$ appears as Theorem 7.1 in [2]).

Let $\mu$ be a $\sigma$-subgaussian  measure. Let $\mu_n$ be the empirical measure formed from an i.i.d. sample of size $n$ from $\mu$. Then
\[
\bbP\big(\abs{d_{W_2}(\mu_n, \mu) - \bbE d_{W_2}(\mu_n, \mu)} \geq t \big) \leq Ce^{-\frac{cnt^2}{\sigma^2}}.
\]

for some absolute constants $C, c > 0$.

Now armed with the above preliminaries we can move onto proving our central results:
As mentioned before we begin by using Lemma <ref> and the hypothesis that $u^\top X$ is $\sigma$-subgaussian to get that for any fixed $u$ (with $\norm{u}=1$) we have
\[P\left(\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq t\right) \leq A'e^{-\frac{cn\epsilon^2}{\sigma^2}}, \]

where the constants $A',c$ are absolute constants

Let $E_{p,\delta}$ be the smallest delta net on the unit sphere in $\bbR^p$, i.e. given any $u \in \bbS_{p-1}$, $\exists$ a $v \in E_{p,\delta}$ such that $\norm{u-v} \leq \delta$. It is known that there exists such a net for any $p$ such that $\abs{E_{p,\delta}} \leq A''\frac{1}{\delta^p}$. Then we have that
\[P\left(\underset{u \in E_{p,\delta}}{\sup} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq \epsilon\right) \leq \left(\frac{1}{\delta}\right)^pAe^{-\frac{cn\epsilon^2}{\sigma^2}}.\]

To go from taking the supremum over the net to that on the entire sphere then we would have to control how small changes in $u$ affect the quantity of interest.  To that end let $u \in \bbS_{p-1}$. Let $v \in \bbE_{p,\delta}$ such that $\norm{u-v} \leq \delta$. Then we have

\begin{align*}
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi \right)} &\leq \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, u\sharp\Psi \right)} \\ &+ \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, u\sharp\Psi \right) - d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi \right)} \\
&\leq \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n} \right)} + \abs{d_{W_2}\left(u\sharp\Psi, v.\Psi \right)}
\end{align*}

Let us then quickly bound both of the two terms above. The second term can be upper bounded as

\begin{align*}
\abs{d_{W_2}\left(u\sharp\Psi, v.\Psi\right)}^2 &= \inf_{\substack{(X,X') : \textit{ the marginals } \\ X \textit{ and } X' \textit{ distributed as } \Psi}} \bbE[(u^\top X - v^\top X')^2] \\
&\leq \bbE_{X \sim \Psi}[(u^\top X - v^\top X)^2] & (\textit{Considering the coupling } X = X') \\
&= (u-v)^T \bbE[XX^T] (u-v) \\
&\leq \delta^2. &(\textit{as } \bbE[XX^T] = \bbI_p)
\end{align*}

Similarly for the first term we note

\begin{align*}
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n} \right)}^2 &\leq \sum_{i=1}^n \frac{(u^\top X_i - v^\top X_i)^2}{n} \\
&= (u-v)^T \left(\frac{1}{n}\sum_{i=1}^n X_iX_i^T \right) (u-v) \\
&\leq \left(\delta \norm{\hat{\Sigma}_n}_2\right)^2,
\end{align*}

where $\hat{\Sigma}_n$ is the sample covariance matrix and $\norm{.}_2$ denotes the operator or $L_2$ norm. We will now use a technical claim that the operator norm of the sample covariance matrix is with high probability smaller than $2 + C'\sigma^2\left(\sqrt{\frac{p}{n}} +\frac{p}{n}\right)$, for some universal constant $C'$. This follows from invoking Proposition <ref>.

Then we can condition on this event as this is true with high probability ($1-e^{-\theta n}$). Thus we have whp

\[\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n} \right)} \leq \delta \left(3 + C'\sigma^2\left(\sqrt{\frac{p}{n}}+\frac{p}{n}\right)\right).\]

Combining everything then we have whp
\[\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi \right)} \leq \delta \left(3 + C'\sigma^2\left(\sqrt{\frac{p}{n}}+\frac{p}{n}\right)\right).\]

We then have that whenever there exists an $u \in \bbS_{p-1}$ with

\[
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq \epsilon + 2\delta \left(3 + C'\sigma^2\left(\sqrt{\frac{p}{n}}+\frac{p}{n}\right)\right),
\]

there exists a $v \in E_{p,\delta}$ whp (with the property that $\norm{u-v}\leq \delta$) such that
$$\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq \epsilon.$$

Thus using the probability bound on the delta net gives us

\begin{align*}&P\left(\underset{u \in \bbS_{p-1}}{\sup} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq \epsilon + 2\delta  \left(3 + C'\sigma^2\sqrt{\frac{p}{n}}+C'\sigma^2\frac{p}{n} \right)\right) \\
&\leq P\left(\underset{v \in E_{p,\delta}}{\sup} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi \right)\right]} \geq \epsilon\right) + e^{-\theta n}\\
&\leq  \left(\frac{1}{\delta}\right)^pAe^{-\frac{cn\epsilon^2}{\sigma^2}} + e^{-\theta n} = Ae^{-\frac{cn\epsilon^2}{\sigma^2} - p\log \delta} + e^{-\theta n},
\end{align*}

where we get the extra $\frac{p}{n}$ term as it is no longer true that $\frac{p}{n} \ll \sqrt{\frac{p}{n}}$ (and from invoking Proposition <ref>)
Using the hypothesis that $\frac{p}{n} \rightarrow \gamma$ and choosing a $\delta$ such that
\[\delta = \frac{\epsilon}{2\left(3 + C'\sigma^2\sqrt{\gamma} + C'\sigma^2\gamma\right)},\]

we get :

\begin{align*}
&P\left(\underset{u \in \bbS_{p-1}}{\sup} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right) - \bbE_{X_1, ..., X_n} \left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]} \geq 2\epsilon\right) \\
&\leq A\exp\left(-\frac{cn\epsilon^2}{\sigma^2} - p \log \frac{\epsilon}{6 + 2C'\sigma^2\sqrt{\gamma} + 2C'\sigma^2\gamma }\right) + e^{-\theta n} \\
&\rightarrow A\exp -n\left(\frac{c\epsilon^2}{\sigma^2} + \gamma \log \frac{\epsilon}{6 + 2C'\sigma^2\sqrt{\gamma} + 2C'\sigma^2\gamma }\right) + e^{-\theta n}
\end{align*}

As there is some constant $\gamma_{\sigma,t}$ such that for all $\gamma \leq \gamma_{\sigma,\epsilon}$, $\frac{c\epsilon^2}{\sigma^2} + \gamma \log \frac{\epsilon}{6 + 2C'\sigma^2\sqrt{\gamma} + 2C'\sigma^2\gamma }$ is positive and lower bounded, the proof follows.

§.§ Proof of Theorem <ref>

Note that Proposition <ref> probabilistically bounds the difference between $d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)$ and the deterministic quantity $\bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right]$ as $u$ varies over all possible unit vectors. Define $\gamma_{\sigma, \epsilon}$ as in Proposition <ref>. Thus when $\gamma \leq \gamma_{\sigma, \epsilon}$ we have with high probability for all unit vectors $u \in \bbR^p$ simultaneously that

\begin{align} \label{eq:pconc-1}
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right) - \bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right]}  \leq \epsilon
\end{align}

Since $u\sharp\Psi$ is $\sigma$-subgaussian we can conclude that all its moments are upper bounded (by a suitable function of $\sigma$). Choose any $s > 4$ and combine the upper bound on $s$'th moment of $u^\top X$ with Proposition <ref> (Corollary 7.17 of [2]).  We then get that there is a constant $C_{\sigma}$, dependent on $\sigma$, such that

\begin{align}\label{eq:pconc-2}
\bbE\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)\right]^2 \leq \bbE\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi \right)^2\right] \leq \frac{C_{\sigma}^2}{\sqrt{n}}.
\end{align}

Combining equations <ref> and <ref> with the triangle inequality we can write with high probability that

\begin{align*}
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right) - d_{W_2}\left(u\sharp\Psi, \Phi\right)} &\leq d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right) \\
&\leq \epsilon + \bbE\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right] \leq \epsilon + \frac{C_{\sigma}^2}{\sqrt{n}}.
\end{align*}

§.§ Proof of Theorem <ref>

Let $u = \alpha v + \sqrt{1-\alpha^2} w$, where $w \in W$, $v \in W^{\perp}$ and $\norm{w} = \norm{v} = 1$. Let $X$ be a random variable distributed as $\Psi$.
Thus $u^\top X = \alpha v^\top X + \sqrt{1 - \alpha^2} w^\top X$, where $v^\top X$ and $w^\top X$ are independent random variables following the distributions $v.\Psi$ and $w\sharp\Psi$ respectively. In other words if $Y_1$ and $Y_2$ are two independent random variables from the distributions $v.\Psi$ and $w\sharp\Psi$ respectively we can write
\begin{align} \label{eq:couplingXYZ}
u^\top X \eqd \alpha Y_1 + \sqrt{1-\alpha^2} Y_2,
\end{align}
where $\eqd$ means equal in distribution. Note then that if $Z$ is a random variable distributed as $\Phi$(i.e.  $N(0,1)$), and $Z_2, Z_3$ are two iid copies distributed as $\Phi$, we can also write
\begin{align} \label{eq:couplingZZZ}
Z \eqd \alpha Z_1 + \sqrt{1-\alpha^2} Z_2.
\end{align}

Let $\Omega$,  be the set of all possible couplings of the distributions of $u\sharp\Psi$ and $\Phi$. Similarly, let $\Omega_1$ (resp. $\Omega_2$) be the set of all possible couplings of the distributions of $v.\Psi$ and $\Phi$ (resp. $w\sharp\Psi$ and $\Phi$). Then there is a natural way to construct a coupling in $\Omega$ given a coupling in $\Omega_1$ and another in $\Omega_2$. That given any joint distribution $\mu$ in $\Omega_1$ whose marginals are $v.\Psi$ and $\Phi$ respectively, define $Y_1$ and $Z_1$ be the corresponding marginal random variables i.e. $(Y_1,Z_1) \sim \mu$. Similarly given any joint distribution $\nu$ in $\Omega_2$, we can define random variables $Y_2, Z_2$ where $(Y_2,Z_2) \sim \nu.$ Note that by construction we can keep the pair $(Y_1,Z_1)$ independent of $(Y_2,Z_2)$. Then equations <ref> and <ref> can be used to define a joint distribution in $\Omega$.
We then derive the following inequality,

\begin{align*}
d_{W_2}\left(u\sharp\Psi, \Phi\right)^2 &=  \underset{(u^\top X,Z) \in \Omega}{inf} \bbE \left[(u^\top X-Z)^2 \right] \\
&\leq \underset{(Y_1,Z_1) \in \Omega_1 \textit{ and } (Y_2,Z_2) \in \Omega_2}{inf} \bbE \left[(\alpha Y_1 + \sqrt{1-\alpha^2} Y_2 - \alpha Z_1 - \sqrt{1-\alpha^2} Z_2)^2 \right] \\
&= \underset{(Y_1,Z_1) \in \Omega_1 \textit{ and } (Y_2,Z_2) \in \Omega_2}{inf} \alpha^2\bbE \left[(Y_1 - Z_1 )^2 \right] + (1-\alpha^2)\bbE \left[(Y_2 -  Z_2)^2 \right] + 2\alpha\sqrt{1-\alpha^2}\bbE \left[(Y_1- Z_1)(Y_2 -  Z_2) \right] \\
&= \underset{(Y_1,Z_1) \in \Omega_1 \textit{ and } (Y_2,Z_2) \in \Omega_2}{inf} \alpha^2\bbE \left[(Y_1 - Z_1 )^2 \right] + (1-\alpha^2)\bbE \left[(Y_2 -  Z_2)^2 \right] \\ %&\textit{as $(Y_1,Z_1)$ independent of $(Y_2,Z_2)$}\\
&= \alpha^2 \underset{(Y_1,Z_1) \in \Omega_1}{inf} \bbE \left[(Y_1 - Z_1)^2 \right]  +  (1-\alpha^2)\underset{(Y_2,Z_2) \in \Omega_2}{inf}\bbE \left[(Y_2 -  Z_2)^2 \right]\\
&= \alpha^2 d_{W_2} (v.\Psi, \Phi)^2 + (1-\alpha^2)d_{W_2} (w\sharp\Psi, \Phi)^2\\
&\leq \alpha^2 d(\Psi, \Phi)^2 + (1-\alpha^2)d(\Psi_{|W}, \Phi)^2.
\end{align*}

Rewriting the above we get

\[\norm{\textit{Proj}_{W^{\perp}}(u)}^2 = \alpha^2 \geq \frac{d_{W_2}\left(u\sharp\Psi, \Phi\right) - d (\Psi_{|W}, \Phi)^2}{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}. \]

Using Theorem <ref> and the hypothesis, we have for an appropriate $\gamma_{\sigma,\epsilon}$ and $C_{\sigma}$ that $d_{W_2}(u\sharp\Psi, \Phi) \geq \sqrt{(1-\delta^2)}d(\Psi,\Phi)$

\[\norm{\textit{Proj}_{W^{\perp}}(u)}^2 \geq \frac{d_{W_2}\left(u\sharp\Psi, \Phi\right) - d (\Psi_{|W}, \Phi)^2}{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2} \geq\frac{(1-\delta^2)d(\Psi,\Phi)^2 - d (\Psi_{|W}, \Phi)^2}{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}.\]

Finally, we get
$\norm{\textit{Proj}_{W}(u)} = \sqrt{1-\norm{\textit{Proj}_{W^{\perp}}(u)}^2} \leq \delta\frac{d(\Psi,\Phi)}{\sqrt{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}}.$

[1]
Peter J Bickel, Gil Kur, and Boaz Nadler.
Projection pursuit in high dimensions.
Proceedings of the National Academy of Sciences,
115(37):9151–9156, 2018.

[2]
Sergey Bobkov and Michel Ledoux.
One-dimensional empirical measures, order statistics, and
Kantorovich transport distances, volume 261.
American Mathematical Society, 2019.

[3]
LJ Cao, Kok Seng Chua, WK Chong, HP Lee, and QM Gu.
A comparison of pca, kpca and ica for dimensionality reduction in
support vector machine.
Neurocomputing, 55(1-2):321–336, 2003.

[4]
Persi Diaconis and David Freedman.
Asymptotics of graphical projection pursuit.
The annals of statistics, pages 793–815, 1984.

[5]
Jerome H Friedman and John W Tukey.
A projection pursuit algorithm for exploratory data analysis.
IEEE Transactions on computers, 100(9):881–890, 1974.

[6]
Mark Girolami and Colin Fyfe.
Negentropy and kurtosis as projection pursuit indices provide
generalised ica algorithms.
In Advances in Neural Information Processing Systems Workshop,
volume 9. Denver, CO, 1996.

[7]
Peter J Huber.
Projection pursuit.
The annals of Statistics, pages 435–475, 1985.

[8]
Ian T Jolliffe and Jorge Cadima.
Principal component analysis: a review and recent developments.
Philosophical transactions of the royal society A: Mathematical,
Physical and Engineering Sciences, 374(2065):20150202, 2016.

[9]
Ian T Jolliffe.
Principal component analysis for special types of data.
Springer, 2002.

[10]
Te-Won Lee and Te-Won Lee.
Independent component analysis.
Springer, 1998.

[11]
Andrea Montanari and Kangjie Zhou.
Overparametrized linear dimensionality reductions: From projection
pursuit to two-layer neural networks.
arXiv preprint arXiv:2206.06526, 2022.

[12]
Michael Novey and Tulay Adali.
Complex ica by negentropy maximization.
IEEE Transactions on Neural Networks, 19(4):596–609, 2008.

[13]
Martin J Wainwright.
High-dimensional statistics: A non-asymptotic viewpoint,
volume 48.
Cambridge university press, 2019.

§ PROOF OF COROLLARY <REF>

Here we restate Corollary <ref> for convenience.

Let $U$ be a $k$-dimensional sub-space of $\bbR^p$, where $k$ is a constant. Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $l < k$ be some integer. Let $v_1 , ..., v_l$ be some vectors in $\bbR^p$. Then there exists some constant $C_{\sigma}$, depending on $\sigma$ such that given $\epsilon > 0$ there exists with high probability a unit vector, $u$ which is orthonormal to all $v_i$ such that

\[
d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right) \geq d_{min}(\Psi_{|U},\Phi) - \epsilon -  \frac{C_{\sigma}}{\sqrt[4]{n}}.
\]

As $l < k$ and $\dim(W) = k$, by rank nullity theorem, there exists a unit vector $u$ in $U$ which is orthogonal to all the vectors $\{v_1, ..., v_l\}$. As $k$ is constant. $k \ll n$ thus we can invoke Proposition <ref> with $p=k$ to get that with high probability

\[
\sup_{v : v \in U \textit{ and } \norm{v} = 1} \abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi_{\Phi}\right) - \bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v^\top X_i}}{n}, v.\Psi_{\Phi}\right)\right]}  < \epsilon.
\]

In particular then

\[
\abs{d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi_{\Phi}\right) - \bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi_{\Phi}\right)\right]}  < \epsilon.
\]

Finally similar to the proof of Theorem <ref> invoking Proposition <ref> gives us

\[
\bbE_{X_1, ..., X_n}\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi_{\Phi}\right)\right] \leq \frac{C_{\sigma}}{\sqrt[4]{n}}.
\]

Therefore we have

\begin{align*}
d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, \Phi\right) &\geq d_{W_2}\left(u\sharp\Psi,\Phi\right) -  d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)  \\
&\geq d_{min}(\Psi_{|U},\Phi) - \epsilon - \bbE\left[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{u^\top X_i}}{n}, u\sharp\Psi\right)\right] \\
&\geq d_{min}(\Psi_{|U},\Phi)  - \epsilon - \frac{C_{\sigma}}{\sqrt[4]{n}}
\end{align*}

§ PROOF OF COROLLARY <REF>

We will need the following simple linear algebraic lemma.
Let $v_1, ..., v_k$ be a set of orthonormal vectors in a vector space $V$. Let $G$ be the subspace spanned by $v_1, ..., v_k$. Let $H$ be a subspace of $V$. Then $\forall g \in G$ and $h \in H$ such that $\norm{g} = \norm{h} = 1$, we have that
\[\left(g^Th\right)^2 \leq \sum_{j=1}^k \norm{\textit{Proj}_{H}(v_j)}^2\]

We first remember from basic linear algebra that for any unit vectors $v \in V$ and $h \in H$, we have $\abs{v.h} \leq \abs{\textit{Proj}_{H}(v)}.$ Then we can write $g = \sum_{j=1}^k \alpha_j v_j$, where $\sum_j \alpha_j^2 = 1$ as $\norm{g} = 1$ and $v_1, \dots, v_k$ form an orthonormal basis of $G$. Combining these we get

\begin{align*}
(g^Th)^2 &= \left(\sum_{j=1}^k  \alpha_j v_j^Th\right)^2 \\
&\leq \left(\sum_{j=1}^k  \alpha_i^2 \right)\left(\sum_{j=1}^k (v_j^Th)^2  \right) &\textit{by Cauchy-Schwarz} \\
&\leq \sum_{j=1}^k \norm{\textit{Proj}_{H}(v_j)}^2 \\
\end{align*}

Armed with this we restate Corollary <ref> and prove it:

Given  integers $m > k+1$, let $\delta$ be a positive real number such that $4\delta^2 < \frac{1}{m} \left(1 - \frac{d(\Psi_{|W}, \Phi)^2}{d(\Psi, \Phi)^2}\right).$ Let $X_1, ..., X_n$ be $n$ data points from $\Psi$. Let $n$,$p$ go to infinity in a way such that $p/n \rightarrow \gamma$ Given $\epsilon > 0$ there is a $\gamma_{\sigma,\epsilon}$, where $\gamma_{\sigma,\epsilon}$ is a constant depending on $\sigma$, $\epsilon$, such that if $\gamma \leq \gamma_{\sigma,\epsilon}$ then with high probability there does not exist a set of $k+1$ orthonormal unit vectors $v_1, ..., v_{k+1}$ such that

\[d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v_j^TX_i}}{n}, \Phi\right) \geq \sqrt{1-4\delta^2} d(\Psi,\Phi) + \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}. \]

To prove this we will use Theorem <ref> along with the trivial linear algebraic Lemma <ref>. We prove by contradiction. Suppose a orthonormal set $v_1, ..., v_{k+1}$ exists satisfying the hypothesis :

\[
d_{W_2}\left(\frac{\Sigma_{i=1}^n \delta_{v_{j}^TX_i}}{n}, \Phi\right) \geq \sqrt{1-4\delta^2} d(\Psi,\Phi) + \epsilon + \frac{C_{\sigma}}{\sqrt[4]{n}}.
\]

We can invoke Theoremm <ref> to get with high probability,
\[\norm{\textit{Proj}_{W}(v_j)} \leq 2\delta \frac{d(\Psi,\Phi)}{\sqrt{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}}.\]

Then let $G = \textit{Span}\{v_1, ..., v_{k+1}\}$. As As $\textit{dim}(G) +\textit{dim}(W) = \textit{dim}(V)+1$,
there exists a non-zero vector $s \in G \cap W$ such that $\norm{s} = 1$.
Invoking Lemma <ref> with $H = W$ and $g=h=s$, we get the contradiction

\[ 1 = \left(s^Ts\right)^2 \leq \sum_{j=1}^{k+1} \norm{\textit{Proj}_{H}(v_j)}^2 \leq \frac{4(k+1)\delta^2 d(\Psi,\Phi)^2}{d(\Psi, \Phi)^2 - d(\Psi_{|W}, \Phi)^2}  < 1. \]

§ PROOF OF LEMMA <REF>

To prove Lemma <ref>, we need a bound on the concentration function of subgaussian random variables. For a Borel set $A$, let $A^r$ denote the $r$-fattening of $A$:

\[
A_r = \{ x : d(x, A) < r\}.
\]

Let $\mu$ be a probability measure on $\R$. Let

\[
\alpha_{\mu}(r) = \{1 - \mu(A^r) : \mu(A) \ge 1/2\}, r > 0,
\]

denote the concentration function of $\mu$.

Let $\mu$ be a $\sigma$-subgaussian probability measure. Then there exist absolute constants $C, c > 0$ such that $\alpha_{\mu}(r) \le C e^{-\frac{cr^2}{\sigma^2}}$ for all $r > 0$.

Without loss of generality, we may assume that $\sigma = 1$. Choose $r_0$ such that $\mu((-r_0, r_0)^c) < \frac{1}{2}$. Then any $A$ such that $\mu(A) > \frac{1}{2}$ must intersect $(-r_0, r_0)$, for otherwise one would get $\mu(A) \le \mu((-r_0, r_0)^c) < \frac{1}{2}$. Take $x_0 \in A \cap (-r_0, r_0)$. Then one must have

\[
(-r, r) \subseteq x_0 + (-(r_0 + r), r_0 + r) \subseteq A^{r_0 + r}.
\]

Now, by subgaussianity, there exist constants $C_1, c_1 > 0$ such that $\mu((-r, r)^c) \le C_1 e^{-c_1 r^2}$ for all $r > 0$. Therefore

\[
1 - \mu(A^{r + r_0}) \le \mu((-r, r)^c) \le C_1 e^{-c_1 r^2} \le C e^{-c(r + r_0)^2},
\]

where the last inequality holds for some constants $C, c > 0$ for all large enough $r$, say $r > r_1$. (For example, one can take $c = \frac{1}{2}c_1, C = C_1 e^{\frac{1}{2}B r_0^2}$ and $r_1 = 2r_0$.) Thus for all $r > r_0 + r_1$, we have that $\alpha_{\mu}(r) \le C e^{-cr^2}$.

We can always increase the constant $C$ so that one has $\sup_{r \in (0, r_0 + r_1]}\alpha_{\mu}(r) \le C e^{-c(r_0 + r_1)^2}$. Then for any $r \le r_0 + r_1$,

\[
\alpha_{\mu}(r) \le \sup_{r \in (0, r_0 + r_1]}\alpha_{\mu}(r) \le Ce^{-c(r_0 + r_1)^2} \le Ce^{-cr^2}.
\]

We conclude that there exist absolute constants $C, c > 0$ such that $\alpha_{\mu}(r) \le Ce^{-cr^2}$ for all $r > 0$.

For clarity's sake we restate Lemma <ref>.

Let $\mu$ be a $\sigma$-subgaussian  measure. Let $\mu_n$ be the empirical measure formed from an i.i.d. sample of size $n$ from $\mu$. Then
\[
\bbP\big(\abs{d_{W_2}(\mu_n, \mu) - \bbE d_{W_2}(\mu_n, \mu)} \geq t \big) \leq Ce^{-\frac{cnt^2}{\sigma^2}}.
\]

for some absolute constants $C, c > 0$.

The proof is the same as the proof of Theorem 7.1 in [2], except that we replace their log-Sobolev assumption on $\mu$ with a subgaussianity assumption, which yields a stronger bound on the concentration function as in Lemma <ref>, which in turns gives us a tail bound of the form

\[
\bbP\big(\abs{d_{W_2}(\mu_n, \mu) - \bbE d_{W_2}(\mu_n, \mu)} \geq t \big) \le Ce^{-\frac{cnt^2}{\sigma^2}}
\]

for some absolute constants $C, c > 0$.


# Universal Quantum Tomography With Deep Neural Networks

Nhan T. Luu, Truong Cong Thang  Department of Computer Science and Engineering
The University of Aizu
Aizuwakamatsu, Japan
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Quantum state tomography is a crucial technique for characterizing the state
of a quantum system, which is essential for many applications in quantum
technologies. In recent years, there has been growing interest in leveraging
neural networks to enhance the efficiency and accuracy of quantum state
tomography. Still, many of them did not include mixed quantum state, since
pure states are arguably less common in practical situations. In this research
paper, we present two neural networks based approach for both pure and mixed
quantum state tomography: Restricted Feature Based Neural Network and Mixed
States Conditional Generative Adversarial Network, evaluate its effectiveness
in comparison to existing neural based methods. We demonstrate that our
proposed methods can achieve state-of-the-art results in reconstructing mixed
quantum states from experimental data. Our work highlights the potential of
neural networks in revolutionizing quantum state tomography and facilitating
the development of quantum technologies.

###### Index Terms:

quantum state tomography, neural network

##  I Introduction

A promising avenue for both research and technological advancements has been
presented thanks to progress in controlling capacity and manipulation of small
quantum systems. These systems offer potential applications in various fields
such as quantum information processing and computation[1, 2, 3], quantum
chemistry simulations[4, 5, 6], secure communication[7, 8], among with others
notable works[4, 9, 10]. In recent breakthroughs, the successful demonstration
of a 53-qubit quantum computer accomplishing a computational task in a
fraction of the time expected on a classical supercomputer[11] underscore the
remarkable speedup achievable through quantum systems. This acceleration is
attributable in part to the exponentially vast state space available for
storing and manipulating information in quantum systems[12, 13, 14].
Nonetheless, while the expansive state space offers opportunities, it also
presents challenges in accurately characterizing and describing these systems.

Quantum State Tomography (QST) [15, 16, 17] aims to ascertain the unknown
quantum state by conducting measurements on a finite set of identical copies
of the system. If the state is characterized by the density matrix $\varrho$,
residing in a $d$-dimensional Hilbert space, approximately $O(d/\varepsilon)$
copies are necessary to achieve an estimate of $\varrho$ with an error
(measured as total variation distance) less than $\varepsilon$[18]. This
underscores the considerable resource demands of QST for large-scale systems.

In a broader context, QST can be viewed as an inverse problem[19, 20].
Consequently, linear inversion[21] is often considered the most
straightforward approach. However, it has drawbacks, such as the potential to
yield a non-physical state and the inability to analytically determine the
mean squared error bound of the estimate. To address these limitations,
several alternative QST methods have been developed, including maximum-
likelihood estimation (MLE)[22, 23], Bayesian tomography[24, 25], compressed
sensing techniques[26, 27], and matrix-product states[28, 29], where MLE
remains the most commonly employed approach.

In recent times, machine learning techniques have been employed in QST,
showing promising outcomes[30, 31, 32]. Specifically, generative models[33,
34, 35], often in the form of restricted Boltzmann machines (RBMs), have
emerged as effective Ansätze with minimal parameters to depict a quantum state
and understand the probability distribution of anticipated outputs[36, 30,
37]. Additionally, deep neural networks have been utilized in QST, enabling
physicists to harness the swift advancements in machine learning
methodologies.

In our research, we introduce QST with 2 type of neural networks (NNs)
architecture: Restricted Feature Based Neural Network (RFB-Net) and Mixed
States Conditional Generative Adversarial Network (MS-CGAN). RFB-Net is a
modified QST-NN[38] with modified regression modules that can classify quantum
states into specific specific class and pointing out its features of the
quantum state within certain constraints, finally using these features and
classes to infer the full state of the system. On the other hand, MS-CGAN is
an architecture derived from QST conditional-GAN (QST-CGAN)[39], which contain
certain modifications from the original work to extend its reconstruction
range to include mixed states. This network is improved for a better
estimation of the density matrix of the system while extend reconstruction
range to multiple type of quantum states. Both of these techniques can
accurately construct many types of given simulated measurement and shown great
promise in practical intermediate-scale quantum systems.

##  II Experimental data

By the assisstance of QuTiP[40], our generated QST dataset include (but not
limited to) 10000 Husimi Q measurements of 32-qubit Coherent states, Fock
states, Thermal states, Cat states, Num states[41, 42], Binomial states[41,
42] and GKP states[41, 43, 44] (number of samples uniformly distributed
between states), along with the original density matrices, labels and
generated values based on certain constraints. Using this diverse set of
states in our dataset, we are able to extract many types of information and
features .

Typically, these states exist within a Hilbert space that has infinite
dimensions. However, we can create a limited representation of these states
with finite dimensions by setting a threshold on their energy levels, which is
why in our experiments a Hilbert-space cutoff of $N_{c}=32$ is used for all
cases. To prevent artefacts resulting from truncation, the maximum photon
number of the states is limited to 16 after the displacements are applied.
Different types of quantum optical states used in the research are defined
based on the original paper of QST-NN[38], including three well-known basic
classes and four states from bosonic codes designed for quantum error
correction.

###  II-A Fock state

Eigenstates of the Fock basis:

$\ket{\psi_{fock}}=\ket{n_{photon}}$ (1)

Figure 1: Samples of Fock state

###  II-B Coherent state

Displaced vacuum states, characterized by the complex displacement amplitude
$\alpha$:

$\begin{split}&\ket{\psi_{coherent}(\alpha)}=\ket{\alpha}=D(\alpha)\ket{0},\\\
&D(\alpha)=exp(\alpha\times a^{{\dagger}}-\alpha\times a)\end{split}$ (2)

$D(\alpha)$ is the displacement operator where $a$ is the annihilation and
$(a^{{\dagger}})$ is the creation operator of the bosonic mode.

Figure 2: Samples of Coherent state

###  II-C Thermal state

Mixed states where the photon number distribution follows super-Poissonian
statistics:

$\rho_{thermal}(n^{th})=\sum_{n=0}^{N_{c}-1}\frac{1}{n^{th}+1}(\frac{n^{th}}{n^{th}+1})^{n}\ket{n}\bra{n}$
(3)

Figure 3: Samples of Thermal state

###  II-D Cat state

Bosonic-code states consisting of superpositions of coherent states up to a
normalization $N$ with projections $\Pi_{r}$ given by:

$\begin{split}&\ket{\psi_{cat}^{\mu}}=\frac{1}{N}\Pi_{(S+1)\mu}\\{\ket{\alpha
e^{i\frac{\pi}{S+1}k}}\\}^{2S+1}_{k=0},\\\
&\Pi_{r}=\sum^{\infty}_{m=0}\ket{2m(S+1)+r}\bra{2m(S+1)+r}\end{split}$ (4)

Figure 4: Samples of Cat state

###  II-E Num state

Specific set of bosonic-code states, consisting of superpositions of a few
Fock states, numerically optimized for quantum error correction, and
characterized by their average photon number $\overline{n}$

Figure 5: Samples of Num state

###  II-F Binomial state

Bosonic-code states constructed from a superposition of Fock states weighted
by the binomial coefficients:

$\ket{\psi_{binomial}^{\mu}}=\frac{1}{\sqrt{2^{N+1}}}\sum^{N+1}_{m=0}(-1)^{\mu
m}\sqrt{\begin{pmatrix}N+1\\\ m\end{pmatrix}}\ket{(S+1)m}$ (5)

Figure 6: Samples of Binomial state

###  II-G GKP state

Finite Gottesmann-Kitaev-Preskill states, limits the lattice and adds a
Gaussian envelope to make the state normalizable:

$\ket{\psi_{GKP}^{\mu}}=\sum_{\alpha\in
K(\mu)}e^{-\delta^{2}\absolutevalue{\alpha}^{2}}e^{-iRe[\alpha]Im[\alpha]}\ket{\alpha},$
---
$K(\mu)=\sqrt{\frac{\pi}{2}}(2n_{1}+\mu)+i\sqrt{\frac{\pi}{2}}n_{2}$
(6)

Figure 7: Samples of GKP state Name |  Constraint
---|---
Fock state |  $1\leq n_{photon}\leq 16$
Coherent state |  $10^{-6}\leq\absolutevalue{\alpha}^{2}\leq 3$
Thermal state |  $1\leq n_{th}\leq 16$
Cat state |  $S\in[0,2],\absolutevalue{\alpha}\in[1,3],r\in[0,2S+1]$
Num state |  $\overline{n}\in\\{1.56,2.67,2.77,4.15,4.34\\}$
Binomial state |  $2\leq N\leq N_{C}/(S+1)-1$
GKP state |  $n_{1},n_{2}\in\\{-20,20\\},\delta\in[0.2,0.5]$
TABLE I: Constraints for states in dataset

##  III Architecture and training settings

###  III-A Restricted Feature Based Neural Network

We utilizes the PyTorch[45] framework to implement a deep convolutional neural
network (CNN) similar to QST-NN. In our proposed architecture, we aim to
extract relevant information from input measurements that are provided in the
form of 32x32 density matrices. To accomplish this, we employ a series of 6
convolutional heads, which are interspersed with Gaussian noise, dropout, and
leaky ReLU activation functions to extract meaningful features from the input.

Subsequently, we feed the extracted features into both a classification tail
and a regression tail. The classification tail is responsible for predicting
the label, while the regression tail computes three features of the state. To
generate the final state prediction, we sum the predicted label output with
the extracted features and utilize the regression tail to predict the three
essential features of the state.

These three features are then fed into a specialized _Reconstructor_ module,
which is designed explicitly to reconstruct the state by taking 3-length value
vectors as input and generating the predicted states.

Figure 8: Architecture of RFB-Net.

We train our model using a linear combination of cross-entropy loss for the
label and the sum of MAE loss for the real and imaginary parts of the
regression output _F_ as a feature vector:

$loss=-\sum_{c=1}^{N}label_{real,N}\log(label_{pred,N})+$
---
$\sum_{i=1}^{N}\absolutevalue{F^{real}_{out}-F^{real}_{pred}}+\sum_{i=1}^{N}\absolutevalue{F^{imag}_{out}-F^{imag}_{pred}}$
(7)

In general, the RFB-Net architecture is capable of not only reconstructing
states from input measurements but also reconstructing the state from input
measurements, which can be viewed as a form of multitasking architecture. In
this architecture, the decoder produces both output labels and reconstructed
operators. The ability to perform multiple tasks simultaneously is a key
feature of multitasking architectures, which are increasingly being used in
machine learning applications. By combining multiple tasks into a single
architecture, multitasking models can improve efficiency and accuracy, while
reducing the computational overhead associated with running separate models
for each task. RFB-Net’s ability to perform both state reconstruction and
label prediction highlights its potential as a versatile and efficient
multitasking architecture for quantum information applications. Furthermore,
the ability to reconstruct the state from input measurements is particularly
useful in quantum information processing, where it is essential to accurately
recover the state of a quantum system from measurements. Overall, RFB-Net’s
multitasking capabilities make it a promising architecture for a wide range of
quantum information tasks.

###  III-B Mixed State Conditional Generative Adversarial Network

Building upon the seminal work of QST-CGAN, we have developed a novel neural
architecture for state prediction, utilizing the powerful and versatile
TensorFlow [46] framework. Our proposed model is designed to accommodate a
batch of measurements in the form of 32x32 density matrices, as well as label
vectors representing the target state. To facilitate efficient processing, we
have implemented a series of pre-processing steps, which include flattening
the measurements to a 1024-length vector and one-hot encoding the label to a
7-length vector.

Subsequently, we feed the label vector through a dense layer, which enables us
to encode it into a 1024-length vector that can be mapped to the same size as
the measurements. By adding the resulting vector and the flattened
measurements, we generate a final feature that captures the essence of the
underlying quantum state. This final feature is then fed through a series of
convolutional layers, batch normalization, and leaky ReLU activation,
resulting in a batch of 2x32x32 tensor, where the two channels correspond to
the real and imaginary components of the quantum state.

Figure 9: Architecture of MS-CGAN.

Following the generation of the complex tensor batch, we proceed to further
process the data by passing it through a _DensityMatrix_ layer, which is
responsible for normalizing the quantum state in accordance with specific
constraints. Subsequently, we pass the normalized state through an
_Expectation_ layer, which is designed to measure the state by calculating the
expected values of various observables.

Unlike the original architecture, which only utilizes the lower half of the
result to obtain a Hermitian matrix, our novel approach simulate a Cholesky
decomposition [22] to the upper half of the result, thereby optimizing the
number of parameters utilized in the process. By subtracting the two obtained
matrices and normalizing the result, we are able to generate a density
operator that encompasses positive semi-definite, negative semi-definite, and
even indefinite outcomes.

This approach is particularly relevant in situations where the creation of an
odd cat state is required, and can be achieved through the subtraction of two
coherent states $\ket{\alpha}$ and $\ket{-\alpha}$:

$\ket{\psi_{cat}^{\mu}}=\frac{1}{\sqrt{2(1-e^{-2\absolutevalue{\alpha}^{2}})}}(\ket{\alpha}-\ket{-\alpha})$
(8)

Figure 10: Demonstration of the improved _DensityMatrix_ process.

Model’s loss function is a linear combination of MAE loss of input measurement
$M_{in}$ and the measurement from created density matrix $M_{pred}$ with the
weighted sum (with a scalar $\alpha=100$) of MAE loss for the real and
imaginary parts ($D^{real}$ and $D^{imag}$) of the density operator:

$loss=\sum_{i=1}^{N}\absolutevalue{M_{in}-M_{pred}}+$
---
$\alpha\times(\sum_{i=1}^{N}\absolutevalue{D^{real}_{true}-D^{real}_{pred}}+$
$\sum_{i=1}^{N}\absolutevalue{D^{imag}_{true}-D^{imag}_{pred}})$
(9)

The model is optimized for 100 epochs on dataset using Adam[47] with $lr=1e-4$
and $\beta=(0.9,0.999)$. Our experimental results show that our proposed
architecture has the potential to cover every quantum state that can be
represented as a density operator, which is a Hermitian matrix with unit
trace. This makes our architecture a more optimized solution for state
prediction problems in various applications, as it can handle a wider range of
quantum states.

###  III-C Accuracy metrics

In both of our experiment, we will be using quantum fidelity to measure the
similarity between two quantum states. Quantum fidelity is a key concept in
quantum information theory and plays a critical role in many quantum
information protocols, including quantum error correction, quantum
teleportation, and quantum cryptography. Given two density matrices $\rho$ and
$\sigma$, the fidelity is defined by :

$F(\rho,\sigma)=(tr(\sqrt{\rho\sigma}))^{2}$ (10)

where if $F(\rho,\sigma)==1$, the two states are identical. If the fidelity is
$0\leq F(\rho,\sigma)<1$ , the two states differ, and the degree of difference
is quantified by the fidelity value.

##  IV Result

Within this section, we shall expound upon the findings of our investigation,
in which we leveraged two distinct neural network architectures for quantum
tomography. Our analysis of these models yielded remarkable outcomes, as both
demonstrated adequate performance in reconstructing the states of quantum
systems.

The experimentation of both models are conducted on a single NVIDIA GeForce
RTX 3080 GPU.

###  IV-A Restricted Feature Based Neural Network

Figure 11: Training loss progression of RFB-Net.

Figure 12: Validation fidelity progression per 100 iterations of RFB-Net.

RFB-Net were initially trained for 100 epochs using Adam[47] optimizer with
$lr=1e-4$ and $\beta=(0.9,0.999)$. Although the model have already converge
majorly post training with this setting, we want to explore how long would it
take for the model to completely converge close to state of the art result.
After throughout examining, we see that it took at least around 1000 epochs
for model to converge close to state of the art result.

Figure 12 demonstrate that the model’s degree of fidelity converges to a
noteworthy value within the initial 100 epochs. However, to achieve an even
higher degree of fidelity, it is necessary to engage in a protracted training
process of over 1300 epochs, which can be considerably time-consuming. The
probable cause behind this gradual convergence is due to the model’s need for
a more extensive fine-tuning process, which would enable it to make accurate
predictions of vectors that contain small values, such as when predicting the
mean photon number of a coherent state. This implies that the model’s ability
to make precise predictions is contingent on the length of the fine-tuning
process, and a lengthier process would enable the model to offer better
predictions.

To attain fidelity of the proposed architecture, the deployment of the
_Reconstruction_ module is imperative to restore the state from the predicted
values. Given the significant computational expense incurred during the
process of reconstruction, a judicious decision was made to evaluate the
architecture’s fidelity after every 100 epochs. This approach allows for a
systematic assessment of the performance of the architecture, while minimizing
the computational burden of reconstructing the state at every epoch.

###  IV-B Mixed State Conditional Generative Adversarial Network

Figure 13: Training loss progression of MS-CGAN.

Figure 14: Validation fidelity progession of MS-CGAN during training process.

The MS-CGAN architecture exhibits a distinct characteristic when compared to
the RFB-Net architecture, in that it is capable of converging to desired
values within a training period consisting of no more than 100 iterations
(Figure 13 and 14), despite requiring an extra label input in comparison to
the RFB-Net architecture. Notably, the MS-CGAN architecture is primarily
focused on reconstruction, but this does not detract from its robustness,
given its swift convergence when trained on relevant data. It is therefore
reasonable to conclude that the MS-CGAN architecture is a viable alternative
to the RFB-Net architecture, given its impressive performance characteristics
in terms of convergence speed and accuracy, despite having distinct
differences in terms of input requirements and model focus.

###  IV-C Compare with existing methods

Figure 15: Validation fidelity progression of QST-CGAN on generated dataset.

Figure 16: Samples of QST-CGAN predictions. Model name |  Fidelity
---|---
RFB-Net |  0.97
MS-CGAN |  0.91
QST-CGAN |  0.19
TABLE II: Validation fidelity of investigated models on generated dataset.

In order to ensure a more robust and comprehensive comparison between our
proposed model and the original model, we have undertaken a meticulous and
rigorous experimental procedure that involves the use of a QST-CGAN model
trained on the same dataset for a prolonged period of 100 epochs, utilizing
the original configuration parameters and settings as described in the seminal
paper. This approach is intended to provide a more nuanced and detailed
analysis of the comparative performance of the two models, by subjecting both
models to identical experimental conditions, and enabling us to assess their
relative strengths and weaknesses in a more systematic and rigorous manner.

Upon conducting this extensive analysis, we have observed that the QST-CGAN
model, while exhibiting notable strengths and advantages in certain areas, has
encountered certain challenges and limitations in the reconstruction of
specific batches of states, resulting in an unstable and erratic convergence
on the dataset. These limitations are reflected in the average fidelity score,
which is found to be hovering around a modest 0.19, indicating a suboptimal
and unsatisfactory level of accuracy and precision. In contrast, our proposed
solution has proven to be more superior and effective in overcoming these
challenges, by demonstrating a higher degree of stability, robustness, and
accuracy, as well as a more consistent and reliable performance.

###  IV-D Benchmarking with noise

The presence of noise, which refers to any unwanted fluctuations or
disturbances in the experimental setup, is an inescapable factor that pervades
most scientific investigations, and consequently poses significant challenges
for accurate and reliable data analysis. Given this ubiquitous nature of
noise, it is imperative that any method employed for state reconstruction is
designed and optimized to be robust and resilient against various types of
noise, to ensure that the final results are valid and trustworthy.

Moreover, it is essential to recognize that noise can infiltrate the system at
multiple stages, leading to a potentially complex and nuanced interplay of
factors that can affect the experimental outcome. Specifically, noise can
manifest itself during the preparation of the state to be reconstructed, where
errors in the preparation procedure can introduce small deviations in the
actual state of interest, resulting in a distorted or altered version of the
target state. Additionally, the measurement protocol itself can be a source of
noise, as any calibration errors or other external influences may cause us to
measure something other than what we intended to measure, leading to further
inaccuracies in the final results. Finally, noise can arise during the data
collection process itself, where amplification errors or photon shot noise can
corrupt the signal, leading to further challenges in data interpretation and
analysis. Therefore, a comprehensive and multifaceted approach is needed to
account for the various sources of noise and ensure that any method for state
reconstruction is adequately robust and resilient against these influences, to
provide accurate and reliable results that can be used to advance scientific
knowledge and understanding.

(a) Mixed state noise

(b) Photon loss noise

(c) Pepper noise

Figure 17: Samples of noisy input and RFB-Net predictions.

(a) Mixed state noise

(b) Photon loss noise

(c) Pepper noise

Figure 18: Samples of noisy input and MS-CGAN predictions.

In this research, our approach involves the utilization of various noise
generating techniques, akin to the employed dataset, taken from QST-NN. In
particular, we have employed three types of noise in our experiment, which
include additive mixed state noise, photon loss noise, and pepper noise.
Rather than introducing the noise during the training phase, whereby the model
would gradually learn and adapt to it, we have opted to evaluate the
performance of two distinct architectures devoid of exposure to noisy data.
This is to ascertain whether our models can generalize well under these
conditions or not.

Upon analyzing the sample input and prediction, it is evident that the RFB-Net
exhibits a commendable performance in noisy environments in comparison to the
MS-CGAN with regards to fidelity, albeit the reconstructed operator showcasing
some dissimilarities from its original form. On the other hand, the MS-CGAN’s
performance is more stable, particularly when dealing with pepper noise added
samples, and shows superior feature extraction capabilities compared to RFB-
Net.

| Noise type
---|---
Model name |  Mixed state noise |  Photon loss noise |  Pepper noise
RFB-Net |  0.80 |  0.51 |  0.30
MS-CGAN |  0.49 |  0.48 |  0.60
TABLE III: Average fidelity of reconstructed noisy data

##  V Conclusion

In summary, this research paper explored the application of neural networks in
quantum tomography, and the findings suggest that they offer a promising
approach to solving the challenges associated with reconstructing quantum
states. The results demonstrate that neural network architectures, when
trained on appropriate datasets, can produce accurate and efficient
reconstructions of quantum states. Furthermore, the approach presented in this
paper is scalable and can be extended to large-scale quantum systems. It is
clear that the combination of neural networks and quantum tomography has great
potential for advancing the field of quantum information processing, and we
anticipate that this research will inspire further exploration and development
of these methods in the future.

## References

  * [1] R. P. Feynman, “Simulating physics with computers,” in _Feynman and computation_. CRC Press, 2018, pp. 133–153.
  * [2] A. Montanaro, “Quantum algorithms: an overview,” _npj Quantum Information_ , vol. 2, no. 1, pp. 1–8, 2016.
  * [3] G. Wendin, “Quantum information processing with superconducting circuits: a review,” _Reports on Progress in Physics_ , vol. 80, no. 10, p. 106001, 2017\.
  * [4] I. M. Georgescu, S. Ashhab, and F. Nori, “Quantum simulation,” _Reviews of Modern Physics_ , vol. 86, no. 1, p. 153, 2014.
  * [5] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” _nature_ , vol. 549, no. 7671, pp. 242–246, 2017.
  * [6] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, “Toward the first quantum simulation with quantum speedup,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 38, pp. 9456–9461, 2018.
  * [7] H. J. Kimble, “The quantum internet,” _Nature_ , vol. 453, no. 7198, pp. 1023–1030, 2008.
  * [8] J. Yin, Y.-H. Li, S.-K. Liao, M. Yang, Y. Cao, L. Zhang, J.-G. Ren, W.-Q. Cai, W.-Y. Liu, S.-L. Li _et al._ , “Entanglement-based secure quantum cryptography over 1,120 kilometres,” _Nature_ , vol. 582, no. 7813, pp. 501–505, 2020.
  * [9] J.-Q. You and F. Nori, “Atomic physics and quantum optics using superconducting circuits,” _Nature_ , vol. 474, no. 7353, pp. 589–597, 2011\.
  * [10] X. Gu, A. F. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, “Microwave photonics with superconducting quantum circuits,” _Physics Reports_ , vol. 718, pp. 1–102, 2017.
  * [11] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell _et al._ , “Quantum supremacy using a programmable superconducting processor,” _Nature_ , vol. 574, no. 7779, pp. 505–510, 2019.
  * [12] M. A. Nielsen and I. L. Chuang, _Quantum computation and quantum information_. Cambridge university press, 2010.
  * [13] C. M. Caves, I. H. Deutsch, and R. Blume-Kohout, “Physical-resource requirements and the power of quantum computation,” _Journal of Optics B: Quantum and Semiclassical Optics_ , vol. 6, no. 8, p. S801, 2004.
  * [14] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, “Supervised learning with quantum-enhanced feature spaces,” _Nature_ , vol. 567, no. 7747, pp. 209–212, 2019.
  * [15] G. M. D’Ariano, M. G. Paris, and M. F. Sacchi, “Quantum tomography,” _Advances in imaging and electron physics_ , vol. 128, no. 10.1016, pp. S1076–5670, 2003.
  * [16] Y.-x. Liu, L. Wei, and F. Nori, “Tomographic measurements on superconducting qubit states,” _Physical Review B_ , vol. 72, no. 1, p. 014547, 2005.
  * [17] A. I. Lvovsky and M. G. Raymer, “Continuous-variable optical quantum-state tomography,” _Reviews of modern physics_ , vol. 81, no. 1, p. 299, 2009.
  * [18] R. O’Donnell and J. Wright, “Efficient quantum tomography,” in _Proceedings of the forty-eighth annual ACM symposium on Theory of Computing_ , 2016, pp. 899–912.
  * [19] A. Tarantola, _Inverse problem theory and methods for model parameter estimation_. SIAM, 2005.
  * [20] R. C. Aster, B. Borchers, and C. H. Thurber, _Parameter estimation and inverse problems_. Elsevier, 2018.
  * [21] B. Qi, Z. Hou, L. Li, D. Dong, G. Xiang, and G. Guo, “Quantum state tomography via linear regression estimation,” _Scientific reports_ , vol. 3, no. 1, p. 3496, 2013.
  * [22] K. Banaszek, G. D’ariano, M. Paris, and M. Sacchi, “Maximum-likelihood estimation of the density matrix,” _Physical Review A_ , vol. 61, no. 1, p. 010304, 1999.
  * [23] A. I. Lvovsky, “Iterative maximum-likelihood reconstruction in quantum homodyne tomography,” _Journal of Optics B: Quantum and Semiclassical Optics_ , vol. 6, no. 6, p. S556, 2004.
  * [24] R. Blume-Kohout, “Optimal, reliable estimation of quantum states,” _New Journal of Physics_ , vol. 12, no. 4, p. 043034, 2010.
  * [25] C. Granade, J. Combes, and D. Cory, “Practical bayesian tomography,” _new Journal of Physics_ , vol. 18, no. 3, p. 033024, 2016.
  * [26] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, “Quantum state tomography via compressed sensing,” _Physical review letters_ , vol. 105, no. 15, p. 150401, 2010.
  * [27] D. Ahn, Y. S. Teo, H. Jeong, F. Bouchard, F. Hufnagel, E. Karimi, D. Koutnỳ, J. Řeháček, Z. Hradil, G. Leuchs _et al._ , “Adaptive compressive tomography with no a priori information,” _Physical Review Letters_ , vol. 122, no. 10, p. 100404, 2019\.
  * [28] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, “Efficient quantum state tomography,” _Nature communications_ , vol. 1, no. 1, p. 149, 2010.
  * [29] B. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. Buyskikh, A. Daley, M. Cramer _et al._ , “Efficient tomography of a quantum many-body system,” _Nature Physics_ , vol. 13, no. 12, pp. 1158–1162, 2017.
  * [30] G. Carleo and M. Troyer, “Solving the quantum many-body problem with artificial neural networks,” _Science_ , vol. 355, no. 6325, pp. 602–606, 2017.
  * [31] Q. Xu and S. Xu, “Neural network state estimation for full quantum state tomography,” _arXiv preprint arXiv:1811.06654_ , 2018.
  * [32] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, “Neural-network quantum state tomography,” _Nature Physics_ , vol. 14, no. 5, pp. 447–450, 2018.
  * [33] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, “Reconstructing quantum states with generative models,” _Nature Machine Intelligence_ , vol. 1, no. 3, pp. 155–161, 2019.
  * [34] S. Lloyd and C. Weedbrook, “Quantum generative adversarial learning,” _Physical review letters_ , vol. 121, no. 4, p. 040502, 2018.
  * [35] M. Kieferová and N. Wiebe, “Tomography and generative training with quantum boltzmann machines,” _Physical Review A_ , vol. 96, no. 6, p. 062327, 2017.
  * [36] I. Glasser, N. Pancotti, M. August, I. D. Rodriguez, and J. I. Cirac, “Neural-network quantum states, string-bond states, and chiral topological states,” _Physical Review X_ , vol. 8, no. 1, p. 011006, 2018.
  * [37] E. S. Tiunov, V. Tiunova, A. E. Ulanov, A. Lvovsky, and A. K. Fedorov, “Experimental quantum homodyne tomography via machine learning,” _Optica_ , vol. 7, no. 5, pp. 448–454, 2020.
  * [38] S. Ahmed, C. Sánchez Muñoz, F. Nori, and A. F. Kockum, “Classification and reconstruction of optical quantum states with deep neural networks,” _Physical Review Research_ , vol. 3, no. 3, p. 033278, 2021.
  * [39] S. Ahmed, C. Sanchez Munoz, F. Nori, and A. F. Kockum, “Quantum state tomography with conditional generative adversarial networks,” _Physical review letters_ , vol. 127, no. 14, p. 140502, 2021.
  * [40] J. R. Johansson, P. D. Nation, and F. Nori, “Qutip: An open-source python framework for the dynamics of open quantum systems,” _Computer Physics Communications_ , vol. 183, no. 8, pp. 1760–1772, 2012.
  * [41] V. V. Albert, K. Noh, K. Duivenvoorden, D. J. Young, R. Brierley, P. Reinhold, C. Vuillot, L. Li, C. Shen, S. M. Girvin _et al._ , “Performance and structure of single-mode bosonic codes,” _Physical Review A_ , vol. 97, no. 3, p. 032346, 2018.
  * [42] M. H. Michael, M. Silveri, R. Brierley, V. V. Albert, J. Salmilehto, L. Jiang, and S. M. Girvin, “New class of quantum error-correcting codes for a bosonic mode,” _Physical Review X_ , vol. 6, no. 3, p. 031006, 2016.
  * [43] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf _et al._ , “Quantum error correction of a qubit encoded in grid states of an oscillator,” _Nature_ , vol. 584, no. 7821, pp. 368–372, 2020.
  * [44] D. Gottesman, A. Kitaev, and J. Preskill, “Encoding a qubit in an oscillator,” _Physical Review A_ , vol. 64, no. 1, p. 012310, 2001.
  * [45] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [46] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “$\\{$TensorFlow$\\}$: a system for $\\{$Large-Scale$\\}$ machine learning,” in _12th USENIX symposium on operating systems design and implementation (OSDI 16)_ , 2016, pp. 265–283.
  * [47] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.