# Exploring Lossy Compressibility through Statistical Correlations of
Scientific Datasets

David Krasowska1, Julie Bessac2, Robert Underwood1, Jon C. Calhoun1, Sheng
Di2, and Franck Cappello2  1 Holcombe Department of Electrical and Computing
Engineering, Clemson University Clemson, USA  3 Mathematics and Computer
Science DivisionArgonne National Laboratory Lemont, USA  3 School of
Computing, Clemson University Clemson, USA

###### Abstract

Lossy compression plays a growing role in scientific simulations where the
cost of storing their output data can span terabytes. Using error bounded
lossy compression reduces the amount of storage for each simulation; however,
there is no known bound for the upper limit on lossy compressibility.
Correlation structures in the data, choice of compressor and error bound are
factors allowing larger compression ratios and improved quality metrics.
Analyzing these three factors provides one direction towards quantifying lossy
compressibility. As a first step, we explore statistical methods to
characterize the correlation structures present in the data and their
relationships, through functional models, to compression ratios. We observed a
relationship between compression ratios and statistics summarizing correlation
structure of the data, which are a first step towards evaluating the
theoretical limits of lossy compressibility used to eventually predict
compression performance and adapt compressors to correlation structures
present in the data.

###### Index Terms:

compression, lossy compression, high performance computing, statistical
correlation analysis

##  I Introduction

Scientific research increasingly uses error-bounded lossy compressors to
achieve greater compression ratios in relation to lossless compressors[1].
This improved performance allows applications to run with larger and more
frequently produced datasets due to faster I/O times and smaller I/O volumes.
The theoretical limit on compressibilty of data using lossless compression is
given by the entropy[2]. The entropy quantifies the information content
present in a symbol from a source sequence based upon its probability of
occurrence. Thus, for a given sequence of symbols the entropy enables
computing the minimum number of bits, on average, needed to represent the
data. For over 70 years, this concept has guided the development and
evaluation of lossless compression algorithms. However, for lossy compression
algorithms, there is currently no known bound for the maximum degree of
compression that can be achieved for some specified point-wise error bound.
Establishing this bound will allow researchers to anticipate compression
performance and adapt compressors to correlation structures of the data
ensuring they get the best compression ratio possible. Similarly, establishing
the limit for lossy compression allows for the maximum efficiency for storing
large scientific datasets.

One possible direction to establish an entropy-like bound for lossy
compression is studying the impact on compressibility of correlation
structures (correlations in space, time, or other dimensions) of the data,
compressor types and error bounds. We refer to “compressibility” as the
maximum compression ratio associated with a given error bound. As many
compressors implicitly or explicitly exploit correlation structures of the
data (e.g. SZ [3], ZFP [4]), analyzing the relationships between correlation
structures and compression performances will allow researchers to anticipate
compression performances. Ultimately, we seek to establish an entropy-like
metric for lossy compression algorithms which can guide the lossy compression
community to optimal development and usage by adapting compressors to
correlation structures present in the data.

In this work, we focus on data correlations and their link to compression
ratios for several compressors. The goal of the research paper is to explore:

  1. 1.

statistical methods to characterize the correlation structures of the data and

  2. 2.

their relationships, through functional models, to compression ratios.

These models will form the first step into evaluating the theoretical limits
of lossy compressibility used to eventually predict compression performance.
In particular, we focus on estimated correlation ranges through variograms and
its effect on compressibility. The variogram is commonly used in geostatistcis
to estimate second-order dependence and more precisely how data are correlated
with distance. Under stationary conditions, the variogram and covariance
function have a direct correspondence. However, in practice the variogram can
be estimated in more relaxed assumptions than the covariance and thus is
preferred by practitioners. To characterize compressibility, we use the
compression ratio, which is an important statistic within lossy compression
due to its positive correlation with storing and processing as much data as
possible.

In this study, we focus on 2D-gridded datasets and through variogram analysis,
perform a characterization of the correlation between grid-points along with a
compression analysis (Section IV-C) with several compressors (Section IV-B).
We consider two types of data (Section IV-A), synthetic datasets consisting of
stochastic correlated Gaussian fields with known correlation structures and
another one consisting of simulations from a hydrodynamical model Miranda [5].

##  II Background

###  II-A Compressors

SZ [3], ZFP [4], and MGARD [6] are some of the leading error bounded lossy
compressors. In this section we explain from an algorithmic perspective how
these compressors exploit correlations in data.

SZ scans through the data block by block: for 2D data the block size is
$16\times 16$. For each block, it makes a prediction of what the data in the
block is using one of two methods: the Lorenzo predictor or the regression
predictor. The Lorenzo predictor uses the values of neighboring points to
estimate the value at the current position. The regression predictor fits a
hyper-plane through the block, and uses the fitted hyper plane to interpolate
the values within the blocks. If these predictions are linearly quantized, and
if the quantized values are sufficiently accurate according to the error
bound, the quantized values are stored. If the predictions are not accurate,
the values are stored exactly. Finally the entire sequence is passed first to
a huffman encoding, and then to the Zstd lossless compressor to exploit
patterns in the quantized sequence. However, since the the predictor does not
observe values outside of its block, it cannot exploit global correlation
structures easily.

ZFP likewise uses local correlations, but uses a completely different
compression principle based on near orthogonal transforms – a similar approach
used by JPEG image compression. ZFP first partitions 2D data into blocks of
size $4\times 4$, then it converts each block of floating point data into a
common fixed point representation, applies the near orthogonal transform,
applies an embedded encoding that orders bits from most significant to least
significant and then applies a truncation to archive a desired tolerance.
Again, because the blocks are compressed independently, the compressor cannot
acquire a global knowledge of the data’s correlations structures.

MGARD however is a newer compressor uses a different approach which can
account for global correlation structures. MGARD relies on the mathematical
theory of multi-grid methods in order to compress the data. MGARD operates by
decomposing the data into multi-level coefficients which represent recursively
defined sub-regions until the block is represented within the allowed
tolerance. These multi-level coefficients are then quantized and then
compressed with either Zlib (in older versions) or Zstd (in the newest
unreleased version). Because these multi-level coefficients can represent
regions of differing sizes and potentially the entire dataset, MGARD can
capture multi-level effects in a way that SZ and ZFP cannot making it an
important comparison for our paper.

###  II-B Variogram

The variogram is a function describing the degree of dependence of a
correlated spatial field [7, 8], it gives a measure of how much two points of
the same field correlate depending on the distance $h$ between those points.
The variogram is a function of the distance $h$ and is typically characterized
by several parameters: nugget (microscale variability), sill (variance of the
studied field), range (distance at which autocorrelation vanishes), see Figure
1.

Figure 1: Example of a variogram as a function of the distance $h$ between
points.

In the following, we focus solely on the variogram range $a$ as it corresponds
to the distance ($h$) where the variogram ($\gamma$) plateaus (at the sill
value) and indicates the distance beyond which the spatial correlation among
grid-points vanishes. Intuitively, the larger the range the stronger the
correlation is across grid-points.

In practice, the empirical semi-variogram is computed on the data via Equation
(1):

$\displaystyle\gamma(h)=\frac{1}{2N(h)}\sum_{|x_{i}-x_{j}|=h}^{N(h)}\left(z(x_{i})-z(x_{j})\right)^{2},$
(1)

where $z$ is the studied field of interest (e.g. velocityx from the Miranda
dataset), $x_{i}$ and $x_{j}$ are grid-points coordinates or indexes, $N(h)$
is the number of grid-points at distance $h$ from each other. The variogram
corresponds to $2\gamma(h)$, in practice the terms semi-variogram and
variogram are used interchangeably. In the following, we will refer to
variogram as $\gamma(h)$. Finally, to estimate the variogram range $a$ we fit
by least-squares a parametric squared-exponential variogram
($\gamma(h)=c_{0}(1-\exp(-h^{2}/a))$) to the empirical variogram range
estimated via Equation (1).

###  II-C Complex correlation patterns

Statistical tools exist beyond the variogram to quantify and extract complex
correlation structures of datasets for instance in order to analyze long-range
dependencies [9], to detect change-point in time series [10], to cluster and
identify correlation regimes [11]. Identifying multiscale components of
scientific datasets mostly relies on eigen or a basis-function decompositions
such as singular value decomposition (SVD) or wavelet decomposition [12, 13].
However due to the often extreme complexity of correlations and dependencies
in data, developing methods to extract for instance spatial and spatio-
temporal heterogeneity or non-stationarity is still on-going research.
Although the detailed use of these techniques is outside the scope of the
current work and left for future efforts, we provide preliminary results and
future work directions with SVD at the end of the Section V.

##  III Related Work

Beyond the classical efforts to estimate compressiblity using entropy, there
has been relatively little attention afforded to the topic of lossy
compressiblity. As a parallel in data reduction techniques, [14] investigated
the determination of thresholds for singular value decomposition of large
matrices based on some optimality loss criteria. [15] identified several
factors that affect compression ratio for SZ and ZFP. They relied on block-
based sampling based approach which considered individual data points tailored
to each considered compressor that they consider and used Gaussian model to
estimate the subsequent compressibility. For SZ, they considered a number of
input predictors such as the number of elements in the dataset, the
quantization interval, information about the Huffman tree constructed, and the
number of points that were unpredictable with SZ’s predictor. For ZFP, they
offered a proof of their sampling methods estimated accuracy and empirically
show it to be 99% accurate across many datasets. In the same vain, [16] based
on their prior work, used deep neural networks to estimate the compression
ratio instead of a Gaussian model. However the built neural network may not
generalize to other applications because it could simply over-learn the
training data and the testing data would not prevent the over-learning [17].
[18] designed an automated methodology to switch between SZ and ZFP based on
which compressor is estimated to have a greater compression ratio. Compression
ratios are estimated in a block-based sampling approach using Shannon entropy
[2] of the sampled quantized blocks to investigate SZ’s behavior. Most of
these works have limited generalizability because of their reliance on
algorithmic details of each studied compressor or reduction technique.

Finally, little effort has been directed to explore explicit links of
correlation structures in the data to reduction and compression techniques and
their performance. [19] has explored the decorrelation efficiency of a
specific reduction method on scientific datasets in order to identify trade-
offs for parameterization. [20] lead an evaluation and comparison of several
lossy reduction techniques based on basis-function decompositions adapted to
temporal and spatial data. [21] developed an adaptive hierarchical geospatial
field data representation (Adaptive-HGFDR) based on blocked hierarchical
tensor decomposition to exploit multidimensional correlations on the data.
However, none of these works systematically investigate the explicit link
between correlation structures and compressibility. Our work goes beyond these
approaches to consider a direction that is compressor independent by looking
only at local spatial relationships in the data. However, unlike these papers
we leave to future work the task of estimating the compression ratio from the
observations about spatial relations we study in this paper.

##  IV Methodology

###  IV-A Datasets

In the following, we refer to datasets as a particular field at a particular
time when in an application (e.g. a single continuous variable in memory). In
particular, we work with a single temporal snapshot of the studied data.

The first dataset consists of synthetic 2D Gaussian fields with a controllable
correlation structure following a squared-exponential correlation model. These
2D fields are 1028$\times$1028 grid-points. We consider these fields as
“ideal” as the correlation range is known and varied to create multiple
correlated fields. Gaussian fields $z$ over a grid defined by indexes $x_{i}$
are generated using the following probability distribution $f$

$f(z(x_{1}),\ldots,z(x_{k}))=\frac{\exp\left(-{\frac{1}{2}}({\mathbf{z}}-{\boldsymbol{\mu}})^{\mathrm{T}}{\boldsymbol{\Sigma}}(\mathbf{x})^{-1}({\mathbf{z}}-{\boldsymbol{\mu}})\right)}{\sqrt{(2\pi)^{k}|{\boldsymbol{\Sigma}}(\mathbf{x})|}}$
(2)

with $\mathbf{z}=(z(x_{1}),\ldots,z(x_{k}))\in\mathbb{R}^{k}$, the mean
${\boldsymbol{\mu}=0}\in\mathbb{R}^{k}$ in this study, and a squared-
exponential correlation
${\boldsymbol{\Sigma}}(x_{i},x_{j})=\sigma^{2}\exp(-|x_{i}-x_{j}|^{2}/a^{2})$,
where the variance $\sigma^{2}$ is set to $1$, $a$ is the correlation range
that is known and varied, and $x_{i}$ are spatial grid-points of the 2D field
images.

Figure 2: Original images of 2D Gaussian fields (top) and Miranda dataset
2D-slices (bottom)

In our evaluation, we consider two types of synthetic Gaussian datasets:
single correlation range fields and multiple correlation range fields. We
consider single correlation range Gaussian fields as a proof of concept that
gives us a high degree of control over correlation structures in the data.
However, because application datasets are likely to exhibit complex
correlation structure as multiple correlation ranges, we consider also
multiple range Gaussian fields. For the multi-range correlation, we generate
Gaussian fields with two distinct correlation ranges contributing equally to
the total field. This provides a controlled case with increased complexity.

The final dataset is generated from the Miranda[5] code, designed for
hydrodynamical large turbulence simulations. These data are more complex than
the Gaussian fields due to multiple correlation ranges and complex
dependencies. These original 3D data dimensions are
256$\times$384$\times$384\. The 3D datasets are split into separate 2D slices
based on equally spaced slices along the first dimension. In this paper, we
use slices of the velocityx variable as shown in Figure 2.

###  IV-B Compressors and Software

