KrMM Recurrences: HPC-Enhanced Estimators for Nonlinear Classical Shadows

KrMM Recurrences: HPC-Enhanced Estimators for Nonlinear Classical Shadows

arXiv:2610.26472v1 [cs.DC] 19 Oct 2026

Matthias Gruber∗‡ , Sabrina Herbst∗‡ , Ivona Brandic∗ , Vincenzo De Maio†∗ ∗ Research Unit of Computational Sustainability

† School of Computing and Mathematical Sciences

TU Wien, Vienna, Austria

University of Leicester, Leicester, United Kingdom

Abstract—While quantum technologies are rapidly evolving, near-term applications heavily rely on hybrid quantum-classical paradigms. Among these, the classical shadows framework has gained a lot of traction, not only for quantum information science, but also for the development and verification of quantum hardware. It allows to efficiently classically estimate many properties of quantum states from randomized measurements alone, a task that is otherwise intractable due to the destructive nature of measurements in quantum systems. There are, however, many properties of interest that are, due to their very problem formulation, not efficiently estimable. In this work, we identify a computational kernel occurring in estimators of polynomial functions, a class of problems that are in general computationally difficult to solve. We show that this structure allows a rowwise partition that, together with a custom SIMD vectorization, enables massive parallelization and can reduce the computational cost significantly. We show applicability and scalability of the approach on different tasks, and our results showcase that the approach allows high sequential throughput and close-to-ideal parallel scaling. The work thus provides a new way of identifying and exploiting structure in classical shadow estimators, resulting in practical advantages on inherently difficult problems. Index Terms—hybrid quantum-classical computing, classical shadows, randomized measurements

I. I NTRODUCTION Hybrid quantum-classical paradigms have emerged as an important means to harness and improve the capabilities of current quantum hardware. An important application is the classical shadows formalism [1], which allows deriving properties of quantum states based on randomized measurements alone, otherwise difficult due to the very nature of destructive measurements. The protocol enables efficient unbiased estimation of many linear properties, such as expectation values of local observables, fidelity estimation of any state against a pure state, as well as entanglement verification for pure states. Nonlinear properties, such as the purity, can be estimated, however, the resulting estimator may become computationally expensive. In particular, polynomial quantities are estimated using the U-statistic formalism [1], i.e., we can estimate Tr(ρk ) by averaging over all size-k subsets of the N measurements. This results in a computational complexity of O(N k ), rendering the evaluation of higher-order moments increasingly expensive. Further, many interesting properties, such as the purity [2], have an exponential sample complexity N in the number of qubits n, making the overall estimation procedure ‡ These authors contributed equally to this work. Corresponding authors: {matthias.gruber,sabrina.herbst}@tuwien.ac.at

inefficient. Nonetheless, such estimators are of practical interest, e.g., for quantifying entanglement [2], thus motivating research for increased scalability. Storing the (estimated) density matrix of a quantum state ρ requires exponential-in-n memory, which can be avoided by utilizing the Kronecker structure imposed by classical shadows. While used extensively for linear properties, for higher-order terms the sheer amount of operations required becomes a major practical bottleneck [3]. To overcome the combinatorial explosion, batch shadows [4] were proposed, where the measurements are placed into distinct batches, and estimators are built by averaging over the snapshots of a batch. This yields an unbiased estimator with improved scalability, at the cost of storing the NB batch shadows as 2n × 2n density matrices. Compared to the vanilla shadows, not all combinations are considered anymore, thus more measurements may be required for the same accuracy. In this work, we identify a recurring computational kernel in polynomial estimators and formulate it as a chain of Kronecker-matrix-matrix (KrMM) recurrences. These recurrences require many small matrix products, allowing independent row-wise computations, with dependencies only between successive steps. This structure can be exploited for more efficient computation, e.g., through the recent FastKron algorithm [5], and a custom SIMD vectorization scheme. We implement our approach on estimators for mixed-state entanglement verification [2] and the Quantum Fisher Information estimation [6]–[8], as well as batch shadows. In summary, our contributions are as follows: (1) we identify a computational kernel underlying classical shadows estimators of polynomial functions; (2) we show that this formulation decomposes into independent row-wise calculations; and (3) we leverage this structure to massively parallelize the calculations, achieving near-ideal parallel scaling, and high throughput. Using state-of-the-art methods of parallel computing, we thus obtain a massive practical speed-up for a problem of increasing interest, that has not been parallelized on this scale so-far. Fig. 1 illustrates our approach. The paper is structured as follows. Sec. II provides background information and Sec. III discusses related work. Sec. IV introduces the KrMM recurrence and execution model. Sec. V discusses our applications and Sec. VI the experimental setup. Finally, we present our results in Sec. VII, discuss limitations in Sec. VIII and conclude the paper in Sec. IX.

$$C_t = \alpha_tA_tF_t + \beta_tC_{t-1}$$ $$q_l^{(k)} = \operatorname{Tr}\left( H\rho^l H \rho^{k+2 -l}\right)$$

Row-wise splitting induced by Kronecker structure

Quantum Fisher Information Common representation

Local Result

$$\{\mathcal{K}_i, \forall i\in R_0\}$$

\(\Gamma^{(0)}_T \)

\(\{\mathcal{K}_i,\forall i \in R_1\}\)

\(\Gamma^{(1)}_T \)

\(\{\mathcal{K}_i,\forall i\in R_p\}\)

\(\Gamma^{(p)}_T \)

CPU 0

$$c_t^{(i)} = \alpha_ta_t^{(i)}F_t + \beta_tc_{t-1}^{(i)}$$

Polynomial classical shadows estimators

Assigned Rows

CPU 1

$$F_t = \bigotimes_{j=1}^n G_{tj}$$

.

Aggregate & Reduce

. .

$$\hat{p}^{(k)} =\binom{N}{k}^{-1}\operatorname{Tr}\left(\sum_{1\le i_1 < \cdots < i_k \le N} \hat{\rho}_{i_1}^{T_B}\dots \hat{\rho}_{i_k}^{T_B}\right)$$

$$\mathcal{K}_i := c_0^{(i)} \overset{\bigotimes_j G_{0j}}{\rightar ow} c_1^{(i)} \overset{\bigotimes_j G_{1j}}{\rightar ow}\dots\overset{\bigotimes_j G_{Tj}}{\rightar ow} c_T^{(i)} $$

Mixed State Entanglement Verification

KrMM Recurrence

$

CPU p

Parallel Execution

Fig. 1: Our contribution. Polynomial classical shadows estimators, including online and batched [4] estimators for mixedstate entanglement verification and the Quantum Fisher Information, (see Sec. V), share a common representation, the KrMM recurrence, which induces a row-wise splitting by the Kronecker structure of Ft . This allows defining a computational kernel Ki , that iteratively computes the row estimators independently of each other, as each computation only depends on the previous timestep, making the computation highly parallelizable. Each worker locally aggregates the results, and communication is only required for the final aggregation and reduction. This results in high-throughput and near-ideal parallel scaling. While M−1 is not a valid quantum channel anymore, it is efficiently computable classically for U , allowing, for singlequbit Clifford unitaries, to obtain an unbiased estimator of ρ

