[SciPy-User] Least squares speed

[josef.pktd at gmail.com]

On Tue, Oct 1, 2013 at 5:24 PM, Nathaniel Smith <njs at pobox.com> wrote:
> Here are several ways to solve the least squares problem XB = Y:
>
> scipy.linalg.lstsq(x, y)
> np.linalg.lstsq(x, y)
> np.dot(scipy.linalg.pinv(x), y)
> np.dot(np.linalg.pinv(x), y)
>
> >From the documentation, I would expect these to be ordered by speed,
> fastest up top and slowest at the bottom. It turns out they *are*
> ordered by speed, except, on my system (linux x86-64, numpy 1.7.1,
> scipy 0.12.0, ubuntu-installed atlas), it's exactly backwards, with
> very substantial differences (more than a factor of 20!):
>
> # Typical smallish problem for me:
> In [40]: x = np.random.randn(780, 2)
>
> In [41]: y = np.random.randn(780, 32 * 275)
>
> In [42]: %timeit scipy.linalg.lstsq(x, y)
> 1 loops, best of 3: 494 ms per loop
>
> In [43]: %timeit np.linalg.lstsq(x, y)
> 1 loops, best of 3: 356 ms per loop
>
> In [44]: %timeit np.dot(scipy.linalg.pinv(x), y)
> 10 loops, best of 3: 62.2 ms per loop
>
> In [45]: %timeit np.dot(np.linalg.pinv(x), y)
> 10 loops, best of 3: 23.2 ms per loop
>
> Is this expected? I'm particularly confused at why scipy's lstsq
> should be almost 40% slower than numpy's. (And shouldn't we have a
> fast one-function way to solve least squares problems?)
>
> Any reason not to use the last option? Is it as numerically stable as lstsq?
>
> Is there any other version that's even faster or even more numerically stable?

If you have very long x, then using normal equation is faster for univariate y.

There are many different ways to calculate pinv
https://github.com/scipy/scipy/pull/289

np.lstsq breaks on rank deficient x IIRC, uses different Lapack
functions than scipy's lstsq

In very badly scaled cases (worst NIST case), scipy's pinv was a bit
more accurate than numpy's, but maybe just different defaults.
numpy's pinv was also faster than scipy's in the cases that I tried.

There is only a single NIST case that can fail using the defaults with
numpy pinv.

(what about using qr or chol_solve ?)

Lots of different ways to solve this and I never figured out a ranking
across different cases, speed, precision, robustness to
near-singularity.

Josef


>
> -n
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user