# MRI Banding Removal via Adversarial Training

Aaron Defazio Tullie Murrell
Facebook AI Research,
Facebook New York
Michael P. Recht
Department of Radiology
NYU Grossman School of Medicine

###### Abstract

MR images reconstructed from sub-sampled Cartesian data using deep learning
techniques show a characteristic banding (sometimes described as streaking),
which is particularly strong in low signal-to-noise regions of the
reconstructed image. These unnatural artifacts have been identified as one of
the largest obstacles preventing the clinical use of machine-learning based
MRI reconstructions. In this work, we propose the use of an adversarial loss
that penalizes banding structures without requiring any human annotation. Our
technique greatly reduces the appearance of banding, without requiring any
additional computation or post-processing at reconstruction time. Our approach
is compatible with any existing reconstruction approach that uses supervised
machine learning, including the current state-of-the-art. We report the
results of a blind comparison against a strong baseline by a group of expert
evaluators (board-certified radiologists), where our approach is ranked
superior at banding removal with no statistically significant loss of detail.
A reference implementation of our method is available in the supplementary
material.

##  1 Introduction

The use of deep-learning approaches for accelerating MR imaging has recently
shown significant promise [7, 19, 1], with learning approaches far out-
performing classical penalized least squares approaches for reconstructing
images from raw subsampled k-space signals. However, existing approaches
produce images with some unnatural structures that prevent radiologists from
accepting the images for clinical use despite the advantages the images have
over classical techniques with respect to reconstruction accuracy metrics such
as the structured similarity metric (SSIM).

Figure 1: A deep-learning MRI reconstruction that shows significant horizontal
banding artifacts

In this work, we describe a method for removing the primary artifact produced
by Cartesian deep-learning reconstruction systems: banding. This banding, as
illustrated in Figure 1, is characterized by a streaking pattern exactly
aligned with the phase-encoding direction (horizontal in the figure). This
banding is anisotropic and non-homogenous across the image, and most visible
in high-noise or low-contrast areas. It is the result of the signal
subsampling process used during Cartesian accelerated MRI, whereby subsampling
occurs in one spatial direction only.

At the 2019 Medical Imaging meets NeurIPS workshop [11], banding artifacts
were shown to occur in reconstructions from each of the top 3 winning teams in
the competition, despite the significant differences between the
reconstruction approaches. These artifacts were identified as a major obstacle
to the use of machine-learning based reconstructions in clinical practice.

Our paper is structured as follows: In Section 2 we formally describe the MRI
reconstruction problem as it applies to current clinical MRI scanners. In
Section 2 we describe how to augment standard deep-learning based MRI
reconstruction methods with our orientation adversary. Section 4 describes the
state-of-the-art reconstruction model that we use in our experiments. In
Section 5 we describe the masking procedure we use, in Section 3 we describe
the classical baseline that we compare against, and in Section 7 we detail how
our model was trained. Finally in Section 8 we detail the results of a blind
evaluation by radiologists.

##  2 Accelerated Parallel 2D MRI

In MR imaging, a spatial image is produced by combining measurements of the
anatomy acquired in the Fourier domain, known as k-space. Let $m$ be the true
greyscale spatial image. Classical approaches produce an image of the anatomy
by acquiring a full cartesian grid of samples from k-space, then applying the
inverse fast Fourier transform. In our notation the estimated greyscale image
$\hat{m}$ is:

$\hat{m}=\mathcal{F}^{-1}(x)$

where $x$ is a $h\times w$ matrix of k-space measurements. In this work, we
focus on 2D MRI images produced by a modern MR imaging system, which contain
two additional complications: parallelization and acceleration.

###  2.1 Parallel Imaging

In a parallel MRI system, more than one receiver coil is used, resulting in a
tensor of acquired k-space images. Instead of each coil imaging the entire
field-of-view, each coil covers a smaller portion of the anatomy. The signal
acquired by coil $i$ of $n_{c}$ coils is given by a Fourier transform
$\mathcal{F}$:

$x_{i}=\mathcal{F}(s_{i}\circ m)+\text{noise},$ (1)

where $s$ is a complex-valued coil-sensitivity map that is applied element-
wise. The $s_{i}$ values can be estimated via well-known auto-calibration
procedures, but this is not necessary for all deep-learning approaches.

The coil signals are commonly combined using the root-sum-squares (RSS)
procedure to produce a spatial image. The RSS estimate at pixel $l,m$ is given
by:

$m_{\text{RSS},lm}=\sqrt{\sum_{i=1}^{n_{c}}\left|m_{i,lm}\right|^{2}},$

where for each coil $m_{i}=\mathcal{F}^{-1}(x_{i})$ is the individual coil
image. The RSS estimate produces images with slightly (often negligibly)
higher noise than more sophisticated approaches but has the advantage of
robustness and simplicity.

###  2.2 Accelerated imaging

Accelerated MRI systems capture a subset of the full k-space system to reduce
scan time. If parallel imaging is used, the system of linear equations given
by Equation 1 is overdetermined even when a subset of the $x_{i}$ pixel values
are known as long as the sub-sampling factor is less than the number of coils,
and so a least-squares solve may be used to produce a spatial image, at the
expense of an increase in noise over non-accelerated MRI [5, 18]. Two-fold
subsampling is widely used in clinical practice, typically using either the
SENSE [18] or GRAPPA [5] classical approaches to reconstruction. Higher
acceleration factors can be achieved using regularized least-squares in the
case of sparse anatomy such as vascular MRI [12], but these approaches produce
poor results for general-purpose MR imaging.

###  2.3 Machine learning approaches to accelerated 2D MRI

In the machine learning approach to 2D MRI reconstruction, a training set of
$n_{\text{data }}$ instances (slices) is gathered, where each instance is a
$k$-space tensor $x^{(j)}:n_{c}\times h\times w$. Then standard parallel
imaging is used to produce spatial images $m^{(j)}$. A scan of a patient
consists of multiple spatially consecutive slices with different scan
modalities of the same anatomy, although for training purposes we treat the
slices independently and sample from the total set of slices across all
patients in the dataset i.i.d.

These training pairs are then used to train a black box predictor $B_{\phi}$,
which maps from raw subsampled k-space tensors (complex valued multichannel
Fourier domain images) to greyscale $h\times w$ spatial images. Given an
image-space loss function $l$, the training loss for datapoint $j$ is:

$L^{(j)}(\phi)=l\left(B_{\phi}\left(M\left(x^{(j)}\right)\right),\,m^{(j)}\right)$

where $M$ is a masking function that zeros out a fraction of the $k$-space
lines. We detail this masking function further in the Section 5. Essentially,
the model $B$ is trained to produce an image as close as possible to the
“ground-truth” fully-sampled parallel MRI as possible, using only a fraction
of the data following the standard empirical risk minimization setup. In our
experiments, we used a combination of SSIM and L1 losses with weighting 0.01
for the L1 component.

##  3 Orientation Adversary

Figure 2: Orientation adversarial training

The primary difficulty with removing banding is the lack of available
annotation; direct supervised machine learning techniques can not be used. Our
insight is to use image orientation during reconstruction as a self-supervised
learning signal. Image banding is aligned with the direction of subsampling in
the mask, in our case horizontally. This signal can be used by an adversarial
training term to identify and penalize banding.

In particular, let $A_{\theta}$ be an adversary model, which maps from spatial
images $h\times w$ to real values $[0,1]$. Its goal is to predict if the given
spatial image contains horizontal $(0$) or vertical $(1)$ banding. We train
this adversary simultaneously with the predictor model, using the same
minibatch to compute stochastic gradient steps simultaneously for both, rather
than in an alternating fashion.

#### Predictor

The training of our prediction model is modified as follows. Before applying
$B$, a random operator is sampled using a Bernoulli variable $r^{(j)}$ with
probability 0.5, either a random flip $R_{1}$ (transpose $h\longleftrightarrow
w$) or the identity operator $R_{0}$. This operator is then applied before and
after the application of $B$:

$\hat{m}^{(j)}=\left(R_{r^{(j)}}\circ B_{\phi}\circ M\circ
R_{r^{(j)}}\right)x^{(j)}.$

The predictor’s loss is then augmented with a term that encourages it to
produce images that fool the adversary:

$L_{B}^{(j)}(\phi)=l\left(\hat{m}^{(j)},\,m^{(j)}\right)+\text{CE}\left(1-r^{(j)},A_{\theta}\left(\hat{m}^{(j)}\right)\right),$

where CE is the binary cross-entropy. Gradients are not propagated to the
adversaries parameters $\theta$ during the predictor update, but are fully
propagated through the adversary from output to input then through the
predictor to its parameters $\phi$. This is accomplished by toggling
requires_grad_ to false for the adversaries parameters in Pytorch during loss
calculation, then toggling it back before the adversary loss calculation.

#### Adversary

The adversary is trained to predict the $r$ variable directly, with the
addition of a regularization term:

$L_{A}^{(j)}(\theta)=\text{CE}\left(r^{(j)},A_{\theta}\left(\hat{m}^{(j)}\right)\right)+\gamma\left\|\nabla_{\hat{m}^{(j)}}A_{\theta}(\hat{m}^{(j)})\right\|^{2}.$

Gradients are not propagated through the predictor during the adversary step,
using Pytorch’s detach operator on $\hat{m}^{(j)}$. We found that
regularization of the adversary is necessary for stable convergence. We use a
simplified gradient penalty as used by [15], which is closely related to WGAN-
GP penalty [6] which is widely used for Generative Adversarial Network (GAN)
training.

##  4 Models

Our technique is applicable using any predictor and adversary model
architecture. We detail the two models we used below. Full source code is
available in the supplementary material.

#### Predictor

We used the state-of-the-art predictor architecture from [20] consisting of a
sequence of 12 U-Net models, interleaved with a soft-projection onto the known
k-space lines. This sequence is followed by an inverse Fourier transform then
a root-sum-squares operation to produce the final image. Each U-Net has 12
channels after the first convolution and has 4 pooling operations (i.e.
convolutions occur at 5 resolutions when the input resolution is included).
This network has only 13 million parameters, but more than 250 convolutional
layers. This model is the current state-of-the-art for MRI reconstruction.

#### Adversary

The shallow ResNet consists of an initial convolution to increase the channel
count to 64, followed by 2 pre-activation ResNet basic blocks, 4x4 max-
pooling, 2 more blocks with 128 channels, then 4x4 max-pooling followed by
average pooling down to a 1x1 image, a ReLU then a linear layer. This
architecture was not heavily optimized, and we expect a better architecture
could give even better results. We found initially that a ResNet-50
architecture took too long to converge, and that it was necessary to include
ResNet blocks that apply at the full-resolution of the image, rather than
after a downsampling operation as performed in the standard ResNet
architecture. Because of the small batch-sizes used for training MRI
reconstruction models, we found that it was necessary to replace BatchNorm
with GroupNorm [21], and we used group-size 32 as our default.

##  5 Masking

In 2D Cartesian MR imaging, the two directions in the Fourier domain are known
as the Frequency and Phase encoding directions. The most common clinical
practice for accelerated MRI is to only capture every second frequency-
encoding line, except for a band of the lowest-frequency lines which are
always included. These lines are typically stored electrically with zero
values filling the unseen values. All coil images for a slice follow the same
mask as the lines are captured simultaneously during acquisition. The
acquisition process can rapidly acquire an entire line in the frequency
direction line at once, so typically no subsampling is performed within the
lines.

For the fastMRI dataset, the phase-encoding direction is always horizontal,
and we use a central band of 16 low-frequency lines. These lines are only
consecutive in k-space if the Fourier domain is viewed with the low-frequency
lines in the center of the image, which is NOT the standard layout returned by
fast-Fourier routines. The low-frequency lines serve multiple purposes [5,
18]:

  * •

They allow estimation of the coil sensitivities through auto-calibration,

  * •

they contain a significant fraction of the total signal energy.

  * •

They allow easier disambiguation of aliasing produced by the mask.

Machine learning models trained using the same number of total lines but
without a central region do not produce images of as high quality. The
disambiguation function is important. A direct IFFT of a masked k-space coil
image produces two images of the anatomy, with one shifted by half the image
width. This is known as aliasing. Without additional k-space lines, the signal
does not contain information indicating which of the two positions in the
image the anatomy is in. The central low-frequency k-space lines contain
information allowing a perfect reconstruction of a low-resolution version of
the full coil image, which the machine learning model can use to disambiguate
the two positions.

For our accelerated reconstruction we performed 4-fold acceleration by
sampling every 4th frequency line, using an offset equispaced mask following
[3], with the central lines as discussed above. This level of acceleration
produces significant artifacts when used with classical accelerated MRI
techniques, but produces very high-quality results when using deep-learning
reconstructions. We also experimented with the random masks preferred by the
theory of compressed sensing; however, we found they performed significantly
worse. This is consistent with prior work on machine learning reconstructions
[20, 3]. We attribute this to the lack of sparsity required by compressed
sensing theory for the Proton Density (PD) and PD Fat Saturated (PDFS) MRI
images used in our experiments.

##  6 Baselines

Figure 3: Baseline (top), with varying levels of blurring, against our
adversarial method (bottom). Columns show different levels of adaptive noise.
Images are best viewed on a high-brightness monitor at full resolution.

As the banding problem is relatively new we are not aware of any existing
baseline methods for banding removal specific to MR imaging. We will instead
apply a classical image processing method for minimizing the visibility of
banding artifacts in images known as dithering [4, 8]. Dithering is simply the
process of adding noise in a specific way to improve the perceptual quality of
an image. We found that adding noise directly to the reconstruction improved
apparent sharpness, but did not by itself reduce the visibility of banding
(Figure 3). However, if an anisotropic blurring filter is used before the
addition of noise, banding was well-suppressed. We used the following
convolutional kernel parametrized by $\alpha$:

$K=\frac{1}{1+2\alpha}\left[\begin{array}[]{ccc}0&\alpha&0\\\ 0&1&0\\\
0&\alpha&0\end{array}\right],$

which blurs in the vertical direction only, the opposite direction from the
banding pattern. Figure 3 shows that some regions are banding free for
$\alpha=1/8$, but banding suppression improves as blurring is further
increased. The addition of noise is crucial, the blurred images in the left
column for $\alpha=1/8$ and $1/4$ still show signs of banding that are much
less apparent after noise is applied (middle and right column). The noise also
increases perceptual sharpness, so the usage of blurring is no longer
immediately apparent. The noise used here is Gaussian noise with variance
equal to the median of the image in an $11\times 11$ pixel region, multiplied
by a constant.

The dithering method has the clear disadvantage of both decreasing the
effective resolution of the image and decreasing the signal-to-noise ratio.
This is clear when examining fine detail in Figure 3. The higher of the two
noise levels shown (0.05) clearly masks banding much more significantly than a
noise level of 0.03, however it was determined to be too much noise to be used
clinically after consultation with a radiologist. We settled on a noise level
of 0.03 as a level which is high but still a clinically relevant baseline.

Another more subtle problem arises from the inhomogeneity of the amount of
banding across the image. Applying a uniform amount of dithered across the
image results in some regions not receiving enough dithering, and other
regions that showed no banding being unnecessarily dithered. Ideally, the
dithering would be applied adaptively across the image, in proportion to the
amount of banding; however, this requires a model or method that outputs the
“banding amount” in a region, which is a non-trivial problem. Our proposed
adversarial method implicitly learns such a model, as part of the
discriminator term in its loss.

Since dithering can improve the perceptual image quality as assessed by a
viewer while at the same time reducing signal-to-noise, resolution, and
measures such as MSE and SSIM compared to the ground truth image, its use must
be carefully considered. We determined that blurring up to $\alpha=1/8$ may be
of use clinically but that $\alpha=1/4$ losses too much fine detail.

##  7 Training

We trained our models on the knee scans from the fastMRI dataset [22] using
the approach detailed above. This is the first large-scale dataset released of
raw full-sampled k-space. For the evaluation, we trained and evaluated on
scans from 1.5 Tesla machines only as banding is less visible in scans from 3
Tesla machines making evaluation more difficult, although our approach works
equally well for 3T scans. Training consisted of 100 pre-training epochs,
where the adversarial term is not used (but the flip operator is included),
using ADAM with learning rate 0.0003, momentum 0.9 and batch-size 8 (1 per GPU
on an 8-gpu system) and no weight decay. The pre-training was followed by 60
epochs of training including the adversarial term, with learning rate 0.0001,
and adversary gradient regularization strength $\gamma=0.1$. We trained with a
factor of 4 acceleration, using 16 central k-space lines, using the train/test
splits distributed as part of the fastMRI dataset.

##  8 Evaluation by Radiologists

We set up our reader study as a three-fold comparison against the adversarial
method, the non-adversarial method, and the dithering baseline. The non-
adversarial method, which uses an identical setup to the adversarial approach
excepting the adversarial terms in the loss, will be denoted the “standard”
approach. Each of the six board-certified radiologists in our study were given
a set of 20 of 40 volumes from the validation set (Each volume is
approximately 25 slice images), so that each volume of the 40 was evaluated
three times independently.

Each radiologist was asked to rank the three methods in terms of the degree of
banding and separately in terms of the retention of fine detail (the
questionnaire given to each radiologist is in the Supplementary Material) for
the volumes provided to them. Equal adaptive noise as detailed in Section 3
was used on all three to avoid the known bias of human evaluators to assess
noisy images as having higher detail levels.

They were given access to a ground-truth corresponding to the non-accelerated
(fully sampled) images. For each volume, the methods were assigned
designations “A”, ”B”, “C”, randomized at the volume level, to ensure that the
evaluators were blind to the reconstruction method used. A two-sided paired
sign-test with Bonferroni multiple comparison correction to the p-values was
used. The tests were performed using the average rank from the 3 radiologist
who ranked each volume. Volume evaluations are considered independent for the
purposes of the test, and the pairing is at the within-volume level. The
statistical analysis was chosen at the experimental design stage to avoid
statistical fishing. The raw experimental results and R script used for our
analysis are given in the Supplementary Material.

p-values | Standard | Dithering | Average rank (higher is better)
---|---|---|---
Adversary | $1.09\times 10^{-11}$ | $2.18\times 10^{-11}$ | 2.83
Dithering | 0.028 | - | 1.74
Standard | - | - | 1.43

Table 1: Banding Removal – results of the two-sided pairwise comparison with
Bonferroni correction for the following question proposed to board-certified
radiologists: “Rank the 3 methods in terms of the amount of visible banding”.
Our proposed method is ranked as better than the two baselines with very high
statistical significance.

p-values | Standard | Dithering | Average rank (higher is better)
---|---|---|---
Adversary | 2.61 | $8.82\times 10^{-4}$ | 2.18
Dithering | $3.25\times 10^{-6}$ | - | 1.5
Standard | - | - | 2.32

Table 2: Detail retention – results of the two-sided pairwise comparison with
Bonferroni correction for the following question proposed to board-certified
radiologists: “Rank the 3 methods in terms of the retention of fine anatomical
detail in comparison to the ground truth”. Our proposed method is not
statistically significantly different in terms of detail retention from the
standard baseline, and highly statistically significantly better than the
dithering approach. P-values range from 0 to 3 due to the Bonferroni
correction.

Adversary | Standard | Dithering
---|---|---
0.725 (0.62-0.82) | 0.967 (0.92-1) | 0.975 (0.90-0.99)

Table 3: Presence of banding – two-sided 95% binomial confidence intervals
with Bonferroni correction for the question proposed to board-certified
radiologists “Is any banding present (yes, no)”. Our approach completely
removes all traces of banding 27.5% of the time. The standard and dithering
approaches are rarely reported as having no banding by radiologists (<4% of
the time).

###  8.1 Evaluation Results

The radiologists ranked our adversarial approach as better than the standard
and dithering approaches with an average rank of 2.83 out of a possible 3.
This result is statistically significantly better than either alternative with
p-values $1.09\times 10^{-11}$ and $2.18\times 10^{-11}$ respectively, and the
adversarial approach was ranked as the best or tied for best in 85.8% of 120
total evaluations (95% CI: 0.78-0.91). The dithering approach is also
statistically significantly better than the standard approach.

We also asked radiologists if banding was present (in any form) in the
reconstructions in each case. This evaluation is highly subjective, as
“banding” is hard to define in a precise enough way to ensure consistency
between evaluators. Considering each radiologist’s evaluation independently,
on average banding is still reported to be present in 72.5% (95% CI:
0.62-0.82) of cases even with the adversarial learning penalty. The
radiologists were not consistent in their rankings; the overall percentages
reported by the six radiologists were 20%, 75%, 75%, 80%, 85%, and 100% for
the adversarial reconstructions. In contrast, for the baseline and dithered
reconstructions, only one radiologist reported less than 100% presence of
banding for each method (80% and 85% presence respectively, from different
radiologists).

Ground Truth Standard Orientation Adversary

Figure 4: A comparison of ground-truth images against a standard accelerated
deep-learning reconstruction and the proposed orientation adversary. Two fat-
suppressed and two non-fat suppressed images are shown, chosen at random from
the validation set. Images are best viewed on a high-brightness monitor.

##  9 Related Work

Machine-learning based reconstruction of MRI images is an ongoing research
direction and not currently in clinical use, so very little discussion of the
downsides to current techniques exists in the literature. We are not are of
any work on banding removal of accelerated MRI images; however, other forms of
banding introduced by balanced steady state free precession (bSSFP) sequences
have been investigated, using both physics-based [2] and supervised machine
learning [10] approaches. These techniques are not directly applicable to the
removal of acceleration banding artifacts.

Adversarial learning has been applied to medical imaging in a number of recent
works, for instance for segmentation [16], Magnetic Resonance Angiography
(MRA) image generation [17], super-resolution [13], and generation of
anonymized datasets [9]. Adversarial training has also been applied directly
for MRI reconstruction [14]. We believe caution should be used when applying
adversarial penalties to directly encourage reconstructions to resemble non-
accelerated ground-truth images, they are prone to introducing phantom
anatomical detail. Our use of an adversarial penalty to detect orientation
does not suffer from this problem as it only ever takes reconstructed images
as input rather than ground-truth images.

## Conclusion

In this work, we have presented an effective technique for producing machine
learning models for accelerated MRI that minimize visible banding artifacts.
Our technique is not specifically tied to any reconstruction approach, it may
be be applied on top of any MRI reconstruction technique that uses supervised
machine learning.

## Acknowledgements

This work was made possible through close collaboration with the fastMRI team
at Facebook AI research, including Nafissa Yakubova, Anuroop Sriram, Jure
Zbontar, Larry Zitnick, Mark Tygert and Suvrat Bhooshan, and our fastMRI
project members at NYU Langone Health [22], with special thanks to Florian
Knoll, Matthew Muckley and Daniel Sodickson.

## Broader Impact Statement

The use of machine learning in medical imaging raises a number of ethical and
societal considerations. The benefits to the use of our method are clear;
images produced using our method have fewer image artifacts than standard
reconstructions. Our method, used in combination with current machine-learning
based reconstruction methods, has the potential to significantly reduce time
spent in MRI scanners, and the associated cost of scanning. This benefits both
patients and medical professionals.

In theory, the removal of the banding artifacts may lead to the occlusion of
fine anatomical detail that would otherwise be present and useful for
diagnosis. Our evaluations by radiologists are strong evidence that the method
does not remove detail, but we believe a larger study directly aimed at
determining differences in diagnosis would be required to establish this. It’s
not a-priori certain that any method of banding removal will reduce detail,
since the ground truth does not have banding, the use of a banding removal
method is just a use of prior knowledge in the Bayesian sense.

Machine learning approaches to MRI reconstruction are potentially prone to
biases in the training data. If anatomy outside of the norms of the training
data is encountered, the reconstruction may not be accurate. When using our
proposed approach, the same considerations are necessary as for any
application of machine learning in medical imaging. At a minimum the model
must be tested for robustness to outlier examples.

## References

  * [1] Hemant Aggarwal, Merry Mani, and Mathews Jacob.  Modl: Model based deep learning architecture for inverse problems.  IEEE Transactions on Medical Imaging, PP, 12 2017.
  * [2] Marcus Björk, Erik Gudmundson, Joëlle K. Barral, and Petre Stoica.  Signal processing algorithms for removing banding artifacts in mri.  Proceedings of the 19th European Signal Processing Conference (EUSIPCO-2011), 2011.
  * [3] Aaron Defazio.  Offset sampling improves deep learning based accelerated mri reconstructions by exploiting symmetry.  Technical report, Facebook, 2019.
  * [4] R. Floyd and L. S. Steinberg.  An adaptive algorithm for spatial gray scale.  In International Symposium Digest of Technical Papers, 1975.
  * [5] MA Griswold, PM Jakob, RM Heidemann, M Nittka, V Jellus, J Wang, B Kiefer, and A Haase.  Generalized autocalibrating partially parallel acquisitions (grappa).  Magnetic Resonance in Medicine, 2002.
  * [6] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron Courville.  Improved training of wasserstein gans.  In Proceedings of the 31st International Conference on Neural Information Processing Systems (NeurIPS 2017), 2017.
  * [7] Kerstin Hammernik, Teresa Klatzer, Erich Kobler, Michael P. Recht, Daniel K. Sodickson, Thomas Pock, and Florian Knoll.  Learning a variational network for reconstruction of accelerated mri data.  Magnetic Resonance in Medicine, 79(6):3055–3071, 2018.
  * [8] J.F. Jarvis, C.N. Judice, and W.H. Ninke.  A survey of techniques for the display of continuous tone pictures on bilevel displays.  Computer Graphics and Image Processing, 1976.
  * [9] Koshino Kazuhiro, Rudolf Werner, Fujio Toriumi, Mehrbod Javadi, Martin Pomper, Lilja Solnes, Franco Verde, Takahiro Higuchi, and Steven Rowe.  Generative adversarial networks for the creation of realistic artificial brain magnetic resonance images.  Tomography (Ann Arbor, Mich.), 12 2018.
  * [10] Ki Hwan Kim and Sung-Hong Park.  Artificial neural network for suppression of banding artifacts in balanced steady-state free precession mri.  Magnetic Resonance Imaging, 37:139 – 146, 2017.
  * [11] Florian Knoll, Tullie Murrell, Anuroop Sriram, Nafissa Yakubova, Jure Zbontar, Michael Rabbat, Aaron Defazio, Matthew J. Muckley, Daniel K. Sodickson, C. Lawrence Zitnick, and Michael P. Recht.  Advancing machine learning for mr image reconstruction with an open competition: Overview of the 2019 fastmri challenge, 2020.
  * [12] Michael Lustig, David Donoho, and John M. Pauly.  Sparse mri: The application of compressed sensing for rapid mr imaging.  Magnetic Resonance in Medicine, 2007.
  * [13] Qing Lyu, Chenyu You, Hongming Shan, Yi Zhang, and Ge Wang.  Super-resolution MRI and CT through GAN-CIRCLE.  In Developments in X-Ray Tomography XII. International Society for Optics and Photonics, 2019.
  * [14] M. Mardani, E. Gong, J. Y. Cheng, S. S. Vasanawala, G. Zaharchuk, L. Xing, and J. M. Pauly.  Deep generative adversarial neural networks for compressive sensing mri.  IEEE Transactions on Medical Imaging, 2019.
  * [15] Lars Mescheder, Andreas Geiger, and Sebastian Nowozin.  Which training methods for GANs do actually converge?  In Proceedings of the 35th International Conference on Machine Learning, 2018.
  * [16] Pim Moeskops, Mitko Veta, Maxime W. Lafarge, Koen A. J. Eppenhof, and Josien P. W. Pluim.  Adversarial training and dilated convolutions for brain mri segmentation.  In Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support, Cham, 2017. Springer International Publishing.
  * [17] Sahin Olut, Yusuf H. Sahin, Ugur Demir, and Gozde Unal.  Generative adversarial training for mra image synthesis using multi-contrast mri.  In PRedictive Intelligence in MEdicine, Cham, 2018. Springer International Publishing.
  * [18] K P Pruessmann, M Weiger, M B Scheidegger, and P Boesiger.  Sense: sensitivity encoding for fast mri.  Magnetic Resonance in Medicine, 1999.
  * [19] J. Schlemper, J. Caballero, J. V. Hajnal, A. N. Price, and D. Rueckert.  A deep cascade of convolutional neural networks for dynamic mr image reconstruction.  IEEE Transactions on Medical Imaging, 2018.
  * [20] Anuroop Sriram, Jure Zbontar, Tullie Murrell, Aaron Defazio, C. Lawrence Zitnick, Nafissa Yakubova, Florian Knoll, and Patricia Johnson.  End-to-end variational networks for accelerated mri reconstruction.  Technical report, Facebook, 2019.
  * [21] Yuxin Wu and Kaiming He.  Group normalization.  International Journal of Computer Vision, 2019.
  * [22] Jure Zbontar, Florian Knoll, Anuroop Sriram, Matthew J. Muckley, Mary Bruno, Aaron Defazio, Marc Parente, Krzysztof J. Geras, Joe Katsnelson, Hersh Chandarana, Zizhao Zhang, Michal Drozdzal, Adriana Romero, Michael Rabbat, Pascal Vincent, James Pinkerton, Duo Wang, Nafissa Yakubova, Erich Owens, C. Lawrence Zitnick, Michael P. Recht, Daniel K. Sodickson, and Yvonne W. Lui.  fastMRI: An open dataset and benchmarks for accelerated MRI.  2018\.