II. BACKGROUND A. Classical Descriptions of Quantum States Near-term quantum computing is often interested in properties of quantum states. However, as quantum measurements are probabilistic and destructive, determining properties of arbitrary quantum states is non-trivial. Therefore, much research has been conducted on the classical reconstruction of quantum states, i.e., quantum state tomography [9]–[11]. In the most general case, it requires, however, exponential-in-n many copies of the state, accompanied by exponential memory complexity and inefficient data processing. Thus, a related line of research has considered whether, rather than reconstructing the full density matrix, it is sufficient to predict only certain properties of interest. This problem, called shadow tomography [12], aims to predict expectation values of a fixed set of observables for a quantum state. The circuits and measurements required by the protocol are, however, challenging to implement on near-term hardware, which motivated the development of classical shadows. B. Classical Shadows The classical shadows formalism [1] considers an unknown fixed state ρ, and an ensemble of unitaries U. The ensemble is commonly chosen to be tensor products of single-qubit Clifford unitaries, as this makes the protocol suitable for nearterm hardware. Next, randomized measurements are taken, i.e., the state is rotated according to a random U ∼ U , and measured in the computational basis. The unitary U and measurement result |b⟩ ∈ {0, 1}n are then stored classically. The average mapping of ρ to the classical description can be seen as a channel, which, when inverted, gives ρ as

ρ̂ =

n  O i=1

 3Ui† |bi ⟩⟨bi |Ui − I .

(2)

A classical shadow is then built as a collection of N independent snapshots. We provide an overview of the formalism in Fig. 2. The sample complexity N is generally problem dependent but can often be shown to be efficient and independent of the system size1 . This allows properties to be predicted via a median-of-means algorithm, which is often possible without reconstructing the whole density matrix. The approach naturally extends to polynomials via Ustatistics formalism, where we average over products of distinct snapshots. For quadratic terms, we would thus consider X 1 Tr(Oρ̂i ⊗ ρ̂j ). (3) ô = N (N − 1) i̸=j

To estimate polynomials of order k, the computational complexity grows as O(N k ). Thus, even for moderate k, the scaling is unfavorable, which led to batch shadows. C. Batch Shadows Batch shadows were originally proposed in the context of estimating operator entanglement [4], aiming to reduce the combinatorial explosion in the estimator by dividing the classical snapshots into NB distinct batches, where NB ≥ k. Denoting by Bi the set of snapshots in the i-th batch, we can then define its batch shadow as X ρ̄i = |Bi |(−1) ρ̂. (4) ρ̂∈Bi

†

−1

Eb,U [U |b⟩ ⟨b| U ] = M(ρ) ⇒ ρ = E[M

†

(U |b⟩ ⟨b| U )]. (1)

1 There are many subtleties involved here, with dependencies on U and the observables of interest. We refer the interested reader to the original source [1].

Snapshot

CPU

$$\times N$$

$$\hat{\rho} = \mathcal{M}^{-1}(U^\dagger\ket{b}\bra{b}U)$$

batch? Yes

No

$$\hat{\mu}_k = \binom{N}{k}^{(-1)}\sum_{1 \le i_1<\dots< i_k \le N} f(\hat{\rho}_{i_1},\dots,\hat{\rho}_{i_k})$$

$$B_i = \{\hat{\rho}_{iN/N_B},\dots,\hat{\rho}_{(i+1)N/N_B}\}$$

$ \mathcal{U}$

$$U\sim\mathcal{U}$$

$ \hat{\mu}_k$

$$\tilde{\rho}^{(1)} = |B_1|^{-1} \sum_{\hat\rho\in B_1} \hat\rho$$ $$(U, \ket{b})$$

$ \rho$

$$U\rho U^\dagger$$

$$\ket{b}$$

$$\hat{\mu}_k = \binom{N_B}{k}^{-1}\sum_{1 \le i_1<\dots< i_k \le N_B} f(\tilde{\rho}^{(i_1)},\dots,\tilde{\rho}^{(i_k)})$$

$ \vdots$

$$\tilde{\rho}^{(N_B)} = |B_{N_B}|^{-1} \sum_{\hat\rho\in B_{N_B}} \hat\rho$$

Randomized Measurements

Classical Postprocessing

Fig. 2: Overview of Plain and Batch Classical Shadows. On the left (orange), quantum measurements are collected by randomly sampling a unitary U from the ensemble U. The process is repeated N times, and the measurements {(Ui , |b⟩i )} are forwarded to the classical computer for post-processing. For the plain shadows (top), the estimator for the k-th order polynomial is then built by iterating over all size-k subsets of the snapshots, scaling as O(N k ). For batch shadows (bottom), the snapshots are split into NB ≥ k distinct batches. Batch shadows are built by averaging over all snapshots of the batch, and the k-th order polynomials are estimated by iterating over all size-k subset of batches. This reduces the number of combinations, at the expense of losing accuracy by not considering all cross products and storing the batch shadow density matrices in memory.

We can then obtain an unbiased estimator for the observable O(k) by averaging over all k-tuples of batches as ô =

 −1 NB k

X

Tr O(k)

1≤b1 <···<bk ≤NB

k O

ρ̄bi .

(5)

i=1

While this formulation reduces the combinatorial terms, it requires keeping the NB exponentially-sized matrices in memory. Further, as not all snapshot-level cross-products are considered anymore, more samples may be required. See [4] for the derivation and discussion on the resulting variance, and [3] for empirical results on sample complexity.

processes [14], [15], which motivated the FastKron algorithm [5]. It obtains advantages over existing approaches [16], [17] in terms of contiguous data access and memory operations on certain matrices (see Sec. IV-C), as well as communication overhead in multi-GPU settings. To the best of our knowledge, it has not yet found applications in quantum information. IV. K ERNEL F ORMULATION AND E XECUTION A. The KrMM Recurrence We define the Kronecker matrix-matrix (KrMM) recurrence with respect to computations over matrices as Ct = αt At Ft + βt Ct−1

III. R ELATED W ORK The computational complexity of shadows applications is primarily governed by the choice of ensemble U and the observables of interest. Many properties remain difficult to estimate, despite being of practical relevance, which has triggered research on pushing the scalability of these estimators. In particular, recently, a formulation of the estimators as online algorithms was proposed [3], where the estimators are updated on-the-fly instead of collecting all measurements first. This leads to a truly hybrid system, where both classical and quantum resources are used simultaneously, while also allowing to monitor convergence to employ early stopping. These practical improvements are particularly important for challenging problems such as operator and mixed-state entanglement [2], [4], the Quantum Fisher Information [6]–[8], or expectation values of multicopy observables [13]. In this work, we exploit the Kronecker structure of the problem. A matrix represented as n 2 × 2 Kronecker factors can be compactly stored using only O(log M ) memory and reduces the cost of dense matrix multiplication from O(M 3 ) to O(M 2 logM ), with M = 2n . This structure often occurs in quantum information, and is ubiquitous in classical shadows by design. However, it also appears in classical Gaussian

t = 1, . . . , T,

(6)

with input data αt , βt ∈ C, At , Ft ∈ CM ×M , and recursion matrix Ct ∈ CM ×M . In this work, we consider n-qubit systems, thus M = 2n . Ft has a Kronecker structure with complex-valued Kronecker factors Gtj ∈ C2×2 , satisfying Ft =

n O

Gtj .

(7)

j=1

Efficiently parallelizing the KrMM recurrence beyond individual matrix products is generally nontrivial, since each iteration depends on all previous ones through Ct−1 and possibly At . This is especially problematic if the matrix dimension M is relatively small but the number of iterations T is large, which is often the case in nonlinear classical shadow applications [18]–[20]. The Kronecker structure, however, allows a row-wise decomposition. Denoting the i-th row of a matrix B as b(i) , we obtain its result through the Kronecker matrix-vector (KrMV) recurrences (i)

(i)

(i)

ct = αt at Ft + βt ct−1 .

(8)

This execution model allows for an M -way parallelism with (i) (i) (i) row kernel Ki := c0 → c1 → . . . → cT , where each iteration depends only on the shared input data αt , βt , Ft , local

(i)

(i)

input data at , and the row-local state ct−1 . Therefore, Ki ⊥ Kj for i ̸= j, allowing to map large independent computations to each worker. Note that each row kernel needs to load the shared input data, thus, practical scalability implies being able to efficiently store Ft as O(log M ) Kronecker factors. B. Complexity and Memory Analysis Dense evaluation of Eq. (6) requires materializing the dense matrices and performing matrix-matrix multiplications, with a per-iteration cost of O(M 3 ). The dominant operation in (i) the KrMV recurrence is at Ft , which can be computed in O(M log M ) by exploiting Ft ’s Kronecker structure. Each of the row kernels performs T KrMV operations, for a total of M · T KrMV products and a complexity of O(T M 2 log M ). Each worker needs O(M ) memory for executing the row kernel Ki . If we are interested in Ct , we still need to store an O(M ) output of M row kernels, resulting in memory cost of O(M 2 ). However, often we are interested in a reduced version of CT , e.g., its trace, requiring to only keep the i-th entry of row i, thus a per-row memory cost of O(1). Thus, we avoid storing the dense matrix and instead only need the working memory for each row kernel, arriving at a reduced O(M ) total memory for sequential execution, while a fully parallel setting with p = M workers requires a dense O(M 2 ) memory cost distributed equally among all workers. C. Sliced-Row Multiply (SRM) Within each row kernel Ki , the dominant operations are the KrMV products. This puts our approach close to state-vector simulation methods [21], [22], where one is also concerned with the action of a Kronecker-structured matrix on a vector. The notable difference is that in state-vector simulation, one aims to distribute large KrMV operations over many workers, while here each worker computes many small KrMV operations in sequence, thus shifting the focus from single KrMV operations to efficiently fusing successive KrMV operations. Therefore, we adopt the Sliced-Row Multiply (SRM) operation defined in the recent FastKron algorithm [5]. Given a Kronecker factor Gtj ∈ C2×2 , SRM partitions the input row into slices of even and odd elements and independently applies Gtj . Subsequently, the complete KrMV product can be evaluated as a sequence of n SRM operations, without constructing the corresponding M × M Kronecker matrix Ft . Moreover, the SRM operation requires only contiguous data access in a highly regular, predictable pattern and is therefore well suited for fusing subsequent operations. This can be seen as a specialization of the approach employed in state-vector simulation methods. Using SRM, we lose the ability to exploit structure within the Kronecker factors, a crucial aspect of statevector simulation, but gain higher regularity and more effective fusion of successive KrMV operations. Note that since the shadow snapshots do not contain identity Kronecker factors, their structure could not be significantly exploited either way. The resulting KrMV evaluation serves as the basic building block of every row kernel. However, SRM cannot leverage

existing linear algebra routines, so we develop a custom SRM implementation with a SIMD vectorization scheme. D. Execution Model A conventional evaluation requires computing each iteration of Eq. (8) independently, therefore risking materializing working data to memory. Thus, we treat each row kernel Ki as a single fused computation, illustrated in Fig. 3. We instantiate vector registers using only consecutive vector loads and perform as many computations on them as possible before writing them back. The state produced is used immediately in the next cycle, avoiding intermediate memory access and reducing memory traffic. This provides contiguous memory accesses, leading to a cache-aware implementation. When computing a set of coupled KrMM recurrences, the state produced by the update cycle is directly used as the input row of the next recursion cycle. Also, the recurrence fusion is independent of the Kronecker acceleration: while the latter reduces the cost of KrMV, the fusion reduces the cost of moving data through the recurrence. The row-wise decomposition naturally exposes parallelism as row kernels share no state, enabling independent distribution across the p workers. Let Rj be the set of rows assigned to worker j. Each worker j evaluates all kernels Ki , i ∈ Rj , and aggregates their results. Since no communication is required during the kernel execution, a single global reduction combines the local results. Further, as all row kernels perform the same operations, a static row partitioning strategy naturally provides load balancing. The absence of communication during kernel execution allows near-ideal parallel scalability until bandwidth limitations become dominant. Whether bandwidth contention limits scalability hinges on the rows fitting into cache, which is usually the case for the system sizes considered in nonlinear classical shadows (n ≲ 12) [18]–[20]. E. Data Layout and SIMD Realization The SRM operation applies the same transformation to each row, therefore the most natural SIMD strategy is vectorizing across rows. However, our execution model assigns each row to different workers and exploits row-level parallelism, making it incompatible with cross-row vectorization. Thus, we implement a column-based SIMD strategy on a custom data layout using architecture-specific permute intrinsics. One SRM cycle with Kronecker factor   a c G := ∈ C2×2 (9) b d G

transforms the row buffers of Fig. 3 as src − → dest. In the following we denote by xl , yl ∈ C, l < M/2 even and odd entries of src in the natural data layout of interleaved real and imaginary parts. Every result element of dest is obtained by summing up an even-odd pair xl , yl , specifically dest[l] = xl · a + yl · b

dest[M/2 + l] = xl · c + yl · d.

(10)

Kronecker factor $$G_{tj}$$

Row state

$$a,b,c,d\in\mathbb{C}$$

$$c^{(i)}$$

Start $$j=n$$ $$t=1$$

broadcast vector registers

$$[x_0^r,y_0^r][x_0^i,y_0^i]$$ $$[x_1^r,y_1^r][x_1^i,y_1^i]$$

$$a_t^{(i)}$$

load $$j=n$$ $$t++$$

shuffle_even 2x $ \mathbf{x}^r$ $ \mathbf{x}^i$ $ \mathbf{y}^r$ $ \mathbf{y}^i$

load

FMA 4x4 $$\mathbf x^r a^r + \mathbf y^r b^r - \mathbf x^i a^i-\mathbf y^i b^i$$ $$\mathbf{x}^r a^i + \mathbf y^r b^i + \mathbf x^i a^r+\mathbf y^i b^r$$ $$\mathbf x^r c^r + \mathbf y^r d^r - \mathbf x^i c^i - \mathbf y^i d^i$$ $$\mathbf{x}^r c^i + \mathbf y^r d^i + \mathbf x^i c^r+\mathbf y^i d^r$$

shuffle_odd 2x

write load

if $$j=1$$

$$\times\beta_t$$

$$\times\alpha_t$$

$$+$$

Aggregate

if $$j > 1$$ $$j--$$

if $$j>1$$ if $$j=1$$

$${src}$$ $${dest}$$ swap after full SRM

Recursion cycle $$\mathcal O(T)$$

write

SRM cycle $$\mathcal O(nMT)$$

Update cycle $$\mathcal O(MT)$$

Fig. 3: Fused row kernel Ki with custom SIMD data layout for a vector size Vs of 2 doubles. Performing one full SRM operation (index j) requires fully covering one row with M/(2Vs ) repetitions of the SRM and update cycle. r indices denote real part, and i indices imaginary part. Green fields correspond to small data that will likely be accessed from cache. Blue fields correspond to large data that will likely be loaded from memory. In the online algorithm, each recursion cycle repeats m times, with the output of the previous update cycle directly serving as the input. Aggregation in the online algorithm corresponds to storing the i-th entry of the row state. In the batch algorithm it corresponds to storing the entire row state in the last iteration. The combined action of shuffle + FMA preserve our data layout, enabling efficient chaining of successive SRM operations. Efficient vectorization requires extracting xl and yl into separate SIMD registers using only consecutive vector loads, while also supporting complex FMA arithmetic. Furthermore, SRM operations are evaluated successively, so the output of one operation must serve as input to the next Kronecker factor without additional reorganization. Thus, we use a data layout where we group Vs real and Vs imaginary floats together, visible in Fig. 3 for a vector size of Vs = 2. This allows extracting the data using shuffle_even and shuffle_odd primitives, which, given two SIMD vectors, fill a new vector with all their even or odd entries. These operations are implemented using architecture-specific permute intrinsics, enabling vectorized complex FMA operations while preserving the layout. The kernel performs a 4-way pattern of 1 load, 1 shuffle, 4 FMA, and 1 write, where the efficiency depends on the architecture support for intrinsics permuting arbitrary elements of two input vectors to one output vector. In practice, the CPU can efficiently hide the cost if native support is given, as in e.g. AVX-5122 . Relying on a composition of multiple intrinsics, e.g., for AVX23 , leads to a noticeable performance loss. The data layout is extendable to all power-of-2 vector sizes VS as long as it fits into one row at least two times, i.e. 2VS ≤ M . V. A PPLICATIONS Beyond describing the formalism for general polynomial functions, we showcase our approach for two specific estimators: mixed state entanglement verification and Quantum 2 For AVX-512 we use the SIMD intrinsic _mm512_permutex2var_pd. 3 For AVX2 we use two _mm256_permute4x64_pd followed by one _mm256_blend_pd.

Fisher Information estimation, which represent practically relevant problems (e.g., quantum-state verification [23] and quantum metrology [24]) and allow for online and batch estimators that can be expressed via the KrMM recurrence. A. Mixed State Entanglement Verification Entanglement of mixed states can be verified using the positive partial transpose [23], [25], [26] condition. In a nutshell, for a state ρ on partitions A and B, negative eigenvalues of the partially transposed density matrix suffice to establish that A and B are entangled, however, this is computationally difficult as it requires access to the full density matrix. An alternative criterion can be found via the moments of the partial transpose, which can be estimated using classical shadows [2], [18], but become increasingly computationally difficult with higherorder moments. We apply the KrMM recurrence on two estimators; the online estimator from [3] and batch shadows. a) Online Algorithm: The k-th moment at shot t ∈ [1, . . . , N ] is defined as  −1 t (k) p̂t = Tr(Ck (t)), (11) k with the matrices Ck (t) as the combinatorial sum of all t snapshots with partial transposition on subsystem B X Ck (t) = ρ̂Ti1B · · · ρ̂TikB . (12) 1≤i1 <···<ik ≤t

It has been shown that the matrices Ck can be computed using a recurrence of Kronecker matrix products [3], which is a natural fit for our KrMM recurrence formalism. We slightly adapt the original formulation, by integrating the normalization

into the running estimate for numerical stability and altering the recursive structure for better kernel fusion, obtaining B Ck (t) = αtk Ck−1 (t)ρ̂Tt+k−1 + βtk Ck (t − 1),

(13)

k t−1 with αtk = t+k−1 and βtk = t+k−1 for k = 1, . . . , m, and initial conditions Ck (0) = 0 and C0 (t) = I. The running moment estimate of all orders k can then obtained in every iteration as p̂kt+k−1 = Tr(Ck (t)). Within each t iteration, the output of iteration k directly serves as the input to k+1, which is exploited in the row kernel to reduce memory traffic. Cheap access to the recursively built intermediate values allows uncertainty estimation via standard formula of the U-statistics variance [27], which we briefly derive in the following. In the first order of t it is given as

k 2 σ12 , (14) Var(p̂kt ) = t which dominates all other terms scaling as t−x with x ≥ 2 for large enough t. The variance of the conditional expectation of the U-statistics kernel, given one input is fixed, is defined as  h i σ12 = Var E Tr(ρ̂Ti1B · · · ρ̂TikB )|ρ̂TikB , (15) which is a by-product of Eq. (13). Subtracting Ck (t − 1) from both sides of the recurrence and taking the trace, together with the definition of the matrices Ck (t), gives

t−1 t h i TB TB TB k k k k

· · · ρ̂ )|ρ̂ p̂ − p̂ = E Tr(ρ̂ . (16) t t−1 ik ik i1 t−1 t−1 k−1

k−1

Computing many of these conditional expectation values at consecutive t allows computing their sample variance, and thus σ12 and Var(p̂kt ). This in turn can be used to check for convergence at a desired accuracy or propagate uncertainty to derived quantities. Extracting all intermediate moment estimates (and uncertainty) costs O(mN 2n ), leaving the total complexity of O(mnN 4n ) unchanged. As a result, the online algorithm allows for a truly hybrid approach, where incoming measurements are immediately processed at O(mn4n ) cost, providing live insight into the quantum system. b) Batch Algorithm: Computing a batch shadow ρ̄i corresponds to averaging over a set of snapshots ρ̂t ∈ Bi . This can be formulated as a KrMM recurrence using a fixed At = I and a running average, giving the defining recurrence as C(t) = αt Iρ̂Tt B + βt C(t − 1),

(17)

with αt = 1t , βt = t−1 t , ρ̂t ∈ Bi , and initial condition C(0) = 0. For NB = k independent batches coveringQall meak surements, this gives the final estimate as p̂kN = Tr( i=1 ρ̄i ). The cost of forming the batch shadows thus scales as O(nN 4n ), whereas their multiplication scales as O(m8n ). However, as N ≫ 2n , the dominant cost in practice is the formation of the batches, thus, we only build the batches once and estimate all m moments up from them. For this we initially build m batches and compute p̂m . We then combine batches m and m − 1, and compute the estimate of p̂m−1 . We repeat this until two batches are left for a total cost of O(m2 8n ) and

combine the batches evenly to retain the most combinations possible. While the KrMM recurrence serves to parallelize the computation of a single batch, the computation can be further parallelized by forming the batches concurrently. With a total cost of O(nN 4n ), the batched algorithm provides a benefit of O(m) over the online algorithm in an offline setting. While batch formation is easily implemented online, extracting moment estimates comes at a cost of O(m2 8n ). For the typically small n used in classical shadows, paying this once is negligible, but in an online setting, this would scale as O(m2 N 8n ), which makes a comparison between the approaches difficult, as they cover different problem domains. B. Quantum Fisher Information The Quantum Fisher Information (QFI) [24], [28] was originally proposed in quantum metrology, however, has since found applications in P quantum information [29]. Consider a quantum state ρ = i λi |i⟩ ⟨i| and a Hermitian operator H. The QFI is then defined as X (λi − λj )2 | ⟨i| H |j⟩ |2 . (18) FQ = 2 λi + λ j (i,j)λi +λj >0

While it cannot be directly expressed as a polynomial function, polynomial lower bounds can be estimated via classical shadows [6], however, the required number of bounds make it prohibitive in practice. Thus, an alternative framing based on the Krylov subspace method was proposed recently as Krylov Shadow Tomography (KST) [8], [20], [30]. Here, increasingly accurate nonpolynomial lower bounds of the QFI are obtained from multiple polynomial quantities. Intriguingly, for a state ρ and operator H, an analytical maximum bound order N ∗ exists, for which the estimate and FQ coincide. Obtaining an unbiased estimate for FQ then reduces to estimating a set of polynomials (k)

ql

= Tr(Hρl Hρk+2−l ),

(19)

for k = 0, . . . , 2N ∗ − 1 and l = 0, . . . , k + 2. This can be computed with batch shadows by inserting the computed batches into Eq. (19), similar to previous works [8], [20]. We can thus directly apply the batch algorithm from Eq. (17). To obtain a fine-grained estimation of the QFI, we can also employ the online algorithm to estimate the polynomials (k) ql . For this, we use N ∗ + 1 independent runs of order (k) m = 2N ∗ + 1, each estimating ql for a fixed value of ∗ l = 0, . . . , N and all k = 1, . . . , m. The remaining values can be obtained by the cyclicity of the trace. To consider the two (k) Hermitians H present in the definition of ql , we adapt the online recurrence from Eq. (13). The leftmost H in Eq. (19) enters as the initial condition C0 (t) = H. For the second H, we split the recurrence into two separate cases as Ck (t) = αtk Ck−1 (t)ρ̂t+k−1 + βtk Ck (t − 1) , k ̸= l + 1 (20) Ck (t) = αtk Ck−1 (t)H ρ̂t+k−1 +βtk Ck (t−1) , k = l+1. (21)

We thus require H as (a sum of) Kronecker products, which is the case in many problems of interest. Here we consider H =

1 2

(j) j=1 σZ

Pn

TABLE I: Hardware configurations.

(j)

, where σZ is the Pauli-Z operator acting on qubit j. This Hamiltonian is useful for entanglement detection and is the same as in the main works on KST [8], [20], [30]. VI. E XPERIMENTAL S ETUP a) Baseline: We compare our implementation to the state-of-the-art Julia package RandomMeas.jl [31], which supports batch shadows out of the box. It partitions all measurements as N = NU ·NM , averaging over NM measurement outcomes grouped under one of the NU unitaries by computing the empirical Born distribution. This changes the dominant cost factor from N to NU . To obtain the full randomized measurement statistics, we use all possible unique unitaries NU = 3n with NM = N/NU outcomes each, and distribute the unitaries equally among all batches. To obtain a fair comparison, we compute all batches in parallel. We extract the moment estimates from the computed batches using the same scheme as in our batch algorithm. b) Quantum States: We use three different states with fixed experimental parameters. Unless mentioned otherwise, results provided for a given state refers to the parameters listed here. We obtain the classical shadows data and ground truth values by numerically simulating the data acquisition phase. Werner states [32] are a family of states, where the entanglement is controlled by a parameter t. We collect N = 6 · 108 measurements for n = 2, 4, 6, 8, 10 from a state with t = 0.95, requiring moments up to order m = 34 to detect entanglement [3]. These parameters were chosen to be experimentally feasible, while demonstrating the scale of our approach; for n = 10, applying the online algorithm costs about 3 Exaflops in total. We use this state for entanglement verification, with the partial transpose on an equal split. To show distributed scalability, we employ 2 and 8 nodes for n = 8 and 10 respectively for the online algorithm. We do not distribute the batch one, as it already outperforms the online algorithm on a single node and relies on the same implementation, thus, would have provided only limited additional insights. The considered quantum state itself does not affect runtime in any way. To analyze performance metrics, we use 1 I), as measurement data can the maximally mixed state ( M be simulated easily. We use 11 independent samples with N = 107 measurements each for n = 2, . . . , 8 qubits with order m = 34 to obtain a meaningful runtime estimate, that is not dominated by start-up costs. For higher n = 9, . . . , 12, we decrease the number of measurements and use 11 samples with N = 5 · 104 each. Results presented using this state correspond to the median out of all 11 samples. The random rank 2, 4, 6 states refer to randomly generated states with fixed rank, used for QFI estimation. We use 11 independent samples with N = 3 · 106 each for n = 2, . . . , 8 qubits. Results presented using this state correspond to the median out of all 11 samples. c) Hardware and Software: We benchmark our approach on three different AMD EPYC processors with the most recent microarchitecture generations Zen 3, 4, and 5, which we compare in Tab. I. The Zen 3 configuration corresponds to

EPYC 7713 EPYC 9354P EPYC 9255 Generation Zen 3 [33] Zen 4 [34] Zen 5 [35] Cores 64 32 24 Threads 128 64 48 Sockets 2 1 1 Base 2.00 GHz 3.25 GHz 3.20 GHz Boost 3.72 GHz 3.80 GHz 4.32 GHz L1 (i/d) 32/32 KB 32/32 KB 32/48 KB L2 512 KB 1 MB 1 MB Memory 512 GB DDR4 503 GB DDR5 755 GB DDR5 Max SIMD ISA AVX2 AVX-512 AVX-512 Vector width 32B 32B∗ 64B FMA pipes 2 2 2 Compiler AOCC 4.1.0 GCC 15.3.0 GCC 15.3.0 MPI OpenMPI 4.1.4 — — ∗ Zen 4 implements AVX-512 via double-pumped 256-bit execution units. 256

Cores Vectorsize

128

2x EPYC 7713 (Zen3)

64 EPYC 9354P (Zen4)

32 EPYC 9255 (Zen5)

16

AVX-512 (Zen4, Zen5)

8

AVX2 (Zen3)

4 SSE4.2

2 1

2

3

4

5

6

7

8

9

10

11

12

#Qubits n

Fig. 4: Used single node parallelism over increasing number of qubits. We use simultaneous multithreading (SMT) once the number of rows exceeds the number of available cores.

a single node of an institutional HPC cluster. For multi-node experiments, we use nodes interconnected with a Mellanox HDR InfiniBand fabric (200 Gbit/s). We compile using the C++17 standard with optimization flags -O3, -march=native, -ffast-math, -ffp-contract=fast, and use OpenMP [36] for shared memory parallelism. To implement the custom SIMD vectorization scheme, we use vector intrinsics provided by the Vector Class Library (VCL), version 2.02.01 [37]. In all experiments, we use double precision floating point numbers. The number of qubits n heavily dictates the parallelism in our partitioning scheme. Unless stated otherwise, for each n we use the maximum core count and vector size. Fig. 4 illustrates this in a single node setting. Our source code and experimental data is accessible at [38]. VII. R ESULTS A. PT-moments estimation a) Accuracy and Convergence: To test our approach for mixed state entanglement verification, we predict PT moments of the Werner state and analyze accuracy by investigating the relative error. Fig. 5 shows convergence over the number of measurements for n = 6, 8, 10 qubits for moment orders 2 and 34 for both algorithms with the KrMM recurrence. The online estimator allows for fine-grained convergence resolution, while the batched estimator only provides estimates

102

102

1

1

Relative Error

10

10

n=6

100

10−1

10−2

10−2

10

−3

10

−3

10

−4

10

−4

106

107

108

101

n=8

100

10−1

105

102

100 10−1 10−2

k=34 Online k=34 Batch k=2 Online k=2 Batch √ O(1/ N )

109

105

Number of Measurements N

106

10−3

n = 10

10−4 107

108

109

105

106

Number of Measurements N

107

108

109

Number of Measurements N

Fig. 5: Relative error of PT moment estimation over number of measurements of the Werner state for moments 2 and 34 for 6, 8, and 10 qubits. Lightened lines and dots depict the relative errors. Thick lines correspond to the probabilistic error bound of the online algorithm obtained through variance estimation.

109

10

10

7

104 10

Runtime [s]

Number of Measurements N

105 1d 8

106

10

N =6 · 108

1h

3

102 1min

O(n2n )

NU =3n N =108

101

5

100 n=10 n=6

104

n=8 n=4 n=2

O(k 2 )

5

10

15

20

25

30

35

Moment Order k

Fig. 6: Number of measurements N needed to reach convergence using the online variance estimation applied to PT moment estimation of the Werner state. We declare convergence once we are 95% confident that the relative error < 10−1 . Lightened lines correspond to N needed for the true relative error < 10−1 for the current, and all following data points.

at discrete evaluation points. Beyond stochastic noise, the relative error of the two algorithms are similar, which coincides with the √ analysis of [4]. For large enough N , they converge as 1/ N , agreeing with the results of [2]. In a real experimental setting, where the true moment values would be unknown, only the probabilistic error bounds (thick lines) are accessible, which only become useful after a critical number of measurements. Before that, the confidence interval includes 0, making the upper error bound diverge. This is not an inherent flaw of the method but simply means that it is not possible to even decide on the sign of the estimate with confidence. This accuracy estimate makes the online algorithm attractive in practice, where quantum measurements arrive one after another, and we wish to stop at a certain accuracy4 . For a more comprehensive convergence analysis, we investigate the number of measurements needed to declare convergence with 95% confidence for an accuracy of 10−1 in Fig. 6. Since we estimate the accuracy with a probabilistic upper bound, the number of measurements to declare convergence 4 Using the Jackknife procedure [4], uncertainty estimates are also available for batch shadows, however, they are not implemented here.

N =107

1s

Online RandomMeas.jl Batch

N =106 10−1 2

4

6

8

10

#Qubits n

Fig. 7: Runtimes of different algorithms applied to the Werner state. For the online algorithm we used 2 nodes at n = 8 and 8 nodes at n = 10. RandomMeas.jl and the batch algorithm were performed on a single node. For RandomMeas.jl we used NM = 107 /NU . Performed on the Zen 3 configuration. is naturally higher than what would be needed to actually reach convergence. The measurements scale as O(k 2 ), which agrees with Eq. (14). The exponentially many measurements N needed to reach convergence make nonlinear classical shadows inherently limited to small quantum system sizes n, since N iterations of the KrMM recurrence are required. Note, however, that this highly depends on the state and observable. b) Performance Benchmark: We compare the runtime of the two algorithms and our baseline RandomMeas.jl in Fig. 7. As the online algorithm essentially corresponds to running the batch algorithm m times, it is generally slower. The core idea of the online algorithm, however, is concurrent execution with the data acquisition itself, limiting the significance of a direct performance comparison. Both the online and batch algorithms initially show a scaling law of O(n2n ), as the complexity is O(n4n /p), with p = 2n . For a single node, we run out of cores to add for n > 7, thus, observe a transition to O(n4n ) scaling. For sufficiently large N , the sampled unitaries eventually correspond to the 3n possible measurements. In this regime, our batch implementation scales linearly with the total number of measurements NU · NM , while our baseline has a fixed cost of NU = 3n . Consequently, for sufficiently large N ,

3.00

Zen 5 Zen 4 Zen 3

2.50

Tauto /TSIMD

0.8

Zen 5 Zen 4 Zen 3

2.75

0.6 Boost

1.0

0.4

Frequency

Per-Core HW Utilization

1.0

0.2

2.00 1.75

Base Zen4 0.8

2.25

Base Zen5

1.50 0.6

Base Zen3 2

4

6

8

10

1.25

12

#Qubits n

2

4

6

8

#Qubits n

10

12

2

4

6

8

10

12

#Qubits n

Fig. 8: (left) (brown) Strong-scaling efficiency of the distributed online algorithm for 12 qubits using the maximally mixed state. To have a sufficiently large workload, we use N = 105 measurements. (pink) Speedup of m times the batch KrMM kernel to the online kernel using the maximally mixed state. Performed on the Zen 3 configuration. (middle) Per-core hardware utilization (FMA) of the online kernel for different microarchitectures (for number of cores, cf. Fig. 4) using the maximally mixed state (single node). The CPU frequency was measured at runtime to faithfully estimate peak performance at each core count. At n = 12 (all cores) we achieve throughputs of 2.28, 1.33, 1.26 TFLOPs for Zen 3, 4, 5 respectively. (right) Speedup of the online kernel with manual SIMD vectorization compared to the auto vectorized one using the maximally mixed state.

RandomMeas.jl becomes asymptotically more favorable than repeatedly evaluating the increasingly large KrMM recurrence. However, this asymptotic regime lies beyond realistic measurement counts, which is limited by the cost of data acquisition [1], [4], [39]. Nevertheless, in the regime accessible on today’s quantum hardware (N ≈ 106 −108 ), the parallelism exposed by the KrMM recurrence yields speed-ups ranging from 1.66× (n = 8, N = 108 ) to 516× (n = 10, N = 106 ). Fig. 8 (left, brown) shows the strong-scaling efficiency of the distributed online algorithm. Our row-wise partition ensures that inter-node communication is limited to a final reduction. Therefore, we observe close-to-ideal strong scaling efficiency; above 0.95 for all multi-node executions with up to 16 nodes. At 2 nodes, we even observe a minor superlinear speed-up, which we attribute to a hardware artifact of the HPC cluster. Note that we only considered power-of-2 node counts P, as it maps perfectly to the row count. Otherwise, we expect slightly reduced efficiency due to non-ideal load balancing. We compare the performance of the single-node online and batch algorithms further in Fig. 8 (left), where the pink line compares the runtime of the online algorithm to m times the batched one. For few qubits, the online algorithm is less efficient than the batched one, likely as a result of higher implementation overhead as, despite efficiently fusing the data access, its implementation still involves much bookkeeping. For small n, the kernel has a tight execution schedule, and the CPU may struggle amortizing these costs. For larger systems the online algorithm becomes slightly more efficient. Having established parallel scalability in a multi-node setting, we turn to single node performance for increasing n and cores, comparing throughput relative to peak performance. Since FMA operations constitute the dominant workload, we take peak FMA throughput as the reference performance freq(n)·#cores(n)·vectorsize(max)·#pipes(FMA)·2, (22)

where the factor 2 comes from one FMA operation performing two flops. To obtain a faithful estimate at varying core counts, we measure the CPU frequency during kernel execution using the perf profiling tool. Thus, at n = 12 (full core usage) we arrive at peak performances of 3.96, 1.77, 2.18 TFLOPs for Zen 3, 4, 5 respectively. The performed flops in the dominant FMA section of the row kernel are given as 14 · N · m · n · 4n , where 14 are the flops to compute α · x + β · y. Fig. 8 (middle) shows the resulting per-core hardware utilization of the online KrMM kernel for the different microarchitectures. Even though more cores are added for increasing n, the per-core hardware utilization does not decrease. Just as in the multi-node setting, all cores work independently, leading to minor parallel overhead. Crucially for the single-node setting, at the considered values of n, the data buffers are small enough to fit into either L1 or L2 cache, leading to no heavy bandwidth contention between threads. Due to the compact Kronecker factor representation, loading ρ̂t constitutes only minor memory traffic per cycle. Both the trivially parallel partition and the high data locality of the row-wise approach thus enable high throughput. Initially, at small n, utilization is poorer, which continuously rises as n increases. The reasons for this are two-fold: (1) we are limited to smaller vector sizes at n = 2, 3, and (2), at small n the SRM cycle is not repeated many times, leading to implementation overhead dominating total cost. From n ≈ 6 qubits onward, we plateau at around 57% for Zen 3, and 75% for Zen 4 and 5. The performance of Zen 5 does, however, decrease to roughly 57% starting at n > 9 qubits. We attribute this performance drop to a transition from L1 to L2 cache, since this coincides exactly with the row buffers leaking L1. Interestingly, this transition manifests neither in the other architectures nor in the auto-vectorized version on Zen 5. Lastly, we show the speed-up of our custom SIMD vectorization to an auto-vectorized implementation in Fig. 8

104

100

10−1

1d

1h

103

Runtime [s]

101

Relative Error

105

n = 8 Online n = 8 Batch n = 4 Online n = 4 Batch

102 10

1min

1

100

1s rank 6 Online rank 4 Online rank 2 Online Batch

10−1 100

101

102

103

104

105

106

Number of Measurements N

3

4

5

6

7

8

#Qubits n

Fig. 9: Relative error of estimating the QFI of the random rank-6 state with n = 4 and n = 8 qubits. Solid lines depict the online algorithm; dots correspond to the batch algorithm.

Fig. 10: Runtimes of QFI using the random rank 2,4,6 states. Batch uses rank 6. For rank 2 we have N ∗ = 3 and for rank 4 N ∗ = 10. Performed on the Zen 4 configuration.

(right). We implemented the auto-vectorized baseline using the same row kernel execution model of Fig. 3, but used a naive data layout of interleaved real and imaginary floats, and performed the arithmetic according to Eq. (10). We can see that our custom vectorization consistently outperforms the auto-vectorized one. It is particularly well suited for AVX-512 architectures (Zen 4, Zen 5), resulting in a rapid increase in speed-up at n = 4 qubits. The lower speedup of Zen 3, coupled with the lower saturated hardware utilization of 57% implies that the hardware support for SIMD intrinsics shuffling data from arbitrary entries of two vectors to a new vector (as in AVX-512) is crucial for the performance of our SIMD scheme.

