[分享] 矩陣對向量運算的函數(像bsxfun(MATLAB))

看板R_Language作者celestialgod (攸藍)時間9年前 (2015/04/16 21:37)推噓1(1推 0噓 5→)

留言6則, 3人參與討論串1/1

[關鍵字]: 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://www.ptt.cc/bbs/R_Language/M.1429191476.A.F82.html

→

Wush978

04/16 22:51, , 1^F

04/16 22:51, 1^F

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