# A Cramér–Wold theorem for elliptical distributions

Ricardo Fraiman  Centro de Matemática, Facultad de Ciencias, Universidad de la
República, Uruguay<EMAIL_ADDRESS>,  Leonardo Moreno  Instituto de
Estadística, Departamento de Métodos Cuantitativos, FCEA, Universidad de la
República, Uruguay<EMAIL_ADDRESS>and  Thomas Ransford  Département de
mathématiques et de statistique, Université Laval, Québec City (Québec),
Canada G1V 0A6<EMAIL_ADDRESS>
(Date: 8 March 2023)

###### Abstract.

According to a well-known theorem of Cramér and Wold, if $P$ and $Q$ are two
Borel probability measures on $\mathbb{R}^{d}$ whose projections $P_{L},Q_{L}$
onto each line $L$ in $\mathbb{R}^{d}$ satisfy $P_{L}=Q_{L}$, then $P=Q$. Our
main result is that, if $P$ and $Q$ are both elliptical distributions, then,
to show that $P=Q$, it suffices merely to check that $P_{L}=Q_{L}$ for a
certain set of $(d^{2}+d)/2$ lines $L$. Moreover $(d^{2}+d)/2$ is optimal. The
class of elliptical distributions contains the Gaussian distributions as well
as many other multivariate distributions of interest. Our theorem contrasts
with other variants of the Cramér–Wold theorem, in that no assumption is made
about the finiteness of moments of $P$ and $Q$. We use our results to derive a
statistical test for equality of elliptical distributions, and carry out a
small simulation study of the test, comparing it with other tests from the
literature. We also give an application to learning (binary classification),
again illustrated with a small simulation.

###### Key words and phrases:

Cramér–Wold, projection, elliptical distribution, Kolmogorov–Smirnov test

###### 2020 Mathematics Subject Classification:

Primary 60B11, Secondary 60E10, 62H15

##  1\. Introduction and statement of main results

Given a Borel probability measure $P$ on $\mathbb{R}^{d}$ and a vector
subspace $H$ of $\mathbb{R}^{d}$, we write $P_{H}$ for the projection of $P$
onto $H$, namely the Borel probability measure on $H$ given by

$P_{H}(B):=P(\pi_{H}^{-1}(B)),$

where $\pi_{H}:\mathbb{R}^{d}\to H$ is the orthogonal projection of
$\mathbb{R}^{d}$ onto $H$. We shall be particularly interested in the case
when $H$ is a line $L$. One can view $P_{L}$ as the marginal distribution of
$P$ along $L$.

According to a well-known theorem of Cramér and Wold [4], if $P,Q$ are two
Borel probability measures on $\mathbb{R}^{d}$, and if $P_{L}=Q_{L}$ for all
lines $L$, then $P=Q$. In other words, a probability measure on
$\mathbb{R}^{d}$ is determined by its complete set of one-dimensional marginal
distributions.

There are several extensions of this theorem, in which one assumes more about
the nature of the measures $P,Q$ and less about the set of lines $L$ for which
$P_{L}=Q_{L}$. For example, if $P$ and $Q$ have moment generating functions
that are finite in a neighbourhood of the origin, and if $P_{L}=Q_{L}$ for all
lines $L$ in a set of positive measure on the unit sphere, then $P=Q$.
Articles on this subject include those of Rényi [24], Gilbert [12],
Bélisle–Massé–Ransford [1] and Cuesta-Albertos–Fraiman–Ransford [6].

If one assumes yet more about $P$ and $Q$, then it is even possible to
differentiate between them using only a finite set of projections. Heppes [18]
showed that, if $P$ and $Q$ are supported on a finite set of cardinality $k$,
and if $H_{1},\dots,H_{k+1}$ are vector subspaces such that $H_{i}^{\perp}\cap
H_{j}^{\perp}=\\{0\\}$ whenever $i\neq j$, then $P_{H_{j}}=Q_{H_{j}}$ for all
$j$ implies that $P=Q$.

Another result of this kind, due to Gröchenig and Jaming [13], is that a
probability measure supported on a quadratic hypersurface in $\mathbb{R}^{d}$
is determined by its projections onto two generic hyperplanes. The proof is
based on the notion of a Heisenberg uniqueness pair, introduced in [17].

Our goal in this note is to establish an analogue of these results for a
certain family of continuous distributions, namely the so-called elliptical
distributions. Here is the definition.

######  Definition 1.

A Borel measure $P$ on $\mathbb{R}^{d}$ is an elliptical distribution on
$\mathbb{R}^{d}$ if its characteristic function has the form

(1)
$\phi_{P}(\xi)=e^{i\mu\cdot\xi}\psi(\xi^{\top}\Sigma\xi),\quad\xi\in\mathbb{R}^{d},$

where $\Sigma$ is a real positive semi-definite $d\times d$ matrix, $\mu$ is a
vector in $\mathbb{R}^{d}$, and $\psi:[0,\infty)\to\mathbb{C}$ is a continuous
function. The measure $P$ is said to be centred if $\mu=0$.

If $P$ is an elliptical distribution with finite second moments, then $\mu$
represents the mean and $\Sigma$ the covariance matrix. However, there are
elliptical distributions even whose first moments are infinite.

The most important elliptical distributions are surely the Gaussian
distributions. They correspond to taking $\psi(t)=e^{-t/2}$ in (1). However,
there are numerous other examples of interest, including multivariate Student
distributions, Cauchy distributions, Bessel distributions, logistic
distributions, stable laws, Kotz-type distributions, and multi-uniform
distributions. For background on the general theory of elliptical
distributions, we refer to [2] and [9]. The latter reference also contains a
more complete list of examples (see [9, Table 3.1]).

To state our results, we need a further definition.

######  Definition 2.

A set $S$ of vectors in $\mathbb{R}^{d}$ is a symmetric-matrix uniqueness set
(or sm-uniqueness set for short) if the only real symmetric $d\times d$ matrix
$A$ satisfying $x^{\top}Ax=0$ for all $x\in S$ is the zero matrix.

For example, if $v_{1},\dots,v_{d}$ are linearly independent vectors in
$\mathbb{R}^{d}$ and if we define $S:=\\{v_{i}+v_{j}:1\leq i\leq j\leq d\\}$,
then $S$ is an sm-uniqueness set. Also, every sm-uniqueness set must contain
at least $(d^{2}+d)/2$ elements, so the example $S$ above is minimal. We shall
justify these statements in §2 below, where sm-uniqueness sets will be
discussed in more detail.

We can now state our first main theorem. Given
$x\in\mathbb{R}^{d}\setminus\\{0\\}$, we denote by $\langle x\rangle$ the one-
dimensional subspace spanned by $x$.

######  Theorem 1.

Let $S\subset\mathbb{R}^{d}$ be a symmetric-matrix uniqueness set. If $P,Q$
are elliptical measures on $\mathbb{R}^{d}$ such that $P_{\langle
x\rangle}=Q_{\langle x\rangle}$ for all $x\in S$, then $P=Q$.

The following corollary is worthy of note. It follows from the fact that, in
$\mathbb{R}^{2}$, any set of three vectors, none of which is a multiple of the
others, forms an sm-uniqueness set.

######  Corollary 1.

An elliptical distribution in $\mathbb{R}^{2}$ is determined by its marginals
along any three distinct lines.

######  Remarks.

(i) Perhaps the main interest of the theorem lies in the fact that, to
differentiate between two elliptical distributions, only finitely many
projections are needed. In that sense, the theorem is an analogue of the
theorems of Heppes and of Gröchenig–Jaming, mentioned earlier. However, there
are also important differences between those results and ours. The results of
Heppes and of Gröchenig–Jaming treat measures supported on certain types of
subsets of $\mathbb{R}^{d}$, whereas Theorem 1 treats continuous ones. Also,
in their theorems, the projections are onto general subspaces $H$, whereas in
our result the projections are onto lines $L$, i.e., only one-dimensional
projections are needed.

(ii) Among the numerous variants of the Cramér–Wold theorem for continuous
distributions, Theorem 1 is unusual (maybe even unique) in that no assumption
is made about the finiteness of moments of $P$ and $Q$. As mentioned earlier,
there are elliptical distributions whose first moments are infinite, for
example the multivariate Cauchy distributions.

(iii) In Theorem 1, it is important to assume that both $P$ and $Q$ are
elliptical distributions. If we suppose merely that one of them is elliptical,
then the result no longer holds. (In the terminology of [1], Theorem 1 is not
a strong-determination result.) Indeed, given any finite set $\mathcal{H}$ of
$(d-1)$-dimensional subspaces of $\mathbb{R}^{d}$, there exist probability
measures $P,Q$ on $\mathbb{R}^{d}$ with $P$ Gaussian and $P_{H}=Q_{H}$ for all
$H\in\mathcal{H}$, but $P\neq Q$. Such examples were constructed by Hamedani
and Tata [16] in the case $d=2$, and by Manjunath and Parthasarathy [21] for
general $d$.

We now turn to the question of sharpness. The following result shows that
Theorem 1 is optimal in a certain sense. In particular, it explains why
symmetric-matrix uniqueness sets enter the picture.

We say that a Borel probability measure $P$ on $\mathbb{R}^{d}$ is non-
degenerate if it is not supported in any hyperplane in $\mathbb{R}^{d}$.

######  Theorem 2.

Let $P$ be a non-degenerate elliptical distribution on $\mathbb{R}^{d}$. Let
$S\subset\mathbb{R}^{d}$, and suppose that $S$ is not a symmetric-matrix
uniqueness set. Then there exists a non-degenerate elliptical distribution $Q$
on $\mathbb{R}^{d}$ such that $P_{\langle x\rangle}=Q_{\langle x\rangle}$ for
all $x\in S$, but $P\neq Q$.

As remarked earlier, a symmetric-matrix uniqueness set in $\mathbb{R}^{d}$
must contain at least $(d^{2}+d)/2$ elements. We thus obtain the following
corollary.

######  Corollary 2.

Let $P$ be a non-degenerate elliptical distribution on $\mathbb{R}^{d}$. Let
$\mathcal{L}$ be a set of lines in $\mathbb{R}^{d}$ containing strictly fewer
than $(d^{2}+d)/2$ lines. Then there exists a non-degenerate elliptical
distribution $Q$ on $\mathbb{R}^{d}$ such that $P_{L}=Q_{L}$ for all
$L\in\mathcal{L}$, but $P\neq Q$.

The rest of the paper is organized as follows. In §2 we discuss in more detail
the notion of symmetric-matrix uniqueness sets. The proofs of our two main
results, Theorems 1 and 2 are presented in §3 and §4 respectively. In §5 we
use Theorem 1 to derive a statistical test for equality of elliptical
distributions, and we carry out a small simulation study of the test. In §6,
we give a further application, to binary classification. Finally, in §7 we
make some concluding remarks and pose a question.

##  2\. Symmetric-matrix uniqueness sets

Recall that a set $S$ of vectors in $\mathbb{R}^{d}$ is a symmetric-matrix
uniqueness set or sm-uniqueness set if the only real symmetric $d\times d$
matrix $A$ satisfying $x^{\top}Ax=0$ for all $x\in S$ is the zero matrix. We
now examine these sets in more detail, beginning with the following simple
result.

######  Proposition 1.

Let $v_{1},\dots,v_{d}$ be linearly independent vectors in $\mathbb{R}^{d}$,
and let

$S:=\\{v_{j}+v_{k}:1\leq j\leq k\leq d\\}.$

Then $S$ is an sm-uniqueness set.

