Direct Measurement of Density Matrices via Dense Dual Bases
Yu Wang Hanru Jiang Yongxiang Liu Keren Li








Beijing Institute of Mathematical Sciences and Applications




Beijing Institute of Mathematical Sciences and Applications



liuyx@pcl.ac.cn
Peng Cheng Laboratory, Shenzhen 518055, China



likr@szu.edu.cn
College of Physics and Optoelectronic Engineering, Shenzhen University, Shenzhen 518060, China
Quantum Science Center of Guangdong-Hong Kong-Macao Greater Bay Area
(Guangdong), Shenzhen 518045. China



             

 
 

§ ABSTRACT  


A pivotal step to learn a quantum system is to select proper observables. 
Pauli observables and mutually unbiased bases are mainstream choices, thus being either utilized in practice or deemed theoretically optimal for state tomography. 
However, for Pauli observables, partially knowing a density matrix cannot be realized within a constant number of observables, and substantial observables are required for full tomography.
For mutually unbiased bases, the existence of informationally complete ones of all dimension is uncertain.
Therefore, alternative observables are necessitated. 
In this work, 2d observables are designed, with their implementations specified, which enabling a complete depiction of any d-dimension quantum state. To emphasize the advantages, we apply them in two scenarios.
Firstly, we show that direct measurement of density matrix elements is feasible without auxiliary systems, and any element is extracted by three selected observables.
Secondly, we show that state tomography of unknown rank-r density matrices, excluding only a negligible set, can be realized with O(r log d) observables. 
The designing observables present a significant advancement in the field, enhancing both the efficiency and feasibility to learn the quantum system.


Efficient understanding of a quantum system fundamentally relies on the selection of observables. Pauli observables and mutually unbiased bases (MUBs) are widely used in practice and often regarded as theoretically optimal for quantum state tomography (QST). However, Pauli observables require a large number of measurements for complete tomography and do not permit direct measurement of density matrix elements with a constant number of observables. For MUBs, the existence of complete sets of d+1 bases in all dimensions remains unresolved, highlighting the need for alternative observables. 
In this work, we introduce a novel set of 2d observables specifically designed to enable the complete characterization of any d-dimensional quantum state. To demonstrate the advantages of these observables, we explore two key applications. First, we show that direct measurement of density matrix elements is feasible without auxiliary systems, with any element extractable using only three selected observables. Second, we demonstrate that QST for unknown rank-r density matrices—excluding only a negligible subset—can be achieved with O(r log d) observables. This significantly reduces the number of unitary operations compared to compressed sensing with Pauli observables, which typically require O(r d log^2 d) operations. Each circuit is iteratively generated and can be efficiently decomposed into at most O(n^4) elementary gates for an n-qubit system. 
The proposed observables represent a substantial advancement in the characterization of quantum systems, enhancing both the efficiency and practicality of quantum state learning and offering a promising alternative to traditional methods.




Direct Measurement of Density Matrices via Dense Dual Bases
    Keren Li
    September 9, 2024
===========================================================




 


§ INTRODUCTION


The density matrix (DM) is a fundamental representation of a quantum state <cit.>, crucial for understanding quantum systems. Determining the DM accurately is a central task in quantum science. Traditionally, this problem is tackled through Quantum State Tomography (QST) <cit.>, which involves performing informationally complete (IC) measurements <cit.> and post-processing the data to estimate the quantum state. In a d-dimensional Hilbert space, a general density matrix has d^2 - 1 independent parameters. An IC positive operator valued measurement should contain at least d^2 projectors <cit.>.

Projective measurements (PMs) onto d + 1 mutually unbiased bases (MUB) are considered optimal IC measurements <cit.>. However, the existence of such d+1 MUBs when d is not a prime power remains an open question in quantum information theory  <cit.>. In an n-qubit system, where d = 2^n, estimating the 4^n Pauli expectation values requires 3^n unitary operations, each followed by projective measurement (PM) on the computational basis <cit.>.  This leads to 6^n different rank-1 projectors and necessitates exponential storage for the data, making QST impractical for large dimensions.



The off-diagonal elements of the DM are crucial as they capture key quantum properties such as entanglement <cit.> and decoherence <cit.>. Determining the entire DM is inefficient when only a subset of elements is needed. To address this, direct measurement protocols (DMPs) have been developed <cit.>, which involve various weak couplings between the main system and ancillary pointers, along with PMs. Although weak measurements slightly disturb the system and are easier to implement, they are inherently biased procedures that introduce unavoidable errors in the reconstructed state. Consequently, DMPs are less precise than QST <cit.>. To improve accuracy and precision, methods similar to weak measurement protocols have been proposed, where the coupling strength can be strong <cit.>.

DMPs can reduce computational effort for specific DM elements but add complexity by coupling ancillary pointers to the main system. For instance, measuring DM element ρ_jk typically involves tailored coupling operations U_j, followed by post-selection of the state |k⟩, and measuring the ancillary pointers <cit.>. For pure states or rank-1 DMs, a single pointer is sufficient <cit.>, but more general DMs require additional pointers <cit.>. Methods like δ-quench measurements have been proposed to eliminate the need for pointers, although they are restricted to specific systems <cit.>. Recent advancements have reduced the number of required pointers to one for multi-qudit DMs <cit.>, and even eliminated all pointers using phase-shifting techniques, though this approach necessitates O(d^2) unitary operations <cit.>. Implementing d^2 independent projectors with O(d^2) unitary operations in QST can also directly measure all DM elements <cit.>, but this method remains more resource-intensive compared to the O(d) coupling operations in two-pointer DMPs. This underscores the ongoing need for more efficient methods that can directly measure DM elements without ancillary systems and with reduced operational complexity.

Thus, we investigate whether O(d) unitary operations and projective measurements (PM) on the computational basis can be used for QST in arbitrary dimension d, eliminating the need for ancillary pointers to directly measure all DM elements. We are particularly interested in whether each DM element can be measured using a constant number of unitary operations from these O(d) operations.

In this work, we achieve this by designing a set of 2d eigenbases for QST, corresponding to 2d unitary operations followed by PM on the computational basis. Unlike the minimal d+1 MUBs, the existence of these observables for any dimension d is ensured through a deterministic construction algorithm. This set of observables also functions as a DMP, enabling each DM element to be measured with just three observables, thereby ensuring accuracy without the need for ancillary systems. Furthermore, these observables are applied to perform QST on quantum states with prior knowledge of rank-r, connecting through DMPs via matrix completion techniques. As a result, we demonstrate that for a rank-r density matrix, only O(r log d) of the 2d observables are required for full-state characterization. This approach dramatically reduces the required unitary operations, showing an exponential decrease compared to the O(rd log^2 d) operations needed when using random Pauli observables from the d^2 set via compressed sensing <cit.>. Finally, we present a unified formula that expresses all these unitary operations, where each can be decomposed into a permutation gate followed by Pauli measurements. The permutation gate itself can be efficiently decomposed into O(n^4) gates on an n-qubit system. This not only enhances the efficiency of quantum state learning but also opens new pathways for practical implementations in high-dimensional quantum systems. 





§ PRELIMINARIES AND DENSE DUAL BASES





Observables and Projective Measurements.– When measuring a DM ρ with an observable O = ∑_k=1^d λ_k |O_k⟩⟨ O_k|, the Born rule states that the measurement outcome λ_k occurs with probability tr(ρ |O_k⟩⟨ O_k|). This corresponds to PM onto the orthonormal eigenbasis {|O_k⟩}_k=1^d. To simulate this PM, we can use the unitary operation U = ∑_k=1^d |O_k⟩⟨ k|, transforming the measurement into the computational basis {|k⟩}_k=1^d. As applying U^† followed by a PM on the computational basis yields the outcome k with probability tr(U^†ρ U |k⟩⟨ k|) = tr(ρ |O_k⟩⟨ O_k|). Our focus is on designing PMs on O(d) eigenbases, implemented with O(d) unitary operations and PM on the computational basis, to directly measure all DM elements. 
 
 

Regardless of directness, 4^n Pauli observables are a common choice to reconstruct an n-qubit unknown state:

    ρ=1/2^n∑_i_1,⋯,i_n=0^3 tr(ρσ_i_1⊗⋯⊗σ_i_n)σ_i_1⊗⋯⊗σ_i_n
 
A total of 4^n expectation values should be measured <cit.>. Experimentally, 3^n Pauli observables, excluding the identity I for each qubit, are sufficient to obtain the 4^n expectation values. However, this requires 3^n unitary operations and computational basis measurements, resulting in 6^n projectors, which is far more than the 4^n needed. For example, the observable Z⊗ Z on a 2-qubit system corresponds to four different eigenstates: |↑↑⟩, |↑↓⟩, |↓↑⟩, and |↓↓⟩. For arbitrary dimension d, the 4^n Pauli matrices can be generalized to d^2 Gell-Mann matrices, Heisenberg-Weyl matrices, and so on. 


PMs onto d+1 MUBs are the minimal and optimal strategy for QST, which produces d(d+1) projectors. 
Two orthonormal bases {|a_j⟩}_j=1^d and {|b_k⟩}_k=1^d are termed as mutually unbiased if |⟨ a_j | b_k ⟩|^2 = 1/d for all j,k. 
For prime power dimension d, d+1 MUBs can be constructed <cit.>. 
In N-qubit systems, all MUB measurements can be efficiently implemented with 2^N+1 decomposed circuits and PMs on the computational basis <cit.>. MUBs are widely used in quantum information theory <cit.>. However, for arbitrary dimension d, the existence of d+1 MUBs remains unknown <cit.>, and strong numerical evidence suggests that only 3 MUBs exist for d=6 <cit.>. 

                  