Each lossy compressor in Table I is run with the following absolute error
bounds: 1E-5, 1E-4, 1E-3, and 1E-2. We chose the software and compressor
versions from the latest available on Spack from a selection of leading error-
bounded lossy compressors. We use the absolute error bound because it is
supported by each of the leading error-bounded lossy compressors.
Additionally, there are formal equivalences between the absolute error bound
mode and other error bounds modes such as the value range relative error bound
mode which are used by compressors such as SZ [3]. Therefore, considering the
compressablity given an absolute error bound is generalizable to other kinds
of bounds. All experiments are run on Clemson’s Palmetto cluster using a node
with two 32 core Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz and 384 GB of RAM.
The OS is Linux CentOS 8 with compiler GCC 8.4.1. Additional software and
packages used in the study are listed in Table I.

Software Version Purpose SZ[3] @2.1.11.1 lossy compressor ZFP[4] @0.5.5 lossy
compressor MGARD[6] @0.1.0 lossy compressor gstat[22] @2.0-7 obtain variogram
range numpy[23] @1.21.1 polyfit function to graph the curves Libpressio[24]
@0.70.0 compress and measure the data

TABLE I: Compressors and software used for the study

###  IV-C Compression Statistics and Statistical Methods

####  IV-C1 Compression Statistics

Compression ratio, the ratio of the uncompressed data size by the compressed
data size, is used to compare the different compressors and their efficiency.
Compression ratio depends on: error bound, compressor used, and correlation
structures within the data. The compression ratio is comparable between
different compressors and error bounds. In the following, this quantity is
computed on the studied datasets from Section IV-A for different compressors
and error bounds, and investigated as function of a measure of several
correlation statistics of data computed through the variogram range described
below.

####  IV-C2 Variogram study

In the following section, we compute the empirical variogram of each 2D data-
slice from the datasets described in Section IV-A and based on the Euclidean
distance between grid-points. The corresponding range $a$ is then estimated
and reported in the following section as estimated variogram range. In
particular, we estimate the variogram ranges on the entire 2D field in order
to assess the overall correlation structure of the fields. However, this is
insufficient to characterize local heterogeneity in datasets, hence we compute
the variogram ranges in windows of a given size that cover the entire 2D field
in a tiled fashion [25]. More specifically, we compute and report the standard
deviation of variogram ranges estimated over the windows covering the entire
field. This statistic provides information on the spatial diversity and spread
of local correlations present in the data.

##  V Experimental Results

###  V-A Compressibility and global correlation

Figure 3: Compression ratios against estimated variogram range for 2D Gaussian
fields with single correlation range (left) and multi-range correlation
(right). Each color is associated with an error bound. Empirical calculations
are depicted with dots and fitted logarithmic regressions are shown with solid
lines. Estimated logarithmic regression coefficients $\alpha$ and $\beta$ are
given in the legend.

In the following figures, the variogram range estimated on entire 2D fields is
referred to as “Estimated global variogram range”, the standard deviation of
locally estimated variogram range on $32\times 32$ windows is denoted as “Std
estimated of local variogram range (H=32)”. Finally, an additional statistic
are considered to illustrate future research directions. It consists in the
standard deviation of locally extracted SVD thresholds on windows of size
$32\times 32$. Thresholds correspond to the number of required singular modes
to recover $99\%$ of the variance of the initial field. This statistic is
referred to as “Std of truncation level of local SVD (H=32)”.

In Figures 3 and 4 , the compression ratios for SZ[3] and ZFP[4] (top panels)
are plotted against the variogram range estimated on entire 2D fields. These
statistics are respectively estimated on the synthetic Gaussian fields (left:
single-range correlation and right: multi-range correlation) and Miranda
datasets. Multiple curves shown on each panel correspond to different error
bounds of the lossy compressor. As the estimated global variogram range
increases indicating stronger dependence between spatial points, dataset
variability decreases yielding smoother, more compressible data. This
increasing relationship between compression ratio and variogram range exhibits
a plateau for highly correlated data (large variogram ranges) suggesting a
limit in compressibility of the data for a given error bound and compressor.
Note that this trend is less visible on the multi-range correlation Gaussian
fields (right column of Figure 3) due to the equal contributions of each
correlation to the total field, preventing any dominant correlation pattern to
prevail and thus to be characterized properly by the variogram range of the
entire field. Compression ratios for MGARD[6] are less sensitive to the
dependencies to correlation ranges present in the data which is likely due to
the global scope of the compressor.

Finally, in order to quantify these relationships and compare them across
different compressors and error bounds, the compression ratios are fitted by
least-squares to logarithmic regressions of the estimated variogram range $a$:
$CR=\alpha+\beta\log(a)+\epsilon$, where $CR$ is the compression ratio, $a$
the estimated variogram range, and $\alpha$ and $\beta$ are estimated
coefficients and reported in the legend box of each panel. The fitted
logarithmic regressions show a good match to the datapoints indicating in most
cases a logarithmic dependence of the compression ratios to the estimated
variogram ranges. Lower compression error bounds exhibit lower variance of
datapoints around the fitted curve and fewer outliers. Regressions fit better
the single-scale correlation Gaussian fields (Figure 3 left column) than
multi-range correlation Gaussian fields and the Miranda data due to their
lower complexity that is captured reasonably by the global variogram range. In
particular, the fitting on datapoints from Miranda data show more dispersion
around the fitted curves but a matching trend.

The variogram range estimated on each entire 2D field represents an average
correlation range observed on the entire field and is not suited to
characterize spatial local heterogeneity nor multiple scales that can be
present in complex scientific datasets, such as in Miranda datasets or multi-
range Gaussian fields. In the following section, we explore local correlations
and their link to compression ratios.

Figure 4: Compression ratios against estimated variogram range for Miranda
velocityx. Empirical calculations are depicted with dots and fitted
logarithmic regressions are shown with solid lines. Each color is associated
with an error bound. Estimated logarithmic regression coefficients $\alpha$
and $\beta$ are given in the legend. Due to the large spread of compression
ratios for SZ (top left panel), the top right panel corresponds to results for
error bounds strictly below 1E-2 in order to ease the reading.

###  V-B Compressibility and local correlation

As the variogram ranges estimated on the entire 2D fields may not provide
enough insights on local heterogeneity, we estimate the variogram ranges on
local windows ($32\times 32$) tiling each entire 2D-field. Additionally,
exploring local heterogeneity is important since many leading error-bounded
lossy compressors exploit some notions of locality in their algorithms as
discussed in Section II-A. In Figure 3, as the global variogram range
increased, the compression ratio fails to have a slope of greater than
absolute value of 1E-2 for mutli-range Gaussian fields (right column). This
illustrates the shortcomings of the global variogram range statistic while
dealing with multi-range correlation datasets. Figure 5 shows the compression
ratios computed on single-range correlation Gaussian fields (left) and multi-
range correlation fields (right) as a function of the local variogram ranges
and the left column of Figure 7 shows the statistics computed on the Miranda
fields. Both figures corroborate and illustrate the benefit of considering
local information for complex datasets to be better explain compression ratios
by local variograms rather than global ones. In particular, we observe that
Figures 3 and 4 depict several close values of compression ratios for close
variogram ranges, whereas this effect is less visible in Figures 5 and 7
indicating a stronger exploratory skill of the local statistic to the
compression ratios. However, results for the single-range correlation Gaussian
fields show a weaker sensitivity of the compression ratios to this local
statistic. This might suggest the need to use several statistics to provide
appropriate explanatory skills. Future work will explore this path.

Figure 5: Compression ratios against standard deviation of the local variogram
range for single range correlation Gaussian fields (right) and for multi-range
correlation Gaussian fields (right). Initial step towards the use of
decomposition techniques to analysis fields with multiscale patterns.
Empirical calculations are depicted with dots and fitted logarithmic
regressions are shown with solid lines. Each color is associated with an error
bound. Estimated logarithmic regression coefficients $\alpha$ and $\beta$ are
given in the legend.

###  V-C Towards multiscale analysis and summary statistics

Figure 6: Compression ratios against standard deviation of the truncation
level of local SVD for single scale Gaussian fields (left) and for multiscale
Gaussian datasets (right). Empirical calculations are depicted with dots and
fitted logarithmic regressions are shown with solid lines. Each color is
associated with an error bound. Estimated logarithmic regression coefficients
$\alpha$ and $\beta$ are given in the legend.

Figure 6 and the right column of Figure 7 illustrate an initial step towards
the use of decomposition techniques to analysis fields with multiscale
patterns. In this section, local SVD are performed on the fields and
summarized via the standard deviation of locally required number of singular
modes to capture $99\%$ of the variance of each local window. This local
statistic, larger values of required singular modes are associated with less
compressible fields so we expect mostly decreasing relationship of compression
ratios to this local statistic. For simplicity, as this section provides
insights to future works, we omit the compressor MGARD as it showed less
sensitivity to the previous statistics. Figure 6 and the right column of
Figure 7 indicate that this local statistic provide a more diverse
representation of the data as seen by the span of unique values over the
x-axis, than the two other statistics based on variogram. This statistic tend
to exhibit several relating trends to compression ratio, highlighting a need
to refine this statistic. Future work will explore transformation of this
statistic or other representations based local SVDs and variograms, as both
methods provided various information to explain compression ratios.

Figure 7: Compression ratios against local statistics, standard deviation of
the local variogram range (left) and standard deviation of the truncation
level of local SVD (left) for Miranda velocityx 2D field. Empirical
calculations are depicted with dots and fitted logarithmic regressions are
shown with solid lines. Each color is associated with an error bound.
Estimated logarithmic regression coefficients $\alpha$ and $\beta$ are given
in the legend. Due to the large spread of compression ratios for SZ, the
central panels correspond to results for error bounds strictly below 1E-2 in
order to ease the reading.

##  VI Conclusions and Future Work

Establishing the limit for lossy compression would allow for the maximum
efficiency for compressing and storing large scientific datasets. Our work
represent a first step toward this goal. We have demonstrated that estimated
global and local variogram ranges can explain compression ratio in a
logarithmic fashion for some compressors and given error bounds. This
hypothesis was illustrated on synthetic Gaussian fields providing a proof of
concept and corroborated by results on a user scientific dataset from Miranda.
With the studied datasets, SZ and ZFP seem to utilize the global and local
spatial correlation ranges in a logarithmic fashion with various coefficients
showing the strength of the dependence. MGARD seems less sensitive to these
trends in its compression ratios.

Since heterogeneous (non-stationary) and multiscale correlations in the data
may be mis-represented by the global spatial variogram, other statistics,
including local variograms and local SVD, have been studied to address these
issues. In continuation of this research, there are a few goals for the
future: i) explore more complex dependent variables (local correlation
combined with multiscale statistics based on decomposition) as candidate
predictors, ii) create more complex synthetic multiscale 2D Gaussian fields
and iii) create a model of compression ratio based on correlation metrics and
error bound. Future work will investigate the effects of correlation
structures on quality metrics of reconstructed data such as PSNR along with a
design of the statistics to a 3D context.

Another aspect for future work is how to quickly assess this metric. The
currently implementation relies on the singular value decomposition which is
slow relative to modern compressors. We plan to leverage a sampling approach
similar to prior work[18, 15]. We are hopeful that increasing levels of
sampling by block can provide an increasingly accurate proxy for our metric.

## Acknowledgments

Clemson University is acknowledged for generous allotment of compute time on
the Palmetto cluster. This material is based upon work supported by the
National Science Foundation under Grant No. SHF-1910197.