###### Proof:.

If $A$ is any symmetric matrix, then, for all $j,k$, we have

$2v_{j}^{\top}Av_{k}=(v_{j}+v_{k})^{\top}A(v_{j}+v_{k})-\frac{1}{4}(v_{j}+v_{j})^{\top}A(v_{j}+v_{j})-\frac{1}{4}(v_{k}+v_{k})^{\top}A(v_{k}+v_{k}).$

Hence, if $x^{\top}Ax=0$ for all $x\in S$, then $v_{j}^{\top}Av_{k}=0$ for all
$j,k$, and so $A=0$. ∎

######  Corollary 3.

If $v_{1},v_{2},v_{3}\in\mathbb{R}^{2}$ and no $v_{j}$ is a multiple of any
other, then $\\{v_{1},v_{2},v_{3}\\}$ is an sm-uniqueness set.

###### Proof:.

Clearly $v_{1},v_{2}$ are linearly independent. Also, replacing them by
suitable non-zero multiples of themselves, we can suppose that
$v_{3}=v_{1}+v_{2}$. The result therefore follows from Proposition 1. ∎

In higher dimensions, criteria for sm-uniqueness sets are more complicated,
though we shall derive some in Proposition 4 and Corollary 5 below. However,
for statistical applications, it suffices merely to have some concrete
examples of sm-uniqueness sets. The next result furnishes a particularly
simple example.

######  Corollary 4.

Let $S$ be the set consisting of those vectors in $\mathbb{R}^{d}$ with either
one or two coordinates equal to $1$ and all the other coordinates equal to
$0$. Then $S$ is an sm-uniqueness set.

Fig. 1 illustrates the set $S$ in dimensions $2$ and $3$.

$(1,0)$$(0,1)$$(1,1)$

$(1,0,0)$$(0,0,1)$$(1,0,1)$$(0,1,0)$$(1,1,0)$$(0,1,1)$

Fig. 1. Visualization of the vectors of the set $S$ from Corollary 4 in
dimensions two and three.

###### Proof of Corollary 4:.

In Proposition 1, take $v_{j}=e_{j}$, where $\\{e_{1},\dots,e_{d}\\}$ is the
standard unit vector basis of $\mathbb{R}^{d}$. By that result, the union of
the two sets $\\{e_{j}+e_{k}:1\leq j<k\leq d\\}$ and $\\{2e_{j}:1\leq j\leq
d\\}$ is an sm-uniqueness set. It clearly makes no difference if, in the
second set, we replace $2e_{j}$ by $e_{j}$. ∎

The sets $S$ in Proposition 1 and Corollary 4 contain $(d^{2}+d)/2$ elements.
The next result shows that this number is minimal.

######  Proposition 2.

An sm-uniqueness set in $\mathbb{R}^{d}$ contains at least $(d^{2}+d)/2$
elements.

###### Proof:.

Let $v_{1},\dots,v_{m}\in\mathbb{R}^{d}$ with $m<(d^{2}+d)/2$. The vector
space $V$ of all $d\times d$ symmetric matrices has dimension $(d^{2}+d)/2$,
so the linear map
$A\mapsto(v_{1}^{\top}Av_{1},\dots,v_{m}^{\top}Av_{m}):V\to\mathbb{R}^{m}$
must have a non-zero kernel. Thus there exists $A\neq 0$ such that
$v_{j}^{\top}Av_{j}=0$ for $j\in\\{1,\dots,m\\}$. Hence
$\\{v_{1},\dots,v_{m}\\}$ is not an sm-uniqueness set. ∎

######  Proposition 3.

If $S$ is an sm-uniqueness set in $\mathbb{R}^{d}$, then $S$ spans
$\mathbb{R}^{d}$.

###### Proof:.

Suppose, on the contrary, that $S$ does not span $\mathbb{R}^{d}$. Then we can
find a unit vector $v\in S^{\perp}$. Let $R$ be a $d\times d$ rotation matrix
such that $Rv=e_{1}$, and let $A$ be the $d\times d$ symmetric matrix given by
$A:=R^{\top}DR$, where $D:=\operatorname{diag}(1,0,\dots,0)$. If $x\in S$,
then $x\perp v$, so $Rx\perp Rv=e_{1}$, so
$Rx\in\operatorname{span}\\{e_{2},\dots,e_{d}\\}$ and hence
$(Rx)^{\top}D(Rx)=0$. In other words, $x^{\top}Ax=0$ for all $x\in S$. On the
other hand, $v^{\top}Av=(Rv)^{\top}DRv=e_{1}^{\top}De_{1}=1$, so $A\neq 0$. We
conclude that $S$ is not an sm-uniqueness set. ∎

The next result is not really needed in what follows, but it provides a
necessary and sufficient condition for a set of vectors in $\mathbb{R}^{d}$ to
form an sm-uniqueness set.

######  Proposition 4.

Given a vector $x\in\mathbb{R}^{d}$, say $x=(x_{1},\dots,x_{d})$, let
$\widehat{x}$ be the upper-triangular $d\times d$ matrix with entries
$\widehat{x}_{ij}:=x_{i}x_{j}~{}(1\leq i\leq j\leq d)$. A set
$S\subset\mathbb{R}^{d}$ is an sm-uniqueness set if and only if
$\\{\widehat{x}:x\in S\\}$ spans the space of upper-triangular $d\times d$
matrices.

###### Sketch of Proof:.

Denote by $\operatorname{Sym}^{2}(\mathbb{R}^{d})$ the symmetric tensor
product of $\mathbb{R}^{d}$ with itself. For each symmetric $d\times d$ matrix
$A$, there is a unique linear functional
$\widetilde{A}:\operatorname{Sym}^{2}(\mathbb{R}^{d})\to\mathbb{R}$ such that
$\widetilde{A}(x\otimes x)=x^{\top}Ax$ for all $x\in\mathbb{R}^{d}$. Thus a
set $S\subset\mathbb{R}^{d}$ is an sm-uniqueness set iff $\\{x\otimes x:x\in
S\\}$ spans $\operatorname{Sym}^{2}(\mathbb{R}^{d})$.

Let $e_{1},\dots,e_{d}$ be the standard unit basis of $\mathbb{R}^{d}$. Then
$\operatorname{Sym}^{2}(\mathbb{R}^{d})$ has a basis $\\{E_{ij}:1\leq i\leq
j\leq d\\}$, where $E_{ij}:=(e_{i}\otimes e_{j}+e_{j}\otimes e_{i})/2$.
Expressing $x\otimes x$ in terms of this basis, we have

$x\otimes x=\sum_{1\leq i\leq j\leq d}x_{i}x_{j}E_{ij}.$

The result follows easily from this. ∎

Since the space of upper-triangular $d\times d$ matrices has dimension
$(d^{2}+d)/2$, we deduce the following corollary.

######  Corollary 5.

Let $S$ be a set of $(d^{2}+d)/2$ vectors in $\mathbb{R}^{d}$. For each $x\in
S$, define $\widehat{x}$ as in Proposition 4. Then $S$ is an sm-uniqueness set
if and only if the matrices $\\{\widehat{x}:x\in S\\}$ are linearly
independent.

######  Remark.

If we view an upper-triangular $d\times d$ matrix as a column vector of length
$D:=(d^{2}+d)/2$, then Corollary 5 becomes a criterion expressed in terms of
the linear independence of $D$ vectors in $\mathbb{R}^{D}$, which can in turn
be reformulated as the non-vanishing of a $D\times D$ determinant. This
furnishes a systematic method of determining whether or not a given set of
$(d^{2}+d)/2$ vectors in $\mathbb{R}^{d}$ is an sm-uniqueness set.

##  3\. Proof of Theorem 1

We break the proof into a series of lemmas. The first of these is fairly
standard.

######  Lemma 1.

Let $P,Q$ be Borel probability measures on $\mathbb{R}^{d}$ and let $L$ be a
line in $\mathbb{R}^{d}$. Then $P_{L}=Q_{L}$ if and only if the characteristic
functions of $P,Q$ satisfy $\phi_{P}(\xi)=\phi_{Q}(\xi)$ for all $\xi\in L$.

###### Proof:.

Let $\xi\in L$. Then $x\cdot\xi=\pi_{L}(x)\cdot\xi$ for all
$x\in\mathbb{R}^{d}$. Consequently,

$\phi_{P}(\xi)=\int_{\mathbb{R}^{d}}e^{ix\cdot\xi}\,dP(x)=\int_{\mathbb{R}^{d}}e^{i\pi_{L}(x)\cdot\xi}\,dP(x)=\int_{L}e^{iy\cdot\xi}\,dP_{L}(y),$

and similarly for $Q$. Hence, if $P_{L}=Q_{L}$, then
$\phi_{P}(\xi)=\phi_{Q}(\xi)$.

Conversely, if $\phi_{P}(\xi)=\phi_{Q}(\xi)$ for all $\xi\in L$, then the
calculation above shows that $P_{L}$ and $Q_{L}$ have the same characteristic
function, and consequently $P_{L}=Q_{L}$. ∎

######  Lemma 2.

Let $P,Q$ be elliptical distributions, with characteristic functions

$\phi_{P}(\xi)=e^{i\mu_{1}\cdot\xi}\psi_{1}(\xi^{\top}\Sigma_{1}\xi)\quad\text{and}\quad\phi_{Q}(\xi)=e^{i\mu_{2}\cdot\xi}\psi_{2}(\xi^{\top}\Sigma_{2}\xi),\quad\xi\in\mathbb{R}^{d}.$

If $P_{\langle x\rangle}=Q_{\langle x\rangle}$ for all $x$ in some spanning
set of $\mathbb{R}^{d}$, then $\mu_{1}=\mu_{2}$.

###### Proof:.

By Lemma 1, if $P_{\langle x\rangle}=Q_{\langle x\rangle}$, then
$\phi_{P}=\phi_{Q}$ on $\langle x\rangle$, in other words
$\phi_{P}(tx)=\phi_{Q}(tx)$ for all $t\in\mathbb{R}$. Recalling the form of
$\phi_{P},\phi_{Q}$, we obtain

(2) $e^{it\mu_{1}\cdot x}\psi_{1}(t^{2}x^{\top}\Sigma_{1}x)=e^{it\mu_{2}\cdot
x}\psi_{2}(t^{2}x^{\top}\Sigma_{2}x),\quad t\in\mathbb{R}.$

Since $\phi_{P},\phi_{Q}$ are characteristic functions, they are both equal to
$1$ at the origin. It follows that $\psi_{1}(0)=\psi_{2}(0)=1$. As
$\psi_{1},\psi_{2}$ are continuous, there exists $\delta>0$ such that both
$\psi_{1}(t^{2}x^{\top}\Sigma_{1}x)\neq 0$ and
$\psi_{2}(t^{2}x^{\top}\Sigma_{2}x)\neq 0$ for all $t\in[-\delta,\delta]$. For
each such $t$, we may divide equation (2) by the same equation in which $t$ is
replaced by $-t$. This gives

$e^{2it\mu_{1}\cdot x}=e^{2it\mu_{2}\cdot x},\quad t\in[-\delta,\delta].$

Differentiating with respect to $t$ and setting $t=0$, we obtain $\mu_{1}\cdot
x=\mu_{2}\cdot x$. Finally, if this holds for all $x$ in a spanning set of
$\mathbb{R}^{d}$, then $\mu_{1}=\mu_{2}$. ∎

######  Lemma 3.

Let $P,Q$ be centred elliptical distributions on $\mathbb{R}^{d}$ such that
$P_{\langle x\rangle}=Q_{\langle x\rangle}$ for all $x$ in some spanning set
of $\mathbb{R}^{d}$. Then there exists a continuous function
$\psi:[0,\infty)\to\mathbb{C}$ and positive semi-definite $d\times d$ matrices
$\Sigma_{1},\Sigma_{2}$ such that

(3)
$\phi_{P}(\xi)=\psi(\xi^{\top}\Sigma_{1}\xi),\quad\phi_{Q}(\xi)=\psi(\xi^{\top}\Sigma_{2}\xi),\quad\xi\in\mathbb{R}^{d}.$

###### Proof:.

By assumption, there exist continuous functions
$\psi_{1},\psi_{2}:[0,\infty)\to\mathbb{C}$ and positive semi-definite
$d\times d$ matrices $\Sigma_{1}^{\prime},\Sigma_{2}^{\prime}$ such that

$\phi_{P}(\xi)=\psi_{1}(\xi^{\top}\Sigma_{1}^{\prime}\xi),\quad\phi_{Q}(\xi)=\psi_{2}(\xi^{\top}\Sigma_{2}^{\prime}\xi),\quad\xi\in\mathbb{R}^{d}.$

The point of the lemma is to show that, adjusting
$\Sigma_{1}^{\prime},\Sigma_{2}^{\prime}$ if necessary, we may take
$\psi_{1}=\psi_{2}$.

If $\Sigma_{1}^{\prime}=\Sigma_{2}^{\prime}=0$, then we may as well take
$\psi_{1}=\psi_{2}\equiv 1$.

Otherwise, at least one of $\Sigma_{1}^{\prime},\Sigma_{2}^{\prime}$ is non-
zero, say $\Sigma_{1}^{\prime}\neq 0$. Since $\Sigma_{1}^{\prime}$ is positive
semi-definite, any spanning set of $\mathbb{R}^{d}$ contains a vector $x$ such
that $x^{\top}\Sigma_{1}^{\prime}x>0$. Hence, there exists
$x_{0}\in\mathbb{R}^{d}$ such that both $P_{\langle x_{0}\rangle}=Q_{\langle
x_{0}\rangle}$ and $x_{0}^{\top}\Sigma_{1}^{\prime}x_{0}>0$. By Lemma 1,

