3:43 p.m.
Hi,
I'm observing an unexpectedly drastic slowdown of mols() when switching
to SVD (with 'set svd on').
My toy benchmark is to regress a 200x10 LHS matrix on a 200x80 RHS
matrix 5000 times. The matrices are not redrawn in order to measure just
the mols speed, not the RNG. On an oldish machine I get:
~ 3 sec without svd
~ 60 sec with svd
For comparison, Python/Numpy's linalg.lstsq takes ~ 8 sec and is said to
use SVD always.
There are at least two possible differences between gretl/svd and numpy:
1) numpy might be using Intel's Math Kernel Library (MKL) in the
background instead of standard Lapack (coming from the conda/Anaconda
distro).
2) The web says that numpy uses Lapack's DGELSD
(http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_g...),
whereas judging from lib/src/gretl_matrix.c in the function
gretl_matrix_SVD_ols (and maybe gretl_matrix_multi_SVD_ols? - don't know
what that is for) it seems that gretl uses DGELSS.
On the Netlib Lapack page it says
(https://www.netlib.org/lapack/lug/node27.html) that "The subroutine
xGELSD is significantly faster than its older counterpart xGELSS", and
it labels DGELSD with "solve LLS using divide-and-conquer SVD".
So, any negative side effects blocking the switch to DGELSD, or is it
something like a free lunch?
thanks
sven
11:06 p.m.
On Fri, 27 Dec 2019, Sven Schreiber wrote:
Hi,
I'm observing an unexpectedly drastic slowdown of mols() when switching
to SVD (with 'set svd on').
My toy benchmark is to regress a 200x10 LHS matrix on a 200x80 RHS
matrix 5000 times. The matrices are not redrawn in order to measure just
the mols speed, not the RNG. On an oldish machine I get:
~ 3 sec without svd
~ 60 sec with svd
For comparison, Python/Numpy's linalg.lstsq takes ~ 8 sec and is said to
use SVD always.
There are at least two possible differences between gretl/svd and numpy:
1) numpy might be using Intel's Math Kernel Library (MKL) in the
background instead of standard Lapack (coming from the conda/Anaconda
distro).
2) The web says that numpy uses Lapack's DGELSD
(http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_g...),
whereas judging from lib/src/gretl_matrix.c in the function
gretl_matrix_SVD_ols (and maybe gretl_matrix_multi_SVD_ols? - don't know
what that is for) it seems that gretl uses DGELSS.
On the Netlib Lapack page it says
(https://www.netlib.org/lapack/lug/node27.html) that "The subroutine
xGELSD is significantly faster than its older counterpart xGELSS", and
it labels DGELSD with "solve LLS using divide-and-conquer SVD".
So, any negative side effects blocking the switch to DGELSD, or is it
something like a free lunch?
It appears to be a free lunch, pretty much. The speed-up is
significant but not huge, something like 30 percent. That's now in
git. But the primary source of difference between the gretl and
numpy times that you quote must be due to the respective Blas/Lapack
implementations.
Here's what I'm seeing on a fairly elderly PC running Fedora with
gretl linked against openblas. I reduced the number of replications
to 2000 (impatience) and created a baseline of accuracy by running
mpols with GRETL_MP_BITS=4096.
# gretl svd on, using DGELSS
gretl (mols): 8.43043 seconds
maxerr (gretl) = 0.000000000000007
python (linalg.lstsq): 5.63282 seconds
maxerr (python) = 0.000000000000007
# gretl svd on, using DGELSD
gretl (mols): 6.00157 seconds
maxerr (gretl) = 0.000000000000007
python (linalg.lstsq): 5.62002 seconds
maxerr (python) = 0.000000000000007
# gretl svd off (Cholesky)
gretl (mols): 0.396789 seconds
maxerr (gretl) = 0.000000000000005
python (linalg.lstsq): 5.62659 seconds
maxerr (python) = 0.000000000000007
From this, three points are apparent: (1) as stated above, DGELSD is
close to 30% faster than DGELSS; (2) numpy is a little faster than
us on SVD; and (3) you get just as accurate results an order of
magnitude faster via Cholesky (the mols default) provided the
regressors (as here) are not horribly collinear.
My test script:
<hansl>
set verbose off
set seed 12345999
n = 200 # observations
k = 80 # regressors
g = 10 # equations
matrix X = mnormal(n,k)
matrix U = mnormal(n,g)
matrix B = ones(k,g)
matrix Y = X*B + U
# get results in multiple precision as baseline;
# I used GRETL_MP_BITS=4096
Bmp = zeros(k,g)
loop j=1..g -q
Bmp[,j] = mpols(Y[,j],X)
endloop
# set svd on # or not
matrix Bgr
set stopwatch
loop 2000 -q
Bgr = mols(Y,X)
endloop
printf "gretl (mols): %#g seconds\n", $stopwatch
printf "maxerr (gretl) = %.15f\n", max(abs(Bgr - Bmp))
mwrite(Y, "Y.mat", 1)
mwrite(X, "X.mat", 1)
foreign language=python
import numpy as np
import time as tm
Y = gretl_loadmat('Y.mat', 1)
X = gretl_loadmat('X.mat', 1)
# numpy OLS
t0 = tm.time()
for i in range(1, 2000):
B = np.linalg.lstsq(X, Y)[0]
print("python (linalg.lstsq): %#g seconds" % (tm.time() - t0))
gretl_export(np.asmatrix(B), 'Bpy.mat', 1)
end foreign
matrix Bpy = mread("Bpy.mat", 1)
printf "maxerr (numpy) = %.15f\n", max(abs(Bpy - Bmp))
</hansl>
Allin
3:17 a.m.
Am 28.12.2019 um 05:06 schrieb Allin Cottrell:
On Fri, 27 Dec 2019, Sven Schreiber wrote:
> So, any negative side effects blocking the switch to DGELSD, or is it
> something like a free lunch?
It appears to be a free lunch, pretty much. The speed-up is
significant but not huge, something like 30 percent. That's now in
git. But the primary source of difference between the gretl and numpy
times that you quote must be due to the respective Blas/Lapack
implementations.
Good, thanks!
From this, three points are apparent: (1) as stated above, DGELSD is
close to 30% faster than DGELSS; (2) numpy is a little faster than us
on SVD; and (3) you get just as accurate results an order of magnitude
faster via Cholesky (the mols default) provided the regressors (as
here) are not horribly collinear.
BTW, I came across this whole question because in
the Python version of
the Tensorflow package (from Google) they have a different
implementation of linalg.lstsq: There you can explicitly switch between
fast (=Cholesky) and not so fast (=SVD I think). (There I could also
vectorize out the loop by using a third dimension, which makes it still
a little bit faster, but that's a different story. Plain numpy also
seems to be moving in this direction starting with 1.8, but it isn't in
numpy.linalg.lstsq yet.)
for i in range(1, 2000):
Strictly speaking this loop will run 1999 times, not 2000. You'd have to
do range(0,2000) or just range(2000). This isn't hansl, hehe!
cheers
sven
5:16 a.m.
Am 28.12.2019 um 05:06 schrieb Allin Cottrell:
>
> git. But the primary source of difference between the gretl and numpy
> times that you quote must be due to the respective Blas/Lapack
> implementations.
One more question about Blas: I have a "Kaby Lake R" mobile CPU here
according to Wikipedia, but gretl's $sysinfo tells me something about
blascore=haswell. If I understand https://github.com/xianyi/OpenBLAS
correctly, there are some more supported features like AVX512 starting
in Skylake (the generation before Kaby Lake). So would it make sense for
gretl to target skylake if available?
thanks
sven
6:46 a.m.
On 28.12.2019 11:16, Sven Schreiber wrote:
One more question about Blas: I have a "Kaby Lake R" mobile
CPU here
according to Wikipedia, but gretl's $sysinfo tells me something about
blascore=haswell. If I understand https://github.com/xianyi/OpenBLAS
correctly, there are some more supported features like AVX512 starting
in Skylake (the generation before Kaby Lake). So would it make sense for
gretl to target skylake if available?
In OpenBLAS there is no such mircokernel: you have HASWELL for haswell,
skylake and most of the kaby lake CPUs and you have SKYLAKEX (with
AVX512) for skylake-x (they have "X" at the and of the name).
Marcin
--
Marcin Błażejowski
9:56 a.m.
Am 28.12.2019 um 12:46 schrieb Marcin Błażejowski:
On 28.12.2019 11:16, Sven Schreiber wrote:
> One more question about Blas: I have a "Kaby Lake R" mobile CPU here
> according to Wikipedia, but gretl's $sysinfo tells me something about
> blascore=haswell. If I understand https://github.com/xianyi/OpenBLAS
> correctly, there are some more supported features like AVX512 starting
> in Skylake (the generation before Kaby Lake). So would it make sense for
> gretl to target skylake if available?
In OpenBLAS there is no such mircokernel: you have HASWELL for haswell,
skylake and most of the kaby lake CPUs and you have SKYLAKEX (with
AVX512) for skylake-x (they have "X" at the and of the name).
Ah OK, good to know. Unfortunate that on the page above they simply
write "Skylake".
thanks
sven
6:30 a.m.
On Fri, 27 Dec 2019, Allin Cottrell wrote:
It appears to be a free lunch, pretty much. The speed-up is
significant but
not huge, something like 30 percent. That's now in git. But the primary
source of difference between the gretl and numpy times that you quote must be
due to the respective Blas/Lapack implementations.
Here's what I'm seeing on a fairly elderly PC running Fedora with gretl
linked against openblas. I reduced the number of replications to 2000
(impatience) and created a baseline of accuracy by running mpols with
GRETL_MP_BITS=4096.
# gretl svd on, using DGELSS
gretl (mols): 8.43043 seconds
maxerr (gretl) = 0.000000000000007
python (linalg.lstsq): 5.63282 seconds
maxerr (python) = 0.000000000000007
# gretl svd on, using DGELSD
gretl (mols): 6.00157 seconds
maxerr (gretl) = 0.000000000000007
python (linalg.lstsq): 5.62002 seconds
maxerr (python) = 0.000000000000007
# gretl svd off (Cholesky)
gretl (mols): 0.396789 seconds
maxerr (gretl) = 0.000000000000005
python (linalg.lstsq): 5.62659 seconds
maxerr (python) = 0.000000000000007
> From this, three points are apparent: (1) as stated above, DGELSD is
close to 30% faster than DGELSS; (2) numpy is a little faster than us on SVD;
and (3) you get just as accurate results an order of magnitude faster via
Cholesky (the mols default) provided the regressors (as here) are not
horribly collinear.
Hm, I'm seeing something weird here: on my home laptop (an 8-core, 4 real
cpus machine) I'm getting results that are fairly consistent with yours.
On the other hand, trying the same on my work PC (a 12-core box) I'm
seeing this:
----------------------------------------------
Old code (with set svd on):
----------------------------------------------
gretl (mols): 7.09208 seconds
maxerr (gretl) = 0.000000000000010
python (linalg.lstsq): 3.92083 seconds
maxerr (numpy) = 0.000000000000008
----------------------------------------------
New code (with set svd on):
----------------------------------------------
gretl (mols): 3.79993 seconds
maxerr (gretl) = 0.000000000000008
python (linalg.lstsq): 3.89591 seconds
----------------------------------------------
New code (with set svd on):
----------------------------------------------
gretl (mols): 120.903 seconds
maxerr (gretl) = 0.000000000000005
python (linalg.lstsq): 5.63783 seconds
It looks as if turning svd off makes mols a good deal slower (by the way:
all 12 cpus go at 100% for the whole time). I have no time now to check
why this happens, but I'll throw a idea in: perhaps, for some reason, I'm
getting the QR decomposition instead of Cholesky?
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:14 a.m.
On Sat, 28 Dec 2019, Riccardo (Jack) Lucchetti wrote:
On Fri, 27 Dec 2019, Allin Cottrell wrote:
> It appears to be a free lunch, pretty much. The speed-up is significant but
> not huge, something like 30 percent. That's now in git. But the primary
> source of difference between the gretl and numpy times that you quote must
> be due to the respective Blas/Lapack implementations.
>
> Here's what I'm seeing on a fairly elderly PC running Fedora with gretl
> linked against openblas. I reduced the number of replications to 2000
> (impatience) and created a baseline of accuracy by running mpols with
> GRETL_MP_BITS=4096.
>
> # gretl svd on, using DGELSS
> gretl (mols): 8.43043 seconds
> maxerr (gretl) = 0.000000000000007
> python (linalg.lstsq): 5.63282 seconds
> maxerr (python) = 0.000000000000007
>
> # gretl svd on, using DGELSD
> gretl (mols): 6.00157 seconds
> maxerr (gretl) = 0.000000000000007
> python (linalg.lstsq): 5.62002 seconds
> maxerr (python) = 0.000000000000007
>
> # gretl svd off (Cholesky)
> gretl (mols): 0.396789 seconds
> maxerr (gretl) = 0.000000000000005
> python (linalg.lstsq): 5.62659 seconds
> maxerr (python) = 0.000000000000007
>
>> From this, three points are apparent: (1) as stated above, DGELSD is
> close to 30% faster than DGELSS; (2) numpy is a little faster than us on
> SVD; and (3) you get just as accurate results an order of magnitude faster
> via Cholesky (the mols default) provided the regressors (as here) are not
> horribly collinear.
Hm, I'm seeing something weird here: on my home laptop (an 8-core, 4 real
cpus machine) I'm getting results that are fairly consistent with yours. On
the other hand, trying the same on my work PC (a 12-core box) I'm seeing
this:
----------------------------------------------
Old code (with set svd on):
----------------------------------------------
gretl (mols): 7.09208 seconds
maxerr (gretl) = 0.000000000000010
python (linalg.lstsq): 3.92083 seconds
maxerr (numpy) = 0.000000000000008
----------------------------------------------
New code (with set svd on):
----------------------------------------------
gretl (mols): 3.79993 seconds
maxerr (gretl) = 0.000000000000008
python (linalg.lstsq): 3.89591 seconds
----------------------------------------------
New code (with set svd on):
Should that be, with set svd off?
----------------------------------------------
gretl (mols): 120.903 seconds
maxerr (gretl) = 0.000000000000005
python (linalg.lstsq): 5.63783 seconds
It looks as if turning svd off makes mols a good deal slower (by the way: all
12 cpus go at 100% for the whole time). I have no time now to check why this
happens, but I'll throw a idea in: perhaps, for some reason, I'm getting the
QR decomposition instead of Cholesky?
If that were the case you should see on stderr,
"gretl_matrix_multi_ols: switching to QR decomp"
Are these 12 "real" cores? The timings I showed were from a quad
core box with hyperthreads and I set OMP_NUM_THREADS=4 to prevent
hyperthreading, which slows things down. (If OpenBLAS doesn't use
OMP, you'd set OPENBLAS_NUM_THREADS if wanted.)
Allin
10:42 a.m.
Am 28.12.2019 um 05:06 schrieb Allin Cottrell:
On Fri, 27 Dec 2019, Sven Schreiber wrote:
> On the Netlib Lapack page it says
> (https://www.netlib.org/lapack/lug/node27.html) that "The subroutine
> xGELSD is significantly faster than its older counterpart xGELSS", and
> it labels DGELSD with "solve LLS using divide-and-conquer SVD".
>
> So, any negative side effects blocking the switch to DGELSD, or is it
> something like a free lunch?
It appears to be a free lunch, pretty much. The speed-up is significant
but not huge, something like 30 percent. That's now in git.
A follow-up about the standalone SVD without a LS problem:
Gretl uses DGESVD, while on
https://www.netlib.org/lapack/lug/node32.html it says: "[DGESDD] is much
faster than [DGESVD] for large matrices, but uses more workspace."
As I guess that nowadays workspace is not an issue anymore except in
extreme cases, this sounds as if this switch should be done as well, no?
thanks
sven
11:40 a.m.
New subject: eigenvalue speed and more (was Re: Re: speed of (m)ols with svd)
Am 28.12.2019 um 16:42 schrieb Sven Schreiber:
Am 28.12.2019 um 05:06 schrieb Allin Cottrell:
> It appears to be a free lunch, pretty much. The speed-up is
significant
> but not huge, something like 30 percent. That's now in git.
A follow-up about the standalone SVD without a LS problem:
OK sorry about this flood of postings, but there might be more.
1) Eigenvalues: An analogous situation to the SVD apparently exists for
eigenvalue problems: Gretl uses DSYEV, but there are also DSYEVD and
DSYEVR. On https://www.netlib.org/lapack/lug/node30.html it says again:
"[DSYEVD] is much faster than [DSYEV] for large matrices, but uses more
workspace."
But then there's also a third competitor: "[RRR / DSYEVR] is the fastest
algorithm of all (except for a few cases), and uses the least
workspace." Sounds almost like a no-brainer.
2) Complex stuff; I looked in lib/src/gretl_cmatrix.c for Lapack
routines used. Short summary:
-) ZHEEV: as with DSYEV; use ZHEEVD or ZHEEVR instead
-) ZGESVD: as with DGESVD; use ZGESDD instead
cheers
sven
4:36 p.m.
New subject: eigenvalue speed and more (was Re: Re: speed of (m)ols with svd)
On Sat, 28 Dec 2019, Sven Schreiber wrote:
Am 28.12.2019 um 16:42 schrieb Sven Schreiber:
> Am 28.12.2019 um 05:06 schrieb Allin Cottrell:
>> It appears to be a free lunch, pretty much. The speed-up is significant
>> but not huge, something like 30 percent. That's now in git.
>
> A follow-up about the standalone SVD without a LS problem:
OK sorry about this flood of postings, but there might be more.
1) Eigenvalues: An analogous situation to the SVD apparently exists for
eigenvalue problems: Gretl uses DSYEV, but there are also DSYEVD and
DSYEVR. On https://www.netlib.org/lapack/lug/node30.html it says again:
"[DSYEVD] is much faster than [DSYEV] for large matrices, but uses more
workspace."
Trying DSYEVD would be fairly straightforward. (Though none of these
alternative functions are simply drop-in replacements for the more
basic versions.)
But then there's also a third competitor: "[RRR / DSYEVR] is
the fastest
algorithm of all (except for a few cases), and uses the least
workspace." Sounds almost like a no-brainer.
Almost... except that its signature is pretty different from the
others, and kinda awkward. It would be a bigger job to make use of
it for our eigensym. This is actually a gripe of mine about the
fancier "alternative" lapack functions. Fair enough to add an
argument or two as needed (as in DGESDD vs DGESVD) but not so nice
to change the semantics of existing arguments, which I think could
probably be avoided.
2) Complex stuff; I looked in lib/src/gretl_cmatrix.c for Lapack
routines used. Short summary:
-) ZHEEV: as with DSYEV; use ZHEEVD or ZHEEVR instead
-) ZGESVD: as with DGESVD; use ZGESDD instead
OK, worth looking at those too.
Allin
12:24 p.m.
On Sat, 28 Dec 2019, Sven Schreiber wrote:
A follow-up about the standalone SVD without a LS problem:
Gretl uses DGESVD, while on
https://www.netlib.org/lapack/lug/node32.html it says: "[DGESDD] is much
faster than [DGESVD] for large matrices, but uses more workspace."
As I guess that nowadays workspace is not an issue anymore except in
extreme cases, this sounds as if this switch should be done as well, no?
OK, that's in git. For now I've kept the old code alongside the new
in case of any problems.
Allin
8:01 a.m.
Am 28.12.19 um 18:24 schrieb Allin Cottrell:
On Sat, 28 Dec 2019, Sven Schreiber wrote:
> A follow-up about the standalone SVD without a LS problem:
>
> Gretl uses DGESVD, while on
> https://www.netlib.org/lapack/lug/node32.html it says: "[DGESDD] is much
> faster than [DGESVD] for large matrices, but uses more workspace."
>
> As I guess that nowadays workspace is not an issue anymore except in
> extreme cases, this sounds as if this switch should be done as well, no?
OK, that's in git. For now I've kept the old code alongside the new in
case of any problems.
Hm, I'm not seeing any speedups (with the ridge package which uses the
svd() function), but admittedly the matrices there are probably small in
the sense above. (Also here I'm comparing to the state of the source
from September or so, I think.) I've also noticed that on stock Linux
Mint the openblas offering seems to have this pthreads thingy that has
caused some confusion in the past. But my impression from the discussion
was that it won't have any real-world impact (?).
thanks
sven
11:02 a.m.
On Sun, 29 Dec 2019, Sven Schreiber wrote:
Am 28.12.19 um 18:24 schrieb Allin Cottrell:
> On Sat, 28 Dec 2019, Sven Schreiber wrote:
>
>> A follow-up about the standalone SVD without a LS problem:
>>
>> Gretl uses DGESVD, while on
>> https://www.netlib.org/lapack/lug/node32.html it says: "[DGESDD] is much
>> faster than [DGESVD] for large matrices, but uses more workspace."
>>
>> As I guess that nowadays workspace is not an issue anymore except in
>> extreme cases, this sounds as if this switch should be done as well, no?
>
> OK, that's in git. For now I've kept the old code alongside the new in
> case of any problems.
>
Hm, I'm not seeing any speedups (with the ridge package which uses the
svd() function) [...]
I'm seeing about 40 percent gain on this:
<hansl>
set verbose off
matrix X = mnormal(200,100)
X = X'X
matrix U
matrix V
set stopwatch
loop i=1..1000 -q
matrix s = svd(X, &U, &V)
endloop
printf "elapsed: %#gs\n", $stopwatch
matrix USV = (U .* s) * V
eval max(abs(X - USV))
</hansl>
For ease of comparison you can force use of the old code by putting
GRETL_OLD_SVD=1 into the environment.
Allin
11:08 a.m.
Am 29.12.19 um 17:02 schrieb Allin Cottrell:
On Sun, 29 Dec 2019, Sven Schreiber wrote:
> Hm, I'm not seeing any speedups (with the ridge package which uses the
> svd() function) [...]
I'm seeing about 40 percent gain on this:
Ah, excellent!
Is this actually expected to have an impact also on the Johansen
procedure, or is some other background code used there?
(I'm not talking about the eigenvalue routines, since these haven't been
changed yet AFAIK. But some SVD appeared there as well I think.)
thanks
sven
1602
days inactive
1604
days old
14 comments
4 participants
participants (4)
-
Allin Cottrell -
Marcin Błażejowski -
Riccardo (Jack) Lucchetti -
Sven Schreiber