Sample complexity of quantum qtate tomography implies by the Aaronson-Rothblum result According to Aaronson-Rothblum arXiv:1904.08747, shadow tomography can accurately predict $m$ two-outcome observables for a quantum system of Hilbert space dimension $d$ by measuring only $O((\log m)^2(\log d)^2/\epsilon^8)$ copies of an unknown quantum state $\rho$ of the system.
Now consider a $n$-qubit system ($d=2^n$). To reconstruct the full density matrix, we only need to know the expectation value of all the $m=4^n-1\simeq d^2$ non-trivial Pauli operators (which are two-outcome observables). Then the Aaronson-Rothblum bound seems to imply that only $O((\log d)^4/\epsilon^8)$ samples are needed to perform a full tomography, which seems to violate the sample complexity bound $O(d^2/\epsilon^2)$ for full tomography by O'Donnell & Wright (2016), Haah et. al. (2017).
How is this paradox resolved?
 A: To resolve the tension, note the results refer to different measures of accuracy.  For instance, in O'Donnell and Wright (2016), Corollary 1.4, the state $\rho \in \mathbb{C}^{d \times d}$ can be estimated with $O(d^2/\epsilon^2)$ copies, where $\epsilon$ refers to the error in trace-distance, i.e. $\frac{1}{2} ||\widetilde{\rho}-\rho|_1 < \epsilon$, using an estimate $\widetilde{\rho}$ for $\rho$.  (They remark this scaling is optimal when $\epsilon$ is fixed, referring to Flammia et al. (2012).)
Meanwhile, the Aaronson-Rothblum result uses $O((\log(M))^2 (\log d)^2 / \epsilon^8)$ copies to estimate the outcomes of $m$ different two-outcome observables to within additive error $\epsilon$ for each observable.
If we use the Aaronson-Rothblum procedure to estimate $\textrm{Tr}(\rho\vec{\sigma})$ for each product of Pauli operators (denoted $\vec{\sigma}$), then the reconstruction that I assume you imagine is
$$ \widetilde{\rho} = \frac{1}{d} \sum_{\vec{\sigma}} c_{\sigma} \vec{\sigma} $$
where $c_\sigma$ is the estimate of $\textrm{Tr}(\rho\vec{\sigma})$ provided by the Aaronsom-Rothblum procedure, with $|c_{\sigma} - \textrm{Tr}(\rho\vec{\sigma})| < \epsilon$. But then a naive bound on the trace distance only yields
$$ ||\widetilde{\rho}-\rho||_1 \leq\sqrt{d}||\widetilde{\rho}-\rho||_2 \leq \sqrt{d} \sqrt{\sum_\vec{\sigma} \frac{1}{d} |c_\sigma-\textrm{Tr}(\rho \vec{\sigma})|^2}  \leq d \epsilon $$
where I used $||X||_1 \leq \sqrt{d} ||X||_2$ and the fact that the sum has $d^2$ terms.  Note the upper bound is $d\epsilon$, not $\epsilon$.  So at least using the naive reconstruction of $\rho$, along with my naive bound on $||\widetilde{\rho}-\rho||_1$, we didn't manage to violate the optimality of the O'Donnell and Wright result.