$\psi_{1}(t^{2}x_{0}^{\top}\Sigma_{1}^{\prime}x_{0})=\psi_{2}(t^{2}x_{0}^{\top}\Sigma_{2}^{\prime}x_{0}),\quad
t\in\mathbb{R}.$

It follows that

$\psi_{1}(s)=\psi_{2}(as),\quad s\geq 0,$

where
$a:=(x_{0}^{\top}\Sigma_{2}^{\prime}x_{0})/(x_{0}^{\top}\Sigma_{1}^{\prime}x_{0})\geq
0$. Feeding this information into the formula for $\phi_{P}$, we obtain

$\phi_{P}(\xi)=\psi_{2}(\xi^{\top}(a\Sigma_{1}^{\prime})\xi),\quad\xi\in\mathbb{R}^{d}.$

Thus (3) holds with $\psi:=\psi_{2}$ and $\Sigma_{1}:=a\Sigma_{1}^{\prime}$
and $\Sigma_{2}:=\Sigma_{2}^{\prime}$. ∎

######  Lemma 4.

Let $P,Q$ be probability measures on $\mathbb{R}^{d}$ with characteristic
functions given by (3). Let $S$ be an sm-uniqueness set, and suppose that
$P_{\langle x\rangle}=Q_{\langle x\rangle}$ for all $x\in S$. Then $P=Q$.

###### Proof:.

We show that either $\Sigma_{1}=\Sigma_{2}$ or $\psi$ is constant. Either way,
this gives that $\phi_{P}=\phi_{Q}$ and hence that $P=Q$.

Suppose then that $\Sigma_{1}\neq\Sigma_{2}$. As $S$ is an sm-uniqueness set,
there exists $x_{0}\in S$ such that
$x_{0}^{\top}(\Sigma_{1}-\Sigma_{2})x_{0}\neq 0$. Exchanging the roles of
$P,Q$ if necessary, we can suppose that

$x_{0}^{\top}\Sigma_{1}x_{0}>x_{0}^{\top}\Sigma_{2}x_{0}\geq 0.$

Since $P_{\langle x_{0}\rangle}=Q_{\langle x_{0}\rangle}$, Lemma 1 gives

$\psi(t^{2}x_{0}^{\top}\Sigma_{1}x_{0})=\psi(t^{2}x_{0}^{\top}\Sigma_{2}x_{0}),\quad
t\in\mathbb{R}.$

It follows that

$\psi(s)=\psi(as),\quad s\geq 0,$

where
$a:=(x_{0}^{\top}\Sigma_{2}x_{0})/(x_{0}^{\top}\Sigma_{1}x_{0})\in[0,1)$.
Iterating this equation, and using the continuity of $\psi$, we obtain

$\psi(s)=\psi(as)=\psi(a^{2}s)=\cdots=\psi(a^{n}s)\underset{n\to\infty}{\longrightarrow}\psi(0),\quad
s\geq 0.$

Therefore $\psi$ is constant, as claimed. ∎

###### Completion of proof of Theorem 1:.

Let $P,Q$ be elliptical distributions on $\mathbb{R}^{d}$, with characteristic
functions

$\phi_{P}(\xi)=e^{i\mu_{1}\cdot\xi}\psi_{1}(\xi^{\top}\Sigma_{1}\xi)\quad\text{and}\quad\phi_{Q}(\xi)=e^{i\mu_{2}\cdot\xi}\psi_{2}(\xi^{\top}\Sigma_{2}\xi).$

Suppose that $P_{\langle x\rangle}=Q_{\langle x\rangle}$ for all $x\in S$,
where $S\subset\mathbb{R}^{d}$ is an sm-uniqueness set. Our goal is to prove
that $P=Q$.

By Proposition 3, $S$ is a spanning set for $\mathbb{R}^{d}$. By Lemma 2, it
follows that $\mu_{1}=\mu_{2}$. Translating both $P,Q$ by $-\mu_{1}$, we may
suppose that both measures are centred. Lemma 3 then shows that $P,Q$ have
characteristic functions given by (3), and finally Lemma 4 gives $P=Q$. ∎

##  4\. Proof of Theorem 2

Once again, we break the proof up into lemmas. The first lemma characterizes
those elliptical distributions on $\mathbb{R}^{d}$ that are non-degenerate
(i.e., not supported on any hyperplane in $\mathbb{R}^{d}$).

######  Lemma 5.

Let $P$ be an elliptical distribution on $\mathbb{R}^{d}$ with characteristic
function

$\phi_{P}(\xi)=e^{i\mu\cdot\xi}\psi(\xi^{\top}\Sigma\xi),\quad\xi\in\mathbb{R}^{d}.$

Then $P$ is non-degenerate if and only if $\psi$ is non-constant and $\Sigma$
is strictly positive definite.

###### Proof:.

Suppose that $P$ is supported on the hyperplane
$H:=\\{x\in\mathbb{R}^{d}:y\cdot x=c\\}$, where
$y\in\mathbb{R}^{d}\setminus\\{0\\}$ and $c\in\mathbb{R}$. Then $P_{\langle
y\rangle}=\delta_{x_{0}}$, where $x_{0}$ is the unique point in $\langle
y\rangle\cap H$. By Lemma 1,

$\phi_{P}(ty)=e^{ity\cdot x_{0}},\quad t\in\mathbb{R}.$

Recalling the form of $\phi_{P}$, we deduce that

$e^{it\mu\cdot y}\psi(t^{2}y^{\top}\Sigma y)=e^{ity\cdot x_{0}},\quad
t\in\mathbb{R}.$

The argument used to prove Lemma 2 shows that $\mu\cdot y=x_{0}\cdot y$, and
hence that

$\psi(t^{2}y^{\top}\Sigma y)=1,\quad t\in\mathbb{R}.$

This in turn implies that either $\psi\equiv 1$ or $y^{\top}\Sigma y=0$. Thus
either $\psi$ is constant or $\Sigma$ is not strictly positive definite.

The converse is proved by running the same argument backwards. ∎

######  Lemma 6.

Let $P$ be a non-degenerate elliptical distribution on $\mathbb{R}^{d}$, with
characteristic function

$\phi_{P}(\xi)=e^{i\mu_{1}.\xi}\psi(\xi^{\top}\Sigma_{1}\xi),\quad\xi\in\mathbb{R}^{d}.$

Then, for each vector $\mu_{2}\in\mathbb{R}^{d}$ and each positive semi-
definite $d\times d$ matrix $\Sigma_{2}$, there exists an elliptical
distribution $Q$ on $\mathbb{R}^{d}$ with characteristic function

$\phi_{Q}(\xi)=e^{i\mu_{2}.\xi}\psi(\xi^{\top}\Sigma_{2}\xi),\quad\xi\in\mathbb{R}^{d}.$

###### Proof:.

Since $\Sigma_{1},\Sigma_{2}$ are positive semi-definite, we can write them as
$\Sigma_{1}=S_{1}^{\top}S_{1}$ and $\Sigma_{2}=S_{2}^{\top}S_{2}$, where
$S_{1},S_{2}$ are $d\times d$ matrices. Further, as $P$ is non-degenerate,
Lemma 5 implies that $\Sigma_{1}$ is strictly positive definite, and so
$S_{1}$ is invertible. Let $A:\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the affine
map defined by

$Ax:=S_{1}^{-1}S_{2}(x-\mu_{1})+\mu_{2},\quad x\in\mathbb{R}^{d},$

and set $Q:=PA^{-1}$. We shall show that $\phi_{Q}$ has the required form.

First of all, we remark that, if $\widetilde{P}$ and $\widetilde{Q}$ denote
the translates of $P,Q$ by $-\mu_{1}$ and $-\mu_{2}$ respectively, then

$\phi_{\widetilde{P}}(\xi)=e^{-i\mu_{1}\cdot\xi}\phi_{P}(\xi)\quad\text{and}\quad\phi_{\widetilde{Q}}(\xi)=e^{-i\mu_{2}\cdot\xi}\phi_{Q}(\xi),\quad\xi\in\mathbb{R}^{d}.$

Further, $\widetilde{Q}=\widetilde{P}T^{-1}$, where
$T:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is the linear map given by

$Tx:=S_{1}^{-1}S_{2}x,\quad x\in\mathbb{R}^{d}.$

By a calculation similar to that in the proof of Lemma 1, we have

$\phi_{\widetilde{P}T^{-1}}(\xi)=\phi_{\widetilde{P}}(T\xi),\quad\xi\in\mathbb{R}^{d}.$

Putting all of this together, we get

$\displaystyle e^{-i\mu_{2}\cdot\xi}\phi_{Q}(\xi)$
$\displaystyle=\phi_{\widetilde{Q}}(\xi)=\phi_{\widetilde{P}T^{-1}}(\xi)=\phi_{\widetilde{P}}(T\xi)$
$\displaystyle=e^{-i\mu_{1}\cdot\xi}\phi_{P}(T\xi)=e^{-i\mu_{1}\cdot\xi}\phi_{P}(S_{1}^{-1}S_{2}\xi).$

Recalling the form of $\phi_{P}$ and the fact that
$\Sigma_{j}=S_{j}^{\top}S_{j}$ for $j\in\\{1,2\\}$, we obtain

$\displaystyle e^{-i\mu_{2}\cdot\xi}\phi_{Q}(\xi)$
$\displaystyle=e^{-i\mu_{1}\cdot\xi}e^{i\mu_{1}\cdot\xi}\psi(\xi^{\top}S_{2}^{\top}S_{1}^{-\top}\Sigma_{1}S_{1}^{-1}S_{2}\xi)$
$\displaystyle=\psi(\xi^{\top}S_{2}^{\top}S_{2}\xi)=\psi(\xi^{\top}\Sigma_{2}\xi).$

Thus $\phi_{Q}$ does indeed have the required form. ∎

###### Completion of the proof of Theorem 2:.

By Lemma 5, since $P$ is a non-degenerate elliptical distribution, its
characteristic function is given by

$\phi_{P}(\xi)=e^{i\mu.\xi}\psi(\xi^{\top}\Sigma_{1}\xi),\quad\xi\in\mathbb{R}^{d},$

where $\psi$ is non-constant and $\Sigma_{1}$ is strictly positive definite.
As $S$ is not an sm-uniqueness set, there exists a symmetric $d\times d$
matrix $A$ such that $x^{\top}Ax=0$ for all $x\in S$ but $A\neq 0$. If
$\epsilon>0$ is small enough, then $\Sigma_{1}+\epsilon A$ is also strictly
positive definite. Fix such an $\epsilon$, and set
$\Sigma_{2}:=\Sigma_{1}+\epsilon A$. By Lemma 6, there exists an elliptical
distribution $Q$ on $\mathbb{R}^{d}$ such that

$\phi_{Q}(\xi)=e^{i\mu.\xi}\psi(\xi^{\top}\Sigma_{2}\xi),\quad\xi\in\mathbb{R}^{d}.$

By Lemma 5 again, $Q$ is non-degenerate, since $\psi$ is non-constant
$\Sigma_{2}$ is positive definite. Also, if $x\in S$, then

$x^{\top}\Sigma_{2}x=x^{\top}\Sigma_{1}x+x^{\top}Ax=x^{\top}\Sigma_{1}x,$

so $\phi_{Q}(tx)=\phi_{P}(tx)$ for all $t$. Hence $P_{\langle
x\rangle}=Q_{\langle x\rangle}$ for all $x\in S$. On the other hand, since
$\Sigma_{1}\neq\Sigma_{2}$ and $\psi$ is non-constant, the argument of the
proof of Lemma 4 shows that $P\neq Q$. ∎

##  5\. Application to testing for equality of multivariate distributions

###  5.1. A Kolmogorov–Smirnov test for elliptical distributions

Several of the variants of the Cramér–Wold theorem mentioned in the
introduction are useful in deriving statistical tests for the equality of
multivariate distributions (see, e.g., [5, 7, 11, 10]). Theorem 1 is no
exception. In this section, we propose a test for the one- and two-sample
problems for elliptical distributions, based on Theorem 1. The problem of
goodness-of-fit testing for elliptical distributions is discussed in several
recent articles (see, e.g., [3, 8, 15]), but always with more restrictive
assumptions than ours (supposing, for example, that the generator function
$\psi$ is known). We also remark that the idea of using projections of
elliptical distributions is mentioned in [22], though that article addresses a
different problem.

Since the one- and two-sample problems are very similar, we describe just the
two-sample problem. Given two samples
$\mathcal{X}_{n}:=\\{X_{1},\ldots,X_{n}\\}\subset\mathbb{R}^{d}$ and
$\mathcal{Y}_{m}:=\\{Y_{1},\ldots,Y_{m}\\}\subset\mathbb{R}^{d}$ of
multivariate elliptical distributions $P$ and $Q$, we consider the testing
problem

$\mathbf{H0}:P=Q\quad\text{vs}\quad P\neq Q,$

based on the samples $\mathcal{X}_{n},\mathcal{Y}_{m}$.

Let $\\{v_{1},\dots v_{D}\\}\subset\mathbb{R}^{d}$ be a fixed symmetric-matrix
uniqueness set, where $D=(d^{2}+d)/2$. Let $F_{n,\langle v_{j}\rangle}$ be the
empirical distribution of $P_{\langle v_{j}\rangle}$ and $G_{m,\langle
v_{j}\rangle}$ be the empirical distribution of $Q_{\langle v_{j}\rangle}$,
for $j\in\\{1,\ldots,D\\}$. We define a random-projection test (which we
abbreviate as RPT) through the statistic,

$\sqrt{\frac{nm}{n+m}}KS_{n,m,D}:=\sqrt{\frac{nm}{n+m}}\max_{j\in\\{1,\ldots,D\\}}\|F_{n,\langle
v_{j}\rangle}-G_{m,\langle v_{j}\rangle}\|_{\infty}.$