## References

  * [1] F. Cappello, S. Di, S. Li, X. Liang, A. M. Gok, D. Tao, C. H. Yoon, X.-C. Wu, Y. Alexeev, and F. T. Chong, “Use cases of lossy compression for floating-point data in scientific data sets,” _The International Journal of High Performance Computing Applications_ , vol. 33, no. 6, pp. 1201–1220, Jul. 2019. [Online]. Available: https://doi.org/10.1177/1094342019853336
  * [2] C. E. Shannon, “A mathematical theory of communication,” _Bell System Technical Journal_ , vol. 27, no. 3, pp. 379–423, Jul. 1948. [Online]. Available: https://doi.org/10.1002/j.1538-7305.1948.tb01338.x
  * [3] X. Liang, S. Di, D. Tao, S. Li, S. Li, H. Guo, Z. Chen, and F. Cappello, “Error-controlled lossy compression optimized for high compression ratios of scientific datasets,” in _2018 IEEE International Conference on Big Data (Big Data)_. IEEE, Dec. 2018. [Online]. Available: https://doi.org/10.1109/bigdata.2018.8622520
  * [4] P. Lindstrom and M. Isenburg, “Fast and efficient compression of floating-point data,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 12, no. 5, pp. 1245–1250, Sep. 2006. [Online]. Available: https://doi.org/10.1109/tvcg.2006.143
  * [5] F. Capello, K. Zhao, S. Di, D. Tao, J. Bessac, and Z. Chen, Jun 2018. [Online]. Available: https://sdrbench.github.io/
  * [6] M. Ainsworth, O. Tugluk, B. Whitney, and S. Klasky, “Multilevel techniques for compression and reduction of scientific data—the multivariate case,” _SIAM Journal on Scientific Computing_ , vol. 41, no. 2, pp. A1278–A1303, Jan. 2019. [Online]. Available: https://doi.org/10.1137/18m1166651
  * [7] G. Matheron, “Principles of geostatistics,” _Economic geology_ , vol. 58, no. 8, pp. 1246–1266, 1963.
  * [8] A. E. Gelfand, P. Diggle, P. Guttorp, and M. Fuentes, _Handbook of spatial statistics_. CRC press, 2010.
  * [9] P. Abry, D. Veitch, and P. Flandrin, “Long-range dependence: Revisiting aggregation with wavelets,” _Journal of Time Series Analysis_ , vol. 19, no. 3, pp. 253–266, 1998.
  * [10] J. Cabrieto, F. Tuerlinckx, P. Kuppens, M. Grassmann, and E. Ceulemans, “Detecting correlation changes in multivariate time series: A comparison of four non-parametric change point detection methods,” _Behavior research methods_ , vol. 49, no. 3, pp. 988–1005, 2017.
  * [11] J. Bessac, P. Ailliot, J. Cattiaux, and V. Monbet, “Comparison of hidden and observed regime-switching autoregressive models for (u,v)-components of wind fields in the Northeast Atlantic,” _Advances in Statistical Climatology, Meteorology and Oceanography_ , vol. 2, no. 1, pp. 1–16, 2016.
  * [12] P. S. Addison, _The illustrated wavelet transform handbook: introductory theory and applications in science, engineering, medicine and finance_. CRC press, 2017.
  * [13] A. Hannachi, I. T. Jolliffe, and D. B. Stephenson, “Empirical orthogonal functions and related techniques in atmospheric science: A review,” _International Journal of Climatology: A Journal of the Royal Meteorological Society_ , vol. 27, no. 9, pp. 1119–1152, 2007.
  * [14] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is $4\backslash\sqrt{3}$,” pp. 5040 – 5053, 2014.
  * [15] T. Lu, Q. Liu, X. He, H. Luo, E. Suchyta, J. Choi, N. Podhorszki, S. Klasky, M. Wolf, T. Liu _et al._ , “Understanding and modeling lossy compression schemes on hpc scientific data,” in _2018 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2018, pp. 348–357.
  * [16] Z. Qin, J. Wang, Q. Liu, J. Chen, D. Pugmire, N. Podhorszki, and S. Klasky, “Estimating lossy compressibility of scientific data using deep neural networks,” _IEEE Letters of the Computer Society_ , vol. 3, no. 1, pp. 5–8, 2020.
  * [17] J. Han, J. Pei, and M. Kamber, _Data Mining: Concepts and Techniques_. Elsevier, 2011.
  * [18] D. Tao, S. Di, X. Liang, Z. Chen, and F. Cappello, “Optimizing Lossy Compression Rate-Distortion from Automatic Online Selection between SZ and ZFP,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 30, no. 8, pp. 1857–1871, Aug. 2019.
  * [19] P. Lindstrom, “Fixed-rate compressed floating-point arrays,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 20, no. 12, pp. 2674–2683, 2014.
  * [20] A. Moon, J. Kim, J. Zhang, and S. W. Son, “Lossy compression on iot big data by exploiting spatiotemporal correlation,” in _2017 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 2017, pp. 1–7.
  * [21] Z. Yu, D. Li, Z. Zhang, W. Luo, Y. Liu, Z. Wang, and L. Yuan, “Lossy compression of earth system model data based on a hierarchical tensor with adaptive-hgfdr (v1. 0),” _Geoscientific Model Development_ , vol. 14, no. 2, pp. 875–887, 2021.
  * [22] E. J. Pebesma, “Multivariable geostatistics in s: the gstat package,” _Computers & Geosciences_, vol. 30, no. 7, pp. 683–691, Aug. 2004. [Online]. Available: https://doi.org/10.1016/j.cageo.2004.03.012
  * [23] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * [24] R. Underwood, S. Di, J. C. Calhoun, and F. Cappello, “Fraz: A generic high-fidelity fixed-ratio lossy compression framework for scientific floating-point data,” _2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2020.
  * [25] J. Bessac, A. H. Monahan, H. M. Christensen, and N. Weitzel, “Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes,” _Monthly Weather Review_ , vol. 147, no. 5, pp. 1447–1469, 2019\. [Online]. Available: https://doi.org/10.1175/MWR-D-18-0384.1


[figure]style=plain,subcapbesideposition=top

# Unveiling NPT bound problem: From Distillability Sets to Inequalities and
Multivariable Insights

Si-Yuan Qi  CAS Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei, Anhui 230026, P. R. China  CAS Center For
Excellence in Quantum Information and Quantum Physics, University of Science
and Technology of China, Hefei, Anhui 230026, P. R. China  Geni Gupur
<EMAIL_ADDRESS>College of Mathematics and Systems Science, Xinjiang
University, Urumqi, Xinjiang 830017, People’s Republic of China  Yu-Chun Wu
<EMAIL_ADDRESS>CAS Key Laboratory of Quantum Information, University of
Science and Technology of China, Hefei, Anhui 230026, P. R. China  CAS Center
For Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei, Anhui 230026, P. R. China  Institute
of Artificial Intelligence, Hefei Comprehensive National Science Center,
Hefei, Anhui 230088, People’s Republic of China  Guo-Ping Guo  CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei, Anhui 230026, P. R. China  CAS Center For Excellence in Quantum
Information and Quantum Physics, University of Science and Technology of
China, Hefei, Anhui 230026, P. R. China  Institute of Artificial Intelligence,
Hefei Comprehensive National Science Center, Hefei, Anhui 230088, People’s
Republic of China  Origin Quantum Computing Hefei, Anhui 230026, P. R. China

(August 31, 2024)

###### Abstract

Equivalence between Positive Partial Transpose (PPT) entanglement and bound
entanglement is a long-standing open problem in quantum information theory. So
far limited progress has been made, even on the seemingly simple case of
Werner states bound entanglement. The primary challenge is to give a concise
mathematical representation of undistillability. To this end, we propose a
decomposition of the N-undistillability verification into $log(N)$ repeated
steps of 1-undistillability verification. For Werner state N-undistillability
verification, a bound for N-undistillability is given, which is independent of
the dimensionality of Werner states. Equivalent forms of inequalities for both
rank one and two matrices are presented, before transforming the two-
undistillability case into a matrix analysis problem. A new perspective is
also attempted by seeing it as a non-convex multi-variable function, proving
its critical points and conjecturing Hessian positivity, which would make them
local minimums.

††preprint: APS/123-QED

##  I Introduction

Entanglement [1] is what distinguished the quantum world from classical
worlds. Unfortunately, our understanding of it remains quite limited to this
day. Entanglement exhibits varying degrees of strength, the highest of which
falls into a subclass called “maximal entanglement”. Bell state is one such
example. No consensus has yet been reached about what qualifies as a “weak
entanglement”. Quantum states with small Schmidt numbers are generally seen as
weakly entangled, since quantum states with Schmidt number 1 are separable.
Undistillable entanglement is another kind of widely accepted weak
entanglement. Relationships between these two kinds of weak entanglements are
discussed in [2].

Entanglement distillation [3] is the process of obtaining pure Bell states
only by Local Operations and Classical Communication (LOCC) from multiple
copies of an entangled bipartite state. If N copies of a given entangled
bipartite state can be transformed into pure Bell states through LOCC, then
the entangled bipartite state is referred to as N-distillable. If not, then
it’s N-undistillable. If a particular entangled bipartite state is
N-undistillable for arbitrary N, it is considered undistillable and is also
referred to as a bound entangled state.

Another property concerned here is “Positive Partial Transpose” (PPT). PPT
refers to a bipartite quantum state whose partial transpose is positive
semidefinite. It is established in [4] that Positive Partial Transpose (PPT)
entanglement implies bound entanglement, and therefore weakly entangled in
terms of distillability. However, the question of whether the converse is true
remains an open problem [5]. This problem is usually referred to as the NPT
bound problem, or the distillability problem. The popular conjecture regarding
this problem is that the converse is not true, and that NPT bound entanglement
exists. A number of special circumstances have been worked out in [6, 7, 8, 9,
10]. Indirect approaches to this problem include expansion of distillation
operations to k-extendible operations [11], catalysis-assisted distillation
[12], and dually non-entangling and PPT-preserving channels [13], linking the
problem with squashed entanglement [14] and linear preserver [15], studying
the transformation of bound entangled states under certain dynamic process
[16], or considering the problem in a broader scenario such as hyperquantum
states [17].

In three separate attempts to directly tackle the distillability problem [18],
[19] and [20], they each had a subset of quantum states singled out, and
proved that for any finite number of copies N, there exist undistillable
states in this subset. However, as N approaches infinity, the subset shrinks
to emptiness. It has been shown in [21] that N-undistillability does not imply
(N+1)-undistillability, by demonstrating a set of distillable-only-by-
arbitrarily-large-number-of-copies state, thus ruling out the easy route of
concentrating on few-copies distillation.

For the Werner states [22], a family of single-parameter states, the existence
of NPT bound states implies existence of NPT bound Werner states [23].
Therefore, it is sufficient to only study Werner states. Due to the inherent
complexity of addressing N-copy Werner state distillability, over the last
decade, direct attacks on the problem have been primarily fixating on the
seemingly straightforward problem of 2-copy Werner state distillability, but
so far only limited progress has been made. In an attempt to attack the 2-copy
Werner state distillability problem, the special case of $4\times 4$ bipartite
states is considered in [24], in which case the distillability condition can
be reduced to a convenient form of “half-property”, and subsequently
reformulated the distillability problem into a matrix analysis problem. For a
special case where the matrices are normal, a proof was given.

In [25] the problem is reformulated and established that for a $d\times d$
bipartite Werner state, the 2-copy Werner state undistillability problem is
equivalent to a certain $2d^{2}\times 2d^{2}$ matrix being positive
semidefinite. Positive definiteness of the block matrices in the upper-left
and lower-right corners has been proved.

In this work, we introduce a method that decomposes bound entanglement
verification into repeated steps of 1-undistillability verification, before
giving a new bound for N-undistillability that is independent of quantum state
dimensions $d$, and therefore different from that given in [18], [19] and
[20]. However, similar to the problem encountered in [18], [19] and [20], our
set of N-undistillable states also shrinks to emptiness for $N\to\infty$.
Also, an equivalent matrix analysis inequality form is presented, which is
applicable to all finite dimensionalities, as opposed to the case in [24] when
only the $d=4$ case is considered. Another equivalent version of inequalities
regarding only rank-one matrices is also included, which takes on a form
reminiscent of Cauchy-Schwartz inequality. Finally, we take on a brand-new
perspective of considering it as a multi-variable function problem, finding
critical points, proving non-convexity and conjecturing Hessian positivity.

We will be using $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ to denote Hilbert
spaces pertaining to A and B subsystems, and $X_{i}$ to denote matrices of
rank less than or equal to $i$. Throughout the paper we will be considering
multi-partite quantum states with identical dimensions, i.e. quantum states on
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, where
$dim(\mathcal{H}_{A})=dim(\mathcal{H}_{B})=d$. This is justified by noticing
that any quantum state of $d_{1}\times d_{2}$ can be equivalently transformed
into a quantum state of identical dimensions $max\\{d_{1},d_{2}\\}\times
max\\{d_{1},d_{2}\\}$, with all unnecessary elements set to zero.

The structure of this paper is as follows: In Section II, we introduce a
matrix transformation that takes N-distillability problem into a higher
dimensional 1-distillability problem. In Section III, the process of
$N$-undistillability verification is transformed into $log(N)$ repeated steps
of 1-undistillability verification. In Section IV, upper bounds on
$N$-undistillability are given, yielding a result similar to that of [18],
[19] and [20], effectively finding shrinking families of N-undistillable
quantum states. In Section V, equivalent partial trace inequalities are
presented before turning them into a matrix analysis problem. The rank-two
matrix analysis problem is also reduced to a rank-one matrix analysis problem,
resulting in a Cauchy-Schwartz inequality look-alike. In Section VI, we take
on a new perspective by treating the problem as a multi-variable function,
proving critical points, non-convexity and conjecturing positive-semidefinite
Hessians.

##  II Statement of the distillability problem

The open problem of equivalence between PPT entanglement and bound
entanglement can be resolved by finding bound entangled NPT states. Since
Peres [26] has already shown that NPT is a sufficient condition for
entanglement, attention can be limited to the state being both “bound”, which
is, N-undistillable for arbitrarily large N, and NPT, which is, having a
negative partial transpose. Compared to ascertaining positive semidefiniteness
of a certain matrix, the major difficulty usually lies in finding a concise
mathematical interpretation for undistillability, and asserting it for two
types of arbitrariness: arbitrary LOCC operation and arbitrarily finitely many
copies of the entangled state. This will be the main focus of this section,
and the analysis of mathematical structure will eventually lead to a
transformation that reduces N-undistillability verification to
1-undistillability verification of a transformed state.

Let’s recall the concept of Schmidt decomposition, Schmidt coefficients and
Schmidt rank [27]. Schmidt decomposition of a pure state $|\psi\rangle$ in
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ is:

$|\psi\rangle=\sum_{i=1}^{d}\alpha_{i}|u_{i}v_{i}\rangle,$ (1)

where $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ are both d-dimensional Hilbert
spaces, with $\\{|u_{1}\rangle,\cdots,|u_{d}\rangle\\}$ and
$\\{|v_{1}\rangle,\cdots,|v_{d}\rangle\\}$ as their respective orthonormal
bases. Schmidt coefficients $\alpha_{i}$ are real, nonnegative and unique up
to re-ordering. The number of non-zero Schmidt coefficients is called Schmidt
rank of pure state $|\psi\rangle$. $SR(|\psi\rangle)$ denotes the Schmidt rank
of $|\psi\rangle$.

A bipartite quantum state is called N-distillable if pure Bell states can be
obtained using LOCC, from N-copies of the state. It has been proven in [28]
that all entangled $2\times 2$ bipartite quantum states are distillable,
therefore, all we need to do is transform the original quantum state into a
$2\times 2$ entangled state by LOCC. This is demonstrated in the following
result: A bipartite quantum state $\rho$ on
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ is N-undistillable iff for any
$d^{N}\times d^{N}$ bipartite pure state $|\psi^{SR\leq 2}\rangle$ on
$\mathcal{H}_{A}^{\otimes N}\otimes\mathcal{H}_{B}^{\otimes N}$ with Schmidt
rank $\leq$ 2, the following holds:

$\langle\psi^{SR\leq 2}|(\rho^{T_{A}})^{\otimes N}|\psi^{SR\leq 2}\rangle\geq
0,$ (2)

where $\rho^{T_{A}}=(T\otimes I)\rho$ is the partial transpose of $\rho$, with
$T$ being the transpose operator. The above expression overlooks the pairing
between different subsystems and so an additional operator $M_{N}$ from
$(\mathcal{H}_{A}\otimes\mathcal{H}_{B})^{\otimes N}$ to
$\mathcal{H}_{A}^{\otimes N}\otimes\mathcal{H}_{B}^{\otimes N}$ is introduced
to make it meaningful:

######  Definition II.1.

Operator $M_{N}:(\mathcal{H}_{A}\otimes\mathcal{H}_{B})^{\otimes
N}\to\mathcal{H}_{A}^{\otimes N}\otimes\mathcal{H}_{B}^{\otimes N}$ is defined
as:

$M_{N}=\sum|i_{A_{1}}\dots i_{A_{N}}i_{B_{1}}\dots i_{B_{N}}\rangle\langle
i_{A_{1}}i_{B_{1}}\dots i_{A_{N}}i_{B_{N}}|,$ (3)

where $\mathcal{H}_{A}^{\otimes
N}=\mathcal{H}_{A_{1}}\otimes\mathcal{H}_{A_{2}}\otimes\cdots\otimes\mathcal{H}_{A_{N}}$
, $\mathcal{H}_{B}^{\otimes
N}=\mathcal{H}_{B_{1}}\otimes\mathcal{H}_{B_{2}}\otimes\cdots\otimes\mathcal{H}_{B_{N}}$,
and $i_{A_{j}},i_{B_{j}}$ stands for the index corresponding to system
$A_{j},B_{j}$.

Operator $M_{N}$ then takes the $A_{1}B_{1}A_{2}B_{2}\cdots A_{N}B_{N}$
structure of $(\rho^{T_{A}})^{\otimes N}$ and merges their A-subsystems and
B-subsystems, resulting in a structure of $A_{1}A_{2}\cdots
A_{N}B_{1}B_{2}\cdots B_{N}$ that suits the pure quantum state, thus
explicitly demonstrating the process of system pairing which is essential to
any further calculations derived from this representation of undistillability.

For better representation of undistillability and further use in following
sections, we use a certain state-operator isomorphism introduced in [24]. We
notice that it is in essence a Choi-Jamiokowski isomorphism followed by a
transposition on the resulting matrix. For convenience and consistency with
previous works, we stick to this “transposed Choi-Jamiokowski isomorphism”
when presenting our findings.

######  Definition II.2 (State-operator isomorphism [24]).

A state-operator isomorphism $\Psi$ is one that takes a pure bipartite state
$|\psi\rangle$ of:

$|\psi\rangle=\sum_{ij}\alpha_{ij}|ij\rangle,$ (4)

to a matrix form of:

$\Psi(|\psi\rangle)=\begin{pmatrix}\alpha_{00}&\alpha_{01}&\cdots&\alpha_{0,d-1}\\\
\alpha_{10}&\alpha_{11}&\cdots&\alpha_{1,d-1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\alpha_{d-1,0}&\alpha_{d-1,1}&\cdots&\alpha_{d-1,d-1}\\\ \end{pmatrix}.$ (5)

Note that Schmidt decomposition of a pure quantum state $|\psi\rangle$ is
related to the singular value decomposition of $\Psi(|\psi\rangle)$ through
state-operator isomorphism. In fact, Schmidt coefficients of $|\psi\rangle$ is
the same as the singular values of $\Psi(|\psi\rangle)$, and Schmidt rank of
$|\psi\rangle$ is equal to the rank of $\Psi(|\psi\rangle)$.

The following is a more concise mathematical representation of
N-undistillability: A bipartite quantum state $\rho$ is N-undistillable iff:

$[\Psi^{-1}(X_{2})]^{\dagger}M_{N}(\rho^{T_{A}})^{\otimes
N}M_{N}^{\dagger}\Psi^{-1}(X_{2})\geq 0,$ (6)

where $X_{2}$ is any $d^{N}\times d^{N}$ matrix of rank $\leq 2$. Note that
$\Psi^{-1}(X_{2})$ is equivalent to $|\psi^{SR\leq 2}\rangle$, a pure
bipartite quantum state with Schmidt rank no larger than 2. For the case of
$N=1$, 1-undistillability is:

$[\Psi^{-1}(X_{2})]^{\dagger}\rho^{T_{A}}\Psi^{-1}(X_{2})=\langle\psi^{SR\leq
2}|\rho^{T_{A}}|\psi^{SR\leq 2}\rangle\geq 0.$ (7)

Effectively, the problem of $\rho$ being N-copy undistillable is equivalent to
the problem of $M_{N}(\rho^{\otimes N})M_{N}^{\dagger}$ being 1-undistillable,
but on a Hilbert space with higher dimensions
$\mathcal{H^{\prime}}_{A}\otimes\mathcal{H^{\prime}}_{B}$,
$dim(\mathcal{H^{\prime}}_{A})=dim(\mathcal{H^{\prime}}_{B})=d^{N}$.

Since discussions of 1-copy distillability problem rarely involves
specification of finite dimensional Hilbert space dimensionality, and most of
the theorems and methods are applicable to arbitrary dimensions, this
transformation poses an advantage. In the next section, a specific symmetry of
$M_{2^{N}}$ transformations is utilized for decomposing the increasingly
complicated structure of $M_{2^{N}}$ into identical steps repeated for $N$
times. In fact, verification of $2^{N}$-undistillability is enough to assert
bound entanglement.

##  III N-undistillability problem as repeated steps

From the observation that (N+1)-undistillability implies N-undistillability,
we can get the following result which relaxes the N-undistillability
requirement for all N to only an infinite sequence of N diverging to infinity:

######  Proposition III.1.

$\\{g_{i}\\}$ is a sequence of positive real numbers diverging to infinity. A
bipartite entangled quantum state is $g_{i}$-undistillable for arbitrary i,
iff the state is bound entangled.

In the 2-copy case, $M_{2}$ transformation is simply:

$E=M_{2}=\sum|i_{A_{1}}i_{A_{2}}i_{B_{1}}i_{B_{2}}\rangle\langle
i_{A_{1}}i_{B_{1}}i_{A_{2}}i_{B_{2}}|,$ (8)

which is merely an exchange of the second and third parts of $B_{1}$ and
$A_{2}$. From now on we will always refer to the 2-copy case of $M_{2}$ as
$E$, emphasizing that it is only a simple exchange as opposed to complicated
$M_{N}$ of larger copy numbers.

According to Eq.6, the quantum state $\rho$ is 2-undistillable iff
$E(\rho\otimes\rho)E^{\dagger}$ is 1-undistillable.

In the 4-copy case, $M_{4}$ can actually be decomposed into two steps:

$\displaystyle A_{1}B_{1}A_{2}B_{2}A_{3}B_{3}A_{4}B_{4}$ (9)
$\displaystyle\to$ $\displaystyle A_{1}A_{2}B_{1}B_{2}A_{3}A_{4}B_{3}B_{4}$
$\displaystyle\to$ $\displaystyle A_{1}A_{2}A_{3}A_{4}B_{1}B_{2}B_{3}B_{4}.$

Step 1: do the 2-copy exchange on subsystems 1,2 ($A_{1}B_{1}A_{2}B_{2}\to
A_{1}A_{2}B_{1}B_{2}$) and subsystems 3,4 ($A_{3}B_{3}A_{4}B_{4}\to
A_{3}A_{4}B_{3}B_{4}$). This step is no different than in the 2-copy case,
except that it is done on two systems simultaneously. Group every neighbouring
2 subsystems(in terms of $A,B$), so that the eight parts are now seen as four
parts ($A_{1}A_{2},B_{1}B_{2},A_{3}A_{4},B_{3}B_{4}$ to $ABAB$). In order to
separate this newly obtained four-part state from 2-copy Werner states, we
denote this state $\rho(e)\otimes\rho(e)$, where
$\rho(e)=E(\rho\otimes\rho)E^{\dagger}$.

Step 2: Do the 2-copy exchange $E$ again, obtaining
$\rho(e^{2})=E(\rho(e)\otimes\rho(e))E^{\dagger}$, where $E$ is taken with
regard to the “merged” parts. This step takes once-transformed Werner states
into twice-transformed Werner states.

The above steps of repeating the $E$ transformation twice, is equivalent to
doing the $M_{4}$ translations directly. The original state $\rho$ is 4-copy
undistillable iff $\rho(e^{2})$ is one-copy undistillable.

Similarly, in the $2^{N}$-copy case, E is decomposed into N steps, with each
step taking state $\rho(e^{i-1})\otimes\rho(e^{i-1})$, producing a state
$\rho(e^{i})=E(\rho(e^{i-1})\otimes\rho(e^{i-1}))E^{\dagger}$, verifying its
positivity on Schmidt rank $\leq 2$ states and effectively ascertaining
$2^{i}$-copy undistillability, before passing two copies of this state on, as
a starting state for the next step. The following theorem demonstrates the
process.

######  Theorem III.2.

Define $\rho(e^{0})=\rho$ to be the initial state, and

$\rho(e^{k})=E(\rho(e^{k-1})\otimes\rho(e^{k-1}))E^{\dagger},$ (10)

then $\rho(e^{k})$ being 1-undistillable is equivalent to $\rho$ being
$2^{k}$-undistillable.

In other words, we are essentially searching for an infinite sequence of
points connected by the operation $E$, with the first point being
$\rho=\rho(e^{0})$, second point being
$\rho(e^{1})=E(\rho(e^{0})\otimes\rho(e^{0}))E^{\dagger}$, and the $(i+2)$-th
point being $\rho(e^{i+1})=E(\rho(e^{i})\otimes\rho(e^{i}))E^{\dagger}$. If
points of this infinite sequence always fall within the set of
1-undistillability, then we can safely say that the starting point $\rho$ is
N-undistillable for arbitrary N, or that it is a bound state.

##  IV Bounds for Werner state N-undistillability

###  IV.1 Werner state 1-undistillability and NPT condition

It has been shown in [23] that it is sufficient to only consider Werner states
for NPT bound problem, that if an NPT bound state exists, there must be an NPT
bound state within the family of Werner states. Necessary and sufficient
conditions of Werner state 1-undistillability and partial transpose negativity
are obtained in this section.

We first introduce a lemma that would assist further attempts concerning inner
product with Schmidt rank $\leq 2$ pure states:

######  Lemma IV.1.

Let $|\psi\rangle$ be an arbitrary pure quantum state and $s_{j}$ are the
Schmidt coefficients of $|\psi\rangle$ arranged in descending order, then

$\;\sum_{j=1}^{k}s_{j}^{2}=max_{|\phi^{SR\leq
k}\rangle}|\langle\psi|\phi^{SR\leq k}\rangle|^{2},$ (11)

where $|\phi^{SR\leq k}\rangle$ is an arbitrary quantum pure state with
Schimidt rank being less than k, namely, $SR(|\phi^{SR\leq k}\rangle)\leq k$.

Proof of this lemma is provided in Appendix A.

A Werner state is as follows:

$\rho_{w}=\frac{1}{d^{2}+\beta d}(I+\beta F),$ (12)

where $F=\sum_{ij}|ij\rangle\langle ji|$, $I$ is the unnormalized identity
matrix, d denote the dimensionality of both parts of the system, and $\beta$
is a parameter characterizing the “portion” of swap operator $F$, ranging in
$-1\leq\beta\leq 1$.

Its partial transpose is:

$\rho_{w}^{T_{A}}=\frac{1}{d^{2}+\beta d}(I+\beta G),$ (13)

where $G=\sum_{ij}|ii\rangle\langle jj|$, which is in fact the unnormalized
density matrix of a pure maximally entangle state $|\Phi\rangle$, namely,

$G=d|\Phi\rangle\langle\Phi|,\;|\Phi\rangle=\sum_{i}\frac{1}{\sqrt{d}}|ii\rangle.$
(14)

When $\beta>0$, $\rho_{w}^{T_{A}}$ is always positive semidefinite, hence we
consider $\beta<0$ case only. Combining Eq.6 and Eq.11, we obtain

$\displaystyle\langle\psi^{SR\leq 2}|\rho_{w}^{T_{A}}|\psi^{SR\leq 2}\rangle$
(15) $\displaystyle=\frac{1}{d^{2}+\beta d}(1+\beta d|\langle\Phi|\psi^{SR\leq
2}\rangle|^{2})\stackrel{{\scriptstyle
L.\ref{fact}}}{{\geq}}\frac{1+2\beta}{d^{2}+\beta d}.$

A Werner state $\rho_{w}$ is 1-undistillable iff $\langle\psi^{SR\leq
2}|\rho_{w}^{T_{A}}|\psi^{SR\leq 2}\rangle\geq 0$ for any $|\psi^{SR\leq
2}\rangle$ pure bipartite quantum state of Schmidt rank no larger than 2.
Therefore, Werner state is 1-undistillable iff $\beta\geq-\frac{1}{2}$.

Furthermore, for arbitrary quantum state $|\psi\rangle$, we have

$\langle\psi|\rho_{w}^{T_{A}}|\psi\rangle=\frac{1}{d^{2}+\beta d}(1+\beta
d|\langle\Phi|\psi\rangle|^{2})\geq\frac{1+\beta d}{d^{2}+\beta d},$ (16)

with the equality achieved when $|\psi\rangle=|\Phi\rangle$. Therefore, Werner
state is NPT iff $\beta<-\frac{1}{d}$. Any NPT 1-undistillable Werner state
can therefore only be found within the parameter range of
$-\frac{1}{2}\leq\beta<-\frac{1}{d}$, provided $d>2$.

###  IV.2 Sufficient conditions for Werner state N-undistillability

We derive new bounds for Werner state N-undistillability in this section. As a
result, for any finite N, a set of N-undistillable states can be found, but
similar to the circumstances encountered in [18], [19] and [20], the set of
N-undistillable states shrinks to emptiness as N approaches infinity.

In the previous section it has been established that Werner state
N-undistillability is equivalent to:

$\langle\psi^{SR\leq 2}|M_{N}(\rho_{w}^{T_{A}})^{\otimes
N}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle\geq 0,$ (17)

where $|\psi^{SR\leq 2}\rangle$ is an arbitrary quantum pure state with
Schmidt rank being no larger than 2.

Take the 2-copy case as an example, recalling operation $E$ defined in Eq.8:

$\displaystyle
E(\rho_{w}^{T_{A}}\otimes\rho_{w}^{T_{A}})E^{\dagger}=(\frac{1}{d^{2}+\beta
d})^{2}(E(I\otimes I)E^{\dagger}$ (18) $\displaystyle+\beta[E(I\otimes
G)E^{\dagger}+E(G\otimes I)E^{\dagger}]+\beta^{2}E(G\otimes G)E^{\dagger}).$

The normalization factor of $(\frac{1}{d^{2}+\beta d})^{2}$ can be ignored
since it doesn’t affect positivity. The first and the last term remains
unchanged in the sense that:

$\displaystyle E(I\otimes I)E^{\dagger}=E(\sum_{ijkl}|ijkl\rangle\langle
ijkl|)E^{\dagger}=\sum_{ijkl}|ikjl\rangle\langle ikjl|$ (19)
$\displaystyle=I(e),$

with $I(e)$ still being the identity matrix, only that it exists on a higher
dimensionality of $d^{4}\times d^{4}$. Similarly,

$\displaystyle E(G\otimes G)E^{\dagger}=E(\sum_{ijkl}|iikk\rangle\langle
jjll|)E^{\dagger}=\sum_{ijkl}|ikik\rangle\langle jljl|$ (20)
$\displaystyle=G(e),$

with $G(e)=|\Phi(e)\rangle\langle\Phi(e)|$ still being the pure state density
matrix of unnormalized $|\Phi(e)\rangle=\sum_{ij}|ijij\rangle$, only that it
exists on a higher dimensionality of $d^{4}\times d^{4}$.

Regarding the positivity of the middle part $E(I\otimes
G)E^{\dagger}+E(G\otimes I)E^{\dagger}$ on $|\psi^{SR\leq 2}\rangle$,

$\displaystyle E(I\otimes G)E^{\dagger}=E(\sum_{ijkl}|ijkk\rangle\langle
ijll|)E^{\dagger}=\sum_{ijkl}|ikjk\rangle\langle iljl|$ (21)
$\displaystyle=d(|\psi_{00}^{s_{2}=s_{4}}\rangle\langle\psi_{00}^{s_{2}=s_{4}}|+\dots+|\psi_{d-1,d-1}^{s_{2}=s_{4}}\rangle\langle\psi_{d-1,d-1}^{s_{2}=s_{4}}|)$
$\displaystyle=d\sum_{ij}|\psi_{ij}^{s_{2}=s_{4}}\rangle\langle\psi_{ij}^{s_{2}=s_{4}}|,$

where$|\psi^{s_{2}=s_{4}}_{ij}\rangle=\frac{1}{\sqrt{d}}\sum_{k=0}^{d-1}|ikjk\rangle$
is the superposition of all multipartite states with identical states on the
second and fourth subsystems, while their first and third subsystems are of
states i and j respectively. $|\psi^{s_{2}=s_{4}}_{ij}\rangle$ is a normalized
pure state of Schmidt rank d, and all Schmidt coefficients are
$\frac{1}{\sqrt{d}}$. Do a decomposition of $|\psi^{SR\leq 2}\rangle$ in the
following way:

$|\psi^{SR\leq 2}\rangle=\sum_{ij}p_{ij}|\psi^{SR\leq 2}_{ij}\rangle,$ (22)

where $|\psi^{SR\leq 2}_{ij}\rangle$ represents the normalized extraction of
all $|i-j-\rangle$ terms from the state of $|\psi^{SR\leq 2}\rangle$. More
specifically, for state $|\psi^{SR\leq 2}\rangle$ of:

$|\psi^{SR\leq 2}\rangle=\sum_{abcd}p_{abcd}|abcd\rangle,$ (23)

the corresponding $|\psi^{SR\leq 2}_{ij}\rangle$ is:

$|\psi^{SR\leq 2}_{ij}\rangle=\frac{1}{N_{p}}\sum_{ab}p_{iajb}|iajb\rangle,$
(24)

where $N_{p}=\sqrt{\sum_{ab}|p_{iajb}|^{2}}$.

It can be proven that the normalized extraction of all $|ixjy\rangle$ terms
must also have Schmidt rank $\leq 2$. Consider state-operator isomorphism in
Def.II.2, the Schmidt rank of $|\psi\rangle$ is equivalent to the rank of the
matrix $\Psi(|\psi\rangle)$. The normalized extraction of all $|i-j-\rangle$
terms composes a new state of $|\psi^{SR\leq 2}_{ij}\rangle$, therefore, its
corresponding matrix is a submatrix of $\Psi(|\psi^{SR\leq 2}\rangle)$,
composed of rows $i\times d+x$ and columns $j\times d+y$, with x,y taking
values from $\\{0,\dots,d-1\\}$. Since the rank of a submatrix is never larger
than the rank of the whole matrix, $|\psi^{SR\leq 2}_{ij}\rangle$ must have
Schmidt rank less than or equal to that of $|\psi^{SR\leq 2}\rangle$.

Using Lem.IV.1 again, it is then clear that:

$|\langle\psi_{ij}^{SR\leq
2}|\psi_{pq}^{s_{2}=s_{4}}\rangle|^{2}=|\langle\psi_{ij}^{SR\leq
2}|\psi_{ij}^{s_{2}=s_{4}}\rangle|^{2}\delta_{ip}\delta_{jq}\stackrel{{\scriptstyle
L.\ref{fact}}}{{\leq}}\frac{2}{d}\delta_{ip}\delta_{jq}.$ (25)

It follows that:

$\displaystyle\langle\psi^{SR\leq 2}|E(I\otimes G)E^{\dagger}|\psi^{SR\leq
2}\rangle=d\sum_{ij}|p_{ij}|^{2}|\langle\psi_{ij}^{SR\leq
2}|\psi_{ij}^{s_{2}=s_{4}}\rangle|^{2}$ (26) $\displaystyle\leq
2\sum_{ij}|p_{ij}|^{2}=2.$

A similar analysis is applied to $E(G\otimes I)E^{\dagger}$, getting

$\displaystyle\langle\psi^{SR\leq 2}|E(G\otimes I)E^{\dagger}|\psi^{SR\leq
2}\rangle=d\sum_{ij}|q_{ij}|^{2}|\langle\phi_{ij}^{SR\leq
2}|\psi_{ij}^{s_{1}=s_{3}}\rangle|^{2}$ (27) $\displaystyle\leq
2\sum_{ij}|q_{ij}|^{2}=2.$

A sufficient condition for Werner state 2-undistillability is:

$\displaystyle\langle\psi^{SR\leq 2}|I(e)+\beta^{2}G(e)+\beta E(I\otimes
G+G\otimes I)E^{\dagger}|\psi^{SR\leq 2}\rangle$ (28) $\displaystyle\geq
1+\beta^{2}|\langle\Psi(e)|\psi^{SR\leq 2}\rangle|^{2}+4\beta\geq 1+4\beta\geq
0.$

Therefore we can conclude that when $\beta\geq-\frac{1}{4}$, $\rho_{w}$ must
be 2-undistillable.

A similar process is followed in the arbitrary N case, yielding a lower bound
regarding all the odd terms in $(1+\beta)^{N}$. Since the odd terms in
$(1+\beta)^{N}$ can be expressed by $\frac{(1+\beta)^{N}-(1-\beta)^{N}}{2}$, a
lower bound can then be obtained for the N-undistillability case:

######  Theorem IV.2.

Let $\beta_{0}$ be the zero point of $1+(1+\beta)^{N}-(1-\beta)^{N}$ within
$[-1,0]$. When $\beta\geq\beta_{0}$, Werner state $\rho_{w}$ is
N-undistillable.

Proof of this is presented in Appendix B. By calculating the zero point a
lower bound is obtained where the corresponding Werner state is
N-undistillable. Note that unlike in previous works, this bound $\beta_{0}$ is
also independent of dimensionality $d$. By subsequently raising the
dimensionality of Werner states within the set, thus raising $-\frac{1}{d}$,
we can always obtain a set of Werner states that falls into
$[\beta_{0},-\frac{1}{d})$, which means that they are both NPT and
N-undistillable for any finite N. However, as N increases, higher
dimensionality are required for undistillability, and so the subset of
N-undistillable NPT states shrinks to zero as N approaches to infinity, which
is similar to the circumstances encountered in [18], [19] and [20]. Therefore,
this method fails to find a state that is undistillable for arbitrarily many
copies, while being NPT at the same time.

##  V The problem as a partial trace inequality

The problem of Werner state N-undistillability can be equivalently written as
inequalities regarding the Frobenius norm of partial traces of a rank-2
matrix.

Take 2-undistillability as an example. $\langle\psi^{SR\leq 2}|E(I\otimes
G+G\otimes I)E^{\dagger}|\psi^{SR\leq 2}\rangle$ can be written with regard to
the Frobenius norms of two partial traces:

$\langle\psi^{SR\leq 2}|E(I\otimes G)E^{\dagger}|\psi^{SR\leq
2}\rangle=||Tr_{2}(X_{2}(\psi))||_{F}^{2},$ (29) $\langle\psi^{SR\leq
2}|E(G\otimes I)E^{\dagger}|\psi^{SR\leq
2}\rangle=||Tr_{1}(X_{2}(\psi))||_{F}^{2},$ (30)

where $||\cdot||_{F}$ denotes Frobenius norm:

$||X||_{F}=\sqrt{\sum_{ij}|X_{ij}|^{2}}=\sqrt{Tr(X^{\dagger}X)}.$ (31)

$X_{2}(\psi)$ is a $d^{2}\times d^{2}$ matrix of rank $\leq$ 2, obtained by
state-operator isomorphism, from $|\psi^{SR\leq 2}\rangle$, a pure bipartite
$d^{2}\times d^{2}$ state with Schmidt rank $\leq$ 2\. For simplicity, we will
be writing it as $X_{2}$ from now on.

The normalization of $|\psi^{SR\leq 2}\rangle$ requires that
$||X_{2}||_{F}^{2}=Tr(X_{2}^{\dagger}X_{2})=1$. Now the problem of Werner
state 2-undistillability is equivalent to: Find a range of $\beta$ that makes
the following inequality always holds for arbitrary $X_{2}$ with rank no
larger than 2 (The requirement of $||X_{2}||_{F}^{2}=1$ can be lifted due to
homogeneity):

$\displaystyle||X_{2}||_{F}^{2}+\beta^{2}|Tr(X_{2})|^{2}$ (32)
$\displaystyle+\beta(||Tr_{1}(X_{2})||_{F}^{2}+||Tr_{2}(X_{2})||_{F}^{2})\geq
0.$

N-undistillability has a similar form:

######  Theorem V.1.

Werner state $\rho_{w}$ is N-undistillable iff the following holds for all
$d^{N}\times d^{N}$ $X_{2}$:

$\displaystyle\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}||Tr_{S}(X_{2})||_{F}^{2}\geq
0,$ (33)

where $Tr_{S}$ takes partial traces of the subsystems in set $S$.

The proof of this is presented in Appendix C. A similar result is obtained by
[29].

The partial trace inequalities are with regard to a matrix $X_{2}$ of rank one
or two. From now on we consider the case where $\beta$ is taken to be
$-\frac{1}{2}$ according to popular guess. For 2-undistillability problem, The
rank one case can be trivially proved by making use of the following lemma:

######  Lemma V.2.

The inequalities:

$||Tr_{1}(X_{1})||_{F}^{2}\leq||X_{1}||_{F}^{2},||Tr_{2}(X_{1})||_{F}^{2}\leq||X_{1}||_{F}^{2},$
(34)

hold for any $d^{2}\times d^{2}$ square matrix $X_{1}$ of rank one.

Proof of this lemma is presented in Appendix E. For $X_{2}$ of rank two, the
inequality of concern when $\beta=-\frac{1}{2}$ is:

$||Tr_{2}(X_{2})||_{F}^{2}+||Tr_{1}(X_{2})||_{F}^{2}-\frac{1}{2}|Tr(X_{2})|^{2}\leq
2||X_{2}||_{F}^{2}=2.$ (35)

For any square matrix $X_{2}$ of rank two, the singular decomposition can be
applied, decomposing the matrix as the sum of two rank one matrices:

$X_{2}=\sigma_{1}X_{2}^{1}+\sigma_{2}X_{2}^{2}=\sigma_{1}\mathbf{u}^{1}\mathbf{v}^{1\dagger}+\sigma_{2}\mathbf{u}^{2}\mathbf{v}^{2\dagger},$
(36)

where $\sigma_{1}^{2}+\sigma_{2}^{2}=1$, and both $\sigma_{1}$ and
$\sigma_{2}$ are positive. We introduce $d\times d$ matrices
$U_{1},U_{2},V_{1},V_{2}$ that are the result of state-operator isomorphisms,
corresponding to $\mathbf{u}^{1},\mathbf{u}^{2},\mathbf{v}^{1},\mathbf{v}^{2}$
respectively.

Eq.35 then becomes:

$\displaystyle\sigma_{1}^{2}(Tr(V_{1}U_{1}^{\dagger}U_{1}V_{1}^{\dagger})+Tr(U_{1}^{\dagger}V_{1}V_{1}^{\dagger}U_{1}))$
(37)
$\displaystyle+\sigma_{2}^{2}(Tr(V_{2}U_{2}^{\dagger}U_{2}V_{2}^{\dagger})+Tr(U_{2}^{\dagger}V_{2}V_{2}^{\dagger}U_{2}))$
$\displaystyle+\sigma_{1}\sigma_{2}2Re[Tr(V_{1}U_{1}^{\dagger}U_{2}V_{2}^{\dagger})+Tr(U_{1}^{\dagger}V_{1}V_{2}^{\dagger}U_{2})]$
$\displaystyle-\frac{|\sigma_{1}Tr(U_{1}V_{1}^{\dagger})+\sigma_{2}Tr(U_{2}V_{2}^{\dagger})|^{2}}{2}\leq
2.$

For simplicity the above is rewritten as:

$\sigma_{1}^{2}P+\sigma_{2}^{2}Q+\sigma_{1}\sigma_{2}R\leq 2,$ (38)

where $P,Q,R$ are:

$P=Tr(U_{1}^{\dagger}V_{1}V_{1}^{\dagger}U_{1})+Tr(V_{1}U_{1}^{\dagger}U_{1}V_{1}^{\dagger})-\frac{|Tr(U_{1}V_{1}^{\dagger})|^{2}}{2},$
(39)
$Q=Tr(U_{2}^{\dagger}V_{2}V_{2}^{\dagger}U_{2})+Tr(V_{2}U_{2}^{\dagger}U_{2}V_{2}^{\dagger})-\frac{|Tr(U_{2}V_{2}^{\dagger})|^{2}}{2},$
(40) $\displaystyle
R=2Re[Tr(V_{1}U_{1}^{\dagger}U_{2}V_{2}^{\dagger})+Tr(U_{1}^{\dagger}V_{1}V_{2}^{\dagger}U_{2})$
(41)
$\displaystyle-\frac{Tr^{*}(U_{1}V_{1}^{\dagger})Tr(U_{2}V_{2}^{\dagger})}{2}].$

Maximization regarding variables $\sigma_{1},\sigma_{2}$ reduces the question
to proving:

$R^{2}\leq 4(2-P)(2-Q),$ (42)

with normalization and orthogonality conditions requiring:

$||U_{1}||_{F}^{2}=||V_{1}||_{F}^{2}=||U_{2}||_{F}^{2}=||V_{2}||_{F}^{2}=1,$
(43)

and

$tr(V_{2}^{\dagger}V_{1})=tr(U_{2}^{\dagger}U_{1})=0.$ (44)

The special case of $U_{1}=V_{1}$, in other words, $X_{2}$ being the sum of a
normal matrix and a rank-1 matrix, has been proved in [29].

In fact, by a slight change of representation, we can also get an equivalent
expression of Werner state N-undistillability inequality, regarding only rank
one matrices:

######  Theorem V.3.

Werner state $\rho_{w}$ is N-undistillable iff

$Re[f_{N}(X_{1},X^{\prime}_{1})]^{2}\leq
f_{N}(X^{\prime}_{1},X^{\prime}_{1})f_{N}(X_{1},X_{1}),$ (45)

where

$\displaystyle
f_{N}(X_{1},X^{\prime}_{1})=\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}Tr[Tr^{\dagger}_{S}(X_{1})Tr_{S}(X^{\prime}_{1})]\geq
0,$ (46)

where
$X_{1}=\mathbf{w}^{\dagger}\mathbf{x},X^{\prime}_{1}=\mathbf{y}^{\dagger}\mathbf{z}$
are $d^{N}\times d^{N}$ rank one matrices, and their component vectors
$\mathbf{w},\mathbf{x},\mathbf{y},\mathbf{z}$ satisfy
$\mathbf{w}\perp\mathbf{y},\mathbf{x}\perp\mathbf{z}$.

Proof of this equivalence is presented in Appendix D.

It is then established that N-undistillability problem for arbitrary N can be
written as inequalities concerning rank one matrices. Although the above form
looks like a Cauchy-Schwartz inequality, $f_{N}(X_{1},X^{\prime}_{1})$ cannot
be seen as an inner product and enable the direct application of Cauchy-
Schwartz inequality. $f_{N}(X_{1},X^{\prime}_{1})$ is obviously conjugate
symmetric, and for rank one matrices $X_{1}$ and $X^{\prime}_{1}$, positivity
of $f_{N}(X_{1},X_{1})$ and $f_{N}(X^{\prime}_{1},X^{\prime}_{1})$ can be
ascertained. However, rank-one matrices do not compose a vector space, since
the sum of two rank-one matrices can be a rank-two matrix. Although arbitrary
finite rank matrices do compose a finite dimensional vector space, proving
positivity of $f_{N}(A,A)$ for arbitrary rank matrix $A$ is both beyond our
ability and our need, since proving positivity of $f_{N}(X_{2},X_{2})$ for
arbitrary rank two matrix $X_{2}$ is, in fact, equivalent to Eq.33, and is
enough for proof of N-undistillability.

##  VI Conversion to a multi-variable function

We make another attempt at this problem by seeing it as a multi-variable
function. For simplicity we consider only the real case of the
2-distillability problem, that is, $X_{1}$, $X^{\prime}_{1}$ are $d^{2}\times
d^{2}$ real matrices. For convenience of expression we use
$f(C,D)=f_{2}(X_{1},X^{\prime}_{1})$ to denote the function concerned. Thus,
the 2-undistillability problem is equivalent to proving the following
inequality for all $d^{2}\times d^{2}$-dimensional rank one matrices $C,D$:

$f^{2}(C,D)\leq f(C,C)f(D,D),$ (47)

where

$\displaystyle f(C,D)=Tr(C^{T}D)$ (48)
$\displaystyle+\beta(Tr(C_{1}^{T}D_{1})+Tr(C_{2}^{T}D_{2}))+\beta^{2}Tr(C^{T})Tr(D),$

where we have used
$Tr_{1}(C)=C_{2},Tr_{2}(C)=C_{1},Tr_{1}(D)=D_{2},Tr_{2}(D)=D_{1}$ for
simplicity. A multi-variable function $g(C)$ is defined if we see $D$ as a
constant $D_{0}$:

$g(C)_{D_{0}}=f(C,C)f(D_{0},D_{0})-f^{2}(C,D_{0}),$ (49)

where $D_{0}$ is an arbitrary $d^{2}\times d^{2}$-dimensional real matrix of
rank 1. Both $C$ and $D_{0}$ can be written as outer products of
$d^{2}$-dimensional real vectors $\mathbf{w}$,
$\mathbf{x},\mathbf{y},\mathbf{z}$:

$C=\mathbf{w}\mathbf{x}^{T},D_{0}=\mathbf{y}\mathbf{z}^{T}.$ (50)

The multi variables of function $g(C)$ are taken to be the components of
vector $\mathbf{w},\mathbf{x}$, represented by $w_{ij},x_{ij}$. Indices i and
j take their values in $[0,1,\cdots,d-1]$, and are combined together to
represent the $(i\times d+j)$-th component of the vector. The problem is
therefore transformed into proving that multi-variable function $g(C)_{D_{0}}$
is always non-negative, for all variables $w_{ij},x_{ij}$ and all possible
parameters $y_{ij},z_{ij}$. It is clear that when $C=D_{0}$,
$g(D_{0})_{D_{0}}=f^{2}(D_{0},D_{0})-f^{2}(D_{0},D_{0})=0$, and we conjecture
that these are global minimums.

We first prove that the gradients at these points are zero regardless of the
parameters $\mathbf{y},\mathbf{z}$ by calculating the Jacobian matrix,
therefore showing that the $C=D_{0}$ points are indeed critical points. Then
the Hessian matrix at this point is presented, in particular its three block
parts, since the Hessian matrix is Hermitian. We conjecture that it is
positive semidefinite. An additional proof of nonconvexity is also given in
Appendix.F, thus eliminating the easy case where any local minimum is also a
global minimum.

###  VI.1 The gradient at $C=D_{0}$ is zero

For simplicity all $g(C)_{D_{0}}$ are written as $g(C)$ in the following text.
Define

$\displaystyle h_{1,ij}(w,x)=2w_{ij}\left(\sum_{pq}x_{pq}^{2}\right)$ (51)
$\displaystyle+2\beta\left[\sum_{pq}x_{pj}w_{iq}x_{pq}+\sum_{pq}x_{iq}w_{pj}x_{pq}\right]$
$\displaystyle+2\beta^{2}x_{ij}\left(\sum_{pq}w_{pq}x_{pq}\right)$
$\displaystyle h_{2,ij}(y,z,x)=y_{ij}\left(\sum_{pq}x_{pq}z_{pq}\right)$ (52)
$\displaystyle+\beta\left[\sum_{pq}x_{pj}y_{iq}z_{pq}+\sum_{pq}x_{iq}y_{pj}z_{pq}\right]$
$\displaystyle+\beta^{2}x_{ij}\left(\sum_{pq}y_{pq}z_{pq}\right)$

The partial derivatives of $f(C,C)$ and $f(C,D_{0})$ are:

$\frac{\partial f(C,C)}{\partial w_{ij}}=h_{1,ij}(w,x),\frac{\partial
f(C,C)}{\partial x_{ij}}=h_{1,ij}(x,w)$ (53) $\frac{\partial
f(C,D_{0})}{\partial w_{ij}}=h_{2,ij}(y,z,x),\frac{\partial
f(C,D_{0})}{\partial x_{ij}}=h_{2,ij}(z,y,w)$ (54)

Notice that $h_{1,ij}(w,x)=2h_{2,ij}(w,x,x),h_{1,ij}(x,w)=2h_{2,ij}(x,w,w)$.

Therefore, at the critical points where $C=D_{0}$,

$\frac{\partial f(C,D_{0})}{\partial
w_{ij}}|_{C=D_{0}}=\frac{1}{2}\frac{\partial f(C,C)}{\partial
w_{ij}}|_{C=D_{0}},$ (55) $\frac{\partial f(C,D_{0})}{\partial
x_{ij}}|_{C=D_{0}}=\frac{1}{2}\frac{\partial f(C,C)}{\partial
x_{ij}}|_{C=D_{0}}.$ (56)

It then follows that the Jacobian at the $C=D_{0}$ point is zero:

$\displaystyle\frac{\partial g(C)}{\partial
w_{ij}}|_{C=D_{0}}=f(D_{0},D_{0})\frac{\partial f(C,C)}{\partial
w_{ij}}|_{C=D_{0}}$ (57) $\displaystyle-2\frac{\partial f(C,D_{0})}{\partial
w_{ij}}|_{C=D_{0}}f(D_{0},D_{0})=0.$ $\displaystyle\frac{\partial
g(C)}{\partial x_{ij}}|_{C=D_{0}}=f(D_{0},D_{0})\frac{\partial
f(C,C)}{\partial x_{ij}}|_{C=D_{0}}$ (58) $\displaystyle-2\frac{\partial
f(C,D_{0})}{\partial x_{ij}}|_{C=D_{0}}f(D_{0},D_{0})=0.$

It is then established that the gradients at the $C=D_{0}$ points are zero,
and the points are indeed critical points.

###  VI.2 Hessian matrix

Hessian matrices at the critical points are presented. Define:

$\displaystyle
h_{3,ijkl}(x)=2\delta_{ik}\delta_{jl}\left(\sum_{pq}x^{2}_{pq}\right)$ (59)
$\displaystyle+2\beta\left[\delta_{ik}\sum_{p}x_{pj}x_{pl}+\delta_{lj}\sum_{q}x_{iq}x_{kq}\right]$
$\displaystyle+2\beta^{2}x_{ij}x_{kl}$ $\displaystyle
h_{4,ijkl}(w,x)=4w_{ij}x_{kl}$ (60)
$\displaystyle+\beta\left[2\delta_{jl}(\sum_{q}w_{iq}x_{kq})+x_{kj}w_{il}+2\delta_{ik}(\sum_{p}w_{pj}x_{pl})+x_{il}w_{kj}\right]$
$\displaystyle+2\beta^{2}\left[\delta_{ik}\delta_{jl}\left(\sum_{pq}w_{pq}x_{pq}\right)+x_{ij}w_{kl}\right]$
$\displaystyle h_{5,ijkl}(y,z)=y_{ij}z_{kl}$ (61)
$\displaystyle+\beta\left[\delta_{jl}(\sum_{q}y_{iq}z_{kq})+\delta_{ik}(\sum_{p}y_{pj}z_{pl})\right]$
$\displaystyle+\beta^{2}\delta_{ik}\delta_{jl}\left(\sum_{pq}y_{pq}z_{pq}\right)$

The Hessian matrices at $C=D_{0}$ can be explicitly written as:

$\begin{pmatrix}(\frac{\partial^{2}g(C)}{\partial w_{ij}\partial
w_{kl}})|_{C=D_{0}}&(\frac{\partial^{2}g(C)}{\partial w_{ij}\partial
x_{kl}})|_{C=D_{0}}\\\ (\frac{\partial^{2}g(C)}{\partial x_{ij}\partial
w_{kl}})|_{C=D_{0}}&(\frac{\partial^{2}g(C)}{\partial x_{ij}\partial
x_{kl}})|_{C=D_{0}}\\\ \end{pmatrix}.$ (62)

where the three independent parts of Hessian matrix are:

