18
$\begingroup$

Edit Greg Hurst answers show how to independently get 3x speed-up using single-threaded MKL implementation and 4x speed-up by using native implementation with multi-threading, can they be combined?

Standard Gaussian sampling is slow in Mathematica. I compared against Python packages and Mathematica is about 1x slower than numpy (OpenBLAS?), 3x slower than PyTorch (MKL?), 12x slower than mtalg (MKL+multithreaded?), 1000x slower than GPU generation (massive parallelism).

Is there something I can do to make these samplers run faster on my MacBook? My first suspicion is that Mathematica implementation is not properly multi-threaded.

m = 1000000;
d = 100;
AbsoluteTiming[RandomVariate[NormalDistribution[], {m, d}];] (* 2.01144 *)

Here's how to compare against mtalg locally.

  1. Install Python

  2. Run pip install mtalg

  3. Run the following code

import time
class timeit:
    def __init__(self, tag=""):
        self.tag = tag

    def __enter__(self):
        self.start = time.perf_counter()
        return self

    def __exit__(self, *args):
        self.end = time.perf_counter()
        interval_ms = 1000 * (self.end - self.start)
        print(f"{interval_ms:8.2f}   {self.tag}")

import mtalg

m=1000000
d=100


with timeit("sample mtalg "):
    a = mtalg.random.standard_normal(size=(m, d))
Improve this question
$\endgroup$
2

4 Answers 4

Reset to default
14
$\begingroup$

We can achieve noticeable speedup with Intel's MKL random number generator, which leverages hardware accelerated code.

m = 1000000;
d = 100;

Default behavior:

AbsoluteTiming[RandomVariate[NormalDistribution[], {m, d}];]

{1.93745, Null}

Single core speed up on Intel machines:

BlockRandom[
  SeedRandom[Method -> "MKL"];
  RandomVariate[NormalDistribution[], {m, d}]; // AbsoluteTiming
]

{0.716416, Null}
Improve this answer
$\endgroup$
2
  • $\begingroup$ Ah, nice. I'm guessing PyTorch is probably also using MKL $\endgroup$
    – Yaroslav Bulatov
    Commented Mar 12, 2023 at 3:18
  • $\begingroup$ background for this problem is given here $\endgroup$
    – Yaroslav Bulatov
    Commented Apr 21, 2023 at 5:52
17
$\begingroup$

Watching Activity Monitor, I don't think RandomVariate is multithreaded.

Luckily random samples from a normal distribution can be compiled. We can exploit Compile’s parallelized Listability to get a speedup.

cf = Compile[{{dummy, _Integer}, {d, _Integer}},
  RandomVariate[NormalDistribution[], d],
  Parallelization -> True,
  RuntimeOptions -> "Speed",
  RuntimeAttributes -> {Listable}
];

Here's a simple comparison on a 2018 6-core Macbook Pro:

m = 1000000;
d = 100;

AbsoluteTiming[RandomVariate[NormalDistribution[], {m, d}];]

{1.93745, Null}

AbsoluteTiming[cf[ConstantArray[1, m], d];]

{0.513919, Null}

Verify the compiled version follows a normal distribution:

Through[{Mean, StandardDeviation}[vals]]

{-0.0000760662, 1.00004}

FindDistribution[vals[[1 ;; -1 ;; 1000]]]

NormalDistribution[-0.00209107, 0.999427]

Histogram[vals]

