Using the GPU

For an introductory discussion of Graphical Processing Units (GPU)and their use for intensive parallel computation purposes, see GPGPU.

One of Theano’s design goals is to specify computations at an abstractlevel, so that the internal function compiler has a lot of flexibilityabout how to carry out those computations. One of the ways we takeadvantage of this flexibility is in carrying out calculations on agraphics card.

There are two ways currently to use a gpu, one of which only supports NVIDIA cards (CUDA backend) and the other, in development, that should support any OpenCL device as well as NVIDIA cards (GpuArray Backend).

CUDA backend

If you have not done so already, you will need to install Nvidia’sGPU-programming toolchain (CUDA) and configure Theano to use it.We provide installation instructions for Linux,MacOS and Windows.

Testing Theano with GPU

To see if your GPU is being used, cut and paste the following program into afile and run it.

  1. from theano import function, config, shared, sandbox
  2. import theano.tensor as T
  3. import numpy
  4. import time
  5.  
  6. vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
  7. iters = 1000
  8.  
  9. rng = numpy.random.RandomState(22)
  10. x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
  11. f = function([], T.exp(x))
  12. print(f.maker.fgraph.toposort())
  13. t0 = time.time()
  14. for i in range(iters):
  15.     r = f()
  16. t1 = time.time()
  17. print("Looping %d times took %f seconds" % (iters, t1 - t0))
  18. print("Result is %s" % (r,))
  19. if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
  20.     print('Used the cpu')
  21. else:
  22.     print('Used the gpu')

The program just computes the exp() of a bunch of random numbers.Note that we use the shared function tomake sure that the input x is stored on the graphics device.

If I run this program (in check1.py) with device=cpu, my computer takes a little over 3 seconds,whereas on the GPU it takes just over 0.64 seconds. The GPU will not always produce the exactsame floating-point numbers as the CPU. As a benchmark, a loop that calls numpy.exp(x.get_value()) takes about 46 seconds.

  1. $ THEANO_FLAGS=mode=FAST_RUN,device=cpu,floatX=float32 python check1.py
  2. [Elemwise{exp,no_inplace}(<TensorType(float32, vector)>)]
  3. Looping 1000 times took 3.06635117531 seconds
  4. Result is [ 1.23178029  1.61879337  1.52278066 ...,  2.20771813  2.29967761
  5.   1.62323284]
  6. Used the cpu
  7.  
  8. $ THEANO_FLAGS=mode=FAST_RUN,device=gpu,floatX=float32 python check1.py
  9. Using gpu device 0: GeForce GTX 580
  10. [GpuElemwise{exp,no_inplace}(<CudaNdarrayType(float32, vector)>), HostFromGpu(GpuElemwise{exp,no_inplace}.0)]
  11. Looping 1000 times took 0.638810873032 seconds
  12. Result is [ 1.23178029  1.61879349  1.52278066 ...,  2.20771813  2.29967761
  13.   1.62323296]
  14. Used the gpu

Note that GPU operations in Theano require for now floatX to be float32 (see also below).

Returning a Handle to Device-Allocated Data

The speedup is not greater in the preceding example because the function isreturning its result as a NumPy ndarray which has already been copied from thedevice to the host for your convenience. This is what makes it so easy to swap in device=gpu, butif you don’t mind less portability, you might gain a bigger speedup by changingthe graph to express a computation with a GPU-stored result. The gpu_from_hostop means “copy the input from the host to the GPU” and it is optimized awayafter the T.exp(x) is replaced by a GPU version of exp().

  1. from theano import function, config, shared, sandbox
  2. import theano.sandbox.cuda.basic_ops
  3. import theano.tensor as T
  4. import numpy
  5. import time
  6.  
  7. vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
  8. iters = 1000
  9.  
  10. rng = numpy.random.RandomState(22)
  11. x = shared(numpy.asarray(rng.rand(vlen), 'float32'))
  12. f = function([], sandbox.cuda.basic_ops.gpu_from_host(T.exp(x)))
  13. print(f.maker.fgraph.toposort())
  14. t0 = time.time()
  15. for i in range(iters):
  16.     r = f()
  17. t1 = time.time()
  18. print("Looping %d times took %f seconds" % (iters, t1 - t0))
  19. print("Result is %s" % (r,))
  20. print("Numpy result is %s" % (numpy.asarray(r),))
  21. if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
  22.     print('Used the cpu')
  23. else:
  24.     print('Used the gpu')

