Characterizing entanglement dimensionality from randomized measurements

We consider the problem of detecting the dimensionality of entanglement with the use of correlations between measurements in randomized directions. First, exploiting the recently derived covariance matrix criterion for the entanglement dimensionality [S. Liu et al., arXiv:2208.04909], we derive an inequality that resembles well-known entanglement criteria, but contains different bounds for the different dimensionalities of entanglement. This criterion is invariant under local changes of $su(d)$ bases and can be used to find regions in the space of moments of randomized correlations, generalizing the results of [S. Imai et al., Phys. Rev. Lett. 126, 150501 (2021)] to the case of entanglement-dimensionality detection. In particular, we find analytical boundary curves for the different entanglement dimensionalities in the space of second- and fourth-order moments of randomized correlations for all dimensions $d_a = d_b = d$ of a bipartite system. We then show how our method works in practice, also considering a finite statistical sample of correlations, and we also show that it can detect more states than other entanglement-dimensionality criteria available in the literature, thus providing a method that is both very powerful and potentially simpler in practical scenarios. We conclude by discussing the partly open problem of the implementation of our method in the multipartite scenario.


I. INTRODUCTION
In recent years an impressive degree of control over a genuinely quantum system has been reached, and quantum states are generated in more and more complex experiments, involving more and more particles. In these experiments it is very important to characterize entanglement, even quantitatively, as a major figure of merit for the nonclassicality of the correlations among the particles. In particular, the socalled entanglement dimensionality plays a prominent role, as it quantifies in some sense the dimension needed for reproducing the states' correlations among a given bipartition [1]. In fact, a high Schmidt number across many bipartitions is also the bottleneck for simulating classically many-body quantum states, for example through algorithms based on tensor networks [2]. Furthermore, high-dimensional entanglement has been proven useful for several practical tasks, ranging from improved security for quantum cryptography [3,4], to noise resistant quantum communication [5,6] and universal quantum computation [7,8]. A variety of methods for bounding the entanglement dimensionality has been developed and successfully implemented in several experiments over the past years [1,9]. Most of these methods require quite demanding measurement capabilities, which in particular cannot always be implemented in very complex experiments, for example with many-body systems.
At the same time, a powerful idea, which has been put for-ward in recent years, is based on randomized measurements, which avoids the need of carefully tuning the measurement directions and bypasses even the requirement of a common reference frame between the particles. In that approach, two parties a and b perform measurements of some observables M a and M b , which are then randomly rotated via local unitaries U a ⊗ U b . Consequently the correlations arising from such measurements are randomly distributed, and the typical theoretical approach is to investigate a sampling from the unitary Haar measure. The idea is then to estimate moments of such randomized correlations and use them to extract global information about the state. Randomized measurements have been proven useful already for a number of tasks, including the violation of Bell inequalities [10][11][12][13][14], the estimation of entanglement measures via Rényi entropies [15][16][17][18], the estimation of fidelities [19,20], topological invariants [21] and the detection of entanglement via negativity of the partial transpose [22,23], or directly with moments of randomized correlations [24][25][26][27][28][29][30][31].
While most of the works considered systems of qubits, entanglement detection for qudits of dimension higher than two has been introduced by making use of second moments of randomized correlators dU a dU b (U a ⊗ U b )M ⊗ M(U a ⊗ U b ) † 2 of certain observables M [25]. Recently, proposals to make also use of fourth-order moments have been explored [28,29], including for qudits of dimension higher than two [32]. At this point, the question arises whether the dimensionality of entanglement can be also witnessed via randomized correlators. One approach that goes in that direction has been explored in the context of average work-extraction protocols [33] and makes use of the so-called r-reduction criterion. Since we want to consider measurements averaged over all local unitaries, we have to look for criteria based on local unitary invariants, such as norms of the correlations tensor [34]. For example, an expression that seems suitable to implement and exploit both second-and fourth-order moment of the randomized correlations is the one-norm of the correlation tensor, which is bounded for separable states [32,35].
In this paper, we first observe that the one-norm of the correlation tensor has different bounds for different entanglement dimensionalities, thus generalizing the criterion in Ref. [35]. We show that this follows as a corollary of the so-called covariance matrix criterion [36][37][38], which has been recently extended to detect the entanglement dimensionality [39]. The condition we derive here can also be related to other well-known entanglement criteria, like the computable cross norm or realignment (CCNR) criterion [40][41][42] to witness the entanglement dimensionality. Second, we present a practical method that uses such local unitary invariant condition to bound the entanglement dimensionality from randomized measurements. As a case study we show how our method works in dimension three and four and discuss the potential limitations for increasing dimensions. Finally, we compare our method with other prominent criteria, like witnesses given by fidelities to target entangled states, or the two-norm of the correlation tensor and discussing the practical applications with a finite sample of randomized correlations. We conclude by summarizing our results and discussing the partly open question of extending this method to the multipartite case, i.e., to the detection of the Schmidt number across every bipartition of a multiparticle state.