For any dimension d, Caves, Fuchs, and Schack constructed d^2 rank-1 projections to directly determine all DM elements <cit.>. 
Denote |ϕ_jk^±⟩≐(|j⟩± |k⟩)/√(2), |ψ_jk^±⟩≐(|j⟩± i|k⟩)/√(2), where j,k∈[d]≐{0,⋯,d-1}. 
The collection of d^2 projected states are 

    𝒜_d={|l⟩,|ϕ_jk^+⟩,|ψ_jk^+⟩: 0≤ j<k≤ d-1;  l∈ [d]}.
 
Every ρ_ij is reconstructed via at most four projectors. 

    ρ_ll=(ρ|l⟩⟨ l|),
          ρ_jk=(ρ(|ϕ_jk^+⟩⟨ϕ_jk^+|-i|ψ_jk^+⟩⟨ψ_jk^+|))-1-i/2(ρ_kk+ρ_jj).
 
The challenges associated with implementing those O(d^2) projectors, as discussed in <cit.>, are highlighted in comparison to using O(d) unitary coupling operations and computational basis measurements in two-pointer DSP. 

We use an efficient algorithm to construct at most 2d eigenbases to cover the states in Eq. (<ref>) for arbitrary dimension d. 

These are referred to as dense dual bases (DDBs). 




 

§ RESULTS AND ANALYSIS


 
 For an arbitrary dimension, all elements of a density matrix can be directly measured using   2d-1 DDBs for even d and 2d DDBs for odd d. 
	