The output from this program is

  1. $ THEANO_FLAGS=mode=FAST_RUN,device=gpu,floatX=float32 python check2.py
  2. Using gpu device 0: GeForce GTX 580
  3. [GpuElemwise{exp,no_inplace}(<CudaNdarrayType(float32, vector)>)]
  4. Looping 1000 times took 0.34898686409 seconds
  5. Result is <CudaNdarray object at 0x6a7a5f0>
  6. Numpy result is [ 1.23178029  1.61879349  1.52278066 ...,  2.20771813  2.29967761
  7.   1.62323296]
  8. Used the gpu

Here we’ve shaved off about 50% of the run-time by simply not copyingthe resulting array back to the host. The object returned by eachfunction call is now not a NumPy array but a “CudaNdarray” which canbe converted to a NumPy ndarray by the normal NumPy casting mechanismusing something like numpy.asarray().

For even more speed you can play with the borrow flag. SeeBorrowing when Constructing Function Objects.

What Can Be Accelerated on the GPU

The performance characteristics will change as we continue to optimize ourimplementations, and vary from device to device, but to give a rough idea ofwhat to expect right now:

  • Only computationswith float32 data-type can be accelerated. Better support for float64 is expected in upcoming hardware butfloat64 computations are still relatively slow (Jan 2010).
  • Matrixmultiplication, convolution, and large element-wise operations can beaccelerated a lot (5-50x) when arguments are large enough to keep 30processors busy.
  • Indexing,dimension-shuffling and constant-time reshaping will be equally fast on GPUas on CPU.
  • Summationover rows/columns of tensors can be a little slower on the GPU than on the CPU.
  • Copyingof large quantities of data to and from a device is relatively slow, andoften cancels most of the advantage of one or two accelerated functions onthat data. Getting GPU performance largely hinges on making data transfer tothe device pay off.

Tips for Improving Performance on GPU

  • Consideradding floatX=float32 to your .theanorc file if you plan to do a lot ofGPU work.
  • Use the Theano flag allow_gc=False. See GPU Async capabilities
  • Preferconstructors like matrix, vector and scalar to dmatrix, dvector anddscalar because the former will give you float32 variables whenfloatX=float32.
  • Ensurethat your output variables have a float32 dtype and not float64. Themore float32 variables are in your graph, the more work the GPU can do foryou.
  • Minimizetranfers to the GPU device by using shared float32 variables to storefrequently-accessed data (see shared()). When usingthe GPU, float32 tensor shared variables are stored on the GPU by default toeliminate transfer time for GPU ops using those variables.
  • If you aren’t happy with the performance you see, try running your script withprofile=True flag. This should print some timing information at programtermination. Is time being used sensibly? If an op or Apply istaking more time than its share, then if you know something about GPUprogramming, have a look at how it’s implemented in theano.sandbox.cuda.Check the line similar to Spent Xs(X%) in cpu op, Xs(X%) in gpu op and Xs(X%) in transfer op.This can tell you if not enough of your graph is on the GPU or if thereis too much memory transfer.
  • Use nvcc options. nvcc supports those options to speed up somecomputations: -ftz=true to flush denormals values tozeros.,–prec-div=false and –prec-sqrt=false options to speed updivision and square root operation by being less precise. You canenable all of them with the nvcc.flags=–use_fast_math Theanoflag or you can enable them individually as in this example:nvcc.flags=-ftz=true –prec-div=false.
  • To investigate whether if all the Ops in the computational graph are running on GPU.It is possible to debug or check your code by providing a value to assert_no_cpu_op_flag, i.e. _warn, for warning raise for raising an error or pdb for putting a breakpointin the computational graph if there is a CPU Op.

