<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.EmailStyle17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:1371800498;
mso-list-type:hybrid;
mso-list-template-ids:534544282 1696358412 67698691 67698693 67698689 67698691 67698693 67698689 67698691 67698693;}
@list l0:level1
{mso-level-number-format:bullet;
mso-level-text:\F0D8;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Wingdings;
mso-fareast-font-family:Calibri;
mso-bidi-font-family:"Times New Roman";}
@list l0:level2
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:"Courier New";}
@list l0:level3
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Wingdings;}
@list l0:level4
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Symbol;}
@list l0:level5
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:"Courier New";}
@list l0:level6
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Wingdings;}
@list l0:level7
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Symbol;}
@list l0:level8
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:"Courier New";}
@list l0:level9
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;
font-family:Wingdings;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"> E <- (exp(A) * (1 - 1 / A) + 1 / A) / (exp(A) - 1)<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"> If A is a matrix by then, isn't exp a very slow (and imprecise:<o:p></o:p></p>
<p class="MsoNormal"> <a href="http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/">
http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/</a>)<o:p></o:p></p>
<p class="MsoNormal"> operation, isn't it? You do it twice on the same matrix.<span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">exp(A) is the element-by-element exponential, so it is not really slow.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">However, when elements of A are small, expm1(A) will be more accurate<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">than exp(A)-1 and you might do better by replacing both exp(A)'s by<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">functions of expm1(A). (Since the limit of E as A->0 is 1/2 you can do<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">better by computing E-1/2.) There is a chance that the increased accuracy<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">means you don't have to do as many iterations.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">Bill Dunlap<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">Spotfire, TIBCO Software<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D">wdunlap tibco.com<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<div style="border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt">
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">From:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> rcpp-devel-bounces@lists.r-forge.r-project.org [mailto:rcpp-devel-bounces@lists.r-forge.r-project.org]
<b>On Behalf Of </b>Paul Johnson<br>
<b>Sent:</b> Saturday, November 09, 2013 5:55 PM<br>
<b>To:</b> Xavier Robin<br>
<b>Cc:</b> Rcpp-devel@lists.R-forge.R-project.org<br>
<b>Subject:</b> Re: [Rcpp-devel] How much speedup for matrix operations?<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Hi<o:p></o:p></p>
</div>
<p class="MsoNormal" style="margin-bottom:12.0pt">If you have the time to test, you will find that some things that you think will be faster are actually slower, while some seemingly unimportant things will give huge acceleration.
<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">I predict that you can learn more about R efficiencies, esp. crossprod & tcrossprod, and after you do that, my guess it can go faster using Rcpp & a BLAS like openBLAS.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I wish you'd try these and see how they affect performance. Then let us know what you find out. If you find anything clear, I'll add to my collection of speedup advices.
<a href="http://pj.freefaculty.org/blog/?p=122">http://pj.freefaculty.org/blog/?p=122</a>.<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<div>
<p class="MsoNormal">1. Avoid repeated costly calculations. Please examine the use of exp(A). and 1/A in this part<br>
<br>
E <- (exp(A) * (1 - 1 / A) + 1 / A) / (exp(A) - 1)<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal" style="margin-bottom:12.0pt">If A is a matrix by then, isn't exp a very slow (and imprecise:
<a href="http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/">
http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/</a>) operation, isn't it?<br>
You do it twice on the same matrix. <o:p></o:p></p>
</div>
<p class="MsoNormal">Did you know that DIVISION is much slower than multiplication in a modern CPU? Surprised me to learn that last year. Can't you re-arrange this so that 1/A is not calculated repeatedly?<o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">2. please examine this usage:<br>
<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt"> delta <- (t(X) %*% E - t(X2) %*% E2)<br>
W <- W + delta<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">This allocates a big bloc "delta" that you don't need to do.
<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt"><br>
W <- W + (t(X) %*% E - t(X2) %*% E2)<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">Please write back what you find out. I'm always eager to have clear "do this, don't do that" examples in the classroom.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">In my blog, look at item 3. That was a big shocker to me.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">pj<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal">On Wed, Nov 6, 2013 at 12:04 PM, Xavier Robin <<a href="mailto:xavier@cbs.dtu.dk" target="_blank">xavier@cbs.dtu.dk</a>> wrote:<o:p></o:p></p>
<div>
<p class="MsoNormal">On 11/6/13 6:38 PM, Romain Francois wrote:<o:p></o:p></p>
<p class="MsoNormal" style="margin-bottom:12.0pt">This very much depends on the code but there is a good chance that RcppArmadillo will generate code making less data copies, etc ...<br>
<br>
Hard to say without seeing the code.<br>
<br>
Romain<o:p></o:p></p>
</div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Most of the code (or at least the slow, highly repeated parts) look like:<o:p></o:p></p>
<p class="MsoNormal"> A <- t(c + t(W) %*% X)<br>
E <- (exp(A) * (1 - 1 / A) + 1 / A) / (exp(A) - 1)<br>
E[abs(A) < sqrt(.Machine$double.eps) * 2 ] <- 0.5<br>
<br>
B <- t(b + W %*% t(E))<br>
X2 <- 1 / (1 + exp(-B))<br>
<br>
A2 <- t(c + t(W) %*% X2)<br>
E2 <- (exp(A2) * (1 - 1 / A2) + 1 / A2) / (exp(A2) - 1)<br>
E2[abs(A2) < sqrt(.Machine$double.eps) * 2 ] <- 0.5<br>
<br>
delta <- (t(X) %*% E - t(X2) %*% E2)<br>
W <- W + delta<o:p></o:p></p>
<p class="MsoNormal"><br>
Where b and c are vectors, W and X matrices. All this is encapsulated in a function, that is called a few thousand times in a for loop, with some sanity checks. (But it didn't appear to have much impact on the speed... if I remove the matrix operations so it
does nothing, it executes nearly instantly). I understand from Dirk and Douglas that it probably isn't going to make a huge difference, though (not by orders).<o:p></o:p></p>
<div>
<p class="MsoNormal"><br>
<br>
Thanks,<br>
Xavier<br>
<br>
-- <br>
Xavier Robin, PhD<br>
Cellular Signal Integration Group (C-SIG) - <a href="http://www.lindinglab.org" target="_blank">
http://www.lindinglab.org</a><o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><br>
-- <br>
Paul E. Johnson<br>
Professor, Political Science Assoc. Director<br>
1541 Lilac Lane, Room 504 Center for Research Methods<br>
University of Kansas University of Kansas<br>
<a href="http://pj.freefaculty.org" target="_blank">http://pj.freefaculty.org</a>
<a href="http://quant.ku.edu" target="_blank">http://quant.ku.edu</a> <o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>