Analysis. When designing each eigenbasis, orthogonality and completeness are the key constraints, precisely, each pair of elements in an eigenbasis should be orthogonal, and there should be d elements in each basis.
Thus we introduce elements {|ϕ_jk^-⟩} and {|ψ_jk^-⟩} into the set 𝒜_d in Eq. (<ref>). 
Consequently, there are 2d^2-d elements. For the initial case d=2, the six elements {|0⟩, |1⟩, |ϕ_01^+⟩, |ϕ_01^-⟩, |ψ_01^+⟩, |ψ_01^-⟩} are the eigenstates of the Pauli observables Z, X, Y. 
Here, the pair {(0, 1)} corresponds to the eigenbasis {|ϕ_01^±⟩} of Pauli X, and `dually', to the eigenbasis {|ψ_01^±⟩} of Pauli Y. For general d, we dually arrange |ϕ_jk^±⟩ and |ψ_jk^±⟩ into different eigenbases. Besides, elements |ϕ_j_1k_1^±⟩ and |ϕ_j_2k_2^±⟩ cannot be grouped into the same eigenbasis if j_1=j_2 or k_1=k_2. The minimal DDBs design should avoid reusing elements in different eigenbasis. 



It can be transformed into finding an optimal game strategy. Suppose we have many bands, each marked with d numbers {0, 1, ⋯, d-1}. The game requires us to cut each band to form combinations (pairs) of two numbers. The sequence is irrelevant. Each band can produce at most ⌊ d/2 ⌋ pairs. When d is odd, there are (d-1)/2 pairs and a single element. The target is to determine how we can use a minimal number of bands to ensure that each pair {(j,k):0 ≤ j < k ≤ d-1} can be formed.

Each partition of the cut band is dually associated with two eigenbases, similar to the case when d=2. When d is even, we need at least C_d^2/(d/2)=d-1 bands to construct the corresponding d-1 partitions, which corresponds to 2(d-1) eigenbases. Together with ℬ_0={|0⟩, ⋯, |d-1⟩}, there are 2d-1 DDBs. 
When d is odd, we require d bands to construct the corresponding d partitions, which corresponds to 2d DDBs. 


For n-qubit system, d=2^n, we illustrate the minimal d-1 partitions construction with n=log d iterations in Fig. (<ref>). 

For example, the seven DDBs corresponding to the three partitions for d=4 are as follows: 
 
    ℬ_0^4    = {|0⟩,|1⟩,|2⟩,|3⟩}, 
    ℬ_1^4    = { |ϕ_01^±⟩,|ϕ_23^±⟩},    𝒞_1^4 = { |ψ_01^±⟩,|ψ_23^±⟩}, 
    ℬ_2^4    = { |ϕ_02^±⟩,|ϕ_13^±⟩},    𝒞_2^4 = { |ψ_02^±⟩,|ψ_13^±⟩}, 
    ℬ_3^4    = { |ϕ_03^±⟩,|ϕ_12^±⟩},    𝒞_3^4 = { |ψ_03^±⟩,|ψ_12^±⟩}.

For n=1, the initial partition is 𝕋^2 = T^2_1 = {(0,1)}. 
For the general n, the partitions consist of two parts. The first part is obtained by merging the partitions of 𝕋^2^n-1 and 𝕋^2^n-1 + 2^n-1, specifically: 

    T^2^n_m_1  = T^2^n-1_m_1∪ (T^2^n-1_m_1+2^n-1),

 where 1≤ m_1≤2^n-1-1. 
 The second part is defined by the intersection of the set {0, …, 2^n-1 - 1} and the set {2^n-1, …, 2^n - 1}. Specifically: 
 
    T^2^n_m_2={(k,2^n/2+[(k+m_2)2^n/2]):k∈ [2^n/2]},

where 2^n-1≤m_2≤ 2^n-1. 

 

It is straightforward to verify that each pair {(j,k) : 0 ≤ j < k ≤ 2^n - 1} is present in the constructed partitions. Using a similar iterative technique, the minimal d-1 partitions for even d, or d partitions for odd d can be constructed. The details are provided in Appendix A. 

Without ancillas, the minimal eigenbases for QST are d+1 MUBs, if they exist, corresponding to d(d+1) projectors <cit.>. In contrast, DDBs can produce at most 2d^2 projectors and are applicable for arbitrary d. Earlier DMPs <cit.> employed four or three Pauli observables on two-pointers, resulting in a total of 16d^2 or 8d^2 projectors, respectively. To the best of our knowledge, 2d (or 2d-1) DDBs require a minimal number of projectors to directly measure all DM  elements, illustrated in Fig. (<ref>). Additionally, three DDBs are sufficient to determine each DM element. These three DDBs can also determine the d diagonal elements and d/2 off-diagonal elements. 

It should be noticed here that the construction of minimal partitions is not unique. We can find many solutions when we use depth-first or breadth-first search algorithms for brute-force search. The above-designed algorithm applies log d iterations, which significantly reduces the time complexity of the construction. 


For d=6 (qubit-qutrit), three MUBs exist, but strong numerical evidence suggests that a fourth MUB does not exist <cit.>. Thus they are not enough for QST. As a comparison, a total of 11 DDBs are sufficient to directly reconstruct all 36 × 36 DM elements. Numerical simulations in Appendix B verify this result.

 



 
 
  

 

The 2d DDBs can determine all DM elements {ρ_jk : j,k ∈ [d]}. When only partial DM elements are known, properties such as entanglement or decoherence can be inferred. Furthermore, the entire DM can be recovered using matrix completion techniques <cit.> with the partial DM elements {ρ_jk : j,k ∈ C}, where C corresponds to the prior knowledge about ρ.

Typically, prior knowledge can greatly reduce the measurement resources for QST in different application scenarios. Examples include matrix product states <cit.>, permutation-invariant states <cit.>, subsets of a fiducial state <cit.>, as well as low-rank states <cit.>, among others. Compressed sensing shows that randomly choosing O(rd log^2 d) out of d^2 Pauli observables suffices to determine a rank-r DM with high probability <cit.>. Here, we consider the recovery of low-rank density matrices using partial DDBs. 
 

 
 

To uniquely determine a rank-r DM of dimension d, O(rlogd/r) DDBs suffice, except on a measure zero set. In this context, the rank r is significantly smaller than the dimension d (r ≪ d).



 

Analysis. When r=1, the rank-1 DMs correspond to pure states. PMs onto 3d-2 states can uniquely determine all pure states, except for a measure-zero set <cit.>. These states can be {|l⟩, |ϕ_jk^+⟩, |ψ_jk^+⟩ : 0 ≤ j < k ≤ d-1, |j-k| ≤ 1; l ∈ [d]}. In another view, the five eigenbases for pure state tomography constructed by Goyeneche et al. <cit.> form the minimal cover of these 3d-2 states. Furthermore, these projectors onto 5d states are proved to be rank-1 strictly complete <cit.>. 
Rank-r strictly complete measurements can uniquely determine any given rank-r state with a probability close to 1. No other physical states of any rank share the same measurement probability distributions, except for a possibly dense set of rank-r states on a measure-zero set.  

 A POVM is rank-r strictly complete if it can determine the DM elements with the following labels:

    C = {(j,k) : |j-k| ≤ r}.

The whole DM can be recovered using the following convex optimization:

    X̂ =    min_Xtr(X|k⟩⟨ j|) - ρ_jk  for  (j,k)∈ C, 
       such that  X ≽ 0, tr(X) = 1.
 
Besides, it is proved that the strictly complete measurements are advantageous due to their compatibility with efficient convex optimization and their robustness to statistical noise and state preparation errors  (see Corollary 1 and 2 in <cit.> for detailed discussions).

These correspond to the 0-th to r-th diagonal lines of the DM, as depicted in Fig. (<ref>). Thus, the task is to select DDBs containing the following states:

    ℬ = {|l⟩, |ϕ_jk^+⟩, |ψ_jk^+⟩ : 0 ≤ j < k ≤ d-1, |j-k| ≤ r; l ∈ [d]}.

Equivalently, this involves finding the number of partitions containing pair (j,k), where |j-k| ≤ r, 0 ≤ j < k ≤ d-1. 





 



For d = 2^n or for a general value of d, we can construct O(r logd/r ) partitions for this task, as detailed in Appendix C. 





An alternative approach to covering the states in Eq. (<ref>) is to use 4r + 1 eigenbases, which is specifically applicable to dimensions of the form d = 2^n <cit.>. In contrast, Result <ref> necessitates a larger number of eigenbases but extends its applicability to arbitrary dimensions d. 

When using the compressed sensing method, random O(rd log^2 d) expectation values of Pauli observables can produce the estimation of the unknown state with high probability <cit.>. They can be implemented with O(rd log^2 d) different styles of unitary operation and PM on the computational basis. Thus PMs onto DDBs provide an exponential decrease in the styles of unitary operations. 

Furthermore, we think the exponential decrease in unitary rotations also leads to a reduction in the required sampling times. 
The strength lies in the different information distraction ability between two-outcome measurements and d-outcome measurements. 

Consider the direct measurement in QST discussed in <cit.>, which involves O(d^2) unitary rotations. They can correspond to PM on computational basis and d^2 - d two-outcome  POVMs: 

    { U_k|0⟩⟨ 0|U_k^†, I - U_k|0⟩⟨ 0|U_k^†},    k = 1, ⋯, d^2 - d.

Here, U_k|0⟩ represents the nontrivial DDB elements in Eq. (<ref>).

Let p_k = tr(ρ_0 U_k|0⟩⟨ 0|U_k^†), where ρ_0 is the unknown DM. To estimate p_k with a confidence level of 1-σ and an error less than ϵ, the required sampling time N_k is given by <cit.>:

    N_k ≈2ln(1/σ) · p_k(1 - p_k)/ϵ^2 = ln(1/σ)/2ϵ^2.



However, the N_k(1-p_k) experiments corresponding to measurement I - U_k|0⟩⟨ 0|U_k^† are essentially wasted, as they do not provide additional information. 

To determine the O(rd) DM elements labeled in Eq. (<ref>), the complexity of sampling times is O(rd ·ln(1/σ)/ϵ^2).

When extending the two-outcome measurement to the d-outcome DDBs measurements, the waste component I - U_k|0⟩⟨ 0|U_k^† is decomposed into d-1 rank-1 projections for other DDBs, which can provide up to d/2-1 pieces of new information. Given that a total of O(rlog d) unitary operations are used, the sample complexity is O(r log d ·ln(1/σ)/ϵ^2). 


Pauli observables can be treated as two-outcome measurements {P_k^+, I - P_k^+}, where P_k^+ is the eigenspace with eigenvalue 1. Consequently, we observe a similar reduction in sampling complexity compared to compressed sensing with randomized Pauli measurements. However, even though the types of unitary operations for DDBs decrease to O(r log(d/r)), the total number of projectors remains O(rd log(d/r)). This indicates that the data required for convex optimization still scales with d, albeit with a slight reduction. 
Besides, the parameter ϵ for each expectation value is intrinsically linked to the dimension d, as it is crucial to minimize the errors for each observable, ensuring that the overall accumulated error remains within a tolerably small threshold. 

To test the results, we process a numerical simulation, which explores the fidelity of reconstructed quantum density matrices under varying numbers of measurements and ranks. 
The simulation begins with random quantum density matrices. For each matrix, multiple reconstructions are performed using different numbers of measurements via DDBs. The fidelity of the reconstructed matrices is then calculated, with the process repeated 20 times to obtain average fidelity values.
Figure <ref>(a) shows how the fidelity changes with the number of measurements for a fixed dimension (d=16) while varying the rank of the density matrix (r=3,4,8,16), demonstrating that increasing the number of measurements generally improves the fidelity of the reconstruction. In addition, we conducted a test on reconstructing density matrices using either DDBs or a Pauli-based compressed sensing method. The compressed sensing involved an initial measurement with varying numbers of random Pauli bases (CSPs), specifically rn, 10rn, drn^2, and 10drn^2, followed by convex optimization to minimize the nuclear norm of the density matrix. In contrast, Result <ref> employed 𝒪(rn) DDBs with semidefinite optimization to minimize the Frobenius distance. As shown in Figure <ref>(b)-(d), corresponding to d=4, 8, 16, our methods demonstrate clear advantages in requiring fewer measurements.

 
 
 

 
The failed set of determinations can be characterized. Consider r× r principal submatrices

    A_k=(                  
    	[         ρ_k,k             ⋯     ρ_k,k+r-1;             ⋮             ⋱             ⋮;     ρ_k+r-1,k             ⋯ ρ_k+r-1,k+r-1 ]),

k=0,…,d-r. The total failure set is when A_i is singular for i = 0, …, d - r - 1 and A_j is singular for j = i ± 1 <cit.>. 
Adaptive strategy can handle a failure case when some diagonal elements are zero. DM is expressed with an ensemble {p_k,|ψ_k⟩}, ρ=∑ p_k|ψ_k⟩⟨ψ_k|.
For example, if ρ_00=0, the first component of all {ψ_k} is zero, implying  ρ_0l=ρ_l0=0 for all l=0,…,d-1.
We can reconstruct the submatrix by erasing the first column and row of ρ. Then the DDBs can be designed for the projected subspace spanned by {|1⟩,…,|d-1⟩}.





Now we consider the circuits implementation for the DDBs in n-qubit case. 

 
 
 
 


In n-qubit systems, each projective measurement onto 2 × 2^n - 1 DDBs can be implemented using a permutational operation followed by a Pauli measurement. Each permutation can be decomposed into O(n^4) elementary gates.  

Precisely, the projective measurements onto computational basis {|0⟩,⋯,|2^n-1⟩} is implemented by Pauli measurement Z^⊗ n. 

In total, 2^n - 1 permutational operations correspond to the whole partitions in Result <ref>. 
For t = 1, the permutation P_1 = I^⊗ n. The subsequent Pauli observables are Z^⊗ (n-1)⊗ X and, dually, Z^⊗ (n-1)⊗ Y.

For t = 2, …, 2^n - 1, find j and k such that t = 2^k + j and j < 2^k. The permutational operation is

    P_t = I^⊗ (n - k) ⊗( |0⟩⟨0| ⊗ I^⊗ (k-1) + |1⟩⟨1| ⊗ (𝒰_k-1)^j ),

where 𝒰_k-1 is a (k-1)-qubit operation and 

    𝒰_k=∑_m=0^2^k-1|m-1 2^k⟩⟨ m|.
 The subsequent Pauli observables are Z^⊗ (n - k)⊗ X ⊗ Z^⊗ (k - 1) and, dually, Z^⊗ (n - k)⊗ Y ⊗ Z^⊗ (k - 1). 
The total involved 2n Pauli observables can be implemented with a Hadarmard gate H (or H̃^†) and computational basis measurement at each qubit. 


 

Analysis. 
It is straightforward that the PM onto computation basis is implemented with Pauli measurement Z^⊗ n.  
The universal expression of permutational operation lies in the iterative constructions of the partitions. 
 

When n = 1, the three DDBs are the eigenbases of the Pauli Z, X, and Y operators. The implementation of a Pauli X (or Y) measurement can be transformed into a sequence of applying the H (or H̃^†) operation followed by a Pauli Z measurement. Here H is the Hadamard gate  and H̃=S· H, S=|0⟩⟨ 0|+i|1⟩⟨ 1|.

When n = 2, the merging partition is T_1^4 = T_1^2 ∪ (T_1^2 + 2) = {(0, 1), (2, 3)}. Expressed in binary form, one DDB is ℬ_1^4 = {(|00⟩± |01⟩) / √(2), (|10⟩± |11⟩) / √(2)}. We can obtain it by acting I ⊗ H on the computational basis. 
The first intersection partition is T_2^4 = {(0, 2), (1, 3)}, corresponding to ℬ_2^4 = {(|00⟩± |10⟩) / √(2), (|01⟩± |11⟩) / √(2)}. The unitary operation U_2^4 = H ⊗ I transforms the computational basis into this form. 
The second intersection partition is T_3^4 = {(0, 3), (1, 2)}, corresponding to the Bell basis ℬ_3^4 = {(|00⟩± |11⟩) / √(2), (|01⟩± |10⟩) / √(2)}. The unitary operation U_3^4 = (|0⟩⟨0| ⊗ I + |1⟩⟨1| ⊗ X) · (H ⊗ I) transforms the computational basis into it. By changing H into H̃, we obtain the circuit for the dual basis.
 
 For general n, when t=1,⋯,2^n-1-1, the partition T_t^2^n is obtained by merging and intersection. Denote U_t as the operation to transform the computational basis into one nontrivial DDB. For the emerging partitions by Eq.(<ref>), U_t=I⊗ U_t^2^n-1, where U_t^2^n-1 represents the t-th circuits at iteration (n-1), starting from I^⊗ (n-2)⊗ H. 
When t=2^n-1, T_2^n-1={(0,2^n-1),(1,2^n-1+1),⋯} by Eq.(<ref>), then U_2^n-1=H⊗ I^⊗ (n-1). When t=2^n-1+j, U_t= (|0⟩⟨0|⊗ I+|1⟩⟨1|⊗ (𝒰^†_n-1)^j])· (H⊗ I^⊗ (n-1)).  

As the operation H (H̃^†) followed by Pauli measurement Z is equal to Pauli measurement X (Y). Thus we only need to perform the permutation I^⊗ n and the ones in Eq.(<ref>) followed by a Pauli measurement to implement PMs onto 2^n-2 DDBs. Together with Pauli measurement Z^⊗ n, all circuits implementation of PMs onto DDBs are depicted. 


In Appendix. D <cit.>, we decompose each permutational operation P_j with at most O(n^4) elementary gates. 
𝒰_k is a k-qubit global quantum gate. Though the classical counterpart is a basic operation, realizing it in the quantum device may be challenging. It can be implemented with O(k^3) elementary gates with the decomposition of incremental gates and Corollary 7.6 for generalized Toffoli gate <cit.>. 
Recent work on high-dimensional quantum systems shows promising techniques to implement similar gates (<cit.>), potentially overcoming this technical barrier. 
  

 






§ CONCLUSION AND DISCUSSION
 
 

 

Effective characterization of quantum states is a fundamental problem in quantum science. We propose an efficient method for constructing and decomposing minimal DDBs to address this issue. In the task of QST, our approach is compared with the commonly used Pauli observables and the optimal d+1 MUBs, as detailed in Table. <ref>. Additionally, this method functions as a direct measurement protocol, allowing for the reconstruction of each DM element using a constant number of unitary operations. The strong von Neumann measurements employed in our scheme ensure accuracy, unlike the commonly used weak measurement protocols, and they eliminate the need for ancillary pointers, resulting in minimal 2d^2 projectors. In response to the challenge presented in <cit.>, our method demonstrates that only O(d) unitary operations, rather than O(d^2), along with PM in the computational basis, are also sufficient to directly reconstruct all DM elements. The sampling of the 2d DDBs could provide an exponential decrease when reconstructing low-rank DM. 
Furthermore, we validate our approach through numerical simulations and cloud-based quantum experiments, demonstrating its practical feasibility. 

 There are several promising directions for future research stemming from this work. Firstly, random DDBs hold significant potential for efficiently predicting properties of unknown quantum states, particularly in n-qubit shadows tomography, where they could enable constant-time post-processing <cit.> for any observable in one experiment—a contrast to the exponential worst-case scenario encountered in random Clifford measurements <cit.>. Secondly, optimizing permutation operations in DDBs, which currently exhibit O(n^4) complexity, presents an intriguing task. Extending to d-level systems, where nontrivial DDBs correspond to unitary matrices with only two nonzero elements per row and column, raises important questions about efficient implementation and identifying physical systems capable of performing these operations seamlessly. Thirdly, integrating DDBs with matrix recovery techniques in QST suggests that with certain prior information, the entire density matrix could be reconstructed using exponentially fewer unitary operations. Our work demonstrates this efficiency with rank-r information, and it is worth exploring whether other types of prior information could similarly simplify the reconstruction process. Lastly, comparing the accuracy of Pauli observables, MUBs, and DDBs in quantum information tasks across different platforms—especially under the influence of errors or when analyzing specific unknown states—could provide valuable guidance in establishing these methods as fundamental tools in characterizing unknown quantum states. These directions underscore the substantial potential of DDBs to advance both theoretical and practical aspects of quantum information science. 

Acknowledgements—
We thank Cheng Qian for the brute-force searches that helped clarify the structure of the DDBs construction. This work was supported by the National Natural Science Foundation of China under Grants No. 62001260, 42330707 (Y.W.), 11701536 (Y.L.), and 11905111 (K.L.); the Beijing Natural Science Foundation under Grant No. Z220002 (Y.W.); and the Major Key Project of PCL (Y.L.). 


Author Contributions—  
Y. W. conceived the idea for this paper, applied the rank-r QST, decomposed the circuits, and drafted the manuscript. H. J. constructed partitions for d = 2^n, while Y. L. handled cases where d/2 is odd. K. L. conducted the error analysis, numerical simulations, and cloud experiments, and provided revisions to the manuscript. All authors reviewed and approved the final version of the manuscript for submission. 





apsrev4-1





§ DETAILED PROOFS IN RESULT 1


In the main text, we have demonstrated how to construct d-1 partitions when d=2^n, such that each pair (j, k) with 0 ≤ j < k ≤ 2^n - 1 is included in one of these partitions. 
Now we give the detailed construction of minimal partitions for arbitrary dimension d. 


    For even d, we can always create d-1 partitions {T^d_1,⋯,T^d_d-1} such that each pair (j, k) with 0 ≤ j < k ≤ d - 1 is in one of the partitions. 



    We will use mathematical induction to prove the statement 𝕊(d), where d represents the even dimension.

Base Case: For d=2, there is only one pair, and its corresponding partition is denoted as 𝕋^2 = {T_1^2} = {(0,1)}. Consequently, 𝕊(2) holds true.

Inductive Step: 
Since d is even, one of the numbers d/2 or d/2+1 is also even. 
To prove 𝕊(d) for a general even d, we assume that either 𝕊(d/2) or 𝕊(d/2+1) holds and then construct the partitions for d. 



Case 1: d/2 is even.

According to the induction hypothesis 𝕊(d/2), there exists a set of partitions {T_t^d/2 : t = 1, 2, ⋯, d/2 - 1} covering each (j, k) with 0 ≤ j < k ≤ d/2 - 1.

Similar to the construction when d is the power of 2, we can construct the target partitions {T_t^d : t = 1, 2, …, d-1} using the following rules:

    T_t^d   =    T_t^d/2∪ (T_t^d/2+d/2),      t=1,2,...,d/2-1, 
    
    	T_t^d    =   {(j,d/2+[(t+j) d/2]): 0≤ j≤ d/2-1},       t=d/2,d/2+1,...,d-1.


Case 2: d/2 is odd.

According to the induction hypothesis 𝕊(d/2 + 1), there exists a set of partitions {T_t^d/2 + 1 : t = 1, 2, ⋯, d/2} covering each (j, k) with 0 ≤ j < k ≤ d/2.

We select pairs {(c_t, d/2)}, where c_t is the neighbor of d/2. We construct the target partitions {T_t^d : k = 1, 2, …, d-1} using the following rules:

 
    T_t^d    =    T_t^d/2+1∪(T_t^d/2+1+d/2)-(c_t,d/2)- (d/2+c_t,d)∪ (c_t,d/2+c_t),   ,   t=1⋯,d/2, 
    
    	T_t^d    =   {(j,d/2+[(t+j) d/2]):0≤ j≤ d/2-1},       t=d/2+1,...,d-1.

We illustrate the construction method for d=6 in Fig.(<ref>). 


Conclusion:
Finally, we prove that this construction covers all pairs (j, k) with 0 ≤ j < k ≤ d - 1, assuming that 𝕊(2⌈d/4⌉) holds.

  Case 1: j and k on the same side of d/2.  If 0 ≤ j < k ≤ d/2 - 1 or d/2 ≤ j < k ≤ d - 1, using the induction hypothesis 𝕊(2⌈d/4⌉), the tuple (j, k) is either in T_t^d/2 or (T_t^d/2 + d/2), covered by the partitions {T_t^d} as shown in Equations (<ref>) and (<ref>).

 Case 2: j and k on different sides of d/2. Assume 0 ≤ j < d/2 ≤ k ≤ d - 1.

      If d/2 is even, the tuple (j, k) can be found in the partitions constructed by Equation (<ref>).
      If d/2 is odd, there are two subcases:
        If k - j < d/2, the tuple (j, k) is in the partitions constructed by Equation (<ref>).
        If k - j = d/2, the tuple (j, k) is in the partitions constructed by Equation (<ref>).

Thus, the set of partitions 𝕋^d = {T_t^d : t = 1, 2, …, d - 1} is generated. Each tuple (j, k) can be found in at least one partition, and therefore 𝕊(d) holds.




    If d is odd, there exist d partitions such that any tuple (j,k) with 0 ≤ j < k ≤ d - 1 can be found in one of these partitions.



By Lemma 1, we can construct d partitions for an even dimension d + 1. Each partition contains a tuple (c_d, d), where c_d is a neighboring number of d. By replacing (c_d, d) with the single number c_d, we obtain d partitions satisfying the requirements for the odd dimension d.




 §.§ Algorithm for arbitrary dimension based on the lemmas







 §.§ Examples of minimal partitions and corresponding DDBs

The deterministic algorithm can generate the minimal partitions required for any dimension d. We will first present the partitions for 1-qubit, 2-qubit, and 3-qubit systems, which have dimensions of 2, 4, and 8, respectively. Subsequently, we will provide the corresponding DDBs for each partition. Then, we will demonstrate the partition and DDBs construction for the odd dimension d = 7.

The case of d=2 is trivial, as there is only one tuple (0,1) with 0≤ j<k ≤ 1. Therefore, 𝕋^2=T^2_1,

    T^2_1={(0,1)}.

The three IC DDBs are:

    ℬ_0^2    = {|0⟩, |1⟩}, 
    ℬ_1^2    = {|0⟩± |1⟩/√(2)},   𝒞_1^2 = {|0⟩± i|1⟩/√(2)}.

For d=4, there are 3 partitions, which are the ones in main text.  𝕋^4={T^4_1, T^4_2, T^4_3}, where

    T^4_1    = {(0,1), (2,3)}, 
    
      T^4_2    = {(0,2), (1,3)},   
      T^4_3 = {(0,3), (1,2)}.

The seven IC DDBs are:

    ℬ_0^4    = {|0⟩, |1⟩, |2⟩, |3⟩}, 
    ℬ_1^4    = {|0⟩± |1⟩/√(2), |2⟩± |3⟩/√(2)},   𝒞_1^4 = {|0⟩± i|1⟩/√(2), |2⟩± i|3⟩/√(2)}, 
    ℬ_2^4    = {|0⟩± |2⟩/√(2), |1⟩± |3⟩/√(2)},   𝒞_2^4 = {|0⟩± i|2⟩/√(2), |1⟩± i|3⟩/√(2)}, 
    ℬ_3^4    = {|0⟩± |3⟩/√(2), |1⟩± |2⟩/√(2)},   𝒞_3^4 = {|0⟩± i|3⟩/√(2), |1⟩± i|2⟩/√(2)}.

For d=8, there are 7 partitions in 𝕋^8. 

    T^8_1    = {(0,1),(2,3),(4,5),(6,7)}, 
    
      T^8_2    = {(0,2),(1,3),(4,6),(5,7)}, 
    
      T^8_3    = {(0,3),(1,2),(4,7),(5,6)}, 
    
      T^8_4    = {(0,4),(1,5),(2,6),(3,7)}, 
    
      T^8_5    = {(0,5),(1,6),(2,7),(3,4)}, 
    
      T^8_6    = {(0,6),(1,7),(2,4),(3,5)}, 
    
      T^8_7    = {(0,7),(1,4),(2,5),(3,6)}.

The fifteen IC DDBs are:

    ℬ_0^8    = {|0⟩, |1⟩, …, |7⟩}, 
    ℬ_1^8    = { |ϕ_01^±⟩, |ϕ_23^±⟩, |ϕ_45^±⟩, |ϕ_67^±⟩},   𝒞_1^8 = { |ψ_01^±⟩, |ψ_23^±⟩, |ψ_45^±⟩, |ψ_67^±⟩}, 
    ℬ_2^8    = { |ϕ_02^±⟩, |ϕ_13^±⟩, |ϕ_46^±⟩, |ϕ_57^±⟩},   𝒞_2^8 = { |ψ_02^±⟩, |ψ_13^±⟩, |ψ_46^±⟩, |ψ_57^±⟩}, 
    ℬ_3^8    = { |ϕ_03^±⟩, |ϕ_12^±⟩, |ϕ_47^±⟩, |ϕ_56^±⟩},   𝒞_3^8 = { |ψ_03^±⟩, |ψ_12^±⟩, |ψ_47^±⟩, |ψ_56^±⟩}, 
    ℬ_4^8    = { |ϕ_04^±⟩, |ϕ_15^±⟩, |ϕ_26^±⟩, |ϕ_37^±⟩},   𝒞_4^8 = { |ψ_04^±⟩, |ψ_15^±⟩, |ψ_26^±⟩, |ψ_37^±⟩}, 
    ℬ_5^8    = { |ϕ_05^±⟩, |ϕ_16^±⟩, |ϕ_27^±⟩, |ϕ_34^±⟩},   𝒞_5^8 = { |ψ_05^±⟩, |ψ_16^±⟩, |ψ_27^±⟩, |ψ_34^±⟩}, 
    ℬ_6^8    = { |ϕ_06^±⟩, |ϕ_17^±⟩, |ϕ_24^±⟩, |ϕ_35^±⟩},   𝒞_6^8 = { |ψ_06^±⟩, |ψ_17^±⟩, |ψ_24^±⟩, |ψ_35^±⟩}, 
    ℬ_7^8    = { |ϕ_07^±⟩, |ϕ_14^±⟩, |ϕ_25^±⟩, |ϕ_36^±⟩},   𝒞_7^8 = { |ψ_07^±⟩, |ψ_14^±⟩, |ψ_25^±⟩, |ψ_36^±⟩}.






For dimension d=7, the seven partitions is constructed by Lemma <ref>, which are  

    T^7_1   =   {(0,1),(2,3),(4,5),6}, 
    
      T^7_2   =   {(0,2),(1,3),(4,6),5}, 
    
      T^7_3   =   {(0,3),(1,2),4,(5,6)}. 
    
      T^7_4   =   {(0,4),(1,5),(2,6),3}, 
    
      T^7_5   =   {(0,5),(1,6),2,(3,4)}, 
    
      T^7_6   =   {(0,6),1,(2,4),(3,5)}, 
    
      T^7_7   =   {0,(1,4),(2,5),(3,6)}.

Then the 14 IC DDBs are the following: 

    ℬ_1^7    = { |ϕ_01^±⟩, |ϕ_23^±⟩, |ϕ_45^±⟩, |6⟩},   𝒞_1^7 = { |ψ_01^±⟩, |ψ_23^±⟩, |ψ_45^±⟩, |6⟩}, 
    ℬ_2^7    = { |ϕ_02^±⟩, |ϕ_13^±⟩, |ϕ_46^±⟩, |5⟩},   𝒞_2^7 = { |ψ_02^±⟩, |ψ_13^±⟩, |ψ_46^±⟩, |5⟩}, 
    ℬ_3^7    = { |ϕ_03^±⟩, |ϕ_12^±⟩, |4⟩, |ϕ_56^±⟩},   𝒞_3^7 = { |ψ_03^±⟩, |ψ_12^±⟩, |4⟩, |ψ_56^±⟩}, 
    ℬ_4^7    = { |ϕ_04^±⟩, |ϕ_15^±⟩, |ϕ_26^±⟩, |3⟩},   𝒞_4^7 = { |ψ_04^±⟩, |ψ_15^±⟩, |ψ_26^±⟩, |3⟩}, 
    ℬ_5^7    = { |ϕ_05^±⟩, |ϕ_16^±⟩, |2⟩, |ϕ_34^±⟩},   𝒞_5^7 = { |ψ_05^±⟩, |ψ_16^±⟩, |2⟩, |ψ_34^±⟩}, 
    ℬ_6^7    = { |ϕ_06^±⟩, |1⟩, |ϕ_24^±⟩, |ϕ_35^±⟩},   𝒞_6^7 = { |ψ_06^±⟩, |1⟩, |ψ_24^±⟩, |ψ_35^±⟩}, 
    ℬ_7^7    = { |0⟩, |ϕ_14^±⟩, |ϕ_25^±⟩, |ϕ_36^±⟩},   𝒞_7^7 = { |0⟩, |ψ_14^±⟩, |ψ_25^±⟩, |ψ_36^±⟩}.




Here the computation basis ℬ_0^7 is saved as the elements {|0⟩,⋯,|6⟩} have been included twice in the other DDBs. 

 


§ NUMERICAL EXPERIMENTS ON DIMENSION SIX


While it is known that the projective measurements onto d+1 MUBs are minimal and optimal QST measurements for a d-dimensional system, their construction for each dimension d is still an open question. The first dimension for which d+1 MUBs have not been constructed is 6, corresponding to a qubit-qutrit system, ℋ_2⊗ℋ_3. 

For d=6, there are 5 partitions in 𝕋^6, 

    T^6_1={(0,1),(2,5),(3,4)}, T^6_2={(0,2),(1,4),(3,5)}, T^3_3={(0,3),(1,2),(4,5)}, 
          T^6_4={(0,4),(1,5),(2,3)}, T^6_5={(0,5),(1,3),(2,4)}.


The corresponding 11 IC DDBs are denoted as 
{ℬ_0^6,ℬ_1^6...,ℬ_5^6,𝒞_1^6,...,𝒞_5^6}. 

    ℬ_0^6    = {|0⟩, …, |5⟩}, 
    ℬ_1^6    = { |ϕ_01^±⟩, |ϕ_25^±⟩, |ϕ_34^±⟩},   𝒞_1^6 = { |ψ_01^±⟩, |ψ_25^±⟩, |ψ_34^±⟩}, 
    ℬ_2^6    = { |ϕ_02^±⟩, |ϕ_14^±⟩, |ϕ_35^±⟩},   𝒞_2^6 = { |ψ_02^±⟩, |ψ_14^±⟩, |ψ_35^±⟩}, 
    ℬ_3^6    = { |ϕ_03^±⟩, |ϕ_12^±⟩, |ϕ_45^±⟩},   𝒞_3^6 = { |ψ_03^±⟩, |ψ_12^±⟩, |ψ_45^±⟩}, 
    ℬ_4^6    = { |ϕ_04^±⟩, |ϕ_15^±⟩, |ϕ_23^±⟩},   𝒞_4^6 = { |ψ_04^±⟩, |ψ_15^±⟩, |ψ_23^±⟩}, 
    ℬ_5^6    = { |ϕ_05^±⟩, |ϕ_13^±⟩, |ϕ_24^±⟩},   𝒞_5^6 = { |ψ_05^±⟩, |ψ_13^±⟩, |ψ_24^±⟩}.



We tested our proposal numerically in a 6-dimensional system. 
We examined four quantum state types: (a) the maximally mixed state I/6, (b) a balanced state 1/6∑_k,j=0^5|k⟩⟨ j|, (c) a separable state, and (d) an entangled state. These states have Hermitian, semi-definite, and unit trace density matrices. States (a) and (b) were directly generated in our simulation, while states (c) and (d) were prepared using local unitary transformations U_2∈ℋ_2 and U_3∈ℋ_3, distributed uniformly according to the Haar measure.
As U_2 and U_3 are local, the entanglement of the resulting state remains unchanged, which can be verified using the Peres-Horodecki criterion. We prepared states (c) and (d) as U_2⊗ U_3 |ϕ⟩, where |ϕ⟩ represents |0⟩ or |1⟩+|2⟩+|3⟩+|5⟩, respectively. 

For the reconstruction of an unknown density matrix ρ, we utilize two methods. The first method is based on semi-definite programming, where ρ̃ represents the estimated form of ρ and is obtained through a parameterized matrix X≥ 0. The mathematical model is formulated as follows:  

    ρ̃=_X∑_i=1^66( (X E_i)-p_i),

where · is a norm function, E_i are from the bases of (ℬ_0,...,ℬ_5,𝒞_1,...,𝒞_5), and
p_i are the measured probabilities on E_i by unknown density matrix ρ.

The second approach is direct reconstruction. This method utilizes a total of 36 probabilities, which is half of the available data.  

Consequently, numerical experiments were conducted 20 times for all tested states. The Monte Carlo method was utilized to simulate p_i with 100×2^ shots. The results of these numerical experiments are illustrated in SFig.(<ref>), where the infidelities are represented by the Frobenius distance. 

    F_f=√(((ρ-ρ̃)· (ρ-ρ̃)^†))).

Error bars come from standard deviations of 20 times repetitions of the simulation.







§ RIGOROUS PROOFS OF RESULT 2




    When d = 2^n, the pair (j, k) with 0 ≤ j < k ≤ 2^n - 1 and |j - k| ≤ r can be found in at most O(r (n - log r)) partitions for minimal DDBs.
 

Analysis. Using n iterations, 2^n - 1 partitions are constructed such that each pair (j, k) with 0 ≤ j < k ≤ 2^n - 1 is covered. When r ≪ 2^n, we can always find m such that 2^m-1 < r ≤ 2^m.

At iteration t = m for dimension 2^m, a total of 2^m - 1 partitions have been constructed by Result <ref>. Thus, at most 2^m - 1 partitions contain (j, k) for 0 ≤ j < k ≤ 2^m - 1. At iteration t = m + 1, new partitions are constructed by Eq. (<ref>) and Eq. (<ref>), resulting in at most 2^m - 1 and r partitions that satisfy |j - k| ≤ r, 0 ≤ j < k ≤ 2^m+1 - 1. Therefore, from t = m to t = log d = n, the number of relevant partitions is less than 2^m - 1 + r(log d - m) < r(log d + 2 - m) < r log(4d/r).  



    For any general dimension d, the pair (j, k) with 0 ≤ j < k ≤ d - 1 and |j - k| ≤ r can be found in at most O(r logd/r)) partitions for minimal DDBs.


Analysis.  Based on the construction in Algorithm 1, we iteratively construct b_k-1 partitions for dimension b_k, where k=L,L-1,⋯,1. 
Here b_L = 2 and b_1 = f(d). 

Firstly, we prove that the number of iterations is exactly L = ⌈log d ⌉. 
For example, consider d = 100. We should construct the partitions iteratively for the dimensions:

    b_7 = 2,   b_6 = 4,   b_5 = 8,   b_4 = 14,   b_3 = 26,   b_2 = 50,   b_1 = 100.
 
With ⌈log 100 ⌉ = 7 iterations, we obtain the partitions for d = 100. 


For any general d, we observe that b_k-1 = 2b_k or b_k-1 = 2b_k - 2, leading to the inequality b_k-1≤ 2b_k. If the iteration count is ⌈log d ⌉ - 1, the maximum value of b_1 is 2^⌈log d ⌉ - 1. 
However, 2^⌈log d ⌉ - 1 < d. This results in a contradiction since b_1 should equal d or d + 1. On the other hand, it is easy to observe that if we start with dimension b_1=f(d) and iterate L = ⌈log d ⌉ times, the final value can always be reduced to 2. This is because that 2^L-1< b_1 ≤ 2^L,  and generally 2^L-k< b_k ≤ 2^L-k+1 for k=2,⋯,L. 

Similar to the case when d=2^n, there is always a number m such that 2^m-1<r ≤ 2^m. Denote L-k+1=m. Then at iteration m, there are at most b_L+1-m-1≤ 2^m-1 partitions for dimension b_L+1-m. 
Repeating the procedure until the iteration L for dimension b_1, the new constructions are based on Eq. (<ref>) and Eq. (<ref>). The number of partitions containing the required elements is less than 2^m-1+r( L-m)< r(2+⌈log d ⌉ -m)<rlog (8d/r), where -m ≤ -log r. Then the number of DDBs required is O(r log(d/r)).  

Each partition corresponds to two eigenbases. Together with ℬ_0, O(r logd/r) DDBs can reconstruct the elements {ρ_jk : j,k ∈ C}. 



§ CIRCUITS ANALYSIS
 

We may as well label the 2^n+1-1 DDBs on n-qubit systems as follows: 

    {ℬ_0^2^n,  ℬ_j^2^n, 𝒞_j^2^n:j=1,⋯,2^n-1}.

The computational basis is ℬ_0^2^n={|0⟩,⋯,|2^n-1⟩}. 
The basis ℬ_j^2^n and 𝒞_j^2^n are dually designed for the same partition. Here we consider the 2^n+1-2 circuits to transform the computational basis to nontrivial DDBs {ℬ_j^2^n, 𝒞_j^2^n:j=1,⋯,2^n-1}. 


1-qubit:  The DDBs ℬ_1^2 and 𝒞_1^2 are in Eq. (<ref>). The DDB circuits to map the computational basis into them are shown in SFig.(<ref>). 
 

2-qubit: The six nontrivial DDBs are in Eq. (<ref>). With the binary form, the three DDBs for the three partitions are the following 

    ℬ_1^4    = {|00⟩± |01⟩/√(2), |10⟩± |11⟩/√(2)}, 
    ℬ_2^4    = {|00⟩± |10⟩/√(2), |01⟩± |11⟩/√(2)}, 
    ℬ_3^4    = {|00⟩± |11⟩/√(2), |01⟩± |10⟩/√(2)}.

The corresponding circuits are:  


n-qubit: 
Denote the 2^n-1 unitary operations as {U_t^2^n:t=1,⋯,2^n-1} for the nontrivial DDBs ℬ_t^2^n. The DDBs {𝒞_t^2^n} is dually obtained by adding i in the basis states. At the last iteration of (n-1)-qubit case, the unitary operations  {U_t^2^n-1 : t = 1, ⋯, 2^n-1 - 1} map the computational basis ℬ_0^2^n-1 into ℬ_t^2^n-1. 


For the n-qubit case, we have T_t^2^n = T_t^2^n-1∪ (T_t^2^n-1 + 2^n-1) for t = 0, ⋯, 2^n-1 - 1, and T_t^2^n = {(j, 2^n-1 + [(j + t)  2^n-1]) : 0 ≤ j ≤ 2^n-1 - 1} for t = 2^n-1, ⋯, 2^n - 1.

Thus, when t = 0, ⋯, 2^n-1 - 1, we have U_t^2^n = I ⊗ U_t^2^n-1. This is because the basis states of ℬ_t^2^n-1 are of the form |k_1⟩± |k_2⟩. The partition T_t^2^n is iteratively constructed by T_t^2^n-1∪ (T_t^2^n-1 + 2^n-1), where t = 0, ⋯, 2^n-1 - 1. Therefore, the basis states of ℬ_k^2^n are in the form |0⟩(|k_1⟩± |k_2⟩) or |1⟩(|k_1⟩± |k_2⟩).

When t = 2^n-1, we have U_2^n-1^2^n = H ⊗ I^⊗ n-1. This follows because the partition T_2^n-1^2^n consists of {(0, 2^n-1), (1, 2^n-1 + 1), ⋯, (2^n-1 - 1, 2^n - 1)}. These numbers can be expressed in binary form, and the corresponding basis states are given by {(|0⟩± |1⟩) ⊗ |j_2, ⋯, j_n⟩ : j_2, ⋯, j_n = 0, 1}. Hence, U_2^n-1^2^n = H ⊗ I^⊗ n-1.

When t = 2^n-1 + 1, ⋯, 2^n - 1, we can express t as 2^n-1 + j. Then, U_k^2^n =  [|0⟩⟨ 0|⊗ I + |1⟩⟨ 1|⊗ (𝒱_n-1)^j] · [H ⊗ I^⊗ (n-1)], where 𝒱_n-1 =𝒰^†_n-1= ∑_m=0^2^n-1 - 1 |m + 1  2^n-1⟩⟨ m|. 

This result is derived as follows. The partition T_t^2^n consists of {(0, 2^n-1 + j), (1, 2^n-1 + (1 + j  2^n-1)), ⋯, (2^n-1 - 1, 2^n-1 + (2^n-1 - 1 + j  2^n-1))}. Expressing these numbers in binary form, the corresponding basis states are given by {|0⟩⊗ |j_2, ⋯, j_n⟩ + |1⟩⊗ |[(j_2, ⋯, j_n) + j]  2^n-1⟩ : j_2, ⋯, j_n = 0, 1}. This basis can be obtained by applying the conditional shift operation |0⟩⟨ 0|⊗ I + |1⟩⟨ 1|⊗ (𝒰^†_n-1)^j on the basis of ℬ_2^n-1^2^n.
  



Remark: The circuits depicted on the middle and right of Fig. (<ref>) are realized during the n-th iteration. At each k-th iteration, where k=2, …, n-1, the newly generated 2^k circuits correspond to those shown in Fig. (<ref>). The operation P_j, as defined in Eq. (9) of the main text, represents the universal form of the conjugate transpose of controlled operations for iterations k=2, …, n. The structure of these circuits provides an alternative interpretation for the total of 2^n-1 partitions, illustrated by the sum 2^0 + 2^1 + … + 2^n-1 = 2^n - 1. 





 §.§ Circuits decomposition of permutational operation


As previously mentioned, the projective measurement (PM) of the state ρ onto the eigenbasis {U|k⟩ : k = 0, ⋯, 2^n - 1} can be implemented by first applying U^†, followed by a PM onto the computational basis ℬ_0^2^n. 
The operation H (H̃^†) followed by 1-qubit computational measurement (Pauli Z measurement) is equivalent to Pauli X (Y) measurement. Thus we just need to decompose the following operations in Fig. (<ref>):

    P_t= |0⟩⟨ 0|⊗ I + |1⟩⟨ 1|⊗ (𝒰_k-1)^j,

where 𝒰_k=∑_m=0^2^k-1|m-1 2^k⟩⟨ m|, k=2,⋯,n-1 and j=1,⋯,2^k-1-1. 
Next, we will demonstrate that all the permutation operations can be efficiently decomposed into polynomially many elementary gates. 
 



For each j = 1, …, 2^k-1, even in exponential time, (𝒰_k)^j can be decomposed into a linear combination of k specific unitary operations, namely (𝒰_k), (𝒰_k)^2^1, …, (𝒰_k)^2^k-1.



Analysis:  For each j ∈{1, …, 2^k-1}, we can express j in its binary form as j⃗= (j_0, j_1, …, j_k-1) where j_i ∈{0,1}. Specifically, j is  written as:


    j = j_0 × 2^0 + j_1 × 2^1 + ⋯ + j_k-1× 2^k-1.


Therefore, we have:


    (𝒰_k)^j = [(𝒰_k)^1]^j_0× [(𝒰_k)^2^1]^j_1×⋯× [(𝒰_k)^2^k-1]^j_k-1.

The circuit is illustrated in Fig. (<ref>). 

For example, if j_0 = 0, then the term [(𝒰_k)^1]^j_0 simplifies to the identity operation I, and can thus be omitted.
Consequently, we only need to identify the nonzero elements in {j_0, j_1, …, j_k-1} and perform at most k operations from the set {(𝒰_k)^1, (𝒰_k)^2^1, …, (𝒰_k)^2^k-1}. The order of these operations can vary since each of them commutes with the others.  




 


Consider the implementation of k unitary operations (𝒰_k)^1, (𝒰_k)^2^1, …, (𝒰_k)^2^k-1. These circuits are equivalent to those for 𝒰_k, 𝒰_k-1, …, 𝒰_1 = X. Specifically, 

    (𝒰_k)^2^l = 𝒰_k-l⊗ I^⊗ l
   for l = 0, …, k-1, as illustrated in Fig. (<ref>). Additionally, the circuit decomposition of 𝒰_l is shown in Fig. (<ref>).  


Analysis. We have (𝒰_k)^2^l=∑_m=0^2^k-1|m -2^l  2^k⟩⟨ m|. 
The binary form of 2^l is the following 

    l⃗=(0, ⋯, 0_k-l-1, 1, 0, ⋯, 0_l).
 
We express the binary form of m as m⃗=(m_0⋯,m_k-1). 
In the binary form, 

    (𝒰_k)^2^l   = ∑_m_0,⋯,m_k-1=0,1 |(m⃗-l⃗)  2^k⟩⟨m⃗| 
       =  ∑_m_0,⋯,m_k-1=0,1  
    |[(m_0⋯ m_k-l-1)-(0, ⋯, 0_k-l-1, 1)] 2^k-l⟩⟨ m_0⋯ m_k-l-1| ⊗  ( |m_k-l⋯ m_k-1⟩⟨ m_k-l⋯ m_k-1|)
       =∑_m=0^2^k-l-1|m-1 2^k-l⟩⟨ m| ⊗ I^⊗ l
       = 𝒰_k-l⊗ I^⊗ l.





 

Now we consider the decomposition of 𝒰_l, where l=1,⋯,k. 
When l=1, 𝒰_1=∑_m=0^1 |m-1  2⟩⟨ m|=|0⟩⟨ 1|+|1⟩⟨ 0|=X.  

For general l, 𝒰_l = ∑_m=0^2^l-1 |m-1  2^l⟩⟨ m| represents a global shift operation. The classical counterpart of this operation, the basic increment by 1 (`+1'), serves as its foundation. 
The circuit implementation of 𝒰_l is also used in pure state (rank-1 density matrix) QST <cit.>, as depicted in Fig. (<ref>). 
 


    Each permutation operation required to perform PM on DDBs can be decomposed into at most O(n^4) elementary 1-qubit and 2-qubit gates.