GPU Async capabilities

Ever since Theano 0.6 we started to use the asynchronous capability ofGPUs. This allows us to be faster but with the possibility that someerrors may be raised later than when they should occur. This can causedifficulties when profiling Theano apply nodes. There is a NVIDIAdriver feature to help with these issues. If you set the environmentvariable CUDA_LAUNCH_BLOCKING=1 then all kernel calls will beautomatically synchronized. This reduces performance but provides goodprofiling and appropriately placed error messages.

This feature interacts with Theano garbage collection of intermediateresults. To get the most of this feature, you need to disable the gcas it inserts synchronization points in the graph. Set the Theano flagallow_gc=False to get even faster speed! This will raise the memoryusage.

Changing the Value of Shared Variables

To change the value of a shared variable, e.g. to provide new data to processes,use shared_variable.set_value(new_value). For a lot more detail about this,see Understanding Memory Aliasing for Speed and Correctness.

Exercise

Consider again the logistic regression:

  1. import numpy
  2. import theano
  3. import theano.tensor as T
  4. rng = numpy.random
  5.  
  6. N = 400
  7. feats = 784
  8. D = (rng.randn(N, feats).astype(theano.config.floatX),
  9. rng.randint(size=N,low=0, high=2).astype(theano.config.floatX))
  10. training_steps = 10000
  11.  
  12. # Declare Theano symbolic variables
  13. x = T.matrix("x")
  14. y = T.vector("y")
  15. w = theano.shared(rng.randn(feats).astype(theano.config.floatX), name="w")
  16. b = theano.shared(numpy.asarray(0., dtype=theano.config.floatX), name="b")
  17. x.tag.test_value = D[0]
  18. y.tag.test_value = D[1]
  19.  
  20. # Construct Theano expression graph
  21. p_1 = 1 / (1 + T.exp(-T.dot(x, w)-b)) # Probability of having a one
  22. prediction = p_1 > 0.5 # The prediction that is done: 0 or 1
  23. xent = -y*T.log(p_1) - (1-y)*T.log(1-p_1) # Cross-entropy
  24. cost = xent.mean() + 0.01*(w**2).sum() # The cost to optimize
  25. gw,gb = T.grad(cost, [w,b])
  26.  
  27. # Compile expressions to functions
  28. train = theano.function(
  29.             inputs=[x,y],
  30.             outputs=[prediction, xent],
  31.             updates=[(w, w-0.01*gw), (b, b-0.01*gb)],
  32.             name = "train")
  33. predict = theano.function(inputs=[x], outputs=prediction,
  34.             name = "predict")
  35.  
  36. if any([x.op.__class__.__name__ in ['Gemv', 'CGemv', 'Gemm', 'CGemm'] for x in
  37.         train.maker.fgraph.toposort()]):
  38.     print('Used the cpu')
  39. elif any([x.op.__class__.__name__ in ['GpuGemm', 'GpuGemv'] for x in
  40.           train.maker.fgraph.toposort()]):
  41.     print('Used the gpu')
  42. else:
  43.     print('ERROR, not able to tell if theano used the cpu or the gpu')
  44.     print(train.maker.fgraph.toposort())
  45.  
  46. for i in range(training_steps):
  47.     pred, err = train(D[0], D[1])
  48.  
  49. print("target values for D")
  50. print(D[1])
  51.  
  52. print("prediction on D")
  53. print(predict(D[0]))

Modify and execute this example to run on GPU with floatX=float32 andtime it using the command line time python file.py. (Of course, you may use some of your answerto the exercise in section Configuration Settings and Compiling Mode.)

Is there an increase in speed from CPU to GPU?

Where does it come from? (Use profile=True flag.)

What can be done to further increase the speed of the GPU version? Put your ideas to test.

