StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p><strong>The code shown below is available here: <a href="http://pastebin.com/4PWWxGhB" rel="noreferrer">http://pastebin.com/4PWWxGhB</a>. Just copy and paste it into a notebook to test it out.</strong></p> <p>I was actually trying to do several functional ways of calculating matrices, since I figured the functional way (which is typically idiomatic in Mathematica) is more efficient.</p> <p>As one example, I had this matrix which was composed of two lists:</p> <pre><code>In: L = 1200; e = Table[..., {2L}]; f = Table[..., {2L}]; h = Table[0, {2L}, {2L}]; Do[h[[i, i]] = e[[i]], {i, 1, L}]; Do[h[[i, i]] = e[[i-L]], {i, L+1, 2L}]; Do[h[[i, j]] = f[[i]]f[[j-L]], {i, 1, L}, {j, L+1, 2L}]; Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}]; </code></pre> <p>My first step was to time everything.</p> <pre><code>In: h = Table[0, {2 L}, {2 L}]; AbsoluteTiming[Do[h[[i, i]] = e[[i]], {i, 1, L}];] AbsoluteTiming[Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];] AbsoluteTiming[ Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];] AbsoluteTiming[Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];] Out: {0.0020001, Null} {0.0030002, Null} {5.0012861, Null} {4.0622324, Null} </code></pre> <p><code>DiagonalMatrix[...]</code> was slower than the do loops, so I decided to just use <code>Do</code> loops on the last step. As you can see, using <code>Outer[Times, f, f]</code> was much faster in this case.</p> <p>I then wrote the equivalent using <code>Outer</code> for the blocks in the upper right and bottom left of the matrix, and <code>DiagonalMatrix</code> for the diagonal:</p> <pre><code>AbsoluteTiming[h1 = ArrayPad[Outer[Times, f, f], {{0, L}, {L, 0}}];] AbsoluteTiming[h1 += Transpose[h1];] AbsoluteTiming[h1 += DiagonalMatrix[Join[e, e]];] Out: {0.9960570, Null} {0.3770216, Null} {0.0160009, Null} </code></pre> <p>The <code>DiagonalMatrix</code> was actually slower. I could replace this with just the <code>Do</code> loops, but I kept it because it was cleaner looking.</p> <p>The current tally is 9.06 seconds for the naive <code>Do</code> loop, and 1.389 seconds for my next version using <code>Outer</code> and <code>DiagonalMatrix</code>. About a 6.5 times speedup, not too bad.</p> <hr> <p>Sounds a lot faster, now doesn't it? Let's try using <code>Compile</code> now.</p> <pre><code>In: cf = Compile[{{L, _Integer}, {e, _Real, 1}, {f, _Real, 1}}, Module[{h}, h = Table[0.0, {2 L}, {2 L}]; Do[h[[i, i]] = e[[i]], {i, 1, L}]; Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}]; Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}]; Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}]; h]]; AbsoluteTiming[cf[L, e, f];] Out: {0.3940225, Null} </code></pre> <p>Now it's running 3.56 times faster than my last version, and 23.23 times faster than the first one. Next version:</p> <pre><code>In: cf = Compile[{{L, _Integer}, {e, _Real, 1}, {f, _Real, 1}}, Module[{h}, h = Table[0.0, {2 L}, {2 L}]; Do[h[[i, i]] = e[[i]], {i, 1, L}]; Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}]; Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}]; Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}]; h], CompilationTarget->"C", RuntimeOptions->"Speed"]; AbsoluteTiming[cf[L, e, f];] Out: {0.1370079, Null} </code></pre> <p>Most of the speed came from <code>CompilationTarget->"C"</code>. Here I got another 2.84 speedup over the fastest version, and 66.13 times speedup over the first version. But all I did was just compile it!</p> <p>Now, this is a very simple example. But this is real code I'm using to solve a problem in condensed matter physics. So don't dismiss it as possibly being a "toy example." </p> <hr> <p>How's about another example of a technique we can use? I have a relatively simple matrix I have to build up. I have a matrix that's composed of nothing but ones from the start to some arbitrary point. The naive way may look something like this:</p> <pre><code>In: k = L; AbsoluteTiming[p = Table[If[i == j && j <= k, 1, 0], {i, 2L}, {j, 2L}];] Out: {5.5393168, Null} </code></pre> <p>Instead, let's build it up using <code>ArrayPad</code> and <code>IdentityMatrix</code>:</p> <pre><code>In: AbsoluteTiming[ArrayPad[IdentityMatrix[k], {{0, 2L-k}, {0, 2L-k}} Out: {0.0140008, Null} </code></pre> <p>This actually doesn't work for k = 0, but you can special case that if you need that. Furthermore, depending on the size of k, this can be faster or slower. It's always faster than the Table[...] version though.</p> <p>You could even write this using <code>SparseArray</code>:</p> <pre><code>In: AbsoluteTiming[SparseArray[{i_, i_} /; i <= k -> 1, {2 L, 2 L}];] Out: {0.0040002, Null} </code></pre> <p>I could go on about some other things, but I'm afraid if I do I'll make this answer unreasonably large. I've accumulated a number of techniques for forming these various matrices and lists in the time I spent trying to optimize some code. The base code I worked with took over 6 days for one calculation to run, and now it takes only 6 hours to do the same thing.</p> <p>I'll see if I can pick out the general techniques I've come up with and just stick them in a notebook to use.</p> <p><strong>TL;DR: It seems like for these cases, the functional way outperforms the procedural way. But when compiled, the procedural code outperforms the functional code.</strong></p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload