<div dir="ltr"><div><div><div><div>There is actually one other command to test, `crossprod`, or in this case `tcrossprod(x, A)`.<br><br></div>Following your code, first run and ensure TRUE is returned.<br><br><div style="margin-left:40px">
<span style="font-family:courier new,monospace">d <- 2000</span><br><span style="font-family:courier new,monospace">A <- matrix(1.0*(1:d^2),nrow=d)</span><br><span style="font-family:courier new,monospace">x <- 1.0*(1:d)</span><br>
<span style="font-family:courier new,monospace">all.equal(A%*%x, tcrossprod(x, A), mvprod1(A,x), mvprod2(A,x))</span><br></div><br></div>Now, using an i7-3740QM 2.7Ghz CPU with 8GB of RAM as a testbed on Windows 7 (with R & Rcpp compiled from source throwing "<span style="color:rgb(0,0,255)"><span style="font-family:arial,helvetica,sans-serif">-march=core-avx-i -O3 -std=gnu++0x --param l1-cache-line-size=64 --param l1-cache-size=32 --param l2-cache-size=256</span></span>" FWIW) I get the following timings:<br>
<br><span style="font-family:courier new,monospace"><font>microbenchmark(A%*%x, tcrossprod(x, A), mvprod1(A, x), mvprod2(A, x), times = 1000L, control = list(order = 'block'))</font><br><br><span class="" style="border-collapse:separate;color:rgb(0,0,0);font-size:13px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:15px;text-indent:0px;text-transform:none;white-space:pre-wrap;word-spacing:0px;background-color:rgb(225,226,229)"><pre tabindex="0" class="" style="font-size:10pt!important;outline-style:none;outline-width:initial;outline-color:initial;border-style:none;border-width:initial;border-color:initial;white-space:pre-wrap!important;word-break:break-all;margin:0px;line-height:1.2">
Unit: milliseconds
expr min lq median uq max neval
A %*% x 13.183992 13.242965 13.269027 13.317917 16.680117 1000
tcrossprod(x, A) 6.967139 6.991869 7.004424 7.020214 9.795921 1000
mvprod1(A, x) 43.089952 43.134848 43.163953 43.489824 45.129644 1000
mvprod2(A, x) 4.564480 4.883693 4.911848 4.928208 7.435495 1000</pre></span></span><br></div>So the crossprod function is decently optimized as it is, but the prod2 is still faster.<br><br></div>Thanks for idea and post!<br>
<br>Avi<br><div><div><br></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Apr 25, 2014 at 12:47 PM, Søren Højsgaard <span dir="ltr"><<a href="mailto:sorenh@math.aau.dk" target="_blank">sorenh@math.aau.dk</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div link="blue" vlink="purple" lang="DA">
<div>
<p class="MsoNormal"><span lang="EN-US">Dear all,<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">When forming the matrix-vector-product one can form the inner products between As rows and x (mvprod1 below) or form a linear combination of As columns (mvprod2 below). The difference in computing time is striking (at
least for me) for large matrices and it really illustrates the “caching issue” (I assume):<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">> microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">Unit: milliseconds<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> expr min lq median uq max neval<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">mvprod1(A, x) 35.783812 36.310956 36.600186 36.991113 38.922889 100<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">mvprod2(A, x) 4.408892 4.436882 4.475368 4.568668 5.477176 100<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> A %*% x 8.514091 8.592230 8.734978 8.811951 9.738653 100<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">Just wanted to share this!<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">Cheers<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">Søren <u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">#include <Rcpp.h><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">using namespace Rcpp;<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">//[[Rcpp::export]]<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">NumericVector mvprod1 (NumericMatrix A, NumericVector x){<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> int nrA=A.nrow(), ncA=A.ncol(), i, j;<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> NumericVector y(nrA);<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> for (i=0; i<nrA; ++i){<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> for (j=0; j<ncA; ++j){<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> y[i] += A(i,j)*x[j];<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> }<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> }<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> return y;<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">}<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">//[[Rcpp::export]]<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">NumericVector mvprod2 (NumericMatrix A, NumericVector x){<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> int nrA=A.nrow(), ncA=A.ncol(), i, j;<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> NumericVector y(nrA);<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> for (j=0; j<ncA; ++j){<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> </span>for (i=0; i<nrA; ++i){<u></u><u></u></p>
<p class="MsoNormal"> y[i] += A(i,j)*x[j];<u></u><u></u></p>
<p class="MsoNormal"> <span lang="EN-US">}<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> }<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"> return y;<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">}<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">/*** R<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">d <- 5<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">A <- matrix(1.0*(1:d^2),nrow=d)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">x <- 1.0*(1:d)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">A%*%x<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">mvprod1(A,x)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">mvprod2(A,x)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">d <- 2000<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">A <- matrix(1.0*(1:d^2),nrow=d)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">x <- 1.0*(1:d)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">library(microbenchmark)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)<u></u><u></u></span></p>
<p class="MsoNormal"><span lang="EN-US">*/<u></u><u></u></span></p>
</div>
</div>
<br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br></div>