Note

  • Only 32 bit floats are currently supported (development is in progress).
  • Shared variables with float32 dtype are by default moved to the GPU memory space.
  • There is a limit of one GPU per process.
  • Use the Theano flag device=gpu to require use of the GPU device.
  • Use device=gpu{0, 1, …} to specify which GPU if you have more than one.
  • Apply the Theano flag floatX=float32 (through theano.config.floatX) in your code.
  • Cast inputs before storing them into a shared variable.
  • Circumvent the automatic cast of int32 with float32 to float64:
    • Insert manual cast in your code or use [u]int{8,16}.
    • Insert manual cast around the mean operator (this involves division by length, which is an int64).
    • Notice that a new casting mechanism is being developed.

Solution

GpuArray Backend

If you have not done so already, you will need to install libgpuarrayas well as at least one computing toolkit. Instructions for doing soare provided at libgpuarray.

While all types of devices are supported if using OpenCL, for theremainder of this section, whatever compute device you are using willbe referred to as GPU.

Warning

While it is fully our intention to support OpenCL, as of May 2014this support is still in its infancy. A lot of very useful opsstill do not support it because they were ported from the oldbackend with minimal change.

Testing Theano with GPU

To see if your GPU is being used, cut and paste the following programinto a file and run it.

  1. from theano import function, config, shared, tensor, sandbox
  2. import numpy
  3. import time
  4.  
  5. vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
  6. iters = 1000
  7.  
  8. rng = numpy.random.RandomState(22)
  9. x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
  10. f = function([], tensor.exp(x))
  11. print(f.maker.fgraph.toposort())
  12. t0 = time.time()
  13. for i in range(iters):
  14.     r = f()
  15. t1 = time.time()
  16. print("Looping %d times took %f seconds" % (iters, t1 - t0))
  17. print("Result is %s" % (r,))
  18. if numpy.any([isinstance(x.op, tensor.Elemwise) and
  19.               ('Gpu' not in type(x.op).__name__)
  20.               for x in f.maker.fgraph.toposort()]):
  21.     print('Used the cpu')
  22. else:
  23.     print('Used the gpu')

The program just compute exp() of a bunch of random numbers. Notethat we use the theano.shared() function to make sure that theinput x is stored on the GPU.

  1. $ THEANO_FLAGS=device=cpu python check1.py
  2. [Elemwise{exp,no_inplace}(<TensorType(float64, vector)>)]
  3. Looping 1000 times took 2.6071999073 seconds
  4. Result is [ 1.23178032  1.61879341  1.52278065 ...,  2.20771815  2.29967753
  5.   1.62323285]
  6. Used the cpu
  7.  
  8. $ THEANO_FLAGS=device=cuda0 python check1.py
  9. Using device cuda0: GeForce GTX 275
  10. [GpuElemwise{exp,no_inplace}(<GpuArray<float64>>), HostFromGpu(gpuarray)(GpuElemwise{exp,no_inplace}.0)]
  11. Looping 1000 times took 2.28562092781 seconds
  12. Result is [ 1.23178032  1.61879341  1.52278065 ...,  2.20771815  2.29967753
  13.   1.62323285]
  14. Used the gpu

Returning a Handle to Device-Allocated Data

By default functions that execute on the GPU still return a standardnumpy ndarray. A transfer operation is inserted just before theresults are returned to ensure a consistent interface with CPU code.This allows changing the deivce some code runs on by only replacingthe value of the device flag without touching the code.

If you don’t mind a loss of flexibility, you can ask theano to returnthe GPU object directly. The following code is modifed to do just that.

  1. from theano import function, config, shared, tensor, sandbox
  2. import numpy
  3. import time
  4.  
  5. vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
  6. iters = 1000
  7.  
  8. rng = numpy.random.RandomState(22)
  9. x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
  10. f = function([], sandbox.gpuarray.basic_ops.gpu_from_host(tensor.exp(x)))
  11. print(f.maker.fgraph.toposort())
  12. t0 = time.time()
  13. for i in range(iters):
  14.     r = f()
  15. t1 = time.time()
  16. print("Looping %d times took %f seconds" % (iters, t1 - t0))
  17. print("Result is %s" % (numpy.asarray(r),))
  18. if numpy.any([isinstance(x.op, tensor.Elemwise) and
  19.               ('Gpu' not in type(x.op).__name__)
  20.               for x in f.maker.fgraph.toposort()]):
  21.     print('Used the cpu')
  22. else:
  23.     print('Used the gpu')