$\displaystyle\frac{\partial^{2}g(C)}{\partial w_{ij}\partial
w_{kl}}|_{C=D_{0}}$ (63)
$\displaystyle=f(D_{0},D_{0})h_{3,ijkl}(x)-2h_{2,ij}(w,x,x)h_{2,kl}(w,x,x),$
$\displaystyle\frac{\partial^{2}g(C)}{\partial x_{ij}\partial
x_{kl}}|_{C=D_{0}}$ (64)
$\displaystyle=f(D_{0},D_{0})h_{3,ijkl}(w)-2h_{2,ij}(x,w,w)h_{2,kl}(x,w,w),$

$\displaystyle\frac{\partial^{2}g(C)}{\partial w_{ij}\partial
x_{kl}}|_{C=D_{0}}=f(D_{0},D_{0})h_{4,ijkl}(w,x)$ (65)
$\displaystyle-2h_{5,ijkl}(w,x)f(D_{0},D_{0})-2h_{2,ij}(w,x,w)h_{2,kl}(x,w,w).$

We conjecture that this Hessian matrix is positive semidefinite, thus making
the critical points local minimums. In Appendix F we prove that the function
is non-convex, which means that further scrutinization is needed for
characterization of this function.

## Conclusion and Discussion

We’ve broken down the process of verifying bound entanglement into iterative
steps of 1-undistillability verifications, by noticing that N-undistillability
for arbitrary N is equivalent to $2^{N}$-undistillability for arbitrary N, and
then utilizing specific symmetric properties. In the Werner state case, new
bounds for N-undistillability of any finite N are presented, a result similar
to that of [18], [19] and [20], but different in the sense that our bounds are
unaffected by the dimensionality of Werner states. Alternative expressions for
inequalities applicable to both rank two and rank one matrices are given.
Subsequently, the problem of two-undistillability is converted into a matrix
analysis problem. The multi-variable function treatment is also attempted, as
well as proving critical points, non-convexity and conjecturing about Hessian
positivity.

Recently, the problem of N-copy Werner state distillability has been
reformulated in [29] into a set of partial trace inequalities, which coincided
with our Theorem V.1. In the regime of 2-copy Werner state distillability, A
special case of the matrix being the sum of a rank-one matrix and a normal
matrix is proved.

We believe the new perspectives presented here will assist further attempts at
this famous open problem.

###### Acknowledgements.

We would like to thank Prof. Fedor Sukochev, Dr. Dmitriy Zanin and Prof. Zhi
Yin for helpful comments and valuable insights about partial trace
inequalities.

## References

  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Huber  _et al._ [2018] M. Huber, L. Lami, C. Lancien, and A. Müller-Hermes, High-dimensional entanglement in states with positive partial transposition, Physical Review Letters 121, 10.1103/physrevlett.121.200503 (2018).
  * Bennett  _et al._ [1996] C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, Purification of noisy entanglement and faithful teleportation via noisy channels, Physical Review Letters 76, 722–725 (1996).
  * Horodecki  _et al._ [1998] M. Horodecki, P. Horodecki, and R. Horodecki, Mixed-state entanglement and distillation: Is there a “bound” entanglement in nature?, Physical Review Letters 80, 5239–5242 (1998).
  * Horodecki  _et al._ [2020] P. Horodecki, Łukasz Rudnicki, and K. Życzkowski, Five open problems in quantum information (2020), arXiv:2002.03233 [quant-ph] .
  * Chen and Doković [2016a] L. Chen and D. u. Doković, Distillability of non-positive-partial-transpose bipartite quantum states of rank four, Physical Review A 94, 10.1103/physreva.94.052318 (2016a).
  * Chen  _et al._ [2023] L. Chen, C. Feng, and H. He, On the distillability problem of two-copy $4\times 4$ Werner states: matrices A, B of at most four nonzero entries, Quant. Inf. Proc. 22, 108 (2023).
  * Chen  _et al._ [2021] L. Chen, H. He, X. Shi, and L.-J. Zhao, Proving the distillability problem of two-copy $4\times 4$ Werner states for monomial matrices, Quant. Inf. Proc. 20, 157 (2021).
  * Qian  _et al._ [2019] L. Qian, L. Chen, D. Chu, and Y. Shen, A matrix inequality related to the entanglement distillation problem (2019), arXiv:1908.02428 [quant-ph] .
  * Ding and Chen [2023] T. Ding and L. Chen, Entanglement distillation in terms of schmidt rank and matrix rank (2023), arXiv:2304.05563 [quant-ph] .
  * Pankowski  _et al._ [2013] L. Pankowski, F. G. S. L. Brandao, M. Horodecki, and G. Smith, Entanglement distillation by extendible maps (2013), arXiv:1109.1779 [quant-ph] .
  * Lami  _et al._ [2023] L. Lami, B. Regula, and A. Streltsov, Catalysis cannot overcome bound entanglement (2023), arXiv:2305.03489 [quant-ph] .
  * Chitambar  _et al._ [2020] E. Chitambar, J. I. de Vicente, M. W. Girard, and G. Gour, Entanglement manipulation beyond local operations and classical communication, Journal of Mathematical Physics 61, 10.1063/1.5124109 (2020).
  * Brandão  _et al._ [2011] F. G. S. L. Brandão, M. Christandl, and J. Yard, Faithful squashed entanglement, Communications in Mathematical Physics 306, 805–830 (2011).
  * Johnston [2011] N. Johnston, Characterizing operations preserving separability measures via linear preserver problems, Linear and Multilinear Algebra 59, 1171–1187 (2011).
  * Sharma and Pandey [2016] K. K. Sharma and S. N. Pandey, Dzyaloshinskii-moriya interaction as an agent to free the bound entangled states (2016), arXiv:1501.00942 [quant-ph] .
  * van der Eyden  _et al._ [2022] M. van der Eyden, T. Netzer, and G. De las Cuevas, Halos and undecidability of tensor stable positive maps, Journal of Physics A: Mathematical and Theoretical 55, 264006 (2022).
  * DiVincenzo  _et al._ [2000] D. P. DiVincenzo, P. W. Shor, J. A. Smolin, B. M. Terhal, and A. V. Thapliyal, Evidence for bound entangled states with negative partial transpose, Physical Review A 61, 10.1103/physreva.61.062312 (2000).
  * Dür  _et al._ [2000] W. Dür, J. I. Cirac, M. Lewenstein, and D. Bruß, Distillability and partial transposition in bipartite systems, Physical Review A 61, 10.1103/physreva.61.062313 (2000).
  * Chen and Doković [2016b] L. Chen and D. u. Doković, Distillability of non-positive-partial-transpose bipartite quantum states of rank four, Physical Review A 94, 10.1103/physreva.94.052318 (2016b).
  * Watrous [2004] J. Watrous, Many copies may be required for entanglement distillation, Physical Review Letters 93, 10.1103/physrevlett.93.010502 (2004).
  * Werner [1989] R. F. Werner, Quantum states with einstein-podolsky-rosen correlations admitting a hidden-variable model, Phys. Rev. A 40, 4277 (1989).
  * Horodecki and Horodecki [1998] M. Horodecki and P. Horodecki, Reduction criterion of separability and limits for a class of protocols of entanglement distillation (1998), arXiv:quant-ph/9708015 [quant-ph] .
  * Pankowski  _et al._ [2010] Ł. Pankowski, M. Piani, M. Horodecki, and P. Horodecki, A few steps more towards npt bound entanglement, IEEE Transactions on Information Theory 56, 4085 (2010).
  * Doković [2016] D. Doković, On two-distillable werner states, Entropy 18, 216 (2016).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Physical Review Letters 77, 1413 (1996).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Horodecki  _et al._ [1997] M. Horodecki, P. Horodecki, and R. Horodecki, Inseparable two spin- $\frac{1}{2}$ density matrices can be distilled to a singlet form, Phys. Rev. Lett. 78, 574 (1997).
  * Rico [2023] P. C. Rico, New partial trace inequalities and distillability of werner states (2023), arXiv:2310.05726 [math-ph] .
  * Horn and Johnson [1991] R. A. Horn and C. R. Johnson,  _Topics in Matrix Analysis_ (Cambridge University Press, 1991).

##  Appendix A Proof of Lemma IV.1

######  Lemma A.1.

Let $|\psi\rangle$ be an arbitrary pure quantum state and $s_{j}$ are the
Schmidt coefficients of $|\psi\rangle$ arranged in descending order, then

$\;\sum_{j=1}^{k}s_{j}^{2}=max_{|\phi^{SR\leq
k}\rangle}|\langle\psi|\phi^{SR\leq k}\rangle|^{2},$ (66)

where $|\phi^{SR\leq k}\rangle$ is an arbitrary quantum pure state with
Schimidt rank being less than k, namely, $SR(|\phi^{SR\leq k}\rangle)\leq k$.

###### Proof.

State-operator isomorphism translates the problem to an equivalent form of
proving:

$|Tr(A^{\dagger}B)|^{2}\leq\sum_{i=1}^{k}s^{2}_{i}(A),$ (67)

where both $A$ and $B$ are $d\times d$ matrices, and as they are the result of
state-operator isomorphism from $|\psi\rangle$ and $|\phi^{SR\leq 2}\rangle$,
are of Frobenius norm 1 due to normalization conditions, with $B$ having rank
k. We use $s_{i}(A)$ and $|\lambda_{i}(A)|$ to denote the singular values and
absolute values of eigenvalues of $A$ arranged in descending order.

$\displaystyle|Tr(A^{\dagger}B)|^{2}\leq(\sum_{i=1}^{d}|\lambda_{i}(A^{\dagger}B)|)^{2}\leq(\sum_{i=1}^{d}s_{i}(A^{\dagger}B))^{2}$
(68)
$\displaystyle\leq(\sum_{i=1}^{k}s_{i}(A^{\dagger})s_{i}(B))^{2}=(\sum_{i=1}^{k}s_{i}(A)s_{i}(B))^{2},$

where the second and third inequality are from Theorem 3.3.13(a) and Theorem
3.3.14(a) in [30], and $d\to k$ change of index holds because the rest of
$B$’s singular values are all zeroes after $d$. respectively. Then by Cauchy-
Schwartz inequality,

$\displaystyle|Tr(A^{\dagger}B)|^{2}\leq(\sum_{i=1}^{k}s_{i}(A)s_{i}(B))^{2}\leq(\sum_{i=1}^{k}s^{2}_{i}(A))(\sum_{j=1}^{k}s^{2}_{j}(B))$
(69) $\displaystyle=\sum_{i=1}^{k}s^{2}_{i}(A),$

where the last equality follows from the normalization condition of $B$. ∎

##  Appendix B Proof of Theorem IV.2

######  Theorem B.1.

Let $\beta_{0}$ be the zero point of $1+(1+\beta)^{N}-(1-\beta)^{N}$ within
$[-1,0]$. When $\beta\geq\beta_{0}$, Werner state $\rho_{w}$ is
N-undistillable.

###### Proof.

We can overlook the $\frac{1}{d^{2}+\beta d}$ factor since it doesn’t affect
positivity.

$(\rho_{w}^{T_{A}})^{\otimes N}=I^{\otimes
N}+\sum_{m=1}^{N}\beta^{m}\sum_{seq(m,N)}G^{i_{1}}\otimes
G^{i_{2}}\otimes\cdots\otimes G^{i_{N}},$ (70)

where $seq(m,N)$ denotes all possible binary $i_{1},\dots,i_{N}$ sequences
with $m$ ones and $N-m$ zeros. For convenience denote
$S_{0}=\\{n|i_{n}=0\\},S_{1}=\\{n|i_{n}=1\\}$.

$\displaystyle M_{N}G^{i_{1}}\otimes G^{i_{2}}\otimes\cdots\otimes
G^{i_{N}}M_{N}^{\dagger}$ (71) $\displaystyle=M_{N}\sum_{j_{1}\dots
j_{2N}}|\cdots j_{2n-1}j_{2n-i_{n}}\cdots\rangle\langle\cdots
j_{2n-1+i_{n}}j_{2n}\cdots|M_{N}^{\dagger}$ $\displaystyle=\sum_{j_{1}\dots
j_{2N}}|\cdots j_{2n-1}\cdots,\cdots j_{2n-i_{n}}\cdots\rangle\langle\cdots
j_{2n-1+i_{n}}\cdots,\cdots j_{2n}\cdots|$
$\displaystyle=d^{m}\sum_{j_{2n-1},j_{2n},n\in
S_{0}}|\psi_{j_{2n-1},j_{2n},n\in
S_{0}}\rangle\langle\psi_{j_{2n-1},j_{2n},n\in S_{0}}|,$

where

$\displaystyle|\psi_{j_{2n-1},j_{2n},n\in S_{0}}\rangle$ (72)
$\displaystyle=\frac{1}{\sqrt{d^{m}}}\sum_{j_{2n-1},n\in S_{1}}|\dots
j_{2n-1}\dots,\dots j_{2n-i_{n}}\dots\rangle.$

Any pure quantum state $|\psi^{SR\leq 2}\rangle$ of Schmidt rank no larger
than 2 can be decomposed into:

$|\psi^{SR\leq 2}\rangle=\sum_{j_{2n-1},j_{2n},n\in
S_{0}}p_{j_{2n-1},j_{2n},n\in S_{0}}|\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq
2}\rangle,$ (73)

where $|\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq 2}\rangle$ is the normalized
extraction of all $|\cdots j_{2n-1}\cdots,\cdots j_{2n}\cdots\rangle,n\in
S_{0}$ terms in $|\psi^{SR\leq 2}\rangle$, namely, if

$|\psi^{SR\leq
2}\rangle=\sum_{k_{1},\cdots,k_{2N}}c_{k_{1},\cdots,k_{2N}}|k_{1},\cdots,k_{2N}\rangle,$
(74)

then,