Since the statistic is not distribution-free, in order to obtain the critical
value for a level-$\alpha$ test, we approximate the distribution using
bootstrap on the original samples $\mathcal{X}_{n},\mathcal{Y}_{m}$ by
generating a large enough number $B$ of values of $KS_{n,m}$, for each
bootstrap sample choosing $n$ vectors from $\mathcal{X}_{n}$ and $m$ vectors
from $\mathcal{Y}_{m}$, with replacement. See, for instance, [14] for the two-
sample bootstrap. We then take as critical value $c^{*}_{\alpha}$, the
$(1-\alpha)$-quantile of the empirical bootstrap sample, i.e., we reject the
null hypothesis when

$\sqrt{\frac{nm}{n+m}}KS_{n,m,D}>c^{*}_{\alpha}.$

The validity of the bootstrap in this case follows from [23, Theorems 3 and
4]. Therefore the proposed test has asymptotic level $\alpha$. Also it is
consistent since, under the alternative hypothesis, $P\neq Q$, so, by Theorem
1, there exists a $j$ such that $P_{\langle v_{j}\rangle}\neq Q_{\langle
v_{j}\rangle}$, so $\max_{j\in\\{1,\ldots,D\\}}\|F_{n,\langle
v_{j}\rangle}-G_{m,\langle v_{j}\rangle}\|_{\infty}>0$, and thus
$\sqrt{\frac{nm}{n+m}}KS_{n,m,D}\to\infty$.

We summarize our conclusions in the following proposition.

######  Proposition 5.

Consider the test proposed above.

  1. (i)

The bootstrap version of the test has asymptotic level $\alpha$.

  2. (ii)

Under the alternative,

$\sqrt{\frac{nm}{n+m}}KS_{n,m,D}\to\infty,$

i.e., the test is consistent.

###  5.2. A small simulation study

In this subsection, we study the power of the proposed test for two elliptical
distributions $F_{1}$ and $F_{2}$ belonging to the multivariate Student
distribution in $\mathbb{R}^{d}$. More precisely, we suppose that

$F_{i}\sim t_{\nu_{i}}(\mu_{i},\Sigma_{i}),\quad i\in\\{1,2\\},$

where $\nu_{i}\in\mathbb{N}$ are the numbers of degrees of freedom,
$\mu_{i}\in\mathbb{R}^{d}$ are the means, and
$\Sigma_{i}\in\mathbb{R}^{d\times d}$ are the covariance matrices. We consider
three different scenarios.

In the first scenario, we fix $\mu_{i}=0$ and $\Sigma_{i}=\mathbf{I}_{d\times
d}$ (the $d\times d$ identity matrix), for $i\in\\{1,2\\}$ and we vary only
the degrees of freedom: $\nu_{1}=1$ and $\nu_{2}\in\\{1,2,3,4\\}$. The power
functions are plotted in Fig. 2 for different sample sizes.

In the second scenario, we vary only the mean value of one of the
distributions, taking $\nu_{i}=2$ and $\Sigma_{i}=\mathbf{I}_{d\times d}$, for
$i\in\\{1,2\\}$, while $\mu_{1}=0$ and $\mu_{2}\in[0,1/2]$. See Fig. 3.

Lastly, in the third scenario, we vary only the covariance matrices. We fix
$\nu_{i}=2$ and $\mu_{i}=0$ for $i\in\\{1,2\\}$, while
$\Sigma_{1}=\mathbf{I}_{d\times d}$ and $\Sigma_{2}=\mathbf{I}_{d\times
d}+\theta\mathbf{1}_{d\times d}$, where $\mathbf{1}_{d\times d}$ denotes the
matrix with value $1$ in all its entries, and $\theta\in[0,1]$. See Fig. 4.

In all cases we generate iid samples of equal sizes $n=n_{1}=n_{2}$, where
$n\in\\{100,500,1000\\}$ and $d\in\\{5,10,50\\}$ for each scenario. The null
distribution is approximated by bootstrap. For the power function, we report
the proportion of rejections in 10000 replicates. The results are quite
encouraging in all the three scenarios considered. Moreover, in the
simulations where the sample size is large, the test size remained around the
prespecified significance level ($\alpha=0.05$).

(a)  Fig. 2. Power function of test proposed in §5.2 for two elliptical
distributions. Scenario 1: changing the degrees of freedom and the sample
sizes. (a)  Fig. 3. Power function of test proposed in §5.2 for two elliptical
distributions. Scenario 2: changing the mean value and the sample sizes. (a)
Fig. 4. Power function of test proposed in §5.2 for two elliptical
distributions. Scenario 3: changing the covariance matrix as a function of
$\theta$ and the sample sizes

###  5.3. Comparison with other tests

The following example compares the performance of the proposed test RPT for
the equality of two distributions, with simulated data from a mixture of
distributions, with some other alternatives known from the literature. The
cases when the mixture has an elliptical distribution (and in particular a
multivariate normal distribution) are considered, as well as the case when the
mixture distribution is not elliptical. The performance of the RPT test is
compared with two other proposals that impose different assumptions about data
distributions. There are many test proposals for equality of multivariate
normal distributions, but most of them are for equality of the mean—with known
or unknown covariance matrix—or for equality of the covariance matrices—with
known or unknown mean vector, but just a few more general. This is the reason
why we have chosen these three particular competitors.

Let $\mu\in\mathbb{R}$ and
$\textbf{1}_{d}=(1,1,\dots,1)^{\top}\in\mathbb{R}^{d}$. Consider three
independent random vectors in $\mathbb{R}^{d}$, with distributions given by

$\displaystyle X_{1,\mu}$
$\displaystyle\sim\mathcal{N}\left(\mu\textbf{1}_{d},\mathbf{I}_{d\times
d}\right),\quad
X_{2,\mu}\sim\mathcal{C}\left(\mu\textbf{1}_{d},\mathbf{I}_{d\times
d}\right),$ $\displaystyle X_{3,\mu}$ $\displaystyle\sim
B\mathcal{N}\left((1+\mu)\textbf{1}_{d},\mathbf{I}_{d\times
d}\right)+(1-B)\mathcal{N}\left(-(1+\mu)\textbf{1}_{,}\mathbf{I}_{d\times
d}\right),$

with $B\sim\mathcal{B}(1/2)$, where $\mathcal{B}$ is the Bernoulli
distribution, and where $\mathcal{N}$ and $\mathcal{C}$ are the Normal and
Cauchy distributions in dimension $d$ respectively. Let $F_{1,\mu},F_{2,\mu}$
and $F_{3,\mu}$ be their respective probability distributions. We define the
convex combination

$F_{\mu}=\alpha_{1}F_{1,\mu}+\alpha_{2}F_{2,\mu}+\alpha_{3}F_{3,\mu},$

with $\alpha_{1},\alpha_{2},\alpha_{3}\geq 0$ and
$\alpha_{1}+\alpha_{2}+\alpha_{3}=1$. Note that, if $\alpha_{2}=\alpha_{3}=0$,
then the distribution is Gaussian; if $\alpha_{3}=0$, then the distribution is
elliptical; and if $\alpha_{3}\neq 0$, then the distribution is not
elliptical.

In all cases we generate iid samples of equal sizes $n_{1}=n_{2}=100$, $d=5$
and $F_{1}\sim F_{\mu_{1}}$ and $F_{2}\sim F_{\mu_{2}}$ with $\mu_{1}=0$ and
$\mu_{2}\in\\{0,0.2,0.4,0.6\\}$. When $\mu_{2}=0$ we are under the null
hypothesis (that the distributions of $F_{1}$ and $F_{2}$ are the same). The
power is studied in different scenarios obtained for different values of
$\alpha=(\alpha_{1},\alpha_{2},\alpha_{3}$).

The power of RPT is compared to:

LRTN:

a simultaneous test for the mean and the covariance matrix for two samples
under the assumption of normality (see [20] for details);

L2Norm:

a more general test involving non-normal distributions (see [19] for details),
and

eDistance:

a non-parametric test developed in [25].

The L2Norm and eDistance tests are implemented in R through the sim2.2018HN
and eqdist.etest functions of the SHT and energy packages respectively. Note
that, as mentioned in [19], the L2Norm test does not include in its
assumptions all elliptical distributions. In the RPT test, the number of the
$k$-nearest neighbors is chosen by cross-validation.

As shown in Fig. 5, if $\alpha_{2}\neq 0$ the tests LRTN and L2Norm perform
very poorly. This is to be expected, as they are not designed for heavy-tailed
distributions. Also, if $\alpha_{2}\neq 0$ and $\alpha_{3}=0$, then the best-
performing test in all cases is RPT. Even if the distributions are not
elliptical, i.e., $\alpha_{3}\neq 0$, RPT is a competitive test with
eDistance.

(a)  Fig. 5. Power functions of the RPT, LRTN, L2Norm, eDistance tests in
§5.3, for $\mu_{2}\in\\{0,0.2,0.4,0.6\\}$ and for different values of the
vector $\alpha$.

##  6\. Application to binary classification

###  6.1. Binary classification of data from elliptical distributions

Linear Discriminant Analysis and Quadratic Discriminant Analysis are two of
the more classical and well-known learning methods, both based on the
Mahalanobis distance. They assume that the distributions of the different
classes are multivariate Gaussian distributions. In this section we consider
the more challenging learning problem where we only assume that the
distributions are elliptical.

We propose the following algorithm for binary classification, based on our
results. Let
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\in\mathbb{R}^{d}\times\\{0,1\\}$ be an iid
training sample, where the distributions of $X_{1}|Y=0$ and $X_{1}|Y=1$ are
unknown elliptical distributions. From the training sample, a proportion
$\omega\in(0,1)$ will be used to assign weights to the different directions.
Let $n_{1}$ be the integer part of $(1-\omega)n$, and let
$\mathcal{X}_{n_{1}}:=\\{(X_{1},Y_{1}),\ldots,(X_{n_{1}},Y_{n_{1}})\\}\in\mathbb{R}^{d}\times\\{0,1\\}$.
We want to classify the new data $X$. Set $D:=(d^{2}+d)/2$.

Algorithm:

  * •

Choose an sm-uniqueness set of $D$ directions in $\mathbb{R}^{d}$.

  * •

We start using the data from the subsample
$\mathcal{X}_{n}\setminus\mathcal{X}_{n_{1}}$ of the training sample to assign
weights $w_{j}$ to each direction as follows.

    1. (i)

For each direction $u_{j},~{}j\in\\{1,\ldots,D\\}$, consider the two
subsamples

(4)
$\\{(Z_{1,j},Y_{1}),\ldots,(Z_{n_{1},j},Y_{n_{1}})\\}\in\mathbb{R}\times\\{0,1\\},$
(5)
$\\{(Z_{{}_{n{{}_{1}}+1,j}},Y_{n{{}_{1}}+1}),\ldots,(Z_{n,j},Y_{n})\\}\in\mathbb{R}\times\\{0,1\\},$

where $Z_{i,j}=\langle u_{j},X_{i}\rangle,~{}i\in\\{1,\ldots,n\\}$.

    2. (ii)

Apply the classical $k$-NN rule based on the subsample (4) to classify
$Z_{i,j},i\in\\{n_{1}+1,\ldots,n\\}$. That is, let the corresponding score be
given by

$\displaystyle S_{X_{i},j}:=\\{$ number of neighbors of the subsample (4)
$\displaystyle\textrm{with label $1$ among the $k$-nearest data to
}Z_{i,j}\\}.$

    3. (iii)

Set $p_{X_{i},j}:=S_{X_{i},j}/k$. If $p_{X_{i},j}>1/2$ then classify the new
data $Z_{i,j}$ as belonging to class $1$ with $i\in\\{n_{1}+1,\ldots,n\\}$
(otherwise classify it as class $0$). Denote by $\widehat{Y}_{i,j}$ the label
assigned to $Z_{i,j}$.

    4. (iv)

We define

$w_{j}:=\frac{1}{n-n_{1}}\sum_{i=n_{1}+1}^{n}\mathcal{I}_{\widehat{Y}_{i,j}=Y_{i}}.$

  * •

Lastly, given a new datum $X$, if

$\sum_{j=1}^{D}w_{j}p_{X,j}\geq\sum_{j=1}^{D}w_{j}(1-p_{X,j}),$

then classify $X$ as belonging to class $1$. Otherwise classify it as class
$0$.

######  Remarks.

(i) One way to choose the set $\mathcal{D}$ is as follows. We order the values
$w_{j},\ j\in\\{1,\ldots D\\}$, and take the $\Delta\%$ largest values of
$w_{j}$, i.e., we let $\\{w^{(1)},\ldots w^{(D)}\\}$ be the order statistics
and $R_{1},\ldots R_{D}$ be the corresponding rank statistics. Then we set

$\mathcal{D}:=\\{j:R_{j}\geq\Delta\times D\\}.$

Alternatively, we can initially take more than $D$ directions, and then choose
the $D$ best ones.

(ii) Using the weights in the last steps, we can take into account that some
directions may not be very good for classification on the projection, see Fig.
6. However, we need to split the training sample into two subsamples.

(a)  Fig. 6. Illustration of Remark (ii) in §6.1: empirical marginal densities
in two random projections of two groups in dimension 2.

###  6.2. An example with simulated data

To illustrate the performance of the proposed algorithm, we consider a binary
classification problem with two classes of size 1000 in dimension 50. The data
in both classes are iid and marginally independent with Cauchy distribution.
The data in group 1 are centred at the origin. The data in group 2 are shifted
in the direction of the vector
$\eta=(\eta_{1},\ldots,\eta_{50})\in\mathbb{R}^{50}$, where $\eta_{i}$ is
chosen uniformly in $[0,3]$ for $i\in\\{1,\ldots,10\\}$, and $\eta_{i}=0$ for
$i\in\\{11,\ldots,50\\}$.