independent runs to estimate all polynomials, so it has a clear runtime disadvantage for higher rank states.

B. Quantum Fisher information a) Accuracy and Convergence: In KST, the rank of the state heavily dictates the effort needed to estimate the QFI. This is formulated in the bound order N ∗ , at which the estimated lower bound and the QFI coincide exactly [8]. In our experiments, we employ a random rank-6 state, where N ∗ = 21. To estimate FQ directly, we need to estimate the polynomials up to order 43, according to Eq. (19). Fig. 9 shows the relative error of estimating FQ of a random rank-6 state with n = 4, 8 qubits for the algorithms. They produce a similar error up to noisy uncertainty, both requiring around 3 · 106 measurements for a relative error of 10−1 . While the estimator for 8 qubits is clearly less accurate than for 4 qubits, the difference is less pronounced than for the PT moments. For the QFI, we refer to the true relative error as it is not as straightforward to estimate uncertainty. While (k) we have access to the variance of ql , propagation to FQ is unclear. The non-polynomial approach lies at the core of KST, however, it also makes uncertainty estimation non-trivial and beyond the scope of this paper. Previously [20], this was addressed by repeating batch shadows 100 times in a bootstrapping scheme over all N measurements. b) Performance Benchmark: Fig. 10 compares the runtime using 3·106 measurements for different system sizes. The batched version shows clear advantages, requiring less than 1 minute even on 8 qubits. The online algorithm requires N ∗ +1

VIII. L IMITATIONS While our approach provides a practical improvement, it does not eliminate the exponential memory and runtime scaling, which is dictated by the exponential sample complexity. Nevertheless, the considered problems are highly relevant, especially in the small- to intermediate-scale systems regime for which classical shadows were originally designed. Our method, together with the online approach and early stopping from [3], extends the range of feasible problem sizes and significantly improves the efficiency of the estimators within the system sizes relevant to current classical shadows applications. IX. C ONCLUSION AND F UTURE W ORK In this paper, we identified a recurring pattern in classical shadow estimators, allowing for significant practical speed-ups on practically relevant but computationally difficult problems. We formulated the Kronecker-matrix-matrix recurrence, resulting in independent row-wise partitions that allow parallelizing the execution. In combination with a custom SIMD vectorization, it enables optimizing the throughput and parallel scaling. Further, we have shown that the structure naturally occurs in estimating polynomial functions in classical shadows applications, with the examples of mixed state entanglement verification and the estimation of the Quantum Fisher Information. We defined online estimators that continuously monitor convergence, allowing for advanced hybrid workflows. In this approach, no dense matrix ever needs to be stored concurrently, leading to reduced memory requirements. More broadly, our framework can be applied to batch shadows, the state-of-the-art estimator for polynomial functions in an offline setting. This enables parallelism beyond the conventional concurrent execution of batch computation, translating into significant speed-ups over the reference implementation. With the establishment of online algorithms for classical shadows it would be an interesting future avenue to develop a pipeline for an actual online setting. We have not considered

this in this work due to (1) limited availability of quantum hardware and (2) limited support to-date from hardware providers for online hybrid algorithms. Further, in light of noisy quantum measurements, it would be interesting to develop mitigation strategies, e.g., through robust shadows [40]. ACKNOWLEDGMENTS For this work, we have used AI to improve reading flow and clarity, however, not for text generation. Further, we used it to aid with visualization. We have reviewed all AI-assisted work by manually verifying correctness and take full responsibility for the final content of this work. The computational results have been achieved using the Austrian Scientific Computing (ASC) infrastructure. SH is a recipient of the Andreas Dieberger-Peter Skalicky scholarship awarded by TU Wien and acknowledges financial support by the netidee Austria funding campaign (Grant no. 7728). This work has been partially funded through the Themis project (Trustworthy and Sustainable Code Offloading), Austrian Science Fund (FWF) (Grant no. 10.55776/PAT1668223), the Standalone Project Transprecise Edge Computing (Triton), Austrian Science Fund (FWF) (Grant no. P 36870-N), BMIMI, BMWET, the State of Upper Austria, and the State of Tyrol within the COMET module Quantum Algorithm Engineering (FFG grant no. 923923) managed by Austrian Research Promotion Agency FFG, and by the Advanced learning paradigms for Photonic quantum processors (ALPAQA) project, funded by the Vienna Science and Technology Fund (WWTF) (Grant no. 10.47379/ICT25062). R EFERENCES [1] H.-Y. Huang, R. Kueng, and J. Preskill, “Predicting many properties of a quantum system from very few measurements,” Nature Physics, vol. 16, pp. 1050–1057, 2020. [Online]. Available: https: //www.doi.org/10.1038/s41567-020-0932-7 [2] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, “Mixed-state entanglement from local randomized measurements,” Phys. Rev. Lett., vol. 125, p. 200501, 2020. [Online]. Available: https://doi.org/10.1103/PhysRevLett.125.200501 [3] M. Marso, S. Herbst, J. Wilkens, V. D. Maio, I. Brandic, and R. Kueng, “An online approach for entanglement verification using classical shadows,” 2026. [Online]. Available: https://arxiv.org/abs/2603.26602 [4] A. Rath, V. Vitale, S. Murciano, M. Votto, J. Dubail, R. Kueng, C. Branciard, P. Calabrese, and B. Vermersch, “Entanglement barrier and its symmetry resolution: Theory and experimental observation,” PRX Quantum, vol. 4, p. 010318, 2023. [Online]. Available: https://link.aps.org/doi/10.1103/PRXQuantum.4.010318 [5] A. Jangda and M. Yadav, “Fast kronecker matrix-matrix multiplication on gpus,” in Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming, ser. PPoPP ’24. Association for Computing Machinery, 2024, p. 390–403. [Online]. Available: https://doi.org/10.1145/3627535.3638489 [6] A. Rath, C. Branciard, A. Minguzzi, and B. Vermersch, “Quantum fisher information from randomized measurements,” Phys. Rev. Lett., vol. 127, p. 260501, Dec 2021. [Online]. Available: https: //link.aps.org/doi/10.1103/PhysRevLett.127.260501 [7] V. Vitale, A. Rath, P. Jurcevic, A. Elben, C. Branciard, and B. Vermersch, “Robust estimation of the quantum fisher information on a quantum processor,” PRX Quantum, vol. 5, p. 030338, Aug 2024. [Online]. Available: https://link.aps.org/doi/10.1103/PRXQuantum.5.030338