A. Schmidt number and its detection
Let us consider a bipartite system with Hilbert space H = C d a ⊗C d b and a pure state |ψ ∈ H. Any arbitrary bipartite pure state can be brought to the normal form under local unitary operations |ψ → U a ⊗ U b |ψ = j λ j | j a j b , where {| j n } with n = a, b are orthonormal basis for the two local Hilbert spaces and λ j are called Schmidt coefficients. This is termed Schmidt decomposition and can be achieved essentially from the singular value decomposition of the coefficient matrix of the state |ψ = kl c kl |k, l . In other words, the transformation that brings the coefficient matrix c kl in its singular value decomposition amounts to a local unitary transformation of |ψ , which in turn provides the Schmidt decomposition. The Schmidt decomposition also gives the reduced density matrix of each party as a = tr b (|ψ ψ|) = j λ j | j a j a |, and similarly for b . Clearly, such a transformation does not map a product state into an entangled state nor vice versa. Moreover, it cannot change the amount of entanglement of the state, since all entanglement monotones are by definition invariant under local unitary transformations. In fact, the number s(|ψ ) of nonzero Schmidt coefficients of a pure state defines the entanglement monotone called Schmidt rank. This monotone can be extended to all bipartite density matrices, i.e., trace-class operators ∈ B(H), via convex-roof construction, this way defining the Schmidt number [43]: where the infimum is taken over all pure-state decomposi- Roughly speaking, this entanglement monotone is connected with the minimal dimension of the Hilbert space needed to reproduce the states' correlations. Note that other notions can be given connected to the dimensionality of entanglement, such as that of genuine multilevel entanglement [44]. Typical methods useful, especially in theory, to bound the Schmidt number from below involve either the evaluation of entanglement witnesses, i.e., hermitian operators W r such that tr(W r r ) ≥ 0 for all states r such that s( r ) ≤ r and tr(W r ) < 0 for some state such that s( ) > r, or the application of so-called r-positive maps, i.e., linear maps L r : C d a ⊗ C d a → C d b ⊗C d b such that id r ⊗L r is positive, where id r is the identity map in the space of r × r matrices [43,45,46].
The most typical Schmidt-number witnesses are operators of the type where |Ψ r = k √ λ k |kk is a given state with Schmidt rank larger than or equal to r and Schmidt coefficients ordered decreasingly [19,45,[47][48][49]. This means that the fidelity with a given target state is always upper bounded by its r greatest Schmidt coefficients for all states with the same Schmidt number, i.e., must hold for all states such that s( ) ≤ r. In particular, for a d-dimensional maximally entangled state the upper bound on the fidelity as in Eq. (3) is given by r/d, since all Schmidt coefficients are equal to each other. The most typical r-positive map used for Schmidt-number detection is given by a generalization of the reduction map [50,51], namely D r (A) = tr(A)1 − 1 r A. When this map is applied locally on a bipartite density matrix one gets and its positivity must hold for all states with Schmidt number smaller than r [43]. Note that this map is dual to the witness with respect to the maximally entangled state, namely, the witness in Eq. (2) with |Ψ r being the r-dimensional maximally entangled state is obtained from the map id a ⊗ D r after applying the Choi-Jamiołkowski isomorphism.
A further important Schmidt-number criterion can be found from the coefficients of expansion of the density matrix in a pair of local bases, namely we have = kl (X ) kl g (a) k ⊗ g (b) l . In particular, considering local bases that are composed of 1/ √ d plus normalized su(d) bases and defining the submatrix X su(d) by omitting the entries 1 (a) it is possible to find a criterion for the entanglement dimensionality that reads [34] Finally, a more general approach can be followed by starting from the (symmetric) covariance matrix associated to a pair of g = (g a , g b ), where g a = (g (a) 1 ⊗ 1, . . . , g (a) d 2 a ⊗ 1) and ) are orthonormal basis of observables for party a and b, respectively. For a general set of observables M = (M 1 , . . . , M K ) and a given state , the (symmetric) covariance matrix has components given by and for a pair of local operator bases it assumes the block form in which the diagonals γ a := Cov a ( g a ) and γ b := Cov b ( g b ) are the symmetric covariance matrices of each party, and the off-diagonal blocks are the cross covariances between the two local observables vectors, and we define the marginals k . Note that we have (v a ) k = l (X ) kl tr(g (b) l ) and similarly for party b. Based on this covariance matrix, a compact criterion that has to hold for all states with a certain Schmidt number r can be found [39] that reads where (omitting the label k for simplicity) Γ r = Cov ψ r ( g) are (positive) covariance matrices associated to pure Schmidtrank-r states ψ r = j λ j | j a j b and {p k } is some probability distribution. Here, each Γ r has in general a block form with the singular values of the blocks satisfying B. Haar-randomized correlators and Bloch-sphere integrals One potential advantage of criteria based on quantities, which are invariant under change of local su(d) bases, is that they can be evaluated from random measurements, if specifically sampled from the local unitary Haar measures. To see that let us consider correlations between generic observables M a and M b , randomized over local unitaries, namely is a local unitary randomly sampled from the Haar measure. The probability distribution of the correlations arising from such randomized measurements can be characterized by its moments The nice feature of such Haar-randomized correlators is that these moments written above in some cases can be calculated by transforming the integration over the unitary Haar measure into an integration over d-dimensional Bloch spheres, namely where σα :=α · σ and σβ :=β · σ with σ = (σ 1 , . . . , σ d 2 −1 ) being an su(d) basis and (α,β) being two real (d 2 − 1)dimensional unit vectors. Furthermore, N is a normalization that is chosen such that S (m) = 1 for pure product states. One important property that makes this possible is that both S (m) and R (m) (M a ⊗ M b ) as defined above are invariant under local unitaries. However, the value of the latter depends on the choice of the observables M a and M b , in particular on their concrete eigenvalues (not on the eigenvectors due to invariance under unitaries). On the other hand, the S (m) are actually invariant over all changes of local su(d) bases, which results in the fact that they depend only on the singular values of the correlation matrix X su(d) .
In fact, after some algebra it has been shown in Ref.

