A modified least squares-based tomography with density matrix perturbation and linear entropy consideration along with performance analysis

Quantum state tomography identifies target quantum states by performing repetitive measurements on identical copies. In this paper, we have two key contributions aimed at improving traditional post-processing computational complexity and sample complexity of quantum tomography protocols. In the first case, we propose a new low-cost positivity constraint method based on density matrix perturbation after the least squares (LS) estimation of the density matrix. In the second case, we propose a new cost function with the maximum linear entropy and LS method to improve the sample average trace distance with reasonably low sample complexity. We call it the LS with the maximum entropy (LSME) method. Our proposed algorithm does not follow the iterative optimization technique, which is true for existing maximum likelihood and entropy-based ones. Performance analysis is conducted for our proposed methods by studying how they compare to the existing techniques for different sample complexities and dimensionalities. Extensive numerical simulations have been conducted to demonstrate the advantages of the proposed tomography algorithms.


Introduction
Quantum computation emerges as a new computing paradigm with the promise of an exponential speed-up over a set of classical algorithms. However, one of the challenges is error-free computation. Quantum computers are more noise prone than their classical counterparts. Therefore, verification of the prepared state, quantum processor, and measurement devices are essential to achieve this fault-tolerant quantum computer.
Quantum state tomography (QST) is one of the key ideas to characterize the target quantum state by measuring multiple copies of the same quantum state [1]. Many reliable methods of QST such as maximum likelihood (ML) based QST [2], linear regression estimation (LRE) with least squares (LS) [3], Bayesian [4], adaptive Bayesian [5,6] and self-guided state estimation [7] have already been exploited. However, reconstruction methods have advantages and disadvantages [8,9], and it is necessary to find the best-suited technique depending on the context of the quantum system. The ML QST suffers a computational disadvantage and slow convergence issues. Recently, several proposals have been approached to counter the problem of slow convergence [10][11][12]. Moreover, the ML estimation returns a point estimate that sometimes may not contain meaningful information like the quantum state's region estimate [13]. Another disadvantage of ML estimation is the spurious zero eigenvalue estimation [14]. On the contrary, LRE-based QST requires sufficient measurement data to reconstruct a quantum state accurately. The reconstructed density matrix may not be physically valid if the sample complexity is also less. Moreover, ML and LRE-based estimation produces systematic biases, which may cause erroneous estimation of important properties like fidelity and entanglement [15]. Though work in [16] uses compressed sensing with incomplete measurements, it is used for low-rank density matrices estimation. On the other hand, Bayesian state tomography excels in performance but demands certain prior knowledge about the quantum state, though randomized bench-marking [17] counters some of this imprecise prior. QST techniques based on machine learning algorithms promise an advantage over the traditional QST algorithms as they can extract information more efficiently and accurately [18,19]. However, they require large data sets for training and the accuracy may depend on this heavily and may not be suitable for low sample complexity QST.
In this paper, we have two key motivations, which exploit the idea of density matrix perturbation and the maximum entropy-based LS method for quantum tomography. The conventional approach to ensure a physically reliable state is to perform the eigenvalue decomposition (EVD) of LS reconstructed density matrix [20], which is computationally expensive for higher dimensions. With the perturbation of the reconstructed density matrix, which avoids EVD, we propose a low-cost solution to mitigate the problem of physically infeasible quantum states that arises from the uncertainty of quantum mechanics.
On the other hand, the maximum entropy approach [21] is useful when the data are not sufficient. Inaccurate statistics may occur because of either information-wise incomplete measurement operators [22] or the unavailability of sufficient identical copies. Quantum measurements inherently introduce uncertainty following Born principle [23], which may induce unwanted errors. For example, quantum tomography is performed over N identical copies with M positive operator-valued measures (POVM). The number of identical copies (N) used for performing QST is known as the sample complexity of the QST. However, the number of POVM element M increases with the Hilbert space dimension for informationally complete tomography. Consequently, we also need an increase in N [24] to obtain an accurate estimate of measurement probability. In conclusion, in higher dimensional scenarios, both N and M increase and the ratio N M could be critical to obtain accurate measurement statistics. For the above-mentioned difficulty, the maximum entropy approach can be crucial in performing optimum quantum state estimation for higher-dimensional QST in low sample complexity. ML with maximum entropy (MLME) based quantum state reconstruction method has been exploited in [25], which relies on the iterative optimization procedure. However, this algorithm is computationally inefficient as it is iterative. In this paper, we propose a LS with maximum entropy (LSME) based QST method to reduce the complexity, which shows little accuracy degradation with respect to its MLME counterpart.
Contributions: With the above background, our contributions to this work are summarized as follows • Existing linear regression-based QST may fail to represent a valid quantum state. This demands classical post-processing to ensure positivity of the density matrix by using EVD, which is computationally expensive. In this work, we propose a modified version of the classical post-processing to ensure positivity of the density matrix. Our approach is based on perturbation of diagonal elements of the reconstructed density matrix instead of performing EVD-based methods to reduce the complexity of classical post-processing. The proposed perturbation uses the diagonal dominance property of a matrix to maintain positivity of the eigenvalues. And this is the key novelty of this proposal. • Quantum entropy quantifies the amount of randomness inside a quantum state. This becomes more prevalent in the case of low sample complexity. LS-based tomography algorithm does not ensure optimized entropy. Therefore, we propose a new cost function, also called LSME, that linearly combines the maximum entropy and LS criteria. The proposed algorithm will simultaneously minimize the LS criterion and maximize the entropy. We obtain a closed-form expression for the density matrix with this method. Our proposed LSME method obtains a better accuracy in state estimation compared to the standard LS one with both low and high sample complexity scenarios. • Evaluation of various analytical aspects of our proposals has been done. We have analyzed the computational complexity advantage for our first proposal to ensure positivity of the density matrix and provided a bound on eigenvalues of the perturbed density matrix. Extensive numerical simulations are performed to study accuracy of the proposed LSME method, which shows very little degradation in performance with respect to the MLME method. The numerical simulations also demonstrate that our proposed method may obtain almost the same accuracy with low sample complexity compared to the EVD-based method. Its accuracy is also comparable in the mixed state case.