Analysis: Denote ∧_m(U) as the generalized Toffoli gate with m+1 input bits, which maps |x_1, …, x_m, y⟩ to |x_1, …, x_m, (∏_k=1^m x_k) ⊕ y⟩. On input (x_1, …, x_k, y), the gate applies U to y if and only if ∏_k=1^m x_k = 1. When m = 1, ∧_1(X) corresponds to the Controlled-NOT operation, expressed as ∧_1(X) = |0⟩⟨ 0| ⊗ I + |1⟩⟨ 1| ⊗ X.

Thus, the operation 𝒰_l in Fig. (<ref>) is a combination of the following gates:

    X, ∧_1(X), …, ∧_l-1(X).


According to Corollary 7.6 in <cit.>, the l-qubit gate ∧_l-1(X) can be decomposed into Θ(l^2) elementary 1-qubit and 2-qubit gates. Consequently, 𝒰_l can be decomposed into O(l^3) elementary gates. For the controlled operation |0⟩⟨ 0| ⊗ I + |1⟩⟨ 1| ⊗𝒰_l, the cost in elementary gates also scales as O(l^3), as it involves a combination of ∧_1(X), …, ∧_l-1(X), ∧_l(X).

Now, consider the worst-case scenario where the permutation operation P_t in Eq. (<ref>) incurs the maximum cost in terms of elementary gates. This occurs when k = n. For each j = 1, …, 2^n-1 - 1, P_t can be implemented with at most n controlled operations using Decomposition <ref>. Therefore, the upper bound for decomposing all permutation operations involved in DDBs is O(n^4). 