A. New criterion based on Bloch invariants
Once the relation between moments of randomized correlations R (m) and expressions like Eq. (14) have been found, one might think to use the former measurements to witness important properties of the quantum state, e.g., entanglement [32]. For this task one needs a further ingredient, namely a suitable entanglement criterion written in terms of the local orthogonal invariants i . Similarly, finding a witness of the entanglement dimensionality in terms of the i would enable the possibility of detecting the former with randomized correlators via similar methods.
In the following, we show that a condition like that can be derived from the covariance matrix criterion in Eq. (9).
Observation 1. For every Schmidt-number-r density matrix the following inequality holds which is also invariant under change of orthonormal pairs of su(d a ) and su(d b ) basis. Here trace norm is defined as tr|A| = tr Proof.-We can rewrite Eq. (9) in the form where we called A := γ a + v a v T a and B := γ b + v b v T b the singleparticle correlation matrices. Since the term we subtract is a projector, Y itself is positive and so is its principal minor where we consider two local bases composed of identity plus su(d) bases and have eliminated the first row and first column from each block. In other words, the above matrix contains only variances and cross correlations relative to the su(d a ) and su(d b ) bases for the two particles.
Multiplying this matrix Y with vectors t µ = αu µ , −βv µ , where u µ and v µ are the left and right singular vectors of X and α, β > 0 and summing over all µ we obtain α 2 tr(A su(d) ) + is the linear entropy of the single-particle reduced density matrix of a generic (optimal) pure Schmidt rank-r state (which is the same for both parties). Here we use the relations tr( . The latter can be directly verified using the generic expression of its singular values from Eq. (11). This way, we obtain (α . Now, minimizing over all α and β we have that the minimum is achieved for , whenever this is positive and results in Eq. (15). This criterion is invariant under changes of orthonormal pairs of su(d a ) and su(d b ) bases, since any such change of bases results in local orthogonal transformations of the correlation matrix, namely In the following we restrict to bipartite spaces with equal dimensions d a = d b = d, in order to discuss the implementation of this criterion in a method to witness entanglement dimensionality from the randomized correlators R (2) and R (4) . For equal dimensions, Eq. (15) reduces to Note that for r = 1 the left-hand side is exactly the same as in the entanglement criterion from Ref. [35], which was also used in Ref. [32] to find entanglement witnesses from randomized the correlators R (2) (M ⊗ M) and R (4) (M ⊗ M). Thus, our condition can be seen as a generalization of such a criterion to the entanglement dimensionality. Moreover, our condition can be compared to Eq. (5), since they are based on basically the same information. The only difference is that Eq. (5) considers a different matrix norm on the left-hand side. Note also that expressing X in its singular value decomposition X = O T a · ξ · O b corresponds to expressing the density matrix as l obtained from the initial canonical bases g (a) l and with ξ k ≥ 0. This is a Schmidt decomposition in operator space, with G a and G b being the Schmidt vectors in operator space. Now, expressing the covariance matrix in these operator Schmidt bases we can also derive the following Schmidt-number-r necessary condition: which in turn reduces to Eq. (18) in the case in which G a and G b are of the form identity plus an su(d) basis. At the same time, it is worth emphasizing that in general the singular value of the full correlations matrix, i.e., the ξ k do not necessarily coincide with the singular values of the su(d) submatrix, i.e., the k . The condition (19) can be shown from the following corollary of Eq. (9), which was presented in Ref. [39]: which is invariant under a full local orthogonal change of bases.
Proof of Eq. (19).-We have the following chain of in- Here, the first is due to |a − b| ≥ |a| − |b|, the second is obtained by substituting the operators Schmidt decomposition into Eq. (20) and using √ ab ≤ 1 2 (a + b). Finally, the last inequality is similarly obtained from a 2 + b 2 ≥ 2 |ab|.
Observe also how the expression in Eq. (19) is similar to the CCNR criterion, which can be put in the form k ξ k ≤ 1 [40][41][42], and in a sense generalizes it to detect different Schmidt numbers. Such a generalization has been already derived in Ref. [54], where the authors actually derived an even stronger criterion based on unitarily invariant norms of the density matrix.

B. Entanglement dimensionality witness from randomized measurements
With all the ingredients discussed in the previous sections we are now in the position of generalizing the method of Ref. [32] to witness the entanglement dimensionality. Thus, we consider the moments R (2) (M ⊗ M) and R (4) (2) and R (4) (M ⊗ M) ∝ S (4) . Then, varying the i and making use of the criterion in Eq. (18) we can calculate the boundary regions of Schmidt-number-r states in the plane (S (2) , S (4) ) by solving the problem namely by maximizing and minimizing the value of S (4) for each fixed value of S (2) with the additional constraints given by Eqs. (14) and (18) (plus positivity of the i ). Such boundaries for d = 3, 4 and the different Schmidt numbers are plotted in Figs. 1 and 2. It can be observed numerically that the depicted regions are fully filled with physical states (cf. Fig. 3 in appendix B). Extremal points are product or maximally entangled states of a given Schmidt rank r, which we label Ψ r + = 1 √ r r−1 j=0 | j j . In particular, the curves that distinguish the different Schmidt numbers are obtained from the minimization problem. These curves can be obtained analytically for all dimensions with the method of Lagrange multipliers and the Karush-Kuhn-Tucker conditions (cf. appendix C). The result is, for each Schmidt number r, a piecewise function with pieces labeled by an integer n = 1, · · · , d 2 − 2: f n (S (2) ), with g n (X) = n(n + 1)X − nB 2 r .

C. Detection with finite statistics and comparison with other criteria
Here we compare our newly derived entanglement dimensionality witnesses with other Schmidt-number criteria. Pre-liminarly, let us notice that the criterion Eq. (19), detects strictly more states than any criterion constructed from the fidelity with respect to a maximally entangled state, namely Ψ d In fact, the latter criterion can also be written as trX ≤ r in some canonical bases, and it is clear that trX ≤ k ξ k for any canonical pair of local bases. Furthermore, there are states that can be detected with a strictly higher Schmidt number with Eq. (19). As an example consider the mixed state W = 1 2 |Ψ 3 + Ψ 3 + | + 1 4 (|23 + |32 )( 23| + 32|), which has Schmidt number s( W ) = 3. Its Schmidt number certified by any fidelity witness (not just with respect to maximally entangled targets) is at most 2 [55], while Eq. (19) detects the actual Schmidt number 3. More in general, one can make statistics with random mixed states and observe that typically corollaries of the covariance matrix criterion like Eq. (20) detect significantly more states than any fidelity witness (3). The same holds for the criterion based on the r-reduction map in Eq. (4).
At the same time, we can observe that states like W , i.e., states that are not detected by any fidelity witness are ideally detected with the correct Schmidt number by our criterion based on randomized measurements. See, for example, Fig. 2. Similarly, we can compare our criterion with Eq. (5). The latter can be simply expressed as and thus corresponds to just vertical lines in the (S (2) , S (4) ) plane, placing per sé no constraints on S (4) . It seems thus reasonable to expect that generically our method would detect more states, given that both of them are saturated by some pure Schmidt-rank-r states. In fact, from Figs. 1 and 2 we can observe that this is the case in d = 3, 4, where there are regions of states (depicted as striped) with Schmidt number larger than what is detected by Eq. (5). From the analytic expression of the boundary curves Eq. (22) we can observe that the same happens actually for every dimension: both criteria provide boundaries that start from the points given by the states Ψ r + .
However, our criterion gives a tighter boundary in the rest of the plane. For example, our criterion has higher tolerance to white noise than Eq. (5) for detecting the maximally entangled state, i.e., the density matrix to have a given Schmidt number up to a higher p. In fact, we have tr|X su Substituting this into Eq. (18) and Eq. (5), we can observe that our criterion gives higher upper bound to p. Notice also that one could use the well-known relation |c i | 2 to derive a bound on tr|X su(d) | from the known bound on X su(d) 2 , which still turns out to be strictly weaker.
So far we have discussed what states are detected based on the exact calculations of the quantities (S (2) , S (4) ), which in principle would correspond to the moments of the randomized correlators, obtained with infinite statistics. Let us now discuss what happens when the randomized moments are evaluated from finite statistics. After a finite sample of N tot random local unitaries, we can estimate the Bloch sphere integrals aŝ where we indicate Denoting with µ(·) and σ(·) the mean and standard deviation of a random variable we have which holds whenever all the individual events are independent and distributed according to the Haar measure. Thus, with these relations, we can evaluate with how much confidence a state is detected assuming that the individual trials are independently drawn from the Haar measure. As a paradigmatic test, we consider d = 3, and Ψ 3 + as the target measured state. For this state, the ideal second-and fourthorder integrals are S (2) Ψ 3 + = 2 and S (4) We take a random numerical sample of N tot = {10 3 , 10 4 , 10 5 } local unitaries from the Haar measure, respectively, and obtain as the estimate for the moments of the randomized correlators the values µ Ŝ (2) When N tot = 10 3 , 10 4 , 10 5 , respectively, taking three-σ intervals as a measure of error we obtain and we can see that the state is detected within the error already with the sample of N tot = 10 3 points, as detailed in appendix D. Similarly, one can calculate the white-noise tolerance of the method, for example with a sample of N tot = 10 4 random unitaries. In this case we get that the state p := (1− p) Ψ 3 + Ψ 3 + + p 9 1 3 ⊗ 1 3 is detected with a Schmidt number equal to 3 with a confidence interval of 2 standard deviations for a noise fraction p = 1/4, while it is detected with a Schmidt number equal to 2 for a noise fraction p = 7/10. See also appendix D for more details.

IV. CONCLUSIONS AND OUTLOOK
In conclusion, we studied corollaries of the covariance matrix criterion for the Schmidt number [39], which in some sense resemble well-known entanglement criteria [35,41,42] and can be seen as a generalization of them to a Schmidt-number witness. Exploiting one of these corollaries, we have generalized the entanglement detection method of Ref. [32] to witness the Schmidt number. This method is based on the measurement of correlators of for a given M = diag(α + , . . . , α + , β, α − , . . . , α − ) and for random local unitaries (U a ⊗ U b ) i sampled from the Haar measure. We have shown how such a method works in practice for bipartite states of equal dimensions d = 3, 4, showing also that it detects strictly more states than the other known criteria for the Schmidt number. We have also observed how our method detects paradigmatic example states after a finite sample of random correlation measurements. Our method can be straightforwardly implemented in higher dimension, as long as the dimensions of the two systems are equal to each other. At the same time, the inequality that we have derived, Eq. (15), also works when the two local dimensions are different, which would enable the potential extension of the method based on randomized correlators to the multipartite case, where the notion of Schmidt number can be extended to that of Schmidt-number vector [56]. As we have observed, our criterion would detect strictly more states than typical Schmidt-number witnesses based on fidelities with respect to maximally entangled states, as well as criteria based on the 2-norm of the correlation matrix, Eq. (5). Its implementation with randomized measurements, however, relies on the results of Ref. [32], which pro-vides concrete matrices M such that the moments R (2) and R (4) in Eq. (12) can be evaluated from a Bloch-sphere integral (13). The question then arises how to generalize these results when the dimensions of the two parties are different, which would then allow for the generalization of our results for the detection of the full Schmidt-number vector in multipartite states. This question is left for future investigation and provides an exciting research direction for implementing our method in the study of entanglement dimensionality in manybody states. Note added: While finishing this paper, we learned that related results were obtained in Ref. [53].
Appendix A: Moments of randomized correlators and Bloch-sphere integrals In this section we summarize the derivation of the expressions in Eq. (14) for the correlators of su(d) generators averaged over the Bloch sphere and their relation with the second and fourth moments of the randomized correlators U a M a U † a ⊗ U b M b U † b . We essentially follow and summarize the results of Ref. [32] and present it also here for making our work more self-contained. Let us expand the averaged correlators as Expanding the m-th power in terms of multinomial coefficients we can express the integrand as and consequently we obtain where a k = l n kl and b l = k n kl . Now, for fixed values of a k and b l the integral on the unit sphere is easily calculated from well-known expression in terms of Γ functions. The only thing that remains is to check all possible combinations of n kl such that kl n kl = m and evaluate the integral for each such combination. In the case of m = 2 there are only two possibilities, namely one value n kl = 2 and all the rest are zero or two values are equal to 1, namely n k 1 l 1 = n k 2 l 2 = 1 and the rest are zero. The second possibility leads to vanishing integrals, so only the first should be considered. In that case one gets where Γ(·) is the Γ function and the last factor comes from the integrations over the unit sphere. Adjusting the normalization such that S (2) = 1 for pure product states we get S (2) = d 2 kl=1 (X su(d) ) 2 kl which is precisely the expression in Eq. (14). In the case of m = 4 there are several more possibilities and the calculation is quite lengthy. Therefore, we refer the interested reader to Ref. [32] for the detailed calculation that leads to the expression in Eq. (14). Note that here we use a different normalization for the correlation matrix X su(d) with respect to Ref. [32], which leads to an additional factor of d m in the expressions of the S (m) .
Let us now briefly summarize the calculations of R (2) (M ⊗ M) and R (4) (M ⊗ M) and the constraints on M such that those are proportional to S (2) and S (4) , respectively. First of all, it has been shown in Ref. [25] that the value of R (2) (M ⊗ M) is independent on the concrete eigenvalues (and eigenvectors) of M, as soon as tr(M) = 0. In particular, choosing the correct normalization tr(M 2 ) = d one obtains R (2) (2) . Concerning the fourth-order moment, one can expand the power inside the integral with the multinomial theorem. Imposing that M is traceless one obtains similarly as for S (m) : where σ k are elements of a su(d) basis. Again, for the case m = 4 one has to check and sum up all possible combinations of n kl such that kl n kl = 4. As a result, one obtains that for a matrix of the form M = diag(α + , . . . , α + , β, α − , . . . , α − ) and for properly chosen eigenvalues α ± and β the fourth-order moment is given by One concrete choice of eigenvalues is given, for the case of odd dimension d, by A similar set of eigenvalues can be provided for the case of even d. See Ref. [32] for the details.
which arises from Eq. (21) in the case of minimization. One can observe that the resulting curves are the nontrivial ones for detecting the Schmidt number (see Figs. 1 and 2 and appendix B). The Lagrangian function for the problem (C1) is where λ i ≥ 0 and ν are Lagrange multipliers. The Karush-Kuhn-Tucker conditions associated to the problem are (dual feasibility) .

(C3)
Let us now calculate the optimal points coming from these conditions. These are extremal points for the problem, not necessarily global minima. However, we see that these conditions are satisfied only by a finite number of points. Therefore, we check all of them and observe which ones provide the global minimum. Let us distinguish two cases: (i) λ 0 > 0 and (ii) λ 0 = 0. (i) Due to complementary slackness, we have i i = r − 1 d . At the same time, if λ i = 0 for some i > 0, then for the corresponding index we have i ≥ 0 and 4 3 i + λ 0 + 2ν i = 0 due to stationarity. This equation has at most two positive solutions, which we call i = s and i = l (for "small" and "large"). If instead λ i > 0, then i = 0 is a solution, and stationarity implies λ i = λ 0 . Thus, in total for the case (i) we have three possible values of each i , namely i = {0, s, l}.
Let us now check within all the possible (finite) combinations of i = {0, s, l}, which one gives the global minimum. Assume the total number of i = s to be n s and the total number of i = l to be n l . The remaining i 's are all zero. We now prove that the minimum is achieved for n s = 1. For this, we consider the following minimization problem: min n s ,n l ∈N s 4 n s + l 4 n l , subject to sn s + ln l = A > 0, First, let us express s and l in terms of the other parameters: s = A n s + n l − −A 2 n l + n s n l B + n 2 l B n s (n s + n l ) 2 , l = A n s + n l + n s n l (−A 2 +(n s +n l )B) n s (n s +n l ) 2 n l . (C5) From s, l > 0 we have Here we test our method in a paradigmatic example with a finite statistical sample of local unitaries. We consider d = 3, and Ψ 3 + as the target measured state and simulate numerically the statistics that would ideally occur with a finite sample of randomized correlators. For our target state, the ideal second-and fourth-order integrals are S (2) Ψ 3 + = 2 and S (4) Ψ 3 + = 5 3 . We take a random numerical sample of, respectively, N tot = {10 3 , 10 4 , 10 5 } local unitaries from the Haar measure and estimate the moments of the randomized correlators aŝ and for independent and Haar-distributed points we get i.e., the estimators are unbiased. Then, we can estimate the error by means of the standard deviation of the estimators, which are given by As shown in the plot in Fig. 4, Schmidt number 3 is certified in all cases. Therefore, a sample size of 10 3 is already enough. Similarly, one can calculate the white-noise tolerance of the method, for example, with a sample of N tot = 10 4 random unitaries. In this case we get that the state p := (1 − p) Ψ 3 + Ψ 3 + + p 1 3 ⊗1 3   9 is detected with a Schmidt number equal to 3 with a confidence interval of 2 standard deviations for a noise fraction p = 1/4, while it is detected with a Schmidt number equal to 2 for a noise fraction p = 7/10, as shown in Fig. 5.