r-当矩阵的值很小时，为什么矩阵乘积较慢

Question

我创建了两个相同维度的矩阵A和B。CCD_ 3包含比CCD_ 4更大的值。矩阵乘法CCD_ 5比CCD_ 6快大约10倍。

为什么

## disable openMP
library(RhpcBLASctl); blas_set_num_threads(1); omp_set_num_threads(1)
A <- exp(-as.matrix(dist(expand.grid(1:60, 1:60))))
summary(c(A))
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.000000 0.000000 0.000000 0.001738 0.000000 1.000000 
B <- exp(-as.matrix(dist(expand.grid(1:60, 1:60)))*10)
summary(c(B))
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.0000000 0.0000000 0.0000000 0.0002778 0.0000000 1.0000000 
identical(dim(A), dim(B))
## [1] TRUE
system.time(A %*% A)
#    user  system elapsed 
#   2.387   0.001   2.389 
system.time(B %*% B)
#    user  system elapsed 
#  21.285   0.020  21.310
sessionInfo()
# R version 3.6.1 (2019-07-05)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Linux Mint 19.2
# Matrix products: default
# BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
# LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

当矩阵包含许多小条目时，这个问题可能与base:：chol((的速度减慢有关。

编辑：有些小数字似乎会减慢计算速度。其他人则不然。

slow <-  6.41135533887904e-164
fast1 <- 6.41135533887904e-150
fast2 <- 6.41135533887904e-170
Mslow <- array(slow, c(1000, 1000)); system.time(Mslow %*% Mslow)
#   user  system elapsed 
# 10.165   0.000  10.168 
Mfast1 <- array(fast1, c(1000, 1000)); system.time(Mfast1 %*% Mfast1)
#   user  system elapsed 
#  0.058   0.000   0.057 
Mfast2 <- array(fast2, c(1000, 1000)); system.time(Mfast2 %*% Mfast2)
#   user  system elapsed 
#  0.056   0.000   0.055 

Answer 1

您很可能希望使用.Machine$double.xmin而不是double.eps。这会将更少的数字设置为零，并具有相同的效果。为了避免低于标准的数字，您可能必须使用将这些数字设置为零的编译器标志重新编译BLAS，而不是引发FP陷阱。

Answer 2

R-devel邮件列表的回复表明，这可能是非正规数字的问题，或者openBLAS处理小数字的速度可能较慢。

发件人https://en.wikipedia.org/wiki/Denormal_number:

在计算机科学中，非正规化数或非正规化数字(现在通常被称为次正规数(填补了年零左右的下溢缺口浮点运算。任何幅值较小的非零数字最小的正规数是"次正规"。[…]极端在这种情况下，涉及非正规操作数的指令可能运行多达100速度慢了几倍。

事实上，B包含非常小的数字：

sum(B<.Machine$double.eps)
[1] 12832980
sort(unique(B[B>0]))[10^(0:3)]
[1] 4.940656e-324 2.280607e-320 6.302966e-295 2.185410e-141

如果将小数字设置为零，则计算具有预期的计算时间：

C <- B; C[abs(C)<.Machine$double.eps] <- 0
system.time(C %*% C)
user  system elapsed 
2.266   0.032   2.298 

是否有办法将.Machine$double.eps以下的值自动设置为零手工检查每个矩阵中的小数字似乎不方便。