$\displaystyle|\psi^{SR\leq 2}_{j_{2n-1},j_{2n},n\in S_{0}}\rangle$ (75)
$\displaystyle=\frac{1}{p_{j_{2n-1},j_{2n},n\in S_{0}}}$
$\displaystyle\sum_{j_{2n-1},j_{2n},n\in S_{1}}c_{\cdots j_{2n-1}\cdots,\cdots
j_{2n}\cdots}|\cdots j_{2n-1}\cdots,\cdots j_{2n}\cdots\rangle,$

where

$p_{j_{2n-1},j_{2n},n\in S_{0}}=\sqrt{\sum_{j_{2n-1},j_{2n},n\in
S_{1}}|c_{\cdots j_{2n-1}\cdots,\cdots j_{2n}\cdots}|^{2}}$ (76)

It can be proved that $|\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq 2}\rangle$
also has Schmidt rank less than or equal to 2. According to state-operator
isomorphism, Schmidt rank of a state is equal to rank of the corresponding
operator. In fact, the operator corresponding to $|\psi_{j_{2n-1},j_{2n},n\in
S_{0}}^{SR\leq 2}\rangle$ is the extraction of $d^{m}$ rows $\cdots
j_{2n-1}\cdots,n\in S_{0},$ and $d^{m}$ columns $\cdots j_{2n}\cdots,n\in
S_{0}$. The row and column numbers $\cdots j_{2n-1}\cdots$ and $\cdots
j_{2n}\cdots$ are written in base-d numeral systems. Therefore, operator
corresponding to $|\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq 2}\rangle$ is a
submatrix of the operator corresponding to $|\psi^{SR\leq 2}\rangle$. The rank
of a submatrix is no larger than the rank of the entire matrix, therefore
$|\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq 2}\rangle$ must have Schmidt rank
no larger than 2. Using Lem.IV.1, it is then clear that:

$|\langle\psi_{j_{2n-1},j_{2n},n\in S_{0}}^{SR\leq
2}|\psi_{j_{2n-1},j_{2n},n\in S_{0}}\rangle|^{2}\leq\frac{2}{d^{m}},$ (77)

and so

$\displaystyle\langle\psi^{SR\leq 2}|M_{N}G^{i_{1}}\otimes
G^{i_{2}}\otimes\cdots\otimes G^{i_{N}}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle$
(78) $\displaystyle\leq d^{m}\frac{2}{d^{m}}\sum_{j_{2n-1},j_{2n},n\in
S_{0}}|p_{j_{2n-1},j_{2n},n\in S_{0}}|^{2}=2.$

Considering the fact that all even terms in $(\rho_{w}^{T_{A}})^{\otimes N}$
are non-negative and using the above inequality on all odd terms, a lower
bound is obtained:

$\displaystyle\langle\psi^{SR\leq 2}|M_{N}(\rho_{w}^{T_{A}})^{\otimes
N}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle$ (79) $\displaystyle\geq
1+\sum_{k=1}^{\lceil\frac{N}{2}\rceil}2C_{N}^{2k-1}\beta^{2k-1}=1+(1+\beta)^{N}-(1-\beta)^{N}$

For $\beta\in(-1,0)$, any $\beta$ greater than the zero point $\beta_{0}$
would make $1+(1+\beta)^{N}-(1-\beta)^{N}$ greater than zero, thus ensuring
the positivity of $\langle\psi^{SR\leq 2}|M_{N}(\rho_{w}^{T_{A}})^{\otimes
N}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle$.

∎

##  Appendix C Proof of Theorem V.1

######  Theorem C.1.

Werner state $\rho_{w}$ is N-undistillable iff the following holds for all
$d^{N}\times d^{N}$ $X_{2}$:

$\displaystyle\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}||Tr_{S}(X_{2})||_{F}^{2}\geq
0,$ (80)

where $Tr_{S}$ takes partial traces of the subsystems in set $S$.

###### Proof.

We’ve already established that N-undistillability is equivalent to:

$\langle\psi^{SR\leq 2}|M_{N}(\rho^{T_{A}})^{\otimes
N}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle\geq 0.$ (81)

In the previous section it has been established that:

$(\rho_{w}^{T_{A}})^{\otimes N}=I^{\otimes
N}+\sum_{m=1}^{N}\beta^{m}\sum_{seq(m,N)}G^{i_{1}}\otimes
G^{i_{2}}\otimes\cdots\otimes G^{i_{N}},$ (82) $\displaystyle
M_{N}G^{i_{1}}\otimes G^{i_{2}}\otimes\cdots\otimes G^{i_{N}}M_{N}^{\dagger}$
(83) $\displaystyle=d^{m}\sum_{j_{2n-1},j_{2n},n\in
S_{0}}|\psi_{j_{2n-1},j_{2n},n\in
S_{0}}\rangle\langle\psi_{j_{2n-1},j_{2n},n\in S_{0}}|,$

where

$\displaystyle|\psi_{j_{2n-1},j_{2n},n\in S_{0}}\rangle$ (84)
$\displaystyle=\frac{1}{\sqrt{d^{m}}}\sum_{j_{2n-1},n\in S_{1}}|\dots
j_{2n-1}\dots,\dots j_{2n-i_{n}}\dots\rangle.$

For a pure quantum state of the form:

$|\psi^{SR\leq 2}\rangle=\sum_{j_{1},\cdots,j_{N}}c_{\cdots
j_{2n-1}\cdots,\cdots j_{2n}\cdots}|\cdots j_{2n-1}\cdots,\cdots
j_{2n}\cdots\rangle,$ (85) $\displaystyle\langle\psi^{SR\leq
2}|M_{N}G^{i_{1}}\otimes G^{i_{2}}\otimes\cdots\otimes
G^{i_{N}}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle$ (86)
$\displaystyle=d^{m}\sum_{j_{2n-1},j_{2n},n\in
S_{0}}|\langle\psi_{j_{2n-1},j_{2n},n\in S_{0}}|\psi^{SR\leq 2}\rangle|^{2}$
$\displaystyle=\sum_{j_{2n-1},j_{2n},n\in S_{0}}|\sum_{j_{2n-1},n\in
S_{1}}c_{\cdots j_{2n-1}\cdots,\cdots j_{2n-i_{n}}\cdots}|^{2}$
$\displaystyle=||Tr_{\\{i_{n}|n\in S_{1}\\}}(X_{2})||_{F}^{2}.$

In the last equality, $X_{2}$ is the result of state-operator isomorphism from
$|\psi^{SR\leq 2}\rangle$. Noticing that

$\langle\psi^{SR\leq 2}|M_{N}I^{\otimes N}M_{N}^{\dagger}|\psi^{SR\leq
2}\rangle=1=\beta^{0}||X_{2}||_{F}^{2},$ (87)

it then follows that

$\langle\psi^{SR\leq 2}|M_{N}(\rho^{T_{A}})^{\otimes
N}M_{N}^{\dagger}|\psi^{SR\leq 2}\rangle\geq 0,$ (88)

is equivalent to

$\displaystyle\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}||Tr_{S}(X_{2})||_{F}^{2}\geq
0,$ (89)

∎

##  Appendix D Proof of Theorem V.3

######  Theorem D.1.

Werner state $\rho_{w}$ is N-undistillable iff

$Re[f_{N}(X_{1},X^{\prime}_{1})]^{2}\leq
f_{N}(X^{\prime}_{1},X^{\prime}_{1})f_{N}(X_{1},X_{1}),$ (90)

where

$\displaystyle
f_{N}(X_{1},X^{\prime}_{1})=\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}Tr[Tr^{\dagger}_{S}(X_{1})Tr_{S}(X^{\prime}_{1})]\geq
0,$ (91)

where
$X_{1}=\mathbf{w}^{\dagger}\mathbf{x},X^{\prime}_{1}=\mathbf{y}^{\dagger}\mathbf{z}$
are $d^{N}\times d^{N}$ rank one matrices, and their component vectors
$\mathbf{w},\mathbf{x},\mathbf{y},\mathbf{z}$ satisfy
$\mathbf{w}\perp\mathbf{y},\mathbf{x}\perp\mathbf{z}$.

###### Proof.

It has been proven that N-undistillability is equivalent to:

$\displaystyle\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}||Tr_{S}(X_{2})||_{F}^{2}\geq
0,$ (92)

Since any rank two matrix can be decomposed via singular value decomposition:

$X_{2}=\sigma_{1}X_{1}+\sigma_{2}X^{\prime}_{1},$ (93)

Eq.92 can be written as:

$\displaystyle\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}||Tr_{S}(X_{2})||_{F}^{2}$
(94)
$\displaystyle=\sum_{S\subset\\{1,\cdots,N\\}}\beta^{|S|}[\sigma^{2}_{1}||Tr_{S}(X_{1})||_{F}^{2}+\sigma_{2}^{2}||Tr_{S}(X^{\prime}_{1})||_{F}^{2}$
$\displaystyle+2\sigma_{1}\sigma_{2}Re(Tr(Tr_{S}(X_{1})^{\dagger}Tr_{S}(X^{\prime}_{1})))]$
$\displaystyle=\sigma_{1}^{2}f_{N}(X_{1},X_{1})+\sigma_{2}^{2}f_{N}(X^{\prime}_{1},X^{\prime}_{1})$
$\displaystyle+2\sigma_{1}\sigma_{2}Re[f_{N}(X_{1},X^{\prime}_{1})]$
$\displaystyle\geq 0.$

The above should hold for all singular values $\sigma_{1},\sigma_{2}$, which
is then equivalent to:

$Re[f_{N}(X_{1},X^{\prime}_{1})]^{2}\leq
f_{N}(X_{1},X_{1})f_{N}(X^{\prime}_{1},X^{\prime}_{1}),$ (95)

always holding when
$X_{1}=\mathbf{w}^{\dagger}\mathbf{x},X^{\prime}_{1}=\mathbf{y}^{\dagger}\mathbf{z}$
are $d^{N}\times d^{N}$ rank one matrices, and their component vectors
$\mathbf{w},\mathbf{x},\mathbf{y},\mathbf{z}$ satisfy
$\mathbf{w}\perp\mathbf{y},\mathbf{x}\perp\mathbf{z}$. ∎

##  Appendix E Proof of Lem. V.2

######  Lemma E.1.

The inequalities:

$||Tr_{1}(X_{1})||_{F}^{2}\leq||X_{1}||_{F}^{2},||Tr_{2}(X_{1})||_{F}^{2}\leq||X_{1}||_{F}^{2},$
(96)

hold for any $d^{2}\times d^{2}$ square matrix $X_{1}$ of rank one.

###### Proof.

A rank one matrix $X_{1}$ can always be written as the outer product of two
vectors: (in this case, bipartite)

$(X_{1})_{ij,kl}=w_{ij}x_{kl}.$ (97)

Its partial traces can be computed accordingly:

$(Tr_{2}(X_{1}))_{i,k}=\sum_{j}w_{ij}x_{kj},(Tr_{1}(X_{1}))_{j,l}=\sum_{i}w_{ij}x_{il}.$
(98)

A direct calculation and application of Cauchy-Schwartz inequality yields the
desired result:

$\displaystyle||Tr_{2}(X_{1})||_{F}^{2}=\sum_{ik}|\sum_{j}w_{ij}x_{kj}|^{2}\leq\sum_{ik}\sum_{j}|w_{ij}|^{2}\sum_{l}|x_{kl}|^{2}$
(99) $\displaystyle=\sum_{ijkl}|w_{ij}|^{2}|x_{kl}|^{2}=||X_{1}||_{F}^{2},$
$\displaystyle||Tr_{1}(X_{1})||_{F}^{2}=\sum_{jl}|\sum_{i}w_{ij}x_{il}|^{2}\leq\sum_{jl}\sum_{i}|w_{ij}|^{2}\sum_{k}|x_{kl}|^{2}$
(100) $\displaystyle=\sum_{ijkl}|w_{ij}|^{2}|x_{kl}|^{2}=||X_{1}||_{F}^{2}.$

∎

##  Appendix F Proving non-convexity

We now prove that the function of $g(C)_{D_{0}}$ is non-convex, by showing
that its local minimum set is not a convex one. For

$C=\mathbf{w}\mathbf{x}^{T},D_{0}=\mathbf{y}\mathbf{z}^{T},$ (101)

we set vectors $\mathbf{y}$ and $\mathbf{z}$ to identical forms of:

$y_{ij}=z_{ij}=\delta_{i0}\delta_{j1}$ (102)

For simplicity we first write the variables $\mathbf{w},\mathbf{x}$ in a
matrix form of:

$\begin{pmatrix}w_{00}&w_{01}&\cdots&w_{d-1,d-1}\\\
x_{00}&x_{01}&\cdots&x_{d-1,d-1}\\\ \end{pmatrix},$

where the two rows correspond to two vectors $\mathbf{w}$, $\mathbf{x}$,
respectively. We take the middle point combination of the following two
points:

$\begin{pmatrix}0&1&\cdots 0\\\ 0&1&\cdots 0\\\
\end{pmatrix},\begin{pmatrix}1&0&\cdots 0\\\ 1&0&\cdots 0\\\ \end{pmatrix},$

both of these points can be proven to have zero Jacobian matrix and positive
definite Hessian, effectively making them local minimums.

For a convex function, its local minimum set must be a convex set, so that
middle point combination of any two points should stay in the set. Their
middle point combination is:

$\begin{pmatrix}1&1&0&\cdots 0\\\ 1&1&0&\cdots 0\\\ \end{pmatrix}.$

The normalization factor doesn’t change whether the Jacobian matrix is nonzero
or not, and therefore is interchangeable and omitted here. At this particular
point, the Jacobian matrix is nonzero and proportional to:

$[1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0],$ (103)

which suggests that the middle point is not a local minimum, thus proving the
non-convexity of the set of local minimums. It then follows that the function
of $g(C,D)$ is non-convex.