Early comments regarding this package mentioned the preference for integration in some manner with existing scipy solvers. It took me a while to see a way to do that. The problem was the lack of a quick method for detecting the likelihood of ill-conditioning (including singularity) which would make a regularized solver preferable to a standard solver like lstsq(). But I was able to develop such a quick method, called "strange()" -- because it detects strange behavior that is usually indicative of an ill-conditioned system. This method takes only about 4n^2 flops, rather than O(n^3) which the regularized solvers seem to take. So: lsmr and arls take a few times as long as lstsq. And lstsq takes a few times as long as "strange()". So inside arls() it is low-cost to call strange() first to see if it is safe to call the faster lstsq(). The special solvers -- arlseq, arlsgt, and arlsnn -- all get the benefit of this change as they depend on arls(). But that this capability applies only to the case of single right-hand-sides, as the benefit of this approach dissolves rapidly when multiple right-hand-sides are used, as all their solutions in arls() share a single call for the SVD. Also updating the in-code user documentation. Ron Jones