This was too long for a comment.
There are two papers that have a nice characterization of positivity in terms of the Bloch vector: quant-ph/0301152 and quant-ph/0302024. You might be able to use dynamic programming to speed up a check of positivity by using the conditions in these papers. Basically, they reduce the positivity constraint to a constraint on the characteristic polynomial of $\rho$. Using the well-known connection between coefficients of characteristic polynomials and elementary symmetric polynomials, they can derive necessary and sufficient conditions for positivity of the eigenvalues.
Then the idea is the following. Store the iterated "star" products of the Bloch vectors (see the paper for the definition... it's sort of like a cross product). There are only $d$ of them. Then use Newton's identity to compute the constraints.
I can't immediately see how to compute the star product in time $o(d^3)$ (little-oh), though, so this probably doesn't work. But I haven't thought about it, either.
Another possible issue is that this might not be numerically stable. Again, I haven't thought about it, that's why this should have been a comment.
UPDATE: Maris has privately told me that his comment below is only meant to imply that a single star product can be computed in time $O(d^2)$. That is insufficient for a fast algorithm. There are $d$ conditions that have to be checked, so the total runtime is still $O(d^3)$ for this approach.
In fact, I highly doubt that there is an algorithm that runs faster than $O(d^\omega)$ (matrix multiplication time), since all you've done is rewrite the matrix in a new matrix basis, and even then you would probably have to exploit structure. I think the best you can hope for is $O(d^3)$, which is an equivalent runtime to just a naive approach.