Improve this answer
$\endgroup$
7
  • $\begingroup$ Ah neat. There's probably another factor of 3x left on the table because from activity monitor PyTorch CPU implementation seems to be single-threaded as well $\endgroup$
    – Yaroslav Bulatov
    Commented Mar 12, 2023 at 2:23
  • 1
    $\begingroup$ Interesting. Changing RandomVariate[NormalDistribution[], d] to RandomReal[{0., 1.}, d] in cf takes the timing from 0.53 to 0.34 on my machine... and so I don't see a way to get even more substantial speed up with this trick. I also tried a custom Box-Muller implementation and no speed up either (as this depends on uniform randoms + other operations). $\endgroup$
    – Greg Hurst
    Commented Mar 12, 2023 at 2:48
  • 1
    $\begingroup$ Do you know how the PyTorch CPU algorithm is implemented? My curiosity is piqued. $\endgroup$
    – Greg Hurst
    Commented Mar 12, 2023 at 2:48
  • 1
    $\begingroup$ It's using some standard library but not sure atm which one. I just tried out mtalg and it was about 4 times faster than PyTorch $\endgroup$
    – Yaroslav Bulatov
    Commented Mar 12, 2023 at 3:03
  • 2
    $\begingroup$ Dang, that’s fast. Is the mtalg algorithm multithreaded? $\endgroup$
    – Greg Hurst
    Commented Mar 12, 2023 at 3:16
12
+200
$\begingroup$

This is my LibraryLink take on this. I use one of the random number generators of the standard library in conjunction with OpenMP.

This will only work on Apple Silicon with OpenMP installed via homebrew.

Needs["CCompilerDriver`"]

Quiet[LibraryFunctionUnload[libf]];
ClearAll[libf];
libf = Module[{lib, code, name},
   name = "libf";
   code = StringJoin["

#include \"WolframLibrary.h\"
#include <algorithm>
#include <random>
#include <omp.h>

// Computes k-th job pointer for job_count equally sized jobs \
distributed on thread_count threads.
template<typename Int, typename Int1, typename Int2>
inline Int JobPointer( const Int job_count, const Int1 thread_count, \
const Int2 k )
{
    return job_count/static_cast<Int>(thread_count)*static_cast<Int>(k) \
+ job_count%static_cast<Int>(thread_count)*static_cast<Int>(k)/static_\
cast<Int>(thread_count);
}

EXTERN_C DLLEXPORT int fun(WolframLibraryData libData, mint Argc, \
MArgument *Args, MArgument Res)
{
    MTensor a_           = MArgument_getMTensor(Args[0]); 
    const mint thread_count = MArgument_getInteger(Args[1]); 
    
    const mint n  = libData->MTensor_getDimensions(a_)[0];
    const mint d  = libData->MTensor_getDimensions(a_)[1];

    //// Create MTensor for the result.
    //MTensor a_;
    //(void)libData->MTensor_new(MType_Real, 1, &n, &a_);
    
    mreal * const a = libData->MTensor_getRealData(a_);


    // Use the potentially slow hardware random number generator for \
seeding.
    // We do this in the sequential code because we cannot rely on \
std::random_device being thread safe.

    std::vector<unsigned int> seeds (4 * thread_count);
    std::random_device r;
    for( mint i = 0; i < 4 * thread_count; ++i )
    {
        seeds[i] = r();
    }

    #pragma omp parallel for num_threads(thread_count) schedule( static \
)
    for( mint thread = 0; thread < thread_count; ++thread )
    {
        std::seed_seq seed { \
seeds[4*thread],seeds[4*thread+1],seeds[4*thread+2],seeds[4*thread+3] \
};
            
        // Create the actual random engine.
        std::mt19937_64 random_engine ( seed );

        std::normal_distribution<mreal> normal_dist {0.,1.};

        const mint i_begin = JobPointer<mint>(n*d,thread_count,thread  );
        const mint i_end   = JobPointer<mint>(n*d,thread_count,thread+1);

        for( mint i = i_begin; i < i_end; ++i )
        {
            a[i] = normal_dist( random_engine );
        }
    }
   
    libData->MTensor_disown(a_);

    return LIBRARY_NO_ERROR;
}"];
   
   lib = CreateLibrary[code, name,
     "Language" -> "C++",
     "CompileOptions" -> {" -Wall", "-Wextra", 
       "-Wno-unused-parameter", "-std=c++11", "-Ofast", "-flto", 
       "-Xpreprocessor -fopenmp", "-lomp"},
     "ShellOutputFunction" -> Print,
     "IncludeDirectories" -> {"/opt/homebrew/opt/libomp/include"},
     "LibraryDirectories" -> {"/opt/homebrew/opt/libomp/lib"}
     ];
   
   LibraryFunctionLoad[lib, "fun", {{Real, 2, "Shared"}, Integer}, 
    "Void"]
   ];

