Joseph Barreto and Daniel Lidar
This repository contains all scripts necessary to generate and plot the data associated with
The latest release of this code is archived at
The data associated with the arXiv submission is archived at
├── alpha_gen.py
├── data_analysis.ipynb
├── data_gen.py
├── sparsity.py
├── environment.yml
├── README.md
├── slurm
│ ├── config.sh.example # copy to config.sh (gitignored) and fill in your credentials
│ ├── data_gen.slurm
│ ├── alpha_gen.slurm
│ └── sparsity.slurm
└── utils
├── cs_utils.py
├── misc_utils.py
├── pauli_utils.py
└── shadow_utils.py
The environment.yml specifies all package and dependency versions. You can recreate this environment by running conda env create -f environment.yml in a conda command prompt. The basic package dependencies are
- numpy
- scipy
- qutip
- scikit-learn
- matplotlib
- pandas
- tqdm
- networkx
- ipykernel
- matplotlib-label-lines
- statsmodels
- threadpoolctl
Run data_gen.py to create the following data files. The folder 2x2_heis_neel is the default test system, but it is generally set to {nx}x{ny}_{ham}_{istate} — or, if --eps is nonzero, {nx}x{ny}_{ham}_{istate}_eps={eps:.1e} (e.g. 2x2_heis_neel_eps=1.0e-01). data_gen.py, alpha_gen.py, and sparsity.py must be run with the same --nx, --ny, --ham, --istate, and --eps to operate on the same directory.
└── 2x2_heis_neel
├── errs.npy
├── ests.npy
├── exacts.npy
├── info.json
└── shadows.npy
| File | Shape | Description |
|---|---|---|
exacts.npy |
(N, NUM_PAULIS) |
Exact expectation values over time |
ests.npy |
(N, NUM_PAULIS, NSNUM) |
Shadow estimates |
errs.npy |
(NUM_PAULIS, NSNUM, 8) |
Error information of shadow estimates |
shadows.npy |
(N, NSMAX, NQ, 2) |
Shadows generated from exact state evolution in QuTiP |
info.json |
Parameters file of physical simulation and data_gen.py |
data_gen.py takes the following arguments
usage: data_gen.py [-h] [--dir DIR] [--nx NX] [--ny NY] [--nb NB] [--ham {tfim,heis}]
[--istate ISTATE] [--n N] [--nsmin NSMIN] [--nsmax NSMAX] [--nsnum NSNUM] [--nw NW]
[--tplimit TPLIMIT] [--eps EPS] [--gamma GAMMA] [--seed SEED]
Command line arguments for CSST
options:
-h, --help show this help message and exit
--dir DIR directory to save data files
--nx NX Number of sites in x direction
--ny NY Number of sites in y direction
--nb NB Max Pauli observable weight
--ham {tfim,heis}
Hamiltonian type
--istate ISTATE Initial state ('ghz','neel','w','r','rp','hr','hrp', or length NQ string of
[0,1,+,-,>,<])
--n N Number of time steps
--nsmin NSMIN Min number of shadows
--nsmax NSMAX Max number of shadows
--nsnum NSNUM Number of shadows to try (logspaced)
--nw NW Number of workers for multiprocessing
--tplimit TPLIMIT Threadpool limit for each worker process
--eps EPS Max percentage of perturbation to Ham coeffs and Lindblad rates
--gamma GAMMA Common decay rate for all collapse operators
--seed SEED Random seed for reproducibility
The Hamiltonian options are
-
'tfim': Transverse-field Ising model$H = J \sum_{\langle a,b\rangle}Z_aZ_b + h \sum_a X_a$ -
'heis': Heisenberg model$H = J\sum_{\langle a,b\rangle}\left(X_aX_b+Y_aY_b+Z_aZ_b\right)$
The initial state options are:
-
'ghz':$\text{GHZ}$ state -
'neel': Néel state, alternating$|+\rangle$ /$|-\rangle$ by sublattice (bipartite 2-coloring of the site graph) — requires a bipartite graph, and is the default--istatefor all three scripts -
'w':$\text{W}$ state -
'r': state whose vector elements are uniformly random complex numbers -
'rp': tensor product of different single-qubit'r'states -
'hr': Haar random state -
'hrp': tensor product of different single-qubit'hr'states -
[0,1,+,-,>,<]: a string such as'0000...'or'>+<-...'indicating which of the 6 Bloch sphere states should be assigned to each qubit-
'0':$|0\rangle$ -
'1':$|1\rangle$ -
'+':$|+\rangle$ -
'-':$|-\rangle$ -
'>':$|+i\rangle$ -
'<':$|-i\rangle$
-
Run alpha_gen.py to create additional data files in the runs directory. Each run of alpha_gen.py will create a new pair of cs_xxxx.npy, params_xxxx.json files. The index xxxx will be the lowest integer unused by any existing pairs.
└── 2x2_heis_neel
├── errs.npy
├── ests.npy
├── exacts.npy
├── info.json
├── runs
│ ├── cs_0001.npy
│ └── params_0001.json
└── shadows.npy
| File | Shape | Description |
|---|---|---|
cs_xxxx.npy |
(NUM_PAULIS,NSNUM,MNUM,ALPHA_NUM,8) |
CSST error data |
params_xxxx.json |
Parameters file for alpha sweeps |
alpha_gen.py takes the following arguments
usage: alpha_gen.py [-h] [--dir DIR] [--nx NX] [--ny NY] [--ham {tfim,heis}] [--istate ISTATE] [--nw NW]
[--tplimit TPLIMIT] [--eps EPS]
[--mmin MMIN] [--mmax MMAX] [--mnum MNUM] [--amin AMIN] [--amax AMAX] [--anum ANUM]
[--snr SNR] [--rescale {0,1}] [--axis {0,1}] [--inverse {0,1}] [--fitint {0,1}]
[--replace {0,1}]
Command line arguments for CSST
options:
-h, --help show this help message and exit
--dir DIR directory to save data files
--nx NX Number of sites in x direction
--ny NY Number of sites in y direction
--ham {tfim,heis} Hamiltonian type
--istate ISTATE Initial state ('ghz','neel','w','r','rp','hr','hrp', or length NQ string of [0,1,+,-,>,<])
--nw NW Number of workers for multiprocessing
--tplimit TPLIMIT Threadpool limit for each worker process
--eps EPS Std dev for perturbation of Ham coeffs
--mmin MMIN Min m
--mmax MMAX Max m
--mnum MNUM Number of ms (logspaced)
--amin AMIN Min alpha power (10^amin)
--amax AMAX Max alpha power (10^amax)
--anum ANUM Number of alphas (logspaced)
--snr SNR Filter out (O_i, N_ST) with SNR < value
--rescale {0,1} Rescale by sqrt(N/m) in CS reconstruction
--axis {0,1} Axis to apply DCT/IDCT to
--inverse {0,1} Use IDCT instead of DCT for basis
--fitint {0,1} Fit intercept in Lasso regression
--replace {0,1} Sample with replacement in time-subsampling
The snr flag specifies a cutoff signal-to-noise ratio (default is cs_xxxx.npy can be a large file (GBs) it may contain majority unfilled entries/NaNs depending on the selected cutoff threshold. This file is read as a np.memmap in the plotting code due to its size.
Run sparsity.py to diagonalize the (perturbed) Liouvillian in the Pauli basis and, for each observable, compute the threshold-sparsity of the elementwise product of the initial state's and that observable's expansion coefficients in this eigenbasis — swept over multiple trials and a range of Hamiltonian perturbation strengths (eps). This writes sparsities.npz into the same {nx}x{ny}_{ham}_{istate} directory used by data_gen.py/alpha_gen.py.
| File | Shape | Description |
|---|---|---|
sparsities.npz |
epss: (EPSNUM,), sparsities: (TRIALS, EPSNUM, NUM_PAULIS) |
Threshold-sparsity of each Pauli observable's expansion coefficients, per trial and per eps |
sparsity.py takes the following arguments
usage: sparsity.py [-h] [--dir DIR] [--nx NX] [--ny NY] [--nb NB] [--ham {tfim,heis}] [--istate ISTATE]
[--trials TRIALS] [--eps EPS] [--epsmin EPSMIN] [--epsmax EPSMAX] [--epsnum EPSNUM]
[--gamma GAMMA] [--nw NW] [--tplimit TPLIMIT] [--seed SEED]
Command line arguments for CSST sparsity sweep
options:
-h, --help show this help message and exit
--dir DIR directory to save sparsity data
--nx NX Number of sites in x direction
--ny NY Number of sites in y direction
--nb NB Max Pauli observable weight
--ham {tfim,heis} Hamiltonian type
--istate ISTATE Initial state ('ghz','neel','w','r','rp','hr','hrp', or length NQ string of [0,1,+,-,>,<])
--trials TRIALS Number of trials per eps value
--eps EPS Reference eps value (kept in the eps sweep, used for labeling)
--epsmin EPSMIN Minimum eps in the sweep
--epsmax EPSMAX Maximum eps in the sweep
--epsnum EPSNUM Number of eps values in the sweep
--gamma GAMMA Common decay rate for all collapse operators
--nw NW Number of workers for multiprocessing
--tplimit TPLIMIT Threadpool limit for each worker process
--seed SEED Random seed for reproducibility
Note that --eps must equal one of the values in the linspace(epsmin, epsmax, epsnum) sweep (it marks the reference eps used elsewhere, e.g. in data_gen.slurm/alpha_gen.slurm); the script asserts this on startup.
All plots are generated by the data_analysis.ipynb notebook, including the eigenbasis sparsity violin plots (plot_sparsities) built from sparsities.npz. The data files can be found on Zenodo for the
Some of the data generating functions in data_gen.py, alpha_gen.py, and sparsity.py use Python's multiprocessing module to speed things up. By default 1 worker is requested, but the code will try to use min(requested cpus, available cpus-1). On a HPC cluster, the number of workers is set to the SLURM environment variable $SLURM_CPUS_PER_TASK. Each worker process additionally caps its internal BLAS/OpenMP threadpool via --tplimit (using threadpoolctl), so NW * TPLIMIT should not exceed the number of CPUs allocated to the job.
The slurm files data_gen.slurm, alpha_gen.slurm, and sparsity.slurm run data_gen.py, alpha_gen.py, and sparsity.py (respectively) in a HPC shell where SLURM commands are recognized. Each script sources slurm/config.sh (gitignored) for cluster-specific settings, so copy the template and fill in your own values:
cp slurm/config.sh.example slurm/config.shThis material is based upon work supported by, or in part by, the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number W911NF2310255.
@misc{barreto2026compressedsensingshadowtomography,
title={Compressed Sensing Shadow Tomography},
author={Joseph Barreto and Daniel Lidar},
year={2026},
eprint={2602.12518},
archivePrefix={arXiv},
primaryClass={quant-ph},
url={https://arxiv.org/abs/2602.12518},
}
@software{joeybarreto_2026_18626487,
author = {joeybarreto},
title = {USCqserver/CSST: v1.0.2},
month = feb,
year = 2026,
publisher = {Zenodo},
version = {v1.0.2},
doi = {10.5281/zenodo.18626487},
url = {https://doi.org/10.5281/zenodo.18626487},
}
@dataset{barreto_2026_18625394,
author = {Barreto, Joseph},
title = {Data repository for CSST},
month = feb,
year = 2026,
publisher = {Zenodo},
doi = {10.5281/zenodo.18625394},
url = {https://doi.org/10.5281/zenodo.18625394},
}