Consider the following general formulation of the standard quantum state tomography problem: given an unknown state $\rho$, a set of (known) observables $\{\mathcal O_k\}_k$ (generally the elements of some POVM), and a corresponding vector of measured probabilities (or more realistically, frequencies) $\mathbf Y$, we want to retrieve a description of $\rho$, that is, the coefficients of the decomposition of $\rho$ with respect to some "canonical" operatorial basis $\{\sigma_k\}$ (the typical example being the basis built with Pauli matrices).
This amounts to solving the linear problem $\mathbf Y = \mathbf X\boldsymbol \theta$ for $\boldsymbol\theta$. Here, $\mathbf Y$ vector of measured frequencies, $\mathbf X$ the matrix whose elements are the coefficients of the decomposition of the observables $\mathcal O_k$ in terms of $\sigma_j$, and $\boldsymbol\theta$ the vector of coefficients obtained decomposing $\rho$ with respect to $\{\sigma_j\}$ (this notation is from (Granade et al. 2017)).
The naive solution to this linear problem is $\boldsymbol\theta=\mathbf X^+\mathbf Y$, with $\mathbf X^+$ the pseudoinverse of $\mathbf X$. However, this method is known to be problematic. For one thing, it is not guaranteed to produce a positive semidefinite estimate for the state. As mentioned by (Granade et al. 2017), possible workarounds include performing constrained least squares, or using a "two-step approach that outputs the closest physical state to a given matrix".
I'm aware of the plethora of alternative approaches to state tomography. However, I'm specifically looking for references discussing the issues with this "naive" linear reconstruction method, and in particular the related problems of numerical (in)stability and lack of positive semidefiniteness of the estimate. The only thing I found was the brief mention in the paper above, the discussion in (Qi et al. 2013), and some discussion of the numerical stability problems in Appendix A of (Opatrný et al. 1997).