Here a usage example. I allocate the array in Mathematica and let it the LibraryFunction libf only fill it. This way you can safe a bit time for the memory allocation in successive calls,

a = ConstantArray[0., {m , d}]; // AbsoluteTiming
libf[a, 8]; // AbsoluteTiming

{0.037611, Null}

{0.176087, Null}

The CompiledFunction cf by Greg Hurst needs about the same total time. So this is not much of an improvement.

The MKL pseudorandom number generators are probably vectorized.

Edit

I played a bit more with the Xoshiro pseudorandom number generators. In conjunction with lowering the precision of generated numbers from double to single precision, this gave a 3-fold speedup compared to the pure C++ implementation above.

One may or may not consider this cheating. For many applications in probability theory one does not need that many of digits.

In order to get this work, you have to download or clone Ryo Suzuki's repository at https://github.com/Reputeless/Xoshiro-cpp. The code below assumes that the file XoshiroCpp.hpp will be located in the subdirectory Xoshiro-cpp in your home directory. However, you can put it anywhere your want and adjust the option "IncludeDirectories" of CreateLibrary.

Needs["CCompilerDriver`"]

Quiet[LibraryFunctionUnload[libf]];
ClearAll[libf];
libf = Module[{lib, code, name}, name = "libf";
   code = StringJoin["

#include \"WolframLibrary.h\"
#include <algorithm>
#include <random>
#include \"XoshiroCpp.hpp\"
#include <omp.h>

// Computes k-th job pointer for job_count equally sized jobs distributed on thread_count threads.
template<typename Int, typename Int1, typename Int2>
inline Int JobPointer( const Int job_count, const Int1 thread_count, const Int2 k )
{
    return job_count/static_cast<Int>(thread_count)*static_cast<Int>(k) + job_count%static_cast<Int>(thread_count)*static_cast<Int>(k)/static_cast<Int>(thread_count);
}

using namespace XoshiroCpp;

EXTERN_C DLLEXPORT int fun(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor a_              = MArgument_getMTensor(Args[0]);
    const mint thread_count = MArgument_getInteger(Args[1]);
    
    const mint n  = libData->MTensor_getDimensions(a_)[0];
    const mint d  = libData->MTensor_getDimensions(a_)[1];
    
    mreal * const a = libData->MTensor_getRealData(a_);


    // Use the potentially slow hardware random number generator for seeding.
    // We do this in the sequential code because we cannot rely on std::random_device being thread safe.

    std::random_device r;
    std::vector<std::uint64_t> seeds ( thread_count);

    for( mint i = 0; i < thread_count; ++i )
    {
        reinterpret_cast<std::uint32_t*>(&seeds[i])[0] = r();
        reinterpret_cast<std::uint32_t*>(&seeds[i])[1] = r();
    }

    #pragma omp parallel for num_threads(thread_count) schedule( static )
    for( mint thread = 0; thread < thread_count; ++thread )
    {
        // Create the actual random engine.
        Xoshiro256Plus random_engine ( seeds[thread] );

        std::normal_distribution<float> dist {0,1};

        const mint i_begin = JobPointer<mint>(n*d,thread_count,thread  );
        const mint i_end   = JobPointer<mint>(n*d,thread_count,thread+1);

        for( mint i = i_begin; i < i_end; ++i )
        {
            a[i] = static_cast<mreal>(dist( random_engine ));
        }
    }
   
    libData->MTensor_disown(a_);

    return LIBRARY_NO_ERROR;
}"];
   lib = CreateLibrary[code, name,
     "Language" -> "C++",
     "ShellOutputFunction" -> Print,
     "CompileOptions" -> {
       " -Wall", "-Wextra", "-Wno-unused-parameter", "-std=c++17", "-Ofast", "-flto", "-Xpreprocessor -fopenmp", "-lomp"
       },
     "IncludeDirectories" -> {
       "/opt/homebrew/opt/libomp/include"(*Put path to omp.h here.*),
       FileNameJoin[{$HomeDirectory,"Xoshiro-cpp"}](*Put any other path here that contains XoshiroCpp.hpp.*)
       },
     "LibraryDirectories" -> {
       "/opt/homebrew/opt/libomp/lib"(*Put path to libomp.dylib here.*)
       }];
   LibraryFunctionLoad[lib, "fun", {{Real, 2, "Shared"}, Integer},
    "Void"]
   ];

Here a usage example that ran on my M1 Max with 8 threads:

m = 1000000;
d = 100;

TRandomVariate = 
 First@RepeatedTiming[a = RandomVariate[NormalDistribution[], {m, d}];]
TAllocation = First@RepeatedTiming[b = ConstantArray[0., {m, d}];]
TXoshiro = First@RepeatedTiming[libf[b, 8];]

(TRandomVariate - TAllocation)/TXoshiro

1.14246

0.0344541

0.0548854

20.1876

If one factors out the allocation time (approach the buffer can be reused with this approach), this is 20 times faster than RandomVariate.

But please keep in mind: These double precision numbers have been converted from pseudorandom single precision numbers; only the leading 23(?) binary digits are truely random; the remaining ones are set to 0.

Moreover, Xoshiro256+ is not cryptographically safe.

Edit 2

I found a couple of inefficiencies in Apple clang's std::normal_distribution<float>. First, it seems to discard half of the entropy generated. Moreover, it uses a rejection sampler to generate random points in the unit disk. This is faster than using Box-Muller, because it safes to evaluate cos and sin. But STL implementation performs a couple of unnecassary floating-point operations in the case of a rejection. Replacing the use of std::normal_distribution<float> by getNormalFloatPair allowed me to shave off another 20% of runtime.

Needs["CCompilerDriver`"];

Quiet[LibraryFunctionUnload[libNormalDistributionInt]];
ClearAll[libNormalDistributionInt];
libNormalDistributionInt = Module[{lib, code, name},
   name = "libNormalDistributionInt";
   code = StringJoin["

#include \"WolframLibrary.h\"
#include <algorithm>
#include <random>
#include \"XoshiroCpp.hpp\"
#include <omp.h>

// Computes k-th job pointer for job_count equally sized jobs distributed on thread_count threads.
template<typename Int, typename Int1, typename Int2>
inline Int JobPointer( const Int job_count, const Int1 thread_count, const Int2 k )
{
    return job_count/static_cast<Int>(thread_count)*static_cast<Int>(k) + job_count%static_cast<Int>(thread_count)*static_cast<Int>(k)/static_cast<Int>(thread_count);
}

using namespace XoshiroCpp;

inline void getNormalFloatPair( Xoshiro256Plus & random_engine, mreal & a, mreal & b )
{
    std::int64_t ix;
    std::int64_t iy;
    std::int64_t is;
    
    constexpr std::int64_t threshold = (std::int64_t(1) << 46);
                
    do
    {
        const std::uint64_t bits = random_engine();
        
        ix = static_cast<std::int64_t>(
            reinterpret_cast<const std::int32_t*>(&bits)[0] >> 8
        );
        
        iy = static_cast<std::int64_t>(
            reinterpret_cast<const std::int32_t*>(&bits)[1] >> 8
        );
    
        is = ix * ix + iy * iy;
    } while ( is > threshold || is == std::int64_t(0) );
    
    const float x = 0x1.0p-23f * static_cast<float>( ix );
    const float y = 0x1.0p-23f * static_cast<float>( iy );
    const float s = 0x1.0p-46f * static_cast<float>( is );
    
    const float r = std::sqrt(- 2.0f * std::log(s) / s );
    a = r * x;
    b = r * y;
}

EXTERN_C DLLEXPORT int fun(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor a_              = MArgument_getMTensor(Args[0]);
    const mint thread_count = MArgument_getInteger(Args[1]);
    
    const mint n  = libData->MTensor_getDimensions(a_)[0];
    const mint d  = libData->MTensor_getDimensions(a_)[1];
    
    mreal * const a = libData->MTensor_getRealData(a_);


    // Use the potentially slow hardware random number generator for seeding.
    // We do this in the sequential code because we cannot rely on std::random_device being thread safe.

    std::random_device r;
    std::vector<std::uint64_t> seeds ( thread_count);


    for( mint i = 0; i < thread_count; ++i )
    {
        reinterpret_cast<std::uint32_t*>(&seeds[i])[0] = r();
        reinterpret_cast<std::uint32_t*>(&seeds[i])[1] = r();
    }

    #pragma omp parallel for num_threads(thread_count) schedule( static )
    for( mint thread = 0; thread < thread_count; ++thread )
    {
        Xoshiro256Plus random_engine ( seeds[thread] );

        const mint i_begin = JobPointer<mint>(n*d,thread_count,thread  );
        const mint i_end   = JobPointer<mint>(n*d,thread_count,thread+1);

        if( i_end > i_begin )
        {
            const size_t i_begin_odd = i_begin % 2;
            const size_t i_end_odd   = i_end   % 2;
            
            mreal x;
            mreal y;
            
            getNormalFloatPair( random_engine, x, y );
            
            a[i_begin] = x;
            
            for( size_t i = i_begin + i_begin_odd; i < i_end - i_end_odd; i+=2 )
            {
                getNormalFloatPair( random_engine, a[i+0], a[i+1] );
            }
            
            a[i_end-1] = y;
        }
    }
   
    libData->MTensor_disown(a_);

    return LIBRARY_NO_ERROR;
}"];
   lib =
    CreateLibrary[code, name, "Language" -> "C++",
     "ShellOutputFunction" -> Print,
     "CompileOptions" -> {" -Wall", "-Wextra",
       "-Wno-unused-parameter", "-std=c++17", "-Ofast", "-flto",
       "-Xpreprocessor -fopenmp", "-lomp"},
     "IncludeDirectories" -> {"/opt/homebrew/opt/libomp/include"(*Put path to omp.h here.*),
       FileNameJoin[{$HomeDirectory,"Xoshiro-cpp"}](*Put any other path here that contains XoshiroCpp.hpp.*)},
     "LibraryDirectories" -> {"/opt/homebrew/opt/libomp/lib"(*Put path to libomp.dylib here.*)}];
   LibraryFunctionLoad[lib, "fun", {{Real, 2, "Shared"}, Integer},
    "Void"]
   ];

Edit 3

I recently started GPU programming with Apple's Metal framework. So I thought it might be a good exercise to write a GPU implementation. Maybe this is also of interest to you?

Rejection sampling is a bad idea on the GPUs because it causes thread divergence. (The SIMD group size of may GPU is 32, so instead of the acceptance probability $p$ on non-SIMD code we would have the actual acceptance rate of only $p^{32}$.) Thus I use the Box-Muller transform to convert from 2 uniformly random variables to 2 normally distributed random variables. On the CPU this was quite expensive due to extra evaluation of trigonometric functions; but the GPU chews through that quite happily.

GPUs require typically quite a lot of boiler-plate code. This is why I created a public github repository.

https://github.com/HenrikSchumacher/Randomizor

As this uses git submodules, you have to clone with

git clone --recurse-submodules [email protected]:HenrikSchumacher/Randomizor.git

(IIRC, you have to create a free github account for that and to set up an ssh key pair.)

This is a demonstation LibaryLink program:

Needs["CCompilerDriver`"]

(*Put your path that contains Randomizor_Metal_Xoshiro.hpp here.*)
dirRandomizor = FileNameJoin[{$HomeDirectory, "github", "Randomizor"}];
(*Put your path to the OpenMP directory here. The place here is chosen by homebrew on Apple Silicon machines.*)
dirOpenMP = "/opt/homebrew/opt/libomp";

Quiet[LibraryFunctionUnload[cRandomizer]];
ClearAll[cRandomizer];
cRandomizer = Module[{lib, code, name},
   name = "cRandomizer";
   code = StringJoin["

#include \"WolframLibrary.h\"

#include \"mathlink.h\"
#include <string>
#include <cstdint>
#include <ostream>
#include <sstream>

namespace mma
{
    // I borrowed this from Szabolcs Horvat's LTemplate

    WolframLibraryData libData;
    
    inline void print(const char *msg)
    {
        if (libData->AbortQ())
        {
            return; // trying to use the MathLink connection during an abort appears to break it
        }

        MLINK link = libData->getMathLink(libData);
        
        MLPutFunction(link, \"EvaluatePacket\", 1);
        
        MLPutFunction(link, \"Print\", 1);
        
        MLPutString(link, msg);
        
        libData->processMathLink(link);
        
        int pkt = MLNextPacket(link);
        
        if (pkt == RETURNPKT)
        {
            MLNewPacket(link);
        }
    }

    // Call _Mathematica_'s `Print[]`, `std::string` argument version.
    inline void print(const std::string &msg)
    {
        print(msg.c_str());
    }
}


extern \"C\" DLLEXPORT mint WolframLibrary_getVersion()
{
    return WolframLibraryVersion;
}

extern \"C\" DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData)
{
    mma::libData = libData;
    return LIBRARY_NO_ERROR;
}

extern \"C\" DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData)
{
    return;
}

#include <iostream>
#include <algorithm>
#include <random>

#define NS_PRIVATE_IMPLEMENTATION
#define MTL_PRIVATE_IMPLEMENTATION
#include <Foundation/Foundation.hpp>
#include <Metal/Metal.hpp>
#include <Accelerate/Accelerate.h>


#define MATHEMATICA
#include \"Randomizor_Metal_Xoshiro.hpp\"


using namespace Randomizor;

EXTERN_C DLLEXPORT int fun(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor a_                  = MArgument_getMTensor(Args[0]);
    const mint GPU_thread_count = MArgument_getInteger(Args[1]);
    const mint GPU_group_size   = MArgument_getInteger(Args[2]);
    const mint OMP_thread_count = MArgument_getInteger(Args[3]);
    
    const mint n_  = libData->MTensor_getDimensions(a_)[0];
    const mint d_  = libData->MTensor_getDimensions(a_)[1];
    
    mreal * const a = libData->MTensor_getRealData(a_);

    // Initialize a GPU device. This simply chooses the first GPU it finds.
    NS::SharedPtr<MTL::Device> device = NS::TransferPtr(
        reinterpret_cast<MTL::Device *>( MTL::CopyAllDevices()->object(0) )
    );

    // Initialize a random number generator on the GPU.
    Randomizor::Randomizor_Metal_Xoshiro gen_Xoshiro ( device, GPU_thread_count, GPU_group_size, OMP_thread_count );

    // The GPU requires its own buffer; it cannot write directly to a because it generates single precision floats, while a contains doubles.
    // This requests an internal buffer of size n_ * d_. Effectively, the buffers size is round up to multiples of GPU_thread_count.
    gen_Xoshiro.RequireReservoir(n_ * d_);

    // Now we fill the internal buffer with normally distributed floats.
    tic(\"Generate floats -- first launch (warm up)\");
    gen_Xoshiro.Fill_Normal();
    toc(\"Generate floats -- first launch (warm up)\");

    print(\"Successive runs take significantly less time.\");

    tic(\"Generate floats\");
    gen_Xoshiro.Fill_Normal();
    toc(\"Generate floats\");

    tic(\"Generate floats\");
    gen_Xoshiro.Fill_Normal();
    toc(\"Generate floats\");

    tic(\"Generate floats\");
    gen_Xoshiro.Fill_Normal();
    toc(\"Generate floats\");

    // TODO: Use bit-fiddling to convert to doubles, so that GPU can write directly to the buffer.

    // copy_buffer uses CPU threads to convert and copy to output buffer.
    // Typically two CPU cores fully saturate the RAM bandwidth; so do not expect good scaling here.

    tic(\"Conversion to doubles; copy to output buffer.\");
    copy_buffer( gen_Xoshiro.Reservoir(), a, n_ * d_, OMP_thread_count );
    toc(\"Conversion to doubles; copy to output buffer.\");

    libData->MTensor_disown(a_);

    return LIBRARY_NO_ERROR;
}"];
   lib = CreateLibrary[code, name, "Language" -> "C++",
     "ShellOutputFunction" -> Print,
     (*"ShellCommandFunction"->Print,*)
     "CompileOptions" -> {
       " -Wall", "-Wextra", "-Wno-unused-parameter", 
       "-mmacosx-version-min=12.0", "-std=c++20", "-Ofast", "-flto",
       "-Xpreprocessor -fopenmp", "-lomp",
       "-framework mathlink", "-framework Foundation", 
       "-framework Accelerate", "-framework Metal"
       },
     "IncludeDirectories" -> {
       FileNameJoin[{dirOpenMP, "include"}],
       dirRandomizor,
       FileNameJoin[{dirRandomizor, "metal-cpp"}]
       },
     "LibraryDirectories" -> {FileNameJoin[{dirOpenMP, "lib"}]}
     ];
   LibraryFunctionLoad[lib, 
    "fun", {{Real, 2, "Shared"}, Integer, Integer, Integer}, "Void"]
   ];

Here is a usage example adapted to my Apple M1 Max with the 32 core GPU:

GPUThreads = 24576 * 4;
GPUThreadsPerThreadGroup = 1024/2;
OpenMPThreads = 8;

m = 1000000;
d = 100;
cRandomizer[b, GPUThreads, GPUThreadsPerThreadGroup, OpenMPThreads];

This is what this program prints on my device:

Generate floats -- first launch (warm up)...

0.058042 s.

Successive runs take significantly less time.

Generate floats...

0.003435 s.

Generate floats...

0.003026 s.

Generate floats...

0.002688 s.

Conversion to doubles; copy to output buffer....

0.032264 s.

You have to fiddle around with the settings for GPUThreads and GPUThreadsPerThreadGroup to get best performance. Here I set GPUThreads to the total number of threads the 32 GPU cores can offer -- according to the online material I found. For some reason I had to reduce GPUThreadsPerThreadGroup from the maximum 1024 to 512; maybe it is because of scarcity of threadgroup memory? Anyways, with the optimal settings, a sample pass after a warm-up is about 700-800 times faster than RandomVariate. Alas, I am not happy because there is some considerable overhead that makes the total program slower than the CPU version:

  1. The warm-up time is way too long. I have yet to find out how to reduce it by offline compilation of the GPU kernels...

  2. The final copying to the output buffer takes 10 times as long as to sample the single precision random variables. Metal cannot work with doubles directly; this is why this copy operation is required. This is annoying because the GPU uses so-called unified memory, so copying data from the device to CPU's RAM is not required, in principle. I could try to implement the conversion to doubles by bit-fiddling and to put that onto the GPU, too...

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ 800 faster is in line with speed-up I'm getting on T4 GPU, impressed that you got it working! I don't understand the issue double conversion, it's slow, so it be skipped it? Machine learning people are doing everything in float16 :) $\endgroup$
    – Yaroslav Bulatov
    Commented Mar 21, 2023 at 16:55
  • $\begingroup$ Yes, the only issue with floats is that LibraryLink does not like to import them, as Mathematica's floating point numbers in convential arrays are double precision. So if you need these numbers within Mathematica, we would have to figure out how to funnel them into a NumericArray (or something similar). $\endgroup$
    – Henrik Schumacher
    Commented Mar 21, 2023 at 22:22
3
$\begingroup$ 
Needs["ExternalEvaluate`"]
session = StartExternalSession["Python"];

pythonCode = "
import numpy as np

m = 1000000
d = 100
random_array = np.random.normal(size=(m, d))
random_array.tolist()
";

AbsoluteTiming[randomArray = ExternalEvaluate[session, pythonCode];]

DeleteObject[session];
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.