当结果已知是对称时加速矩阵乘法

Question

我知道矩阵乘法的结果是对称的。是否有一个R包或一些标准方法，我可以通过只计算下半部/上半部三角形然后将结果复制到另一半来加速我的计算。

我知道tcrossprod受益于这个事实，当只有一个参数提供，但我想提供两个矩阵。

这里就是结果是对称的一个例子：

n <- 100 
m <- 200 
s<-matrix(runif(n^2),n,n) 
s[lower.tri(s)] <- t(s)[lower.tri(s)] 
x <- matrix(runif(m*n), m, n) 
x %*% s %*% t(x) 

tcrossprod似乎并没有成为解决方案：

library(microbenchmark) 
microbenchmark(x %*% s %*% t(x), tcrossprod(x %*% s, x)) 

我试图使用RCPP，甚至没有复制一步，这是比R的乘法慢（虽然我坦率地承认我是一个初学者C++/Rcpp用户）：

w <- s %*% t(x) 
mm = Rcpp::cppFunction(
'NumericMatrix mmult(NumericMatrix m , NumericMatrix v) 
{ 
    NumericMatrix out(m.nrow(), v.ncol()); 

    for (int i = 0; i < m.nrow(); i++) 
    { 
    for (int j = 0; j < i + 1; j++) 
    { 
     for(int k = 0; k < m.ncol(); k++){ 
     out(i,j) += m(i,k) * v(k,j) ; 
     } 
    } 
    } 
    return out; 
}' 
) 

microbenchmark(mm(x, w), x %*% w) 

我认为如果.Internal functiondo_matprod中的sym变量被暴露并且可以被用户设置为真，这将被解决。不过，我真的不希望惹这样的事情......

Answer 1

它似乎并不像matrix包采取andvantage对称性：

> n <- 100 
> x <- s <- matrix(runif(n^2),n,n) 
> s[lower.tri(s)] <- t(s)[lower.tri(s)] 
> 
> library(Matrix) 
> s_sym <- Matrix(forceSymmetric(s)) 
> class(s_sym) # has the symmetric class 
[1] "dsyMatrix" 
attr(,"package") 
[1] "Matrix" 
> 
> library(microbenchmark) 
> microbenchmark(x %*% x, s %*% s, s_sym %*% s_sym) 
Unit: microseconds 
      expr min lq mean median uq max neval 
     x %*% x 461 496 571 528 625 1008 100 
     s %*% s 461 499 560 532 572 986 100 
s_sym %*% s_sym 553 568 667 624 701 1117 100 

没有任何迹象表明，它应在帮助文件：

基本矩阵产品，%*%实现我们所有的矩阵和 也为sparseVector类，完全类似的r基地矩阵 和矢量对象。功能crossprod和tcrossprod是矩阵 产品或“交叉产品”，理想情况下有效实施，不需要 计算t(.)。当易于检测时，例如，在crossprod(m), 一个参数情况下，它们也返回分类矩阵。 tcrossprod()取矩阵的转置矩阵的交叉乘积。 tcrossprod(x)正式相当于，但 快，呼吁x %*% t(x)，所以tcrossprod(x, y)而不是 x %*% t(y)。

用于您的解决方案是让使用Rcpp包装功能和R_ext/BLAS.h提供的BLAS功能。你可以做到这一点，如下所示：做一个func.cpp像这样的：

// added to get $(BLAS_LIBS) in compile flags 
//[[Rcpp::depends(RcppArmadillo)]] 
#include <Rcpp.h> 
#include <R_ext/BLAS.h> 

/* 
    Wrapper for BLAS dsymm. See dsymm http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_ga253c8edb8b21d1b5b1783725c2a6b692.html#ga253c8edb8b21d1b5b1783725c2a6b692 
    Only works with side = 'R' 
    Note intput is by refernce with & 
*/ 
// [[Rcpp::export]] 
Rcpp::NumericMatrix blas_dsymm(
    char uplo, int m, int n, double alpha, 
    const Rcpp::NumericMatrix &A, const Rcpp::NumericMatrix &B){ 
    // set lda, ldb and ldc 
    int lda = n, ldb = m, ldc = m; 

    // make new matrix with dim(m, n) 
    Rcpp::NumericMatrix C(m, n); // default values are zero 
    double beta = 0; 

    F77_NAME(dsymm)(
    "R" /* side */, &uplo, &m, &n, &alpha, 
    A.begin(), &lda, B.begin(), &ldb, &beta, C.begin(), &ldc); 

    return(C); 
} 

然后运行下列R-脚本：

> n <- 100 
> m <- 200 
> s<-matrix(runif(n^2),n,n) 
> s[lower.tri(s)] <- t(s)[lower.tri(s)] 
> x <- matrix(runif(m*n), m, n) 
> 
> library("Rcpp") 
> sourceCpp("func.cpp") 
> 
> out <- x %*% s 
> out_blas <- blas_dsymm(
+ uplo = "U", m = nrow(x), n = ncol(x), 
+ alpha = 1, A = s, B = x) 
> 
> all.equal(out, out_blas) 
[1] TRUE 
> 
> library(microbenchmark) 
> microbenchmark(
+ dense = x %*% s, 
+ BLAS = blas_dsymm(
+  uplo = "U", m = nrow(x), n = ncol(x), 
+  alpha = 1, A = s, B = x)) 
Unit: microseconds 
    expr  min  lq  mean median  uq  max neval 
dense 880.989 950.3225 1114.744 1066.866 1159.311 2783.213 100 
    BLAS 858.866 938.6680 1169.839 1016.495 1225.286 3261.633 100 

它似乎并没有在这里做任何区别。请注意，您需要安装RcppArmadillo和Rcpp软件包。

Answer 2

不要用for循环重新编码矩阵乘法。 线性代数库对此进行了高度优化，您可能会慢10倍（或更糟糕）。

对于矩阵计算，您不会通过使用RcppArmadillo或RcppEigen获得太多（或松散）。

如果你想获得，你可以改变你正在使用的数学库，例如使用带有Microsoft R Open的MKL。