Here the theano.sandbox.gpuarray.basic.gpu_from_host() callmeans “copy input to the GPU”. However during the optimization phase,since the result will already be on th gpu, it will be removed. It isused here to tell theano that we want the result on the GPU.

The output is

  1. $ THEANO_FLAGS=device=cuda0 python check2.py
  2. Using device cuda0: GeForce GTX 275
  3. [GpuElemwise{exp,no_inplace}(<GpuArray<float64>>)]
  4. Looping 1000 times took 0.455810785294 seconds
  5. Result is [ 1.23178032  1.61879341  1.52278065 ...,  2.20771815  2.29967753
  6.   1.62323285]
  7. Used the gpu

While the time per call appears to be much lower than the two previousinvocations (and should indeed be lower, since we avoid a transfer)the massive speedup we obtained is in part due to asynchronous natureof execution on GPUs, meaning that the work isn’t completed yet, just‘launched’. We’ll talk about that later.

The object returned is a GpuArray from pygpu. It mostly acts as anumpy ndarray with some exceptions due to its data being on the GPU.You can copy it to the host and convert it to a regular ndarray byusing usual numpy casting such as numpy.asarray().

For even more speed, you can play with the borrow flag. SeeBorrowing when Constructing Function Objects.

What Can be Accelerated on the GPU

The performance characteristics will of course vary from device todevice, and also as we refine our implementation.

This backend supports all regular theano data types (float32, float64,int, …) however GPU support varies and some units can’t deal withdouble (float64) or small (less than 32 bits like int16) data types.You will get an error at compile time or runtime if this is the case.

By default all inputs will get transferred to GPU. You can prevent aninput from getting transferred by setting its tag.target attribute to‘cpu’.

Complex support is untested and most likely completely broken.

In general, large operations like matrix multiplication, orelement-wise operations with large inputs, will be significatlyfaster.

GPU Async Capabilities

By default, all operations on the GPU are run asynchronously. Thismeans that they are only scheduled to run and the function returns.This is made somewhat transparently by the underlying libgpuarray.

A forced synchronization point is introduced when doing memorytransfers between device and host.

It is possible to force synchronization for a particular GpuArray bycalling its sync() method. This is useful to get accurate timingswhen doing benchmarks.

Software for Directly Programming a GPU

Leaving aside Theano which is a meta-programmer, there are:

  • CUDA: GPU programming API by NVIDIA based on extension to C (CUDA C)

    • Vendor-specific
    • Numeric libraries (BLAS, RNG, FFT) are maturing.

  • OpenCL: multi-vendor version of CUDA

    • More general, standardized.
    • Fewer libraries, lesser spread.

  • PyCUDA: Python bindings to CUDA driver interface allow to access Nvidia’s CUDA parallelcomputation API from Python

    • Convenience:

Makes it easy to do GPU meta-programming from within Python.

Abstractions to compile low-level CUDA code from Python (pycuda.driver.SourceModule).

GPU memory buffer (pycuda.gpuarray.GPUArray).

Helpful documentation.

  • Completeness: Binding to all of CUDA’s driver API.

  • Automatic error checking: All CUDA errors are automatically translated into Python exceptions.

  • Speed: PyCUDA’s base layer is written in C++.

  • Good memory management of GPU objects:

Object cleanup tied to lifetime of objects (RAII, ‘Resource Acquisition Is Initialization’).

Makes it much easier to write correct, leak- and crash-free code.

PyCUDA knows about dependencies (e.g. it won’t detach from a context before all memoryallocated in it is also freed).

(This is adapted from PyCUDA’s documentationand Andreas Kloeckner’s website on PyCUDA.)

  • PyOpenCL: PyCUDA for OpenCL

Learning to Program with PyCUDA

If you already enjoy a good proficiency with the C programming language, youmay easily leverage your knowledge by learning, first, to program a GPU with theCUDA extension to C (CUDA C) and, second, to use PyCUDA to access the CUDAAPI with a Python wrapper.

