{"id":6502,"date":"2018-11-20T13:24:48","date_gmt":"2018-11-20T12:24:48","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=6502"},"modified":"2018-11-20T13:24:48","modified_gmt":"2018-11-20T12:24:48","slug":"accelerated-simd-mersenne-twister-random-number-generator-in-matlab","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=6502","title":{"rendered":"Accelerated SIMD Mersenne Twister random number generator in MATLAB"},"content":{"rendered":"

In a recent blog post, Daniel Lemire<\/a> explicitly demonstrated that vectorising random number generators<\/a> using SIMD instructions could give useful speed-ups.\u00a0 This reminded me of the one of the first times I played <\/a>with the Julia language<\/a> where I learned that Julia’s random number generator used a SIMD-accelerated implementation of Mersenne Twister<\/a> called dSFMT to generate random numbers much faster than MATLAB’s Mersenne Twister implementation.<\/p>\n

Just recently, I learned that MATLAB<\/a> now has its own SIMD accelerated version of Mersenne Twister which can be activated like so:<\/p>\n

\r\nseed=1;\r\nrng(seed,'simdTwister')\r\n<\/pre>\n
This new Mersenne Twister implementation gives different random variates to the original implementation (which I demonstrated is the same as Numpy’s implementation in an older post<\/a>) as you might expect<\/p>\n
\r\n>> rng(1,'Twister')\r\n>> rand(1,5)\r\nans =\r\n0.4170 0.7203 0.0001 0.3023 0.1468\r\n>> rng(1,'simdTwister')\r\n>> rand(1,5)\r\nans =\r\n0.1194 0.9124 0.5032 0.8713 0.5324\r\n<\/pre>\n
So it’s clearly a different algorithm and, on CPUs that support the relevant instructions, it’s about twice as fast!\u00a0 Using my very unscientific test code:<\/p>\n
format compact\r\nnumber_of_randoms = 10000000\r\ndisp('Timing standard twister')\r\nrng(1,'Twister')\r\ntic;rand(1,number_of_randoms);toc\r\ndisp('Timing SIMD twister')\r\nrng(1,'simdTwister')\r\ntic;rand(1,number_of_randoms);toc\r\n<\/pre>\n
gives the following results for a typical run on my Dell XPS 15 9560<\/a> which supports AVX instructions<\/a><\/p>\n
number_of_randoms =\r\n    10000000\r\nTiming standard twister\r\nElapsed time is 0.101307 seconds.\r\nTiming SIMD twister\r\nElapsed time is 0.057441 seconds\r\n<\/pre>\n
The MATLAB documentation does not tell us which algorithm their implementation is based on but it seems to be different from Julia’s. In Julia, if we set the seed to 1 as we did for MATLAB and ask for 5 random numbers, we get something different from MATLAB:<\/p>\n
julia> using Random\r\njulia> Random.seed!(1);\r\njulia> rand(1,5)\r\n1\u00d75 Array{Float64,2}:\r\n 0.236033  0.346517  0.312707  0.00790928  0.48861\r\n<\/pre>\n
The performance of MATLAB’s new generator is on-par with Julia’s although I’ll repeat that these timing tests are far far from rigorous.<\/p>\n
\r\njulia> Random.seed!(1);\r\n\r\njulia> @time rand(1,10000000);\r\n  0.052981 seconds (6 allocations: 76.294 MiB, 29.40% gc time)\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
In a recent blog post, Daniel Lemire explicitly demonstrated that vectorising random number generators using SIMD instructions could give useful speed-ups.\u00a0 This reminded me of the one of the first times I played with the Julia language where I learned that Julia’s random number generator used a SIMD-accelerated implementation of Mersenne Twister called dSFMT to […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[53,11,63,67],"tags":[],"class_list":["post-6502","post","type-post","status-publish","format-standard","hentry","category-making-matlab-faster","category-matlab","category-random-numbers","category-scientific-software"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1GS","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6502","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=6502"}],"version-history":[{"count":3,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6502\/revisions"}],"predecessor-version":[{"id":6506,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6502\/revisions\/6506"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=6502"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=6502"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=6502"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}