Using l-1 ancilla qubits, ∧_l-1(X) can be decomposed into O(l) elementary gates <cit.>. Consequently, each permutation operation in the DDB circuit can be decomposed into at most O(n^3) elementary gates. 
It is noteworthy that, with the following decomposition and strategy, the cost of gates could be further reduced. 

    The operations (𝒰_l)^j and (𝒰_l)^2^l-j can be implemented using the same number of gates, where j = 1, …, 2^l - 1.


Analysis: We have 𝒰_l = ∑_m=0^2^l-1 |m-1  2^l⟩⟨ m|. Thus, (𝒰_l)^j · (𝒰_l)^2^l-j = I. Therefore, if we perform the conjugate transpose circuit of (𝒰_l)^j, we obtain the circuit for (𝒰_l)^2^l-j. 

Strategy: By combining the analysis from Decompositions <ref> and <ref>, the circuit decomposition of (𝒰_l)^j can be further simplified compared to the binary expression. For instance, when j = 2^l - 1, we should integrate the circuits for 𝒰_l, (𝒰_l)^2^1, …, and (𝒰_l)^2^l-1 as described in Decomposition <ref>. However, according to Decomposition <ref>, it suffices to implement a single circuit for (𝒰_l)^†. As a result, the circuit components for (𝒰_l)^2^1, …, and (𝒰_l)^2^l-1 can be omitted, leading to a more efficient implementation.

In general, we can define a finite set of integers:

    S = {± 1, ± 2^1, …, ± 2^l-1}.

For any j ∈{1, …, 2^l - 1}, we can identify the minimal elements in the set S such that their sum equals j. We then decompose 𝒰_l^j according to these minimal elements in S, rather than just using the binary form corresponding to S' = {1, 2^1, …, 2^l-1}.  




 

 §.§ Three circuits for arbitrary DM element


 Now we consider the three circuits of PMs onto DDBs to directly reconstruct arbitrary unknown DM element ρ_jk, 0≤ j <k≤ d-1. 
 The basis for determining the diagonal element is the computational basis.  Thus the measurement is Pauli measurement Z^⊗ n.
 The other two circuits are constructed in the following way. 


    
  * Write the binary representation of j and k, which are j_1j_2⋯ j_n and k_1k_2⋯ k_n respectively. 
    
  * Find the first different qubit of |j_1j_2⋯ j_n⟩ and |k_1k_2⋯ k_n⟩. 
    We may as well denote it as q_s. Namely, 
    
    {
            j   =j_1,⋯,j_s-1,j_s=0,j_s+1,⋯,j_n
    
            k   =k_1,⋯,k_s-1,k_s=1,k_s+1,⋯,k_n
        .

      Denote the difference between the binary numbers j_s+1,⋯,j_n and k_s+1,⋯,k_n as 
      
    l=∑_m=s+1^n (k_m-j_m)× 2^n-m.

    
  * The permutational operation for ρ_jk is defined by 
    
    P_j,k=I^⊗ s-1⊗ [|0⟩⟨ 0|⊗ I+|1⟩⟨ 1|⊗ (𝒰_n-s)^l].



The circuits for the PMs onto the nontrivial DDBs are depicted in Fig. (<ref>). When we go through all ρ_jk, the required circuit types are O(2^n) instead of O(4^n). We can verify the function of the conjugate transpose of the circuits. 

After applying the operation H, the state |j⟩ = |j_1 ⋯ j_n⟩ evolves to 

    |j_1⋯ j_s-1⟩|0⟩ + |1⟩/√(2) |j_s+1⋯ j_n⟩.

Following the conditional permutation operation, the final state becomes 

    |j_1⋯ j_n⟩ + |k_1⋯ k_n⟩/√(2) = |ϕ^+_jk⟩.

Thus, in the left circuit of Fig. (<ref>), if the measurement result is j_1,⋯,j_n, it corresponds to the projected state |ϕ^+_jk⟩. Similarly, if the result is k_1,⋯,k_n, it corresponds to |ϕ^-_jk⟩. 
In the right circuit, if the measurement result is j_1,⋯,j_n (or k_1,⋯,k_n), it corresponds to the projected state |ψ^+_jk⟩ (or |ψ^-_jk⟩). 

    



 





















































§ CLOUD EXPERIMENTS


To test the performance of our strategy, real experiments were carried out on two quantum computers, superconducting qubits on IBM quantum Lab and nuclear spins on Spinq. 


For quantum chip on ibmq-manila, it is with a one dimensional structure and 32 quantum volume, shown in SFig.(<ref>). Table.(<ref>) has its detailed parameters.
Only qubit 0 and 1 are used, with frequencies at 4.963 Ghz and 4.838 Ghz, with anharmonicity -0.34335 Ghz and -0.34621 Ghz, respectively. The error of CNOT gate from 0 on 1 is 6.437e-3 which costs 277.333ns, while the error of CNOT gate from 1 on 0 is 6.437e-3 which costs 312.889ns.









For entire experiments, states such as 

    |ψ_1i⟩   =   α_i|0⟩^⊗ 2+β_i|1⟩^⊗ 2
    |ψ_2i⟩   =   (α_i|0⟩+β_i|1⟩)^⊗ 2

are tested, with i=1,...,21 and α_i=cos(θ_i/2), β_i=sin(θ_i/2), θ_i= (i-1) π/20. 
The prepared circuits are depicted in SFig.(<ref>).

With the construction method, 7 measurement circuits are generated, which is shown in SFig.(<ref>).




For entire experiments, two sets were conducted. The first is through our strategy with data analysis of direct calculation and SDP. The second is the standard tomography protocol. The main results such as Frobenius's distance and fidelities are calcultaed in the manuscript. Here we list the density matrix for each experiment-prepared state as supplementary. 
SFig.(<ref>) and SFig.(<ref>) are from Eq.(<ref>) and Eq.(<ref>), respectively. Although 21 experiments were conducted, only 11 density matrices are listed from 0 to π.
SFig.(<ref>) (SFig.(<ref>)) is divided into two-lines, the first line is the real parts of the density matrix and the second line is the image parts.
Meanwhile, transparency parts are theoretical values, and solid parts are from experiments. 
 


  















As for the cloud quantum computer of SpinQ, it is a liquid NMR based architecture, that uses crotonic acid as their qubit system. 
As it is shown in SFig.(<ref>). 4 carbon nuclei are denoted as 4 qubits, where related parameters are listed in table.
With an external programmable radio-frequency pulses as control fields, almost 4-qubit quantum logic gates can be achieved. 
In the table, diagonal elements are frequencies, while off-diagonal elements are J-couplings, which are all measured at room temperature. 





In order to demonstrate our strategy, 30 measurement circuits are generated, which is shown in SFig.(<ref>). ℬ_0 is ignored here as it is a trivial one as conventional tomography strategy. However, we do not realize the following circuits, since some entangled gates are out of the current device's capability. Decoherence time doesn't allow us more quantum gates. 
Thus, we have to simulate circuits in SFig.(<ref>) by decomposing each measurement basis into Pauli and summing them up in the end, indirectly completing the proposal.  



Therefore, For entire experiments, states such as 

    |ψ_1i⟩   =   α_i|0⟩^⊗ 4+β_i|1⟩^⊗ 4
    |ψ_2i⟩   =   (α_i|0⟩+β_i|1⟩)^⊗ 4

are tested, with i=1,...,11 and α_i=cos(θ_i/2), β_i=sin(θ_i/2), θ_i= (i-1) π/10. 
The prepared circuits are depicted in SFig.(<ref>).

Similarly, results are presented with direct calculation and SDP. As a comparison, standard tomography protocol was also conducted. The main results of Frobenius's distance and fidelities are calculated in the manuscript. Here we only list the density matrix for each experiment-prepared state.


SFig.(<ref>) and SFig.(<ref>) are from Eq.(<ref>) and Eq.(<ref>), respectively. Although 11 experiments were conducted, only 6 density matrices are listed, where θ_i=0, π/5, 2π/2, 3π/5, 4π/5, π.
SFig.(<ref>) (SFig.(<ref>)) is divided into two lines, the first line is the real parts of the density matrix and the second line is the image parts.
Meanwhile, transparency parts are theoretical values, and solid parts are from experiments. As GHZ-like states were prepared through 3-CNOT gates, which cost around 100ms, the decoherence affects the states heavily. 


  



 §.§ Table of fidelity for entire experiments

At the end of this section, we list the table for the entire experiments in Table.(<ref>).




§ ERROR ANALYSIS

A real situation is that measurement bases {ℬ_0, ℬ_k, 𝒞_k, 1≤ k≤ d } cannot be perfectly realized. In most cases, we could create very close ones, which cause slight differences when obtaining the measured probability. That is to say, possibly we project an unknown target ρ onto an approximate basis state |̃ϕ̃_̃ĩ⟩̃, instead of the exact one, |ϕ_i⟩. 
For certain ideal and realized basis operators, the i-th basis states (eigenstates) have such relation,  

    |̃ϕ̃_̃ĩ⟩̃=|ϕ_i⟩+ϵ_i|e_i⟩,

where |ϕ_i⟩ is the exact i-th basis state, ϵ_i is a constant amplitude, and |e_i⟩ is one random state of haar measure. 
Obviously, vast repeated measurements produce averaged effects,  

    ∫|e_i⟩ d e_i =0,   ∫|e_i⟩⟨e_i| d e =I/d.

Therefore, probabilities measured on ρ is with such disturbance,

    (ρ|̃ϕ̃_̃ĩ⟩̃⟨̃ϕ̃_̃ĩ|̃)=1/1+ϵ_i^2(ρ|ϕ_i⟩⟨ϕ_i|)+ϵ_i^2/1+ϵ_i^21/d.


As for the procedure to reconstruct a certain density matrix d-dimension ρ, the trace distance is employed to show the performance of protocols under the above error assumption. Specifically, 

    ||ρ-σ||_2=[(ρ-σ)· (ρ-σ)^†]=∑_i,jξ_ijξ_ij^⋆,

where σ is the reconstructed matrix and ξ_ij=ρ_ij-σ_ij and i,j ∈[1, d].
Additionally, ϵ_i are assumed to be at the same level, i.e, ϵ.

Therefore, for diagonal elements, |ρ_ii| and its measured deviation,

    ∑_iξ_iiξ_ii^⋆   =   ∑_i|-ϵ^2/1+ϵ^2ρ_ii+ϵ^2/1+ϵ^21/d|^2 
       ≤    2 ∑_i|ϵ^2/1+ϵ^2ρ_ii|^2+ 2∑_i|ϵ^2/1+ϵ^21/d|^2 
       ∼   𝒪(ϵ^4)

Here, the coefficients ahead are ignored as we assumed ∑_i |ρ_ii|^2 is bounded.
For off-diagonal elements, 

    ρ_jk   =   (ρ|k⟩⟨ j|)
       =   (ρ|ϕ_jk^+⟩⟨ϕ_jk^+|)-i(ρ|ψ_jk^+⟩⟨ψ_jk^+|)-1-i/2(ρ_kk+ρ_jj),

where notations in Eq.(<ref>) are listed in main text.
As with multivariable derivative formula, |ξ_jk| is bounded, where 

    |ξ_jk|^2   ≤   Δ^2 (ρ|ϕ_jk^+⟩⟨ϕ_jk^+|)+ Δ^2(ρ|ψ_jk^+⟩⟨ψ_jk^+|)+|ξ_kk|^2+|ξ_jj|^2,

and 

    Δ^2 (ρ|ϕ_jk^+⟩⟨ϕ_jk^+|)≤ 2× (ϵ^2/1+ϵ^2(ρ|ϕ_jk^+⟩⟨ϕ_jk^+|))^2+2× (ϵ^2/1+ϵ^21/d)^2.

Accordingly, an approximate error is evaluated, 

    ∑_j≠ k|ξ_jk|^2 ≤𝒪(ϵ^4){∑_j≠ k(ρ|ϕ_jk^+⟩⟨ϕ_jk^+|)^2+∑_j≠ k(ρ|ψ_jk^+⟩⟨ψ_jk^+|)^2}+𝒪(ϵ^4)

For normalization condition,  ∑_j≠ k(ρ|ϕ_jk^+⟩⟨ϕ_jk^+|∼𝒪(1). In summary, 

    ∑_j≠ k|ξ_jk|^2 ∼𝒪(ϵ^4).

Specifically, under a random error assumption with the same error strength, the error of the protocol, which is expressed as a distance of measurement,  is in a higher-order formation as with respect to individual measurement devices. 

The above analysis is under the assumption that basis states have a discrepancy ϵ. However, we didn't consider the size of the device. The more qubits involved, the basis states would be less accurate. As 𝒪(poly(n)) quantum gates are required to implement specific measurement operators. With a reasonable assumption that each element gate is with a bound error ε, ϵ∼𝒪(poly(n))ε). Accordingly, the total error caused by measurement setups of our protocol is in polynomials with respect to each individual quantum gate and the size of the system. 




 


 




  