We consider $75\%$ of the data as the training sample and the remaining data
as the test sample. We repeat the experiment 100 times under the same
conditions. At each replication, the classifier accuracy is calculated.

Fig. 7 shows the accuracy of different classifiers: random forest (RF),
support vector machine (SVM), generalized linear model (GLM), and our
algorithm using random projections (RP). In the RP algorithm, we take
$\Delta=50\%$ and $\omega=0.25$. It is observed that our algorithm performs
better than the other competitors.

(a)  Fig. 7. Accuracy in the test sample of the RF, SVM, GLM and RP
classifiers in the example in §6.2.

##  7\. Conclusion

In this article we have studied the problem of which sets $\mathcal{L}$ of
lines in $\mathbb{R}^{d}$ determine elliptical distributions in the sense
that, if $P,Q$ are elliptical distributions whose projections satisfy
$P_{L}=Q_{L}$ for all $L\in\mathcal{L}$, then $P=Q$. Combining our two main
results, Theorems 1 and 2, we have a precise characterization of such
$\mathcal{L}$. In particular, there exist such sets $\mathcal{L}$ of
cardinality $(d^{2}+d)/2$, and this number is best possible.

We have applied our results to derive a Kolmogorov–Smirnov type test for
equality of elliptical distributions, and we have carried out a small
simulation study of this test, comparing it with other tests from the
literature. The results are quite encouraging. As a further application of our
results, we have proposed an algorithm for binary classification of data
arising from elliptical distributions, again supported by a small simulation.

Let us conclude with a question. As mentioned in the introduction, our main
result, Theorem 1, is, in a certain sense, a continuous analogue of Heppes’
result on finitely supported distributions. Recently, the authors obtained the
following quantitative form of Heppes’ theorem, expressed in terms of the
total variation metric:

$d_{TV}(P,Q):=\sup\\{|P(B)-Q(B)|:B\text{~{}Borel}\\}.$

######  Theorem 3 ([10, Theorem 2.1]).

Let $Q$ be a probability measure on $\mathbb{R}^{d}$ whose support contains at
most $k$ points. Let $H_{1},\dots,H_{k+1}$ be subspaces of $\mathbb{R}^{d}$
such that $H_{i}^{\perp}\cap H_{j}^{\perp}=\\{0\\}$ whenever $i\neq j$. Then,
for every Borel probability measure $P$ on $\mathbb{R}^{d}$, we have
$d_{TV}(P,Q)\leq\sum_{j=1}^{k+1}d_{TV}(P_{H_{j}},Q_{H_{j}})$.

Is there likewise a quantitative version of Theorem 1?

## Acknowledgments

Fraiman and Moreno were supported by grant FCE-1-2019-1-156054, Agencia
Nacional de Investigación e Innovación, Uruguay. Ransford was supported by
grants from NSERC and the Canada Research Chairs program.

### Supplementary material

All the code in R is available at:

https://github.com/mrleomr/CRAMER-WOLD-ELLIPTICAL

## References

  * [1] C. Bélisle, J.-C. Massé, and T. Ransford, _When is a probability measure determined by infinitely many projections?_ , Ann. Probab. 25 (1997), no. 2, 767–786. MR 1434125
  * [2] S. Cambanis, S. Huang, and G. Simons, _On the theory of elliptically contoured distributions_ , J. Multivariate Anal. 11 (1981), no. 3, 368–385. MR 629795
  * [3] F. Chen, M. D. Jiménez-Gamero, S. Meintanis, and L. Zhu, _A general Monte Carlo method for multivariate goodness–of–fit testing applied to elliptical families_ , Comput. Statist. Data Anal. (2022), 107548.
  * [4] H. Cramér and H. Wold, _Some theorems on distribution functions_ , J. London Math. Soc. 11 (1936), no. 4, 290–294. MR 1574927
  * [5] J. A. Cuesta-Albertos, R. Fraiman, and T. Ransford, _Random projections and goodness-of-fit tests in infinite-dimensional spaces_ , Bull. Braz. Math. Soc. (N.S.) 37 (2006), no. 4, 477–501. MR 2284883
  * [6] by same author, _A sharp form of the Cramér-Wold theorem_ , J. Theoret. Probab. 20 (2007), no. 2, 201–209. MR 2324526
  * [7] A. Cuevas and R. Fraiman, _On depth measures and dual statistics. A methodology for dealing with general data_ , J. Multivariate Anal. 100 (2009), no. 4, 753–766. MR 2478196
  * [8] G. R. Ducharme and P. Lafaye de Micheaux, _A goodness-of-fit test for elliptical distributions with diagnostic capabilities_ , J. Multivariate Anal. 178 (2020), 104602, 13. MR 4076341
  * [9] K. Fang, S. Kotz, and K. Ng, _Symmetric multivariate and related distributions_ , Monographs on Statistics and Applied Probability, vol. 36, Chapman and Hall, Ltd., London, 1990. MR 1071174
  * [10] R. Fraiman, L. Moreno, and T. Ransford, _A quantitative Heppes theorem and multivariate Bernoulli distributions_ , J. R. Stat. Soc. Ser. B. Stat. Methodol. (2023), qkad003.
  * [11] by same author, _Application of the Cramér–Wold theorem to testing for invariance under group actions_ , Preprint, 2022.
  * [12] W. M. Gilbert, _Projections of probability distributions_ , Acta Math. Acad. Sci. Hungar. 6 (1955), 195–198. MR 70868
  * [13] K. Gröchenig and P. Jaming, _The Cramér-Wold theorem on quadratic surfaces and Heisenberg uniqueness pairs_ , J. Inst. Math. Jussieu 19 (2020), no. 1, 117–135. MR 4045081
  * [14] P. Hall and M. Martin, _On the bootstrap and two-sample problems_ , Austral. J. Statist. 30A (1988), no. 1, 179–192.
  * [15] M. Hallin, G. Mordant, and J. Segers, _Multivariate goodness-of-fit tests based on Wasserstein distance_ , Electron. J. Stat. 15 (2021), no. 1, 1328–1371. MR 4255302
  * [16] G. G. Hamedani and M. N. Tata, _On the determination of the bivariate normal distribution from distributions of linear combinations of the variables_ , Amer. Math. Monthly 82 (1975), no. 9, 913–915. MR 383643
  * [17] H. Hedenmalm and A. Montes-Rodríguez, _Heisenberg uniqueness pairs and the Klein-Gordon equation_ , Ann. of Math. (2) 173 (2011), no. 3, 1507–1527. MR 2800719
  * [18] A. Heppes, _On the determination of probability distributions of more dimensions by their projections_ , Acta Math. Acad. Sci. Hungar. 7 (1956), 403–410. MR 85646
  * [19] M. Hyodo and T. Nishiyama, _Simultaneous testing of the mean vector and covariance matrix among $k$ populations for high-dimensional data_, Comm. Statist. Theory Methods 50 (2021), no. 3, 663–684. MR 4204361
  * [20] J. Lim, E. Li, and S.-J. Lee, _Likelihood ratio tests of correlated multivariate samples_ , J. Multivariate Anal. 101 (2010), no. 3, 541–554. MR 2575403
  * [21] B. G. Manjunath and K. R. Parthasarathy, _A note on Gaussian distributions in $\mathbb{R}^{n}$_, Proc. Indian Acad. Sci. Math. Sci. 122 (2012), no. 4, 635–644. MR 3016817
  * [22] J. P. Nolan, _Multivariate elliptically contoured stable distributions: theory and estimation_ , Comput. Statist. 28 (2013), no. 5, 2067–2089. MR 3107292
  * [23] J. T. Praestgaard, _Permutation and bootstrap Kolmogorov-Smirnov tests for the equality of two distributions_ , Scand. J. Statist. 22 (1995), no. 3, 305–322. MR 1363215
  * [24] A. Rényi, _On projections of probability distributions_ , Acta Math. Acad. Sci. Hungar. 3 (1952), 131–142. MR 53422
  * [25] G. J Székely, M. L. Rizzo, et al., _Testing for equal distributions in high dimension_ , InterStat 5 (2004), no. 16.10, 1249–1272.


&=\int e^{-x^2/\Delta^2} p(s) ds-\int_{-\infty}^{-x}e^{-s^2/\Delta^2}p(s) ds-\int _x^\infty e^{-s^2/\Delta^2} p(s) ds\\
&\geq 1-\eta -e^{-x^2/\Delta^2}\left(\int_{-\infty}^{-x} p(s) ds+\int _x^\infty p(s)ds\right)\\
&\geq 1-\eta-e^{-x^2/\Delta^2}\int p(s) ds\\
&\geq 1-\eta-e^{-x^2/\Delta^2}.
\end{align}
\end{subequations}
\end{preuve}

\subsection{Proof of proposition \ref{pro: concentration g delta}}\label{proof: concentration g delta}
\begin{preuve}
We start with the proof of (\ref{eq: concentration g delta V1}). With the notation defined just before the formulation of the proposition, the hypothesis $O^\Delta\geq 1-\eta$ translate to $\int ds e^{-s^2/\Delta^2} p(s)\geq 1-\eta$. We also denote $q(s)=d\P[\hat g_i+\delta_j=s]=\frac{1}{\mu\sqrt{\pi}}\int dg p(g) e^{-(s-g)^2/\mu^2}$. We thus have:

\begin{subequations}
\begin{align}
\P[\abs{\hat g_i+\delta_j}\leq x]&=\int _{-x}^x q(s)ds\\
&\geq\int_{-x}^x e^{-\frac{x^2}{\Delta^2-\mu^2}}q(s)ds\\
&=\int e^{-\frac{x^2}{\Delta^2-\mu^2}}q(s)ds-\int_{-\infty}^{-x}e^{-\frac{x^2}{\Delta^2-\mu^2}}q(s)ds-\int_x^\infty e^{-\frac{x^2}{\Delta^2-\mu^2}}q(s)ds\\
&\geq \frac{1}{\mu\sqrt{\pi}}\int e^{-\frac{s^2}{\Delta^2-\mu^2}}e^{-(s-g)^2/\mu^2} p(g)dsdg-e^{-\frac{x^2}{\Delta^2-\mu^2}}\left(\int_{-\infty}^{-x} q(s)ds+\int_x^\infty q(s)ds\right)\\
&\geq \sqrt{1-\left(\frac{\mu}{\Delta}\right)^2}\int e^{-g^2/\Delta^2}p(g)dg -e^{-\frac{x^2}{\Delta^2-\mu^2}}\int q(s) ds\\
&\geq (1-\eta)\sqrt{1-\left(\frac{\mu}{\Delta}\right)^2}-e^{-\frac{-x^2}{\Delta^2-\mu^2}}.
\end{align}
\end{subequations}
For the second inequality (\ref{eq: concentration g delta V2}), we write:
\begin{subequations}
\begin{align}
\P[\abs{\hat g_i+\delta_j}\leq x ]&\geq \P[\abs{\hat g_i}\leq x/2~\&~\abs{\delta_j}\leq x/2]\\
&=\P[\abs{\hat g_i}\leq x/2]\P[\abs{\delta_j}\leq x/2]\\
&\geq (1-\eta -e^{-x^2/4\Delta^2})(1-\frac{2\mu}{x}e^{-x^2/4\mu^2})\\
&\geq 1-\eta-e^{-x^2/4\Delta^2}-\frac{2\mu}{x}e^{-x^2/4\mu^2},
\end{align}
\end{subequations}
using the bound on gaussian distribution
\begin{subequations}
\begin{align}
\P[\abs{\delta_j}\leq x]&=1-2\int_x^\infty e^{-s^2/\mu^2}ds\\
&=1-2\mu\int_{x/\mu}^\infty e^{-t^2}dt\\
&\geq 1-\frac{\mu}{x}\int _{x/\mu}^\infty  2te^{-t^2}\\
\end{align}
\end{subequations}
The last one (\ref{eq: concentration g delta V3}) is obtained by taking a generic parameter $t$ and computing
\begin{subequations}
\begin{align}
\P[\abs{\hat g_i+\delta_j}\leq x ]&\geq \P[\abs{\hat g_i}\leq t~\&~\abs{\delta_j}\leq (1-t)]\\
&=\P[\abs{\hat g_i}\leq t]\P[\abs{\delta_j}\leq (1-t)]\\
&\geq (1-\eta -e^{-x^2t^2/\Delta^2})(1-\frac{\mu}{(1-t)x}e^{-x^2(1-t)^2/\mu^2}).
\end{align}
\end{subequations}
Now, if we set $t=\frac{\Delta}{\mu +\Delta}$ such that $1-t=\frac{\mu}{\Delta+\mu}$, we get

\begin{equation}
\P[\abs{\hat g_i+\delta_j}\leq x ]\geq 1-\eta-\left(1+\frac{\mu+\Delta}{x}\right)\exp(-\frac{x^2}{(\mu+\Delta)^2}).
\end{equation}

\end{preuve}

\end{document}

%Removed content, kept at the end for safekeeping

\section{Things to remove?}
\begin{Pro}
The set of translated graph states:
\begin{equation}
\Big\{\prod_{i=1}^n \hat Z_i(s_i)\ket{G}\Big|\forall i,s_i\in \R\Big\},
\end{equation}
constitute a (singular) basis of the Hilbert space.
\end{Pro}

\begin{Pro}
If one execute a nullifier check by $\hat g_i$ on a translated graph state $\prod_{j=1}^n Z_j(s_j)\ket{G}$, then we get:
\begin{equation}
\hat g_i \prod_{j=1}^n Z_j(s_j)\ket{G}=s_i \prod_{j=1}^n Z_j(s_j)\ket{G},
\end{equation}
\end{Pro}