[8] Y.-H. Wang and D.-J. Zhang, “Superiority of krylov shadow tomography in estimating quantum fisher information: from bounds to exactness,” npj Quantum Information, vol. 12, no. 1, p. 74, 2026. [Online]. Available: https://doi.org/10.1038/s41534-026-01216-z [9] K. Banaszek, M. Cramer, and D. Gross, “Focus on quantum tomography,” New Journal of Physics, vol. 15, no. 12, p. 125020, dec 2013. [Online]. Available: https://doi.org/10.1088/1367-2630/15/12/ 125020 [10] R. Blume-Kohout, “Optimal, reliable estimation of quantum states,” New Journal of Physics, vol. 12, no. 4, p. 043034, apr 2010. [Online]. Available: https://doi.org/10.1088/1367-2630/12/4/043034 [11] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, “Quantum state tomography via compressed sensing,” Phys. Rev. Lett., vol. 105, p. 150401, Oct 2010. [Online]. Available: https: //link.aps.org/doi/10.1103/PhysRevLett.105.150401 [12] S. Aaronson, “Shadow tomography of quantum states,” in Proceedings of 50th Annual ACM SIGACT Symposium on the Theory of Computing (STOC’18), ser. STOC 2018. New York, NY, USA: Association for Computing Machinery, 2018, p. 325–338. [Online]. Available: https://doi.org/10.1145/3188745.3188802 [13] B. Vermersch, A. Rath, B. Sundar, C. Branciard, J. Preskill, and A. Elben, “Enhanced estimation of quantum properties with common randomized measurements,” PRX Quantum, vol. 5, p. 010352, Mar 2024. [Online]. Available: https://link.aps.org/doi/10.1103/ PRXQuantum.5.010352 [14] C. Rasmussen and C. Williams, Gaussian processes for machine learning. United States: MIT Press, 2006, vol. 2. [Online]. Available: https://doi.org/10.7551/mitpress/3206.001.0001 [15] A. Wilson and H. Nickisch, “Kernel interpolation for scalable structured gaussian processes (kiss-gp),” in Proceedings of the 32nd International Conference on Machine Learning, ser. Proceedings of Machine Learning Research, F. Bach and D. Blei, Eds., vol. 37. Lille, France: PMLR, 07–09 Jul 2015, pp. 1775–1784. [Online]. Available: https://proceedings.mlr.press/v37/wilson15.html [16] M. Davio, “Kronecker products and shuffle algebra,” IEEE Transactions on Computers, vol. C-30, no. 2, pp. 116–125, 1981. [Online]. Available: https://www.doi.org/10.1109/TC.1981.6312174 [17] A. N. Langville and W. J. Stewart, “The kronecker product and stochastic automata networks,” Journal of Computational and Applied Mathematics, vol. 167, no. 2, pp. 429–447, 2004. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0377042703009312 [18] A. Neven, J. Carrasco, V. Vitale, C. Kokail, A. Elben, M. Dalmonte, P. Calabrese, P. Zoller, B. Vermersch, R. Kueng, and B. Kraus, “Symmetry-resolved entanglement detection using partial transpose moments,” npj Quantum Information, vol. 7, no. 1, p. 152, 2021. [Online]. Available: https://doi.org/10.1038/s41534-021-00487-y [19] Y. Zhou and Z. Liu, “A hybrid framework for estimating nonlinear functions of quantum states,” npj Quantum Information, vol. 10, no. 1, p. 62, Jun. 2024. [Online]. Available: https://doi.org/10.1038/ s41534-024-00846-5 [20] Q.-X. Zhang, Q.-M. Ding, H. Chen, S. Ling, Y.-L. Mao, D.-J. Zhang, X. Yuan, and Z.-D. Li, “Gap-closing experimental estimation of quantum fisher information,” Science Advances, vol. 12, no. 37, p. eaef4342, 2026. [Online]. Available: https://www.science.org/doi/abs/ 10.1126/sciadv.aef4342 [21] T. Häner and D. S. Steiger, “0.5 petabyte simulation of a 45-qubit quantum circuit,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’17. New York, NY, USA: Association for Computing Machinery, 2017. [Online]. Available: https://doi.org/10.1145/3126908.3126947 [22] R. Shi, G. Schieffer, P.-H. Lin, M. Gokhale, A. Herten, and I. Peng, “High-performance vector-length agnostic quantum circuit simulations on arm processors,” in 2026 IEEE International Parallel and Distributed Processing Symposium (IPDPS), 2026, pp. 928–941. [Online]. Available: https://www.doi.org/10.1109/IPDPS65963.2026.00080 [23] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, “Quantum entanglement,” Rev. Mod. Phys., vol. 81, pp. 865–942, Jun 2009. [Online]. Available: https://link.aps.org/doi/10.1103/RevModPhys. 81.865 [24] S. L. Braunstein and C. M. Caves, “Statistical distance and the geometry of quantum states,” Phys. Rev. Lett., vol. 72, pp. 3439–3443, May 1994. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett. 72.3439

[25] R. Horodecki and M. Horodecki, “Information-theoretic aspects of inseparability of mixed states,” Phys. Rev. A, vol. 54, pp. 1838– 1843, Sep 1996. [Online]. Available: https://link.aps.org/doi/10.1103/ PhysRevA.54.1838 [26] A. Peres, “Separability criterion for density matrices,” Phys. Rev. Lett., vol. 77, pp. 1413–1415, Aug 1996. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.77.1413 [27] W. Hoeffding, “A Class of Statistics with Asymptotically Normal Distribution,” The Annals of Mathematical Statistics, vol. 19, no. 3, pp. 293 – 325, 1948. [Online]. Available: https://doi.org/10.1214/aoms/ 1177730196 [28] S. L. Braunstein, C. M. Caves, and G. Milburn, “Generalized uncertainty relations: Theory, examples, and lorentz invariance,” Annals of Physics, vol. 247, no. 1, pp. 135–173, 1996. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0003491696900408 [29] L. Pezzé and A. Smerzi, “Entanglement, nonlinear dynamics, and the heisenberg limit,” Phys. Rev. Lett., vol. 102, p. 100401, Mar 2009. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett. 102.100401 [30] D.-J. Zhang and D. M. Tong, “Krylov shadow tomography: Efficient estimation of quantum fisher information,” Phys. Rev. Lett., vol. 134, p. 110802, Mar 2025. [Online]. Available: https: //link.aps.org/doi/10.1103/PhysRevLett.134.110802 [31] A. Elben and B. Vermersch, “RandomMeas.jl: A Julia Package for Randomized Measurements in Quantum Devices,” Quantum, vol. 10, p. 2086, 2026. [Online]. Available: https://doi.org/10.22331/ q-2026-04-28-2086 [32] R. F. Werner, “Quantum states with einstein-podolsky-rosen correlations admitting a hidden-variable model,” Phys. Rev. A, vol. 40, pp. 4277– 4281, Oct 1989. [Online]. Available: https://link.aps.org/doi/10.1103/ PhysRevA.40.4277 [33] Advanced Micro Devices, Inc., “Software optimization guide for AMD EPYC™ 7003 processors,” Advanced Micro Devices, Inc., Tech. Rep. 56665, 2020, rev. 3.00. [34] ——, “Software optimization guide for the AMD Zen4 microarchitecture,” Advanced Micro Devices, Inc., Tech. Rep. 57647, 2023, rev. 1.01. [35] ——, “Software optimization guide for the AMD Zen5 microarchitecture,” Advanced Micro Devices, Inc., Tech. Rep. 58455, 2024, rev. 1.00. [36] L. Dagum and R. Menon, “Openmp: an industry standard api for shared-memory programming,” IEEE Computational Science and Engineering, vol. 5, no. 1, pp. 46–55, 1998. [Online]. Available: https://www.doi.org/10.1109/99.660313 [37] A. Fog, C++ Vector Class Library, 2023, version 2.02.01, Apache License 2.0. https://github.com/vectorclass/version2. [38] M. Gruber, S. Herbst, I. Brandic, and V. D. Maio, “KrMMshadows,” 2026, Software Repository. [Online]. Available: https: //github.com/matthias2911/KrMM-shadows [39] A. Elben, S. T. Flammia, H.-Y. Huang et al., “The randomized measurement toolbox,” Nature Reviews Physics, vol. 5, pp. 9–24, 2023. [Online]. Available: https://doi.org/10.1038/s42254-022-00535-2 [40] S. Chen, W. Yu, P. Zeng, and S. T. Flammia, “Robust shadow estimation,” PRX Quantum, vol. 2, p. 030348, 2021. [Online]. Available: https://link.aps.org/doi/10.1103/PRXQuantum.2.030348

Related documents

• Deep thermalization and Hilbert space ergodicity#1078319
• Optimal classical shadow estimation of unitary channels at Heisenberg limit#1141568
• Fast surrogate modelling of EIT in atomic quantum systems using LSTM neural networks#1035322
• Hyperbolic Restricted Boltzmann Machine Neural Quantum State#1061937
• Finite-Data Non-Identification of Effective Quantum Coarse Dynamics#996446