The following resources will assist you in this learning process:

The following examples give a foretaste of programming a GPU with PyCUDA. Onceyou feel competent enough, you may try yourself on the corresponding exercises.

Example: PyCUDA

  1. # (from PyCUDA's documentation)
  2. import pycuda.autoinit
  3. import pycuda.driver as drv
  4. import numpy
  5.  
  6. from pycuda.compiler import SourceModule
  7. mod = SourceModule("""
  8. __global__ void multiply_them(float *dest, float *a, float *b)
  9. {
  10.   const int i = threadIdx.x;
  11.   dest[i] = a[i] * b[i];
  12. }
  13. """)
  14.  
  15. multiply_them = mod.get_function("multiply_them")
  16.  
  17. a = numpy.random.randn(400).astype(numpy.float32)
  18. b = numpy.random.randn(400).astype(numpy.float32)
  19.  
  20. dest = numpy.zeros_like(a)
  21. multiply_them(
  22.         drv.Out(dest), drv.In(a), drv.In(b),
  23.         block=(400,1,1), grid=(1,1))
  24.  
  25. assert numpy.allclose(dest, a*b)
  26. print(dest)

Exercise

Run the preceding example.

Modify and execute to work for a matrix of shape (20, 10).

Example: Theano + PyCUDA

  1. import numpy, theano
  2. import theano.misc.pycuda_init
  3. from pycuda.compiler import SourceModule
  4. import theano.sandbox.cuda as cuda
  5.  
  6. class PyCUDADoubleOp(theano.Op):
  7.  
  8.     __props__ = ()
  9.  
  10.     def make_node(self, inp):
  11.         inp = cuda.basic_ops.gpu_contiguous(
  12.            cuda.basic_ops.as_cuda_ndarray_variable(inp))
  13.         assert inp.dtype == "float32"
  14.         return theano.Apply(self, [inp], [inp.type()])
  15.  
  16.     def make_thunk(self, node, storage_map, _, _2):
  17.         mod = SourceModule("""
  18.     __global__ void my_fct(float * i0, float * o0, int size) {
  19.     int i = blockIdx.x*blockDim.x + threadIdx.x;
  20.     if(i<size){
  21.         o0[i] = i0[i]*2;
  22.     }
  23.   }""")
  24.         pycuda_fct = mod.get_function("my_fct")
  25.         inputs = [storage_map[v] for v in node.inputs]
  26.         outputs = [storage_map[v] for v in node.outputs]
  27.  
  28.         def thunk():
  29.             z = outputs[0]
  30.             if z[0] is None or z[0].shape != inputs[0][0].shape:
  31.                 z[0] = cuda.CudaNdarray.zeros(inputs[0][0].shape)
  32.             grid = (int(numpy.ceil(inputs[0][0].size / 512.)), 1)
  33.             pycuda_fct(inputs[0][0], z[0], numpy.intc(inputs[0][0].size),
  34.                        block=(512, 1, 1), grid=grid)
  35.         return thunk

Use this code to test it:

  1. >>> x = theano.tensor.fmatrix()
  2. >>> f = theano.function([x], PyCUDADoubleOp()(x))  
  3. >>> xv = numpy.ones((4, 5), dtype="float32")
  4. >>> assert numpy.allclose(f(xv), xv*2)  
  5. >>> print(numpy.asarray(f(xv)))

Exercise

Run the preceding example.

Modify and execute to multiply two matrices: x * y.

Modify and execute to return two outputs: x + y and x - y.

(Notice that Theano’s current elemwise fusion optimization isonly applicable to computations involving a single output. Hence, to gainefficiency over the basic solution that is asked here, the two operations wouldhave to be jointly optimized explicitly in the code.)

Modify and execute to support stride (i.e. to avoid constraining the input to be C-contiguous).

Note

  • See Other Implementations to know how to handle random numberson the GPU.
  • The mode FAST_COMPILE disables C code, so also disables the GPU. Youcan use the Theano flag optimizer=’fast_compile’ to speed upcompilation and keep the GPU.