\begin{Pro}
The wave function of a vacuum squeezed stat $\ket{\psi_\sigma}$ of squeezing parameter $\zeta=\ln(\sigma)$ is:
\begin{equation}
\psi_\sigma(x)=\frac{1}{(\pi\sigma^2)^{1/4}}\exp[-\frac{x^2}{2\sigma^2}],
\end{equation}
so that, $\ket{\psi_\sigma}=\int dx \psi_\sigma(x)\ket{x}$.
\end{Pro}

If we take, $σ=∞$, then $|ψ_∞⟩=|0⟩^p$ is the state of zero momentum defined previously.
Computation of the two different POVM
\begin{subequations}
\begin{align}
\hat M_g&=\frac{1}{\mu\nu \pi}\int dp~dx~dp'~dx'~\exp[-\frac{(p'-p)^2}{\nu^2}-\frac{(x'-x)^2}{\mu^2}]\delta(g-p+x)\ketbra{p'}\otimes\ketbra{x'}\\
&=\frac{1}{\mu\nu \pi}\int dp~dp'~dx'~\exp[-\frac{(p'-p)^2}{\nu^2}-\frac{(x'-p+g)^2}{\mu^2}]\ketbra{p'}\otimes\ketbra{x'}\\
&=\frac{1}{\mu\nu \pi}\frac{\sqrt{\pi}}{\sqrt{\frac{1}{\nu^2}+\frac{1}{\mu^2}}}\int dp'~dx'~\exp[-\frac{(p'-x'-g)^2}{\nu^2+\mu^2}]\ketbra{p'}\otimes\ketbra{x'}\\
&=\frac{1}{\sqrt{\pi}\sqrt{\nu^2+\mu^2}}\int dp'~dx'~dg'~ \delta(g'-p'+x')\exp[-\frac{(g'-g)^2}{\nu^2+\mu^2}]\ketbra{p'}\otimes\ketbra{x'}\\
&=\frac{1}{\delta\sqrt{\pi}}\int dg'~\exp[-\frac{(g-g')^2}{\delta^2}]\hat P_{g'}\\
&=\hat E_g,
\end{align}
\end{subequations}
if $δ^2=ν^2+μ^2$. This shows that the two point of view coincide.

The full theorem 4
\begin{Thm}
{\bf (Soundness)} For any initial state $\hat \rho$ and integer $k\leq N_r:= (\mu-n)N_\text{test}$ we have:
\begin{equation}\label{eq: soundness with prior}
\Tr(\hat \Pi_\text{acc}\otimes[\1-\hat\Pi_1\times\cdots\times\hat \Pi_k]\hat \rho)\leq \Tr(\hat\Pi_\text{acc}\hat\rho)\left(1-\frac{\binom{N}{k}}{\binom{N_r}{k}}\right)+n\exp(-\frac{\nu^2 N_\text{test}}{1+\frac{1}{\mu-n}})\frac{\binom{N}{k}}{\binom{N_r}{k}},
\end{equation}
with as before $N=nN_\text{test}(\mu/n-1-\mu f-\mu\nu)$ and under the constraint that $\mu/n\geq 1+\mu f+\mu\nu$. Using the simple inequality $\Tr(\hat \Pi_\text{acc}\hat \rho)\leq 1$, we get a bound the bound
\begin{equation}
\Tr(\hat \Pi_\text{acc}\otimes[\1-\hat\Pi_1\times\cdots\times\hat \Pi_k]\hat \rho)\leq  1-\frac{\binom{N}{k}}{\binom{N_r}{k}}\left [1-n\exp(-\frac{\nu^2 N_\text{test}}{1+\frac{1}{\mu-n}})\right].
\end{equation}
We can also lower bound the quotient of binomial coefficient to get a simpler expression
\begin{equation}
\Tr(\hat \Pi_\text{acc}\otimes[\1-\hat\Pi_1\times\cdots\times\hat \Pi_k]\hat \rho)\leq  1-\left(1-\frac{knN_\text{test}\mu(f+\nu)}{(\mu-n)N_\text{test}-k+1}\right)\left [1-n\exp(-\frac{\nu^2 N_\text{test}}{1+\frac{1}{\mu-n}})\right].
\end{equation}
Starting back from equation (\ref{eq: soundness with prior}), we can get the conditional statement
\begin{equation}
\Tr(\hat\Pi_1\times\cdots\times\hat\Pi_k \hat \rho_r)\geq \frac{\binom{N}{k}}{\binom{N_r}{k}}\left(1-\frac{n\exp(-\frac{\nu^2 N_\text{test}}{1+\frac{1}{\mu-n}})}{\Tr(\hat\Pi_ \text{acc}\hat\rho)}\right),
\end{equation}
and with the lower bound on the binomial coefficients
\begin{equation}
\Tr(\hat\Pi_1\times\cdots\times\hat\Pi_k \hat \rho_r)\geq \left(1-\frac{knN_\text{test}\mu(f+\nu)}{(\mu-n)N_\text{test}-k+1}\right)\left(1-\frac{n\exp(-\frac{\nu^2 N_\text{test}}{1+\frac{1}{\mu-n}})}{\Tr(\hat\Pi_ \text{acc}\hat\rho)}\right),
\end{equation}
where $\hat \rho_r$ is the state over the remaining register after a successful protocol.
\end{Thm}

The right hand side of the eq. (\ref{eq: soundness original}) and (\ref{eq: soundness lambda}) are the product of two terms. The right most one, has a dependence in $(Π̂_accρ̂)$, which is the {\it a priori} probability to accept the the state. However, in practice this probability, which is link to how much Bob is trusted by Alice, is not known. Luckily, we can modify a little bit the choice of the parameters by introducing a parameter $α$,
\begin{align}
\mu=2n && N_\text{test}=\left\lceil (2\ln n+2\alpha\ln \lambda)\lambda^2\right\rceil && \nu=1/\lambda && f=1/\lambda,
\end{align}
we get
\begin{equation}
F^{\sqrt{\epsilon^2+\delta^2}}\geq \left(1-\frac{4n}{\lambda}\right)\left(1-\frac{1}{\lambda^\alpha\Tr(\hat\Pi_\text{acc}\hat\rho)}\right).
\end{equation}

By doing so, we get a different exponent in the scaling of the right most term. This mean that, this term, is more quickly suppressed than the other one, and the lack of knowledge on the term $(Π̂_accρ̂)$ is less relevant.\\

\begin{equation}
O_k^{\sqrt{\epsilon^2+\delta^2}}\geq \left(1-\frac{4n}{\lambda}\frac{k}{1-\frac{k-1}{nN_\text{test}}}\right)\left(1-\frac{1}{\lambda\Tr(\hat\Pi_\text{acc}\hat\rho)}\right).
\end{equation}
or a simpler one if we are a bit crude with the inequalities

\subsection{Quantum steering}

When considering gaussian states, it has been shown that whether steering with gaussian measurement is possible, can simply be checked by looking at the correlation matrix. Since, if a state is steerable by gaussian measurement then it is steerable with any measurement, this give us a sufficient condition to verify if a Gaussian state is steerable. Notice however that this result seems to be applicable on the gaussian state.\\

To define the correlation matrix, denote $R=(x̂_1,p̂_1,x̂_2,p̂_2)$. Then the correlation matrix $M$ is given by
\begin{equation}
M_{i,j}=\frac{1}{2}\expval{\hat R_i\hat R_j+\hat R_j\hat R_i}-\expval{\hat R_i}\expval{\hat R_j}.
\end{equation}
The result is then that the state is steerable by the second party with gaussain measurement if and only if
\begin{equation}
M+\frac{i}{2}(\omega\oplus 0)\geq 0,
\end{equation}
as a matrix inequality and where $ω=
0   1
-1   0

If we consider a finitely squeezed version of a two mode graph state with wave function
\begin{equation}
f(x_1,x_2)= \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{x_1^2+x_2^2}{4\sigma^2}+ix_1x_2\Omega),
\end{equation}
we can compute its correlation matrix and get
\begin{equation}
\sigma^2 & 0 & 0 & \sigma^2\Omega\\
0 & \frac{1}{4\sigma^2}+\sigma^2\Omega^2 & \sigma^2\Omega & 0\\
0 & \sigma^2\Omega & \sigma^2 & 0 \\
\sigma^2\Omega & 0 & 0 & \frac{1}{4\sigma^2} +\sigma^2\Omega^2.
\end{pmatrix}
\end{equation}
As a sanity check, we can verify that the variance of the nullifier is small
\begin{subequations}
\begin{align}
\operatorname{Var}(\hat g_1)&=\expval{(\hat p_1-\Omega\hat x_2)^2}\\
&=\Omega^2\expval{\hat x_2^2}+\expval{\hat p_1^2}-\Omega\expval{\hat x_2\hat p_1+\hat p_1\hat x_2}\\
\end{align}
\end{subequations}
which indeed goes to zero with infinite squeezing ($σ→∞$). To check the steerability condition, we simply need to verify when the following matrix is positive:
\begin{equation}
\begin{pmatrix}
\sigma^2 & \frac{i}{2} & 0 & \sigma^2\Omega\\
-\frac{i}{2} & \frac{1}{4\sigma^2}+\sigma^2\Omega^2 & \sigma^2\Omega & 0\\
0 & \sigma^2\Omega & \sigma^2 & 0 \\
\sigma^2\Omega & 0 & 0 & \frac{1}{4\sigma^2} +\sigma^2\Omega^2.
\end{pmatrix}
\end{equation}
By feeding it in a software, we see that whenever $Ω≠0$, it always has at least one negative eigenvalue. This means that, whatever the squeezing and the strength of the correlation, these states display steering capabilities.\\

It turns out that this condition is also sufficient for general states. So it suffice to study the correlation matrix to show if a state is steerable.

\subsection{Comparison}
Our result is two fold: a bound on the probability to falsely accepted a state, and a lower bound on overlap for states that have been accepted. As stated before the first bound reads
\begin{equation}
\Tr((\hat\Pi_\text{acc}\otimes[\1-\hat\Pi_\text{good}^N])\hat\rho)\leq \epsilon:=n \exp(-\frac{\nu^2N_\text{test}}{2}),
\end{equation}
Looking at it directly it seems that we provide an exponential scaling, with $N_total=2μ/ν^2ln(n/ϵ)$. However, the condition we state is formulated with the POVM $Π_good^N$ which count whether there are at least $N$ `good' state among the remaining one. Being accepted by this POVM is thus a less strict condition that the one formulated on eq. (\ref{eq: soundess wu}). In order for this bound to be relevant, we need in the asymptotic limit that $N/N_r→1$.\\

If we put the choice of the parameters defined on eq. (\ref{eq: choice with lambda}) we have $ϵ=1/λ$ and $N=(1-4n/λ)N_r$, for $N_test=⌈(2lnn+2lnλ)λ^2⌉$. Thus the number of copies scales as
\begin{equation}
N_\text{total}=\mathcal O\left(\frac{\ln(\frac{1}{\epsilon})}{\epsilon^2}\right),
\end{equation}
which is still substantially better than the previous result and can be compared to the i.i.d case, which seems unreasonable. However, when discussing the specialized result (section \ref{sec: special cases}), we gave another choice, which leads to a substantially better scaling. For this choice, we had $N_test=⌈(2lnn+2αlnλ)λ^2⌉$ and $ϵ=1/λ^α$. Playing with the relations, we get that in this case
\begin{equation}
N_\text{total}=\mathcal O\left(\frac{\ln(\frac{1}{\epsilon})}{\epsilon^{2/\alpha}}\right),
\end{equation}
which is a even better scaling as $α$ can be any arbitrarily large. Even worse, we can make the choice $N_test=⌈(2lnn+2 λ)λ^3⌉$, to get $ϵ=e^-λ$
\begin{equation}
N_\text{total}=\mathcal O\left(\operatorname{Poly}(\ln(\frac{1}{\epsilon}))\right),
\end{equation}
which is even crazier. However, we have to keep in mind that although the scaling is much better, the statement it self in weaker, in the sens that whenever $ν≠0$, we have $N<N_r$ and thus we are only certifying a fraction of the remaining states and not all of them. This can be directly understood, when looking at the second bound  (Thm. \ref{thm: fidelity bound general})
\begin{equation}
F^{\sqrt{\epsilon^2+\delta^2}}\geq \left(1-\frac{n\mu(f+\nu)}{\mu-n}\right)\left(1-\frac{n\exp(-\nu^2 N_\text{test}/2)}{\Tr(\hat\Pi_\text{acc}\hat\rho)}\right).
\end{equation}
We see that $F^√(ϵ^2+δ^2)$ is bounded from bellow by a product of terms that are close to $1$. In order to understand the scaling, one simply has to study the scaling of each term and consider the worst one. In all the previous discussion, the choices for the parameters $μ=2n$, $ν=f=1/λ$ where the same. So the first term of the product is $1-4n/λ$ in all these three situations (notice that this expression is independent of the definition of $N_test$). The scaling of this term is clearly worse than the one we get for the second term when inputting the different choice of the parameters. Meaning that in all situations, the scaling  of the number of registers so that $O^Δ≥ϵ$, is
\begin{equation}
N_\text{total}=\mathcal O\left(\frac{\ln(\frac{1}{\epsilon})}{\epsilon^2}\right).
\end{equation}
This means that in order to do comparison with other protocol, one need to try considering overlap bound.
and for $k=1$.
\begin{align}
\Tr(\hat \Pi_\text{acc}\otimes[\1-\hat\Pi_1]\hat \rho)\leq\frac{4n+1}{\lambda}.
\end{align}

\end{document}