Result changes depending on whether matrix is sparse or not

Why do you think it is not a good idea?
Computationally wise it might be better use direct methods to find the solution.
But accuracy wise and certainly if one needs the Least Norm solution of the Least Squares there is nothing wrong about it.

Could it be that the Myth about “Don’t compute the inverse” is overblown?
See How Accurate is inv(A)*b?

By the way, MATLAB has the function lsqminnorm() which uses the complete orthogonal decomposition (COD) to compute the minimum norm least squares solution.
If it is available in Julia one could use it as well as it should be faster than pinv() which uses the SVD.