作者celestialgod (攸藍)
看板R_Language
標題[分享] 矩陣對向量運算的函數(像bsxfun(MATLAB))
時間Thu Apr 16 21:37:53 2015
[關鍵字]: bsxfun in R
[出處]:
http://stackoverflow.com/questions/3643555/multiply-rows-of-matrix-by-vector
[重點摘要]:
pastebin link:
http://pastebin.com/jGtHeYuj
M = 3
N = 5
(mat = replicate(N, 1:M))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 1 1 1 1
# [2,] 2 2 2 2 2
# [3,] 3 3 3 3 3
(vec = 1:N)
# [1] 1 2 3 4 5
(sweep(mat, MARGIN=2, vec,`+`))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 2 3 4 5 6
# [2,] 3 4 5 6 7
# [3,] 4 5 6 7 8
(sweep(mat, MARGIN=2, vec,`*`))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 2 3 4 5
# [2,] 2 4 6 8 10
# [3,] 3 6 9 12 15
(sweep(mat, MARGIN=2, vec,`/`))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 0.5 0.3333333 0.25 0.2
# [2,] 2 1.0 0.6666667 0.50 0.4
# [3,] 3 1.5 1.0000000 0.75 0.6
(sweep(mat, MARGIN=2, vec,`-`))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 0 -1 -2 -3 -4
# [2,] 1 0 -1 -2 -3
# [3,] 2 1 0 -1 -2
(sweep(mat, MARGIN=2, vec, pmax))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 2 3 4 5
# [2,] 2 2 3 4 5
# [3,] 3 3 3 4 5
(sweep(mat, MARGIN=2, vec, pmin))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 1 1 1 1
# [2,] 1 2 2 2 2
# [3,] 1 2 3 3 3
## usage in scaling
N = 10; p = 3
dat = matrix(rnorm(N*p), N)
dat_standardized = sweep(sweep(dat, 2, colMeans(dat)), 2,
apply(dat, 2, sd), '/')
dat_normalized = sweep(sweep(dat, 2, apply(dat, 2, min)), 2,
apply(dat, 2, function(x) diff(range(x))), '/')
## usage in computing kernel
kernelMatrix_f = function(X, center, sigma){
exp(sweep(sweep(X %*% t(center), 1,
rowSums(X**2)/2), 2, rowSums(center**2)/2) / (sigma**2))
}
## compare with Rcpp and kernlab
library(kernlab)
library(Rcpp)
library(RcppArmadillo)
sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
NumericMatrix kernelMatrix_cpp(NumericMatrix Xr,
NumericMatrix Centerr, double sigma) {
uword n = Xr.nrow(), b = Centerr.nrow(), row_index, col_index;
mat X(Xr.begin(), n, Xr.ncol(), false),
Center(Centerr.begin(), b, Centerr.ncol(), false),
KerX(X*Center.t());
colvec X_sq = sum(square(X), 1) / 2;
rowvec Center_sq = (sum(square(Center), 1)).t() / 2;
KerX.each_row() -= Center_sq;
KerX.each_col() -= X_sq;
KerX *= 1 / (sigma * sigma);
KerX = exp(KerX);
return wrap(KerX);
}')
N = 8000
p = 500
b = 1000
X = matrix(rnorm(N*p), ncol = p)
center = X[sample(1:N, b),]
sigma = 3
all.equal(kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center)@.Data,
kernelMatrix_cpp(X, center, sigma))
# TRUE
all.equal(kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center)@.Data,
kernelMatrix_f(X, center, sigma))
# TRUE
library(microbenchmark)
microbenchmark(Rfun = kernelMatrix_f(X, center, sigma),
kernlab = kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center),
Rcpp = kernelMatrix_cpp(X, center, sigma),
times = 20L
)
# Unit: milliseconds
# expr min lq mean median uq max neval
# Rfun 872.8854 981.1995 1036.2209 1074.0526 1098.9185 1132.9356 20
# kernlab 773.7489 841.9028 1059.3098 862.7476 883.3874 2979.2541 20
# Rcpp 490.2462 501.2993 520.1283 522.5060 532.5426 571.0949 20
速度不輸package: kernlab。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 36.225.212.167
※ 文章網址: https://webptt.com/m.aspx?n=bbs/R_Language/M.1429191476.A.F82.html
1F：→ Wush978: 這不知道能不能用BLAS算 04/16 22:51
DGER可以做 A = alpha*x*y' + A (+, -的case)
alpha: constant, x: m by 1, y: 1 by n, A: m by n
x 對應到上面的vec
y 是全部都1的向量
A 是mat
alpha = 1
http://tinyurl.com/mqwv2hc
2F：→ Wush978: 就我的理解，你應該直接#include "cblas.h"然後和R 04/17 01:50
3F：→ Wush978: linking就可以了。因為R應該內帶blas 04/17 01:50
4F：→ Wush978: 另外沒意外的話，用blas可能會是最快的算法 04/17 01:51
不過現在看起來BLAS能做的case只有+/-，*或是/或是min, max都無法
5F：推 penolove5566: 謝分享 04/17 02:20
6F：→ clickhere: 那是fortran函式,用 include "R_ext/BLAS.h" 04/17 02:53
也試過了
※ 編輯: celestialgod (36.225.212.167), 04/17/2015 14:43:28