That may mean that your problem is badly conditioned — have you considered using a Tikhonov (ridge regression) regularization? That corresponds to setting a damping parameter \lambda > 0 in IterativeSolvers.lsmr. This should improve the conditioning and convergence rate, and robustify the solutions to noise, at the expense of residual — you’ll have to experiment to find an appropriate regularization strength for your application.
The corresponding “direct”/non-iterative algorithm (ala the code above) is a variant of the normal equations:
B = reshape(cholesky(Hermitian(Wᵀ * Wᵀ' + λ * I)) \ (Wᵀ * Y), m, n)
(which is dominated by the O(N(mn)^2) cost of computing W^T W).