Proposed linear regression-based QST estimation with eigenvalue perturbation
The first proposal discusses the method to reduce complexity of classical post-processing for linear regression-based QST. We illustrate a mathematical framework of the proposed method as follows.

Background
Let us assume that the target quantum state ρ has dimension d and it can be parameterized completely in terms of (d 2 − 1) traceless Hermitian operators, σ i for i = 1, 2, . . . d 2 − 1. Each σ i has the following properties [26] Tr(σ i σ j ) = 2δ ij (1) where i ∈ [1, d 2 − 1] denotes the indices of the operator, δ is the Kronecker delta function and Tr(.) is the trace operation of a matrix. A quantum density matrix can be decomposed by a set of traceless Hermitian matrices (σ i ) as follows [27] where the scalar coefficient θ i is defined by Tr(ρσ i ) and I is the d × d identity matrix.
Standard full quantum tomography is performed by POVM that has M elements {M j } M j =1 which can be decomposed by σ i as follows where Tr(M j σ i ) and f 0j = Tr(M j × I) = Tr(M j ) is the coefficient corresponding to the identity operator for the jth POVM element. The probability of observing a POVM outcome corresponding to the jth element M j is defined as Following equation (5), we can find the linear model of the measurement probability as follows Let us assume thatp j represents the observed value of p j after measurement andp = [p 1 ,p 2 , . . . ,p M ] T . Linear regression-based quantum state tomography is the method of estimating parameters θ i for i = 1, 2, . . . , d 2 − 1 and reconstructing the quantum density matrix by solving a set of linear equations described in equation (3). It is represented as where f 0 ∈ R M×1 is a column vector containing the f 0j coefficients, θ is the coherence vector [θ 1 , θ 2 . . . , θ d 2 −1 ] T . Also z is the observation noise vector, whose jth element is defined as z j =p j − p j . The matrix H contains all the f ij coefficients and can be represented as follows Following equation (7), the LS estimate ofθ is defined as follows [3] where y = (p − f 0 )d d − 1 and H T is the transpose matrix of H. We define ρ ls as the reconstructed density matrix using the LS method.
One of the drawbacks of the linear regression-based method is that it does not guarantee to return a non-negative density matrix [26]. A suitable method for ensuring the positivity is based on solving and thresholding of eigenvalues [20] given as reference in algorithm 1. However, the optimization becomes a computationally complex problem for large dimension having a complexity O(d 3

Proposed density matrix perturbation based linear regression
We propose a matrix perturbation-based method to ensure positivity of the LS-based reconstructed density matrix without any EVD.

Lemma 1. Suppose that A = [a ij
] is an n × n Hermitian and strictly diagonally dominant matrix, where a ij denotes an element of A corresponding to its ith row and jth column. If a ii ⩾ 0 for all i = 1, 2, . . . , n, then A is positive semi-definite [28].
A is said to be a diagonally dominant matrix if the following condition is satisfied: Moreover, A can be called strictly diagonally dominant when We propose that the positivity of the density matrix is obtained by perturbing it with a diagonal matrix D based on lemma 1. The task is to find a real d × d diagonal matrix D = diag(ϵ 11 , ϵ 22 , . . . ϵ dd ) such that where ρ po is the proposed reconstructed density matrix using eigenvalue perturbation. Our objective is to ensure that the estimated matrix ρ po is positive semi-definite. Therefore, we propose to ensure that ρ po is a diagonally dominant matrix with non-negative diagonal elements by adding a perturbation matrix D to ρ ls . Then, the solution of (12) can be obtained using an optimization framework as follows where i and j are the row and column indices, respectively. In equation (13), the constraint ensures the diagonally dominant criterion. Moreover, minimization of |ϵ ii | ensures that it does not perturb the diagonal element when the constraint is satisfied for ρ ls ii in equation (13).

Proposed optimization algorithm
Our objective is to find a suitable parameter ϵ ii for each row of the reconstructed density matrix such that the corresponding perturbed density matrix becomes diagonally dominant with the constraint that every diagonal element remains non-negative. Let us define T h i as follows We consider the values of ρ ls ii for three different sub-cases and analyse the optimum value of ϵ ii for each of the following.

If ρ ls
ii < 0, then T h i < 0 is followed. As a consequence, the minimum value of |ϵ ii | occurs for (13) is satisfied for ϵ ii ⩾ 0. For this scenario, the value of T h i will be T h i ⩾ 0. Hence, the minimum value for |ϵ ii | occurs for ϵ ii = 0.
The steps are summarized in algorithm 2.
The proposed algorithm promises computational advantage over algorithm 1 in [20] at ensuring positivity of the density matrix by avoiding eigenvalue computation. We provide an analytical bound on the perturbed eigenvalues of the reconstructed positive semi-definite density matrix and analyse the computational advantage in section 4. However, this computational advantage may come as a trade-off between computational complexity and accuracy.

Proposed LSME method
Quantum entropy quantifies the amount of randomness in a quantum state. The entropy of an ensemble of quantum states can be identified from the spectral decomposition of the density matrix. In this work, we propose a new cost function with a combination of least squares and maximum entropy.

Background
One of the key metrics of quantum entropy is the von Neumann entropy [29], which is the analog of classical Shannon entropy in quantum information defined by where e i is the ith eigenvalue of the density matrix. From equation (15), we can conclude that von Neumann entropy is a concave function of eigenvalues of the density matrix [30]. For clarity, a function f c (x) becomes concave on an interval, if, for any two points x 1 and x 2 within the interval and for any value of a parameter 't' within 0 ⩽ t ⩽ 1, the following inequality holds [31] The purity of the quantum ensemble can be quantified in terms of the coherence vector θ. However, according to the definition of equation (15), estimation of eigenvalues of density matrix is necessary to obtain the von Neumann entropy. Finding the eigenvalue requires diagonalization of this matrix, a numerically difficult process for higher-dimensional quantum states having the computational complexity of O(d 3 ).
The critical aspect of our proposal is to utilize a linear entropy of the quantum state and decompose it in terms of the coherence vector θ. Quantum linear entropy (E l ) can be obtained by approximating the natural logarithm (ln ρ) with the first order term (ρ − 1) of Mercator series in equation (15). Therefore, the E l is defined as follows where the 2-norm of θ is defined as The linear entropy requires less computational complexity compared to the von Neumann one, as it does not involve eigenvalue estimation or diagonalization of the density matrix. We propose to use E l to compute a quantum state's entropy for higher-dimensional Hilbert space. From the definition of ∥θ∥ 2 2 given above, it needs (d 2 − 1) multiplications and almost equal number of additions. Hence. the computational complexity of E l is in O(d 2 ) after the estimation of θ.

Proposed modified cost function
The existing MLME approach [25] maximizes both the likelihood function and the von Neumann entropy simultaneously. However, the optimization procedure in [25] is an iterative one and it requires to calculate the logarithm of the density matrix. This approach becomes computationally non-viable as the dimension of this matrix increases. With this background, we define a QST algorithm that combines the maximum entropy approach with LRE-based cost function. For this QST protocol, the modified cost function with the optimization framework is defined as arg θ min ∥y − Hθ∥ where γ is a tuning parameter. However, the process to find the optimal values for γ has been discussed in [25]. The optimal γ depends on the number of identical copies available for measurement. In section 5, we discuss the dependency on γ with N in our simulation. The eigenvalues of the density matrix for two-dimensional Hilbert space are For a higher dimensional density matrix, the decomposition of e i in terms of θ becomes a computationally complex and iterative process. To alleviate this issue, we propose to solve the optimization problem with the linear entropy-based cost function with respect to θ using (17) as followŝ Proposition 1. The closed-form solution of the LSME method for the density matrix, defined as ρ le , can be expressed as (20) is convex over θ for γ ⩾ 0. To show this, we consider the second-order derivative (Hessian) of C(θ) with respect to θ. The Hessian matrix is as follows For a non-negative value of γ, the Hessian matrix is always a positive semi-definite matrix as it is a sum of two different positive semi-definite matrices [31].
As the optimization framework is destined to be a convex optimization problem, the optimum value of θ le can be obtained by solving the equation ∇C(θ le ) = 0. Therefore, we obtain as From (23), we obtain theθ le as followsθ From (24), we can reconstruct the density matrix from the optimized value ofθ le .
The LSME method does not ensure positivity of the density matrix due to measurement uncertainty, as discussed previously. To ensure the positivity constraint of the density matrix, we may augment our proposed low complexity framework as described in section 2. The steps of the proposed algorithm are summarized in algorithm 3.

Bound on eigenvalues of the perturbed positive semi-definite density matrix
In the first contribution of this paper, we propose a low-complexity method to ensure the positivity constraint for the reconstructed density matrix. The framework based on (13) is reiterated as follows where T ρt is the normalization factor defined by . Dividing by T ρt in equation (25) is done to ensure the unit trace property of the density matrix. The diagonal element perturbation (prior to the normalization step) is done as follows where ρ po N = T ρt ρ po is the non-normalized density matrix and ij denotes the row and column index. Let us assume that we need to change the k diagonal elements of the reconstructed density matrix to satisfy the positivity constraint. For this purpose, we perturb ρ ls with a diagonal matrix D having k non-zero diagonal elements. Here, k is an integer random variable taking value in [0, d]. From (25), T ρt can be calculated as follows where d t is the trace of the diagonal matrix D. Let us assume that the perturbation value in the ith row of matrix D is ϵ ii and the number of perturbations occurs in k diagonal elements. For the analysis, we start with a simple model without normalization of the density matrix and we can write equation (25) as follows From equation (28), we can write the following inequality using the Bauer-Fike theorem [32].
where λ j is the jth eigenvalue of ρ po and e is any eigenvalue (without index) of ρ ls and ||D|| 2 is the 2-norm of matrix D. The matrix 2-norm ||D|| 2 can be defined as follows [32] for any vector v. Following this definition of matrix 2-norm, let us define the condition number of any matrix X as follows [32] where X is a unitary matrix such that X −1 ρ ls X = diag(e ρ1 , e ρ2 , . . . e ρ d ). To be more specific, ρ ls is a normal matrix which can be diagonalized by the unitary matrix X. From this, we can conclude K 2 (X) = 1. We can obtain eigenvalue bound from equation (29) as follows min λ j ∈λ(ρ po ) As a consequence, we can conclude that the deviation of the eigenvalues of ρ po from ρ ls depends on the operator norm of the perturbation matrix. Assuming E(.) as an expectation operation, we can calculate the expected trace of the matrix D following Bayes theorem [33] as follows where E(ϵ ii |k) is the conditional expected shift in each row of the matrix ρ ls for a particular value of k and P(k) is the probability mass function of k.
Note: E(ϵ ii |k) and P(k) may not be obtained using any standard distributions. However, they can be estimated from a simulation model when tomography is performed using real-world measurement data. Standard non-parametric estimation techniques can be used including direct arithmetic averaging or kernel density estimation (KDE) technique [34,35].

Overall error matrix with 2-norm bound
As our proposal is based on the perturbation of the matrix, the minimum difference of the eigenvalue will be a basic check for accuracy deviation using equation (32). However, it needs to capture the complete perturbation of the matrix sufficiently. We now attempt to constraint the operator norm of the error matrix arising due to the proposed algorithm. We can write it after the normalization as follows In step (a), we have used the following matrix 2-norm property for two matrices A and B as follows However, ρ ls is a Hermitian matrix and hence a normal one as well. Therefore, its eigenvectors are orthonormal. The perturbation elements ϵ ii are all positive. Assume that e ls i is the ith eigenvalue of ρ ls . Following the matrix 2-norm properties of a normal matrix [32], we can write From (34), we can write the inequality as follows Note that this bound is very loose, but it still gives a correlation with the maximum perturbation.

Bound on the eigenvalues of the perturbed density matrix with the consideration of measurement uncertainty and large dimensional scenario 4.2.1. Special case: large matrix scenario
In the previous section, we have modeled the perturbation process based on the implicit method as described in equation (25). In this section, we represent the perturbation as a linear model as follows where P t is a general perturbation matrix that ensures the positivity of ρ ls . As ρ ls and ρ po are Hermitian matrices with unit trace property, it ensures that the matrix P t is also Hermitian and traceless. If the original quantum state is ρ, then the difference with LS reconstructed state can be modeled following equation (3) as follows where θ corresponds to the original state, θ ls corresponds to the LS reconstructed state and ∆θ = θ − θ ls . Following equation (40), we can modify equation (39) as follows where ∆θ i σ i is the perturbation due to measurement uncertainty. From the definition, it is to be noted that ρ ∆θ is a Hermitian matrix. As mentioned earlier in (39), P t is also Hermitian. The elements of a Hermitian matrix may not be Gaussian and independent and identically distributed (IID). However, for large dimension scenario, the perturbation matrices ρ ∆θ and P t follow the same statistical result of a Wigner matrix, specifically the distribution of the eigenvalues [36]. There are many weaker definitions of this Wigner matrix that give the same distribution of the eigenvalues. Let us define the Gamma function as Γ(α, x) =´∞ x e −t t α−1 dt. (41) for the large dimension perturbation matrix can be stated as follows

Proposition 2. The expected value of the bound in
where σ 2 l is the variance of matrix P t with the following relation [36] E Tr( and λ, e represent any eigenvalue (without indexing) of the matrices ρ ls and ρ po in large dimensional scenario, respectively [36]. I 1 , I 2 can be calculated as follows Proof. According to the Bauer-Fike [32] theorem, from (39) we can write the following Hence, from (29), we obtain the expected value of the bound for large dimension matrix after removing the index of λ j [36] as follows where p(ν) is the probability density function (pdf) of maximum eigenvalue of a large dimensional P t matrix.
In step (a), we have used the matrix p-norm property from (35). In step (b), we have replaced the matrix 2norm with the maximum eigenvalues of the corresponding matrices. In step (c), we have used the definition of expectation E(.) from [33]. For Wigner ensemble, it follows the Tracy-Widom distribution [37,38], which can be approximated as follows p(λ) has another component around mean 2σ l defined in a fragile region within 2σ l ± δ p , with δ p being an extremely small positive-valued quantity. The probability will be negligible in this region, and we ignore this term. We now calculate´λ 2 p(λ) d 2 dλ over both the tails of the distribution. The integration with the left tail can be computed as where we use´t m exp(−βt n )dt = − Γ(γt,βt n ) Similarly, for the right tail, we obtain the integration as follows Therefore, the bound can be calculated by combining (50) and (51) in (48).
The average measurement error for the jth POVM outcome can also be expressed in terms of E|∆θ i | as follows where in step (a), we have used the representation of M j from equation (4) and used the property from equation (1). Therefore, we can conclude that measurement error is manifested through E|∆θ i | and the overall error in the reconstructed density matrix is linearly affected by this measurement noise as well.

Bound on the trace distance metric for large matrix scenario
We analyze the bound on the trace distance metric for a large matrix scenario in our error analysis. The trace distance between two quantum density matrices ρ 1 , ρ 2 is defined as follows [39] T(ρ 1 , where ν i is the ith eigenvalue of the error matrix (ρ 1 − ρ 2 ).

Proposition 3.
The expected trace distance, which is the error metric of our proposed QST algorithm, is calculated for a large matrix scenario as follows where σ l is the variance of the error matrix element as defined earlier as well.
Proof. In our model, we obtain the estimate ρ po as outlined in (39). The trace distance depends on the error matrix. Let us assume that ν represents any eigenvalue of the large matrix (ρ ls − ρ po ) where indexing (with i) is omitted because all eigenvalues follow the same empirical distribution defined as p(ν) = 4σ 2 l − ν 2 2πσ 2 l [36]. Therefore, we model the trace distance per dimension of the matrix as follows where in step (a), the distribution of p(ν) is used. For a large matrix scenario, the eigenvalues of the error matrix follow the distribution of those of the Wigner matrix. With a change in variable u ≜ 4σ 2 l − ν 2 , we can calculate the integral as follows Therefore, the expected value of the trace distance per dimension is bounded by 4σ l 3π .

Complexity analysis for positivity constraint
The proposed algorithm 2 of matrix perturbation ensures computational complexity advantage for classical post-processing, which ensures positivity of the density matrix.

Computational complexity of the proposed algorithm 2
We give a brief about the computational complexity of the existing algorithm 1 first. The eigenvalue thresholding-based method, described in algorithm 1 has an initial complexity of O(d 3 ) for performing eigenvalue decomposition of the d × d density matrix, which needs higher computational complexity. Considering all the steps, the existing algorithm 1 has complexity of 2O(d 3 ) + 2O(d) . In the proposed method, the complexity advantage comes from avoiding eigenvalue decomposition. The analysis of computational complexity is given step-by-step.
1. We need to check the d independent inequalities defined by (13). Initially, T h i from (14) is computed for each row for reconstructing the density matrix. As a comparison to the algorithm 1, we do not need to compute the eigenvalues of the density matrix. Therefore, our proposed algorithm has a complexity of O(d 2 ). In conclusion, the proposed algorithm obtains the computational complexity advantage due to the avoidance of eigenvalue decomposition. In table 1, we provide a comparison of computational complexity between existing and our proposed algorithm 2.

Complexity comparison between MLME and proposed LSME
We discuss the computational complexity of MLME and the proposed LSME-based algorithms below.

Complexity of MLME
MLME is based on an iterative optimization algorithm. MLME algorithm contains two crucial operations [25] at each iteration. The first operation is to obtain a matrix, which requires computing the logarithm of the density matrix. For this, we need to perform EVD of the density matrix having complexity O(d 3 ). The second operation is to update the density matrix based on the matrix obtained in the previous step. This step has the complexity in O(2d 3 ) as it requires computing matrix multiplication of three d × d matrices. So the total complexity of the algorithm is in O(d 3 ). In conclusion, if we need to run k e iterations to obtain the optimum solution, the total complexity becomes in O(k e d 3 ).

Complexity of LSME
As discussed in algorithm 3, the proposed LSME only needs to perform matrix multiplication and inversion (step 5). This procedure has computational complexity in O(d 3 ). If we follow the proposed algorithm 2 to ensure the positivity of the density matrix, it requires additional complexity of O(d 2 ) as discussed in the previous analysis. So the total complexity of the LSME algorithm is in O(d 3 ). The existing MLME has a complexity of O(d 3 ) for each iteration. With the increase in the iterations k e , the complexity of MLME increases. The overall complexity comparison is given in table 2.

Results and discussion
The novelty of QST algorithms can be quantified by the complexity of the algorithm and the deviation of the reconstructed state from the target quantum state. We have considered average trace distance as the key error metric [8] in the simulation. We have also considered the Frobenius norm of the error matrix as another error metric. We performed the simulation in MATLAB 2022b platform. The random density matrices are generated following algorithm 4. Some key commands are given. For example, rand(d, d) is used to generate a d × d matrix with every element drawn from the uniform distribution in the interval (0,1) and qr(X) is used as the QR decomposition of a matrix 'X' .

Algorithm 4.
Generating random quantum states using MATLAB microcode.

1: Input:
We reconstructed a random set of quantum states without prior knowledge and evaluated the accuracy of both of our proposed algorithms against the existing ones. We have used symmetric informationally complete measure (SIC-POVM) with the minimal d 2 POVM elements to perform the tomographic measurement on the identical copies of quantum states. A fiducial vector can generate a particular set of SIC-POVM elements as described in [40]. We have extended our result to d = 8 for the comparative analysis with the different dimensional scenarios. In our observation, we use the sample average trace distance error metric. The term sample average indicates that we perform QST for random quantum states and then evaluate the average deviation for each reconstructed quantum state. We have generated non-identical representation of target quantum states (ρ t1 , ρ t2 , . . . ρ tw ) by following algorithm 4. Next, we perform the proposed QST techniques on each of these prepared states, and the corresponding set of reconstructed states are (ρ t1 ,ρ t2 . . .ρ tw ). Finally, we investigated the trace distance between the prepared and reconstructed states through our simulation platform. For such purpose, we define the sample average trace distance (S a ) as follows S a = T(ρ t1 ,ρ t1 ) + T(ρ t2 ,ρ t2 ) + · · · + T(ρ tw ,ρ tw ) w . (57)

Accuracy of the proposed algorithm 2
For the proposed algorithm 2 on the positivity of the reconstructed density matrix, we achieve an advantage over computational complexity in the classical post-processing stage as analyzed in section 4. However, ensuring that it has been achieved with negligible degradation of error metrics is equally important. In the simulation for a qutrit state reconstruction, we have carried out a comparative analysis of the sample average trace distance for different sample complexities (N). From figure 1, we observe that our proposed method obtains nearly comparable accuracy concerning the existing algorithm 1 [20] at low sample complexity, while its performance is slightly degraded at higher sample complexity. However, with increased sample complexity, the sample average trace distance metric decreases with similar characteristics as the existing method. For example, with N ⩽ 10 3 , the S a is lower for the proposed method. As the sample complexity increases, the sample average trace distance metric follows a similar trend for the proposed algorithm 2 compared to the existing EVD-based LS method. However, the proposed algorithm performs slightly worse than the existing one with the larger sample complexity. This is evident for N ⩾ 10 3 in figure 1. This section aims to identify how the error sample points are scattered around the mean S a . As S a is the sample average of the trace distance errors, then each T(ρ t i ,ρ t i ) is the error sample point over which the average is calculated as defined in equation (57). For clarity, the trace distance T(ρ t i ,ρ t i ) is the ith error sample point between the ith target state ρ t i and its reconstructed versionρ t i via the proposed methods. For  this purpose, we use the 'swarmchart' function in Matlab platform [41], which helps to visualize the kernel density estimate of error sample points for multiple groups of data through a violin plot [42]. The kernel density estimate is one procedure for estimating the probability density function from sample points [43]. The violin's width represents the frequency of the error sample points. In the violin plots, the x-axis represents various groups under which the y-axis data are generated. It also corresponds to a group name, not necessarily a numerical value. In the present work, multiple groups correspond to Hilbert space dimensions (d) or sample complexity(N) as depicted in the corresponding figures. In figure 2, the violin plot compares the kernel density estimate of error sample points for sample complexity N = 100, 400, 800. The green dot in figure 2 corresponds to the sample average trace distance for the proposed algorithm 2, whereas the blue dot corresponds to the sample average trace distance for the EVD-based method. This notation has been followed for the subsequent section. It is observed in figure 2 that at N = 100, the green dot is around S a = 0.21, while the blue dot is around 0.24. However, at N = 800, both of them are almost identical. This asserts that for low sample complexity, our proposed algorithm 2 excels in performance compared to the EVD-based LS method, while, in the higher complexity case, algorithm 2 performs similarly or exhibits a slight degradation. We have already observed this in figure 1.
A similar study has been performed for higher dimensional scenarios in figure 3 with d = 6 for N = 800, 1200, 1600. For a lower sample complexity scenario (N = 100 for d = 3 and N = 800 for d = 6), the errors across random sample points have relatively equal kernel density around the mean S a for both the existing and proposed method. However, with higher sample complexity (N = 800 for d = 3 and N = 1600 for d = 6), error samples are more scattered with higher trace distance metrics for the proposed method. It is observed that at N = 800, the S a value is around 0.2 for our proposed algorithm 2 (green dot), while it is 0.24 for the EVD-based LS method (blue dot). For the larger N value, the difference between the green and the blue dot is negligible.
Discussion: With low sample complexity, the measurement statistics are inaccurate. The existing LS estimator may return a density matrix with negative eigenvalues for such a scenario. The existing EVD-based post-processing step sets these negative eigenvalues to zero to obtain the positivity constraint [20]. To preserve the unit trace property, it perturbs the small, non-negative eigenvalues and sometimes sets them to zero as well. Therefore, the eigenvalues of the reconstructed density matrix will be erroneous. This affects the trace distance error metric. However, our proposed algorithm directly manipulates the elements of the estimated density matrix to ensure positivity and does not necessarily set the negative ones to zero. This results in a performance advantage for extremely low sample complexity. This is also reflected in figures 2 and 3. In contrast, a higher sample complexity scenario ensures more accurate measurement statistics. Consequently, such a scenario decreases the chance of obtaining the negative eigenvalues for ρ ls . However, a few cases may occur where the LS-based reconstructed density matrix is positive semi-definite but not strictly diagonally dominant. In these few exceptional cases, our proposed method still perturbs the diagonal elements following algorithm 2, leading to state estimation inaccuracy. This is reflected in figure 2. The error sample points are slightly more inaccurate for the proposed method than the existing method. However, the trace distance performance difference is minimal.
We further evaluate the accuracy of our proposed algorithm 2 across different Hilbert space dimensions. In figure 4, we have plotted the trace distance of error sample points with the variation of d. To evaluate the accuracy of the estimate, we have kept the sample complexity N = 1000 for each of the different Hilbert space dimensions. At d = 3, S a for our proposed algorithm 2 has a higher numerical value than the existing method. The density of the error sample points is also scattered towards higher trace distance metrics. However, our proposed method obtains a much better accuracy than the existing one at d = 7. The key objective is to explore the novelty of our proposed algorithm 2 with the increase in dimension while neutralizing the effect of sample complexity over dimension. For such reason, we have plotted a similar figure while keeping N M ratio constant for each of the dimensions in figure 5, where M is the number of POVM elements. The difference between the S a values of the existing method and of algorithm 2 is minimal for different Hilbert space dimensions with a constant N M = 70. This is reflected in figure 5. However, error sample points are scattered towards the higher trace distance metric for our proposed method.
Discussion: This indicates that our proposed algorithm 2 returns a slightly inaccurate estimate as the Hilbert space dimension increases. From the observations, we can conclude the following.
1. N = 1000 is relatively high sample complexity for d = 3 scenario. Consequently, S a for our proposed method is more than the existing EVD-based one, as depicted in figure 4. 2. If we increase the Hilbert space dimension, the number of POVM elements increases. If we keep N constant, the N M ratio degrades. This leads to inaccurate measurement statistics. In such a scenario, our proposed method obtains a better estimate. This is reflected at d = 7 in figure 4. 3. The loss of accuracy due to perturbation in diagonal elements increases with Hilbert space dimension. A higher dimension increases the chances of more such perturbations due to an increase in the number of matrix elements. This results in more scattered sample points, as depicted in figure 5. However, our proposed method has a computational complexity advantage, as explained in the previous section.
We have explained the variation of the sample average trace distance metric with N and d in our simulation model. However, it is equally important to evaluate the accuracy of our proposed algorithm 2 for different eigen-spectrum scenarios. In our simulation model of d = 8, we have performed the post-processing part of the QST based on the proposed algorithm 2 and existing algorithms. The key motivation for the proposed matrix perturbation-based method is to achieve the positive semi-definiteness of the density matrix with a low-cost approach. However, this may result in accuracy loss in some cases compared to the existing methods. We have analyzed the accuracy with the Frobenius norm [32] of the error matrix. The Frobenius norm of a matrix B m×n is defined as follows where b ij is the element of B corresponding to ith row and jth column. In fact, we notice that for a nearly maximally mixed state, our proposed method (algorithm 2) obtains a better accuracy for the reconstructed density matrix in terms of the Frobenius norm of the error matrix ∥ρ −ρ∥ F as shown in figure 6. However, as the purity of the matrix increases, the error matrix Frobenius norm increases with the proposed solution, as shown in figure 7. Discussion: The proposed algorithm is based on diagonal element perturbation followed by normalization. The diagonal perturbation replaces zero eigenvalues with a small positive value. In the case of more purity in the state, one of the eigenvalues is very large. Therefore, the considerably small amount of additive deviation (due to the positiveness preservation of the density matrix) in the lower eigenvalues is taken from the pure state's highest eigenvalue. This leads to the more considerable Frobenius norm difference of the error matrix as shown in figure 7. However, for a nearly maximally mixed state, where most of the eigenvalues are significant, a small amount of additive deviation, even at the smaller eigenvalues, may lead to only a smaller error norm deviation compared to the existing solutions.
We appraise the accuracy of the proposed algorithm 2 through our simulation model for d = 8 with eigen-spectrum variation of the density matrix. At first, we consider a nearly maximally mixed state scenario and evaluate its accuracy. All the eigenvalues of the target quantum state ρ are almost equal for a nearly  maximally mixed state. The proposed algorithm does not suffer accuracy degradation with respect to the existing method, which is depicted in figure 8. However, we extended the simulation for different eigen-spectrum scenarios. In figure 10, ρ has five equal eigenvalues, whereas the remaining three are almost zero. A similar numerical simulation is performed where the eigenvalues are randomly distributed in figure 9. From both figures, we observe that algorithm 2 obtains a better sample average error metric. However, it faces an accuracy degradation for pure state estimation scenarios where every eigenvalue is near zero except a single one. This is depicted in figure 11.
Discussion: The existing method directly thresholds the eigenvalues. In contrast, our proposed method distributes small positive values to the zero eigenvalue indices by leveraging it from the highest one. In figures 9 and 10, although few of the individual eigenvalues may deteriorate for the existing method, we still achieve an overall advantage in trace distance metric. However, for pure state scenarios, only one eigenvalue is non-zero and the proposed algorithm will perturb most of them, which increases the trace distance metric and error of individual eigenvalues as depicted in figure 11.

Comparison with LS
As mentioned earlier, the maximum entropy approach is suitable while either number of sample copies is insufficient or measurement statistics give an error-prone estimate of outcome probability. Our proposal is based on estimating the best coherence vector θ with measurement outcomes while maximizing system's entropy jointly by tuning the parameter γ in equation (20). For the value of tuning parameter γ = 0, the optimization framework in (20) is the existing linear regression-based quantum tomography.
In figure 12, we have demonstrated the mean and standard deviation of sample trace distance error metric for different γ with qutrit (d = 3), and N = 200. The plot gives an idea of the optimal γ, which is at 0.1 in this experiment setup. It also shows how the error bar and mean increase as we move away from the optimal γ. It obtains the lowest S a at this γ value. It is also observed that the standard deviation is the lowest at this optimal γ. Figure 13 plots the S a variations for different γ values with N = 500. It is observed that S a for N = 500 shows an improvement compared to N = 200 as shown in figure 12. Our numerical result reflects that the optimum γ value occurs at γ = 0.052.

Discussion:
1. From figure 12, it is observed that the proposed algorithm alg:3 (γ > 0) obtains a better accuracy than existing linear regression-based QST (γ = 0). We also observe that there exists an optimal γ, i.e. the tuning parameter for algorithm 3 at 0.1. 2. In figure 13, the trace distance metric decreases as the sample complexity increases 2.5 times compared to figure 12. The higher sample complexity obtains an accurate estimate of measurement statistics, and the LS-based estimation obtains a good estimate. Therefore, the optimum γ decreases with the increase in N. This has been depicted in figure 13 as the optimum γ value decreases to 0.052 from 0.1. This is because as the N increases, the inaccuracy in measurement decreases, and it needs a smaller value of γ. 3. It is also noted that the complexity for performing EVD of a d × d matrix is O(d 3 ), which is independent of sample complexity N. Therefore, the cost of EVD for N = 500 and N = 200 is same.
We extend our study to observe the variation of the error sample points around the mean S a for algorithm 3 through violin plots in figure 14 for LS (γ = 0) and LSME. We have chosen LS with algorithm 2. For LSME, γ = 0.06 and γ = 0.12 were chosen for comparison, with N = 200 and d = 3. The mean values of error sample points, i.e. S a , are shown. It is seen that the difference in S a between LS and LSME with optimal     γ is almost 0.028. It is also observed that the error sample points are more widely scattered around the mean for LS than for LSME.
In figure 15, we plot the variation of sample average trace distance concerning sample complexity and γ simultaneously. We observe that the accuracy of the QST algorithm degrades, i.e. the error metric increases with the decrease in sample complexity. We also observe that there will be an optimal value of γ for a given N.
In figure 16, we demonstrate the variation of optimum γ values with the Hilbert space dimension and sample complexity variation. We observe that for a fixed Hilbert space dimension at lower sample complexity, a higher γ value is required. Another observation is, with the constant sample complexity, the optimum γ increases with the increase in Hilbert space dimension.
Discussion: From these observations, we can conclude the following 1. The measurement statistics become inaccurate with the decrease in sample complexity, resulting in the accuracy degradation of QST as observed in figure 5. As a consequence, the accuracy of the QST algorithm degrades. However, the proposed LSME algorithm achieves better accuracy for such scenarios than the existing linear regression-based method, as reflected in figure 14. 2. The proposed LSME's performance rate depends on the choice of the optimal values of γ. The optimal choice of γ depends on the sample complexity for a particular d−dimensional Hilbert space. As the sample complexity increases, probabilistic measurement values converge towards the theoretical ones as shown in figure 15. The optimal value of γ becomes very small as the sample complexity increases. Similarly, the measurement statistics also become erroneous if we increase the Hilbert space dimension while keeping the sample complexity constant. This is because an increase in Hilbert space dimension demands more sample copies for obtaining sufficient statistics. From figure 16, it is observed that the optimal γ value increases with the increase in d for a fixed N, while it decreases with the increase in N for a fixed d. Figure 17. Comparison of sample average trace distance error metric between our proposed LSME method and existing MLME method for γ = 0.006 for d = 3. The accuracy of our proposed method is comparable to the existing methods. Figure 18. Comparison of sample average trace distance error metric between the proposed method and existing methods for N = 300 with variation in γ for d = 3. Optimum γ can be different for various algorithms.

Comparison with MLME
We now compare the accuracy of MLME and that of LSME in our simulation. It is noted that LSME alone does not guarantee to return a positive semi-definite estimate of the density matrix. Figures 17 and 18 show the comparison of trace distance error metric between the existing MLME and the proposed LSME methods with the variation in N and γ, respectively. To ensure positivity, we have used the proposed LSME method with EVD and algorithm 2 separately. However, the computational complexity of LSME with EVD is higher than algorithm 3. LSME with EVD requires a post-processing complexity of O(d 3 ) whereas algorithm 3 has a post-processing complexity of O(d 2 ) following table 1. Discussion and observation: 1. From figure 17, it is observed that the LSME with EVD almost matches the MLME as the sample complexity increases. It is also observed that the LSME with algorithm 2 attains minimal accuracy degradation compared to the MLME. 2. We observe a similar performance from figure 18, with a constant N = 300. In the figure, the optimal value of γ for MLME occurs at γ = 0.01, whereas the optimal γ values are 0.08 and 0.06 for the proposed algorithm 3 and LSME with EVD, respectively as prominently marked. MLME-based estimation still obtains a better estimate than the proposed algorithm 3. However, our algorithm has a smaller computational complexity.  Our proposed low-complexity LSME method does not degrade the accuracy metric excessively compared to the existing algorithms.

Observation of various components of LSME
Our proposed LSME algorithm aims to minimize the least squares and enhance the entropy part. This can be visualized through the variation of the least squares cost function and linear entropy cost function with respect to different γ values. First, we define the reconstructed measurement probability corresponding to the M j POVM element as p re j ≜ Tr(ρ re M j ), where ρ re is the reconstructed density matrix. Let us define p re as a d 2 × 1 vector consisting of p re j values, where j = 1, 2, . . . , d 2 . The difference between the measured probability vectorp (in equation (7)) and reconstructed probability vector p re is closely related to the linear regression-based cost function. However, we combine another cost function that maximizes the linear entropy in LSME. The variation of both the cost functions with the tuning factor γ is depicted in figure 19 for a qutrit scenario. In figure 19, we observe that for γ = 0, the difference in probability statistics is minimum; however, the entropy part is also minimum. With the increase in γ, the difference between p rc and p increases, and the entropy of the quantum system also increases. For example, if the γ value is substantial, it will maximize the linear entropy quantity and shift the reconstructed state to a maximally mixed state. Following this, the |p − p re | 2 cost function will increase as the measured probability vector will have a significant difference with the reconstructed probability vector p re .
We extend our simulation result for different Hilbert space dimensions for a fixed sample complexity of N = 800. From figure 20, we observe that the proposed LSME obtains a better estimate than the existing LS algorithm in higher dimensions. The observed probability statistics are inaccurate for higher dimensions for a fixed sample complexity, and in such a scenario, our proposed algorithm obtains better accuracy as depicted in figure 20. In contrast, for lower dimensions, the measurement statistics are more accurate. As a consequence, the sample average trace distance has little difference.
Discussion: For any QST algorithm, the computational complexity of QST is also an essential factor. We have estimated the average time (over many states) required for classical post-processing for different Hilbert space dimensions. The average time is calculated over different samples generated through algorithm 4 for N = 800. From figure 21, it is observed that the proposed LSME is faster than the existing iterative accelerated projected gradient ML (APG-ML) [10] method for various Hilbert space dimensions. This reflects our claim of computational complexity advantage, as discussed in our previous complexity analysis.

Conclusion
This study aims to explore QST algorithms that improve computational complexity or perform better in lower sample complexity. The proposed method to ensure the positivity of the quantum density matrix promises improved computational complexity over the existing one in the classical post-processing stage. The second proposal combines the maximum entropy approach with LS-based criteria. As the dimension of the quantum state space increases, the number of POVM elements increases exponentially, requiring an exponentially large number of identical copies. Our proposed maximum linear entropy-based QST method achieves better accuracy in such a scenario. From the simulation results, a key observation is that the proposed algorithms may reduce the possibility of obtaining zero eigenvalues. Our future research aims to determine whether this could avoid the zero eigenvalue estimation problem. We would also like to explore the algorithm accuracy as a function of dimension and purity. Another research direction is to explore the maximum entropy-based QST to estimate the entropy, fidelity, and entanglement of quantum states rather than performing the full tomography.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.