Parallel Forall Mon, 18 Jul 2016 18:01:09 +0000 en-US hourly 1 GPU-Accelerated Black Hole Simulations Tue, 12 Jul 2016 05:17:11 +0000 In February 2016, scientists from the Laser Interferometry Gravitational-Wave Observatory (LIGO) announced the first ever observation of gravitational waves—ripples in the fabric of spacetime. This detection directly confirms for the first time the existence of gravitational radiation, a cornerstone prediction of Albert Einstein’s general theory of relativity. Furthermore, it suggests that black holes may occur in nature considerably more commonly than anticipated, and pioneers an entirely new form of astronomy based on direct observations of gravity rather than of light.

Artist’s illustration of the merging binary black hole systems GW150191 (left) and GW151226 (right). Image Courtesy Caltech/MIT/LIGO Laboratory. Ccredit: LIGO/A. Simonnet.
Artist’s illustration of two merging binary black hole systems GW150191 (left) and GW151226 (right). While they are portrayed together here, these two merging systems were detected at different times and in different parts of the sky. Image courtesy Caltech/MIT/LIGO Laboratory. Ccredit: LIGO/A. Simonnet.
The LIGO Laboratory operates two detector sites, one near Hanford in eastern Washington, and another near Livingston, Louisiana. This photo shows the Livingston detector site. Image Courtesy Caltech/MIT/LIGO Laboratory.

The gravitational wave experiment LIGO has now reported the first few detections (of many to come) of the gravitational radiation from binary black holes (BBH: a pair of closely orbiting black hole). These detections are made by firing lasers at mirrors at the ends of two widely separated perpendicular tunnels, called “arms”: upon return from the mirrors, the interference pattern generated by the lasers measures the relative length of the arms. Since gravitational waves contract and expand the arms at different rates, the waves can be measured by studying the fluctuations in the interference pattern.

Unfortunately many other things can expand and contract the arms besides gravitational radiation, such as thermal fluctuations in the mirrors, traffic, and in one (perhaps apocryphal) incident, rifle fire directed at the experiment by bored rural Louisianans. The original detection GW150194 was an exception: the typical gravitational wave signal is expected to be far too much weaker than this noise to be visible within the LIGO data stream to the naked eye.

LIGO scientists instead extract the signal using a technique called matched filtering. First, a template bank of predicted gravitational waveforms from binaries of different system parameters—different masses and spins—is generated. The correlation between each template in the bank and each segment of the actual experimental data stream is then computed. Roughly, if a particular template’s correlation exceeds a threshold value, that template can be said to lie within the data, and a binary system has been detected.

For this to work, we need very accurate predictions of the gravitational radiation from binaries of differing system parameters to serve as templates. Newtonian gravity is too inaccurate for this purpose (in particular, it predicts no gravitational radiation). The system needs instead to be modelled using the full machinery of general relativity (GR). But the problem is far too complicated for simple pen-and paper techniques. We thus must integrate the full nonlinear Einstein equations numerically.

Here are the Einstein equations in their four-dimensional presentation:


The objects on the left-hand side encode the curvature of spacetime. In “flat” spacetime (with zero curvature) two small rocks released at a small initial separation will follow a straight line through spacetime: time will progress forward, but their position in space will not change. In “curved” spacetime (with non-zero curvature) the rocks will instead follow a curved path, so their position in space will change with time. That’s a fancy way of saying that “the rocks will fall”. Thus curvature encodes the effects of gravity.

The right-hand side encodes the matter content of spacetime. One way for the curvature to be non-zero is thus for the right-hand side to also be non-zero: for spacetime to contain matter. In that case the spacetime will curve in such a way as to bend the little pebbles towards the matter. Thus, matter gravitates.

The Einstein equations in this form describe the entire spacetime geometry across both space and time. This is elegant, but not particularly useful for numerical purposes.  We need to instead cast the equations into a 3+1 formulation which allows them to be solved as an initial value problem. To do this, we divide 4-dimensional spacetime into a stack of 3-dimensional spatial slices, each corresponding to a particular “moment” in time, as Figure 1 shows. There are an infinite number of ways to cut spacetime up, some of which will behave better numerically than others, but we need to pick one to do the simulation.

Figure 1: Schematics of the 3 + 1 formulation
Figure 1: Schematics of the 3+1 formulation. Left: spacetime (represented as the three-dimensional figure) is divided into spatial slices (the two-dimensional blue sheets). Each slice is labelled by a particular moment in time. Alternative slicings are possible: this reflects the freedom in general relativity to choose which events in spacetime are “simultaneous” with one another. The dot at t=t3 shows a particular event, and the double cone encloses the other events causally connected to that one. The dotted line shows the trajectory of a particular observer through spacetime. Right: a black hole’s event horizon (the grey surface) appears as a lower-dimensional object on each slice. Image credit: Prayush Kumar.

Post-slicing, the Einstein equations break into a system of four (groups of) PDEs. There are many equivalent such systems; here is one of them (the so called “ADM equations”: physically, but not numerically, equivalent to the ones we actually solve):

R + K^2 - K_{ij}K^{ij} = 16\pi e
D_jK^j_i - D_iK = 8\pi j_i
\partial_t K_{ij} = -D_i D_j \alpha + \alpha(R_{ij} - 2K_{ik}(R_{ij} - \frac{1}{2}\gamma_{ij}(S - e))K^{kj} + KK_{ij} - 8\pi\alpha + \mathcal{L}_\beta K_{ij}
\partial_t \gamma_{ij} = -2\alpha K_{ij} + \mathcal{L}_\beta \gamma_{ij}

It’s beyond the scope of this post to define all these symbols, but it is worth pointing out that the equations come in two groups. The first two do not depend directly on time and are mathematically elliptic PDEs. They are known as the constraint equations because they constrain the gravitational field of each time slice. If they are violated at any moment, the field at that moment is not allowed by general relativity. Therefore, by solving the constraints, we can obtain initial conditions for our simulation. Figure 2 shows an example.

Figure 2: Initial data for a BBH simulation.
Figure 2: Initial data for a BBH simulation at mass ratio 10 (one black hole 10 times more massive than the other). The colours on the simulation volume show a measure of the local curvature of space, while those on the holes themselves show a measure of the horizon’s curvature. A slice of the simulation grid is overlaid. Note the difference in fineness of the grid near the smaller and larger holes. This is one of the reasons why simulations involving large mass ratios are difficult. Image credit: Harald Pfeiffer.

The second two ADM equations involve time directly and are mathematically hyperbolic. We call them the evolution equations because they, given a gravitational field which satisfies the constraint equations at some time, determine how that field evolves for all times. Thus, by integrating these equations we can find the future state of the field and model the black hole merger.

This sets the essential program. First, we choose masses, spins, separations, and orbital speeds for the black hole binary. Next we solve the constraint equations to find a valid initial BBH configuration with those parameters. Finally, we integrate the evolution equations forward to find how the BBH appears in the future. We can then extract the gravitational waveform from the simulation and compare it to experimental data.

The following three videos show visualizations of simulations meant to model GW150914 (first, second), the first inspiral detected by LIGO between one 36- and one 29-solar-mass black hole, and the second detection GW151226 (third), which was a merger between one 14.2 and one 7.5 solar mass black hole. The first visualization uses a ray tracing algorithm developed by the SXS collaboration to realistically compute the actual appearance of the black holes in space. The second shows the apparent horizons of the black holes atop a coloured plane visualizing the curvature of spacetime and a vector field showing the trajectories of nearby particles. The third, whichi is of GW151226, shows the frequency and strain of the gravitational radiation from the inspiral as viewed on Earth.

While watching the simulations you may notice the black holes evolve through four essential phases. During inspiral they orbit around one another, emitting gravitational energy which causes their orbit to decay. Eventually the dissipation becomes so fast that the motion is more of a fall than an orbit, and the holes are said to plunge. Quickly following plunge is the merger, when the holes join together into a third, larger black hole which may be highly asymmetric. That hole then radiates its asymmetries away during ringdown, leaving behind an isolated, symmetric black hole remnant.

BBH simulations within a certain parameter range are now fairly routine. There is, unfortunately, an essential problem with using numerical integrations to create a template bank in practice: the simulations are very, very expensive, consuming tens of thousands of FLOPS per grid point per time slice. Even the most computationally friendly regions of parameter space can take days of wallclock time to run, and more difficult simulations take months.

The computational complexity depends importantly on two system parameters. A simulation with a higher mass ratio involves one black hole which is larger than another. Therefore, the simulation must simultaneously resolve two distance scales; this shrinks the simulation time step. Simulations with higher spins, in turn, involve more small-scale fluctuations in the gravitational field; this also shrinks the time step.

Unfortunately, quickly spinning high-mass-ratio systems are among the most interesting to simulate: these systems show the strongest relativistic effects, and furnish easy comparison with analytic theory (see Figure 3). We therefore seek to accelerate the simulations.

Figure 3: Applicability of direct numerical integration (orange box) compared to “pen-and-paper” approximation techniques (red and blue boxes) for two-body mergers. The y scale shows the compactness of the merging objects (denser objects are more compact), while the x shows the mass ratio. To compare with theory we need to perform the very expensive simulations in the overlapping regions. Image credit: Marten van de Meent.

GPU Acceleration

The computations involved in the BBH integration are both highly compute-bound and highly parallelizable. This makes them (technologically) an ideal candidate for GPU acceleration. In practice, however, existing black hole codes are quite difficult to port, simply because the codes are very large and complicated.

We have chosen as target a leading black hole code, the Sppectral EEinstein Code “SpEC”. SpEC employs so-called spectral methods to solve the Einstein equations, which gives favourable algorithmic scaling at the cost of fragility and code complexity. Indeed the code is simply too large to either rewrite from scratch or port line by line.

Thus, we have taken a two-pronged strategy. By profiling the code we first identify “hot spots”: a small number of modules which consume greater than 20% or so of the overall run time. A few of these turn out to be simple enough to write equivalent GPU kernels by hand. Some, however, are not. For these we exploit the fact that our functions can be viewed as linear operations on data. Thus, for each function there exists a matrix \mathbf{M} such that multiplication of the input data by M gives equivalent output.

When the simulation begins we create caches of the matrices \mathbf{M}. When one of the hot spot modules is run for the first time, we trace out the matrix by feeding successive delta functions through the module. Afterwards we can perform the same effect as the module by a simple matrix multiplication, using CUDA-accelerated BLAS operations (cuBLAS). This also improves portability, since CPU BLAS can easily be substituted.

Figure 4: Speedup of GPU-accelerated vs single-CPU version of two “hotspot” modules of SpEC (purple and goldenrod) using the matrix-multiplication approach described above. The black line shows the speedup of the entire code simulating an isolated black hole.In the accelerated version used here, this speedup is limited by kernel launch time of small operations (see below). The comparison to a single CPU is appropriate in our case, since due to scaling constraints and low-level design choices in the code SpEC would replace each single CPU core with an entire GPU. This and all subsequent benchmarks compare a single NVIDIA Tesla M2090 GPU with one core of an Intel Xeon (Sandybridge) E5-2520 2.0GHz CPU.

Tensor Expressions

Having ported the hot spots using hand-written or matrix multiplication kernels we are still left with a large number of smaller operations to be performed on the CPU. This drastically impairs the net speedup since a full CPU-GPU data transfer is thus required at each hot spot. We therefore need some way to keep all the data on the GPU, except at the end of each time step when inter-node communication is required.

Mathematically, the operations we are performing take the form of tensor expressions. Given a vector space, a tensor T is an assignment of a multidimensional array to that space’s basis \mathbf{f}. Should the vector space undergo a change of basis given by a transformation matrix R, a tensor transforms roughly via a multiplication of each component of that array by either R or its inverse. Tensors on more general spaces such as manifolds can then be defined in terms of the tangent spaces to the manifold. For numerical purposes, however, we will view tensors as simply fields of multidimensional arrays, one for each point of the simulation grid.

For example, consider the internal forces—the stress—felt by a small three-dimensional object. Each of three planes through the object (xx, xy, and xz) feels a force which itself has three components (x, y, and z). Thus to fully describe the stress on the object we need nine numbers: the x-direction force on the xx plane, the y-direction force on the xy plane, etc. It is convenient to bundle these objects together into a single “stress tensor” \sigma_ij.

i and j here are called the indices of the tensor. Each index has a dimension specifying the possible values it may take: in this case the dimension of both indices is 3, since each may take three possible values. The total number of indices in the tensor is called the tensor’s rank. Thus, \sigma_ij is a rank-2 tensor. It may happen to be true that the same number is always obtained when the tensor indices are swapped: \sigma_ij=\sigma_ji . In this case we say the tensor has a symmetry upon the relevant indices.

Tensors of the same rank can be added, multiplied, or divided component-wise. They can also be contracted over an index to produce a new tensor of lower rank. Suppose, for example, that we wish to compute the stress experienced by each plane of an object in some particular direction. We can compute this stress—which will have three components and thus be a rank-1 tensor (a vector) T_j—by contracting \sigma_ij with a vector n_i that points in a particular direction:

T_j = \sigma_{ij}n_i \equiv \sum_i \sigma_{ij} n_i.

We express contraction by writing two tensors beside one another with a repeated index: this always implies a sum over that index.

Most operations in a physics code—and certainly in a general relativity code—can be viewed as operations on tensors. Traditionally in an HPC environment we would encode such operations using nested for loops. For example, consider the “index raising” operation, ubiquitous to any general relativity code:

c_{kj} = a_{ik}*b_{ij}, i,k,j \in 0...N.     (1)

The above operation might be implemented in C using code like the following.

  for (int k=0; k<N; ++k) {     \\*
    double sum=0;
      sum += a[i][k]*b[i][j];

These loops are tedious and hard to maintain. For example, suppose that we wished to express the same operation upon a b_ij with a symmetry (so b_ij = b_ji): this would amount to replacing the line of code marked \\* above with the following line.

for (int k=0; k<j; ++k) {

This is a very subtle change which can be easily missed while debugging. Such errors become difficult to spot as the complexity and number of expressions increases, and worse, they often compile normally. Creating and maintaining a CUDA port of such code greatly magnifies the difficulty.

On the other hand, given an arbitrary tensor expression, the appropriate code block in either C or CUDA is actually sufficiently predictable that its generation can feasibly be automated. In our original approach to this problem we used C++ expression templates to automatically port the array operation within the for loop to the GPU. Each operation instantiates a templated class BinaryOp&lt;LHS, OP, RHS&gt;, which in the above case would take the following form.

BinaryOp<Array, Assign, BinaryOp<Array, Multiplication, Array>>

Such expressions can be automatically ported to the GPU by compiling the code into a library, linking to an executable that loops over each unique BinaryOp in the library and generates a CUDA kernel from each, then linking those kernels back to a second compilation. Because all computations can via this approach be performed on the GPU, it is no longer necessary to perform repeated CPU-GPU synchronizations. This creates, however, a different problem, illustrated in Figure 5: the sheer number of CUDA calls becomes so large that call latency impairs the overall runtime dramatically.

Figure 5: Wall clock time in seconds of SpEC simulating an isolated black hole on a single CPU without (left) and with (right) GPU. The total GPU-accelerated runtime is overlaid upon the CPU plot (black line marked t_GPU). The stacked coloured regions show the contribution to the overall run time of selected “hot spot” modules. The green region shows the contribution from the rest of the code.

Notice in particular how the relative size of the shaded green region has increased in GPU-SpEC (right). Also notice that the scaling of this region is roughly constant. This region is dominated by kernel launch overhead of the many small ported operations described above: each iteration of the loop suffers about 20 microseconds of GPU latency, which ends up dominating the overall run time.

We thus need to bundle our operations together into larger units. To this end we have developed a C++ package, “TenFor”, enabling expressions syntactically very similar to (1) to be written directly into C++ source code. In this way we can avoid the for loops entirely. Using TenFor, (1) is expressed directly in source code as

c(k_,j_) = Sum(a(i_,k_)*b(i_,j_));

While its symmetric equivalent would be

c(Sym<0,1>, k_, j_) = Sum(a(Sym<0,1>, i_, k_)*b(i_,j_));

These expressions are presented to the compiler as a recursive structure of variadic templates, labelling the operations to be performed, the index structures of the tensors, and any relevant symmetries such that mistakes yield errors at worst at link time. In this case, the expression structure would be

TBinaryOp<Tensor, k, j, Assign, TBinaryOp<Null, Sum, <TBinaryOp<Tensor, i, j, Multiply, Tensor<Array<i, j>>>>>

Using a Makefile option, TenFor can be configured to implement expressions in one of three ways. In Nonaccel mode expressions are compiled directly and evaluated using a nested series of expression templates, providing similar functionality to existing tensor packages.

Novel functionality is provided by the other two modes, Accel_CPU and Accel_CUDA. In these modes, user code containing tensor expressions is first compiled into a library that links TenFor. TenFor adds to that library a registry of each unique tensor expression in the source code. The user library is then linked to the TenFor executable CodeWriter, which outputs a C++ (CUDA) function (kernel) from each. These functions are then written to disk as a text file. A second compilation of the user code links the new functions, which perform the operations.

This allows tensor expressions to be written in source code in a compact and maintainable way that directly expresses their mathematical purpose. The expression can be trivially ported from CPU to GPU using CodeWriter, allowing the entire code base to remain on the GPU. In general, the automatically ported code achieves peak efficiency for small, memory-bound tensor operations, with efficiency improving as the operation grows more complex.

Figure 6: Speedup vs 1 CPU core of the operation (1) using Black: automatic porting of individual arrays; Red: TenFor porting, sending pointers to the tensor components to the GPU at each kernel launch time; Blue: TenFor porting, sending pointers to the tensor components once at construction (currently implemented).


While the full port of the production BBH code is not yet complete, our acceleration of the single black hole test case shows promising results. The rapid porting afforded by the TenFor C++ package has been a big help. We are currently in the process of creating an open-source port of TenFor which does not depend directly upon SpEC. This will be hosted at Bitbucket:

For more information on gravitational waves and black holes, I recommend you visit the visit the highly informative Simulating Extreme Spacetimes project website, as well as the LIGO homepage.

Thumbnail image courtesy Caltech/MIT/LIGO Laboratory.

]]> 0
Meet Laconia: Michigan State’s New Supercomputer Thu, 07 Jul 2016 05:41:30 +0000 Laconia supercomputer ribbon cutting.
MSU President Lou Anna K. Simon cutting the ceremonial ribbon on Laconia.


The new supercomputer will be used by Michigan State University researchers for a variety of projects, from running molecular simulations for new drug designs to developing technology for human recognition.

To the Point: A Spartan Tribute

iCER named Laconia after a region of Greece that was home to the original Spartans, the mascot of Michigan State University. It’s also the origin of the English word laconic, which means concise or succinct in speech. iCER thought this was fitting because when you’re working hard, you don’t have time for idle chat. People who have been described as laconic include Socrates, Abraham Lincoln and Albert Einstein: pretty good company!

High Performance in High Demand at MSU

Michigan State University (MSU) is a Research 1 university with a burgeoning computational sciences focus, which has greatly increased demand on its existing computational resources. Researchers at MSU tackle biological, chemical materials, engineering, digital humanities, etc. types of problems buttressing the need for general computational resources. The new machine doubles MSU’s existing capabilities and on the first days was running at 99% utilization, demonstrating the high level of demand for computing at MSU.

The new machine cost nearly $3M with funds coming from university and federal research funds as appropriate. In the end, a large portion of the machine will be sold in a buy-in program making it largely supported by the user community. The machine is part of the Institute of Cyber-enabled Research (iCER) which has a staff of about 20 people that will support the machine via system administration or via support of the researchers themselves.

Laconia has 320 nodes in seven racks, each equipped with two Intel Xeon E5-2680v4 2.4GHz 14-core CPUs, 128 gigabytes of RAM. The cluster also has 50 Lenovo NeXtScale nx360 server nodes with four NVIDIA K80 GPUs each, for a total of 200 dual-GPU accelerators in Laconia. Users have access to 240 GB of solid-state disk local scratch storage and the nodes are connected with InfiniBand EDR/FDRnetworking.

The machine delivers 535.9 teraFLOP/s of performance (Linpack Rmax) and it ranked #205 of the TOP500 supercomputer list (announced on June 20th 2016 at ISC2016). It is also expected to rank #30 on the Green500 list.

Applications of the Laconia Supercomputer

Laconia will be used mostly by local researchers at MSU.  The following sections detail some example uses.

Using GPUs to Detect Neutrinos at IceCube Neutrino Laboratory in Antarctica

Tyce DeYoung, Associate Professor in the Michigan State University Department of Physics and Astronomy, is deputy spokesperson for the IceCube Collaboration, an interdisciplinary project touching on both astrophysics and fundamental particle physics.

Figure 1: A cross section of the neutrino detectors, set up deep within the ice of Antarctica at IceCube Neutrino Laboratory.

On the astrophysical side, it has been known for over a century that the Earth’s atmosphere is constantly bombarded by cosmic rays: protons and atomic nuclei coming from somewhere out in space. We still don’t know their sources, but these cosmic ray particles are somehow accelerated to energies more than ten million times higher than the world’s most powerful particle accelerator—the Large Hadron Collider in Geneva—can manage. Researchers expect that neutrinos detectable by IceCube should also be produced in the sources of the cosmic rays, whatever they are, and they hope the neutrinos IceCube detects will help identify these sources and probe the physics taking place in those environments, which are far more powerful than anything that can be replicated here on Earth.

IceCube is also trying to understand neutrinos themselves better. Neutrinos are fundamental particles, like electrons or quarks, but their properties are very unlike anything else we know. They are more than one million times less massive than any other particle—it’s unknown exactly how much lighter they are because nobody has yet been able to measure a mass that small. They don’t interact with other matter via either electromagnetic or strong nuclear forces, which is why they’ve been called “ghost particles”—most of them will happily pass right through a planet without interacting with it at all. We don’t understand why neutrinos have such odd properties, or whether they might be related to other things we can’t explain, such as why the universe is made of matter rather than antimatter, why all the particles we know about come in triplets, and how fundamental forces like electromagnetism and gravity are related.

The fundamental challenge in neutrino physics is detecting enough of them. They usually are passing through matter without interacting with it, and this includes the IceCube’s detectors. Huge detectors are needed to detect the rare neutrinos that do interact with other matter. IceCube takes that to the limit by monitoring a cubic kilometer of the Antarctic ice cap—one billion tons of ultrapure ice—to watch for neutrinos interacting.

IceCube is designed to watch as much ice as possible with the 5,100 sensors it has available, which means that the sensors are spaced tens or hundreds of meters apart. That means IceCube doesn’t record very much information about any single event, so the project has to make the absolute best use of every bit of data recorded. IceCube must model every one of the tens of millions of photons emitted during a neutrino event since it doesn’t know which hundred or thousand of them it might detect. This modeling must be done very quickly, since there are several thousand other particles passing through the detector every second, or about one million times more than the neutrinos of interest. This isn’t done in real-time, but it gives an idea of the scale of the problem.

Michigan State University is one of 48 U.S. institutions participating in IceCube, and its science relies heavily on GPU-accelerated software running on high-performance computing clusters at MSU and around the world. IceCube’s photon tracking software is all GPU-accelerated; making it almost two orders of magnitude faster than it would be running on CPUs. Using GPUs is the only way researchers have been able to do the tracking directly, rather than resorting to numerical approximations that would degrade ability to understand the detected neutrino events.

IceCube also uses GPUs to model the passage of low-energy neutrinos through the Earth, where some rather complicated quantum mechanical phase effects can cause the neutrinos to change type (called the neutrinos’ “flavor”). Measuring the details of these flavor changes can provide information about neutrino properties like the differences between their masses, in a way that’s analogous to using a tuning fork to measure how far out of tune a piano string is. This is another significant computational challenge; the code runs about 75 times faster on GPUs than it does on CPUs.

In 2013, the IceCube Collaboration announced the discovery of a flux of very high-energy neutrinos, presumably associated with the sources of high-energy cosmic rays. Researchers are currently working on identifying the specific sources of these neutrinos. They have ruled out one class of objects—exploding giant stars—that many people had thought produced the highest energy cosmic rays. The data is also starting to disfavor another popular theory—supermassive black holes at the centers of galaxies—but the jury’s still out on that one.

IceCube is also making progress studying neutrino properties, with a recent paper submission reporting that they don’t see evidence of the long-conjectured “sterile” neutrinos, and the MSU group hopes it will have another major paper measuring some parameters related to neutrino masses and flavor oscillations later this year.

According to Professor DeYoung, “This research will give us a bit more about the universe we live in — and historically, the better we have understood our universe, the more able we have been to provide a comfortable and healthy life for ourselves. Ben Franklin didn’t envision electric lights, much less computers, when he was flying his kite, and Albert Einstein wasn’t thinking about fiber-optic communications when we was working out the quantum mechanics underlying lasers. I suspect that almost every element of modern technology can be traced back to this sort of purely curiosity-based inquiry eventually.  A hundred years from now, what we learn from research like this may play a key role in some sort of technology we can’t even imagine today.”

Using GPUs in Biometrics: Developing the Technology for Human Recognition

Automated human recognition through distinctive anatomical and behavioral traits is a growing area of research called biometrics. This technology can be used to reduce identity theft, improve homeland security, identify suspects, and prevent financial fraud.

Figure 2: An example of a large-scale face search problem.

Emerging applications in face and fingerprint recognition, such as supporting identification tasks involving tens of millions of identities (homeland security, national ID, law enforcement applications) and more difficult capture scenarios (surveillance imagery) have motivated the application of deep learning methods. The success of deep learning methods in general has been driven by the availability of large training datasets, and the computational power afforded by modern GPUs.

University Distinguished Professor Anil Jain in the Michigan State University Department of Computer Science & Engineering, along with his Ph.D. student, Charles Otto, uses deep learning methods for developing face representations for large-scale search and clustering applications. Using the GPUs at MSU’s High Performance Computing Center, a 10-convolutional layer network with millions of parameters can be trained in just two days, resulting in an effective representation for challenging face recognition applications. In a face search experiment involving photos of the Tsarnaev brothers, convicted of the Boston Marathon bombing, the proposed cascade face search system could find the photo of Dzhokhar Tsarnaev, the younger brother, at rank-1 in one second on a 5-million-image gallery and at rank-8 in seven seconds on an 80-million-image gallery.

Using GPUs to Improve Drug Efficiency

Figure 3: The inhibitor TPPU emerging from a deeply buried binding pocket of the drug target sEH (soluble epoxide hydrolase), with selected water molecules shown in blue. The mean residence time of TPPU in sEH is about 20 minutes, and the Dickson Lab develops computational techniques that can sample even these extremely rare transitions.

From anti-cancer therapeutics to novel antibiotics, the world is in constant need of new drugs to treat common and rare diseases. These drugs can dramatically improve the quality of life around the globe. A key piece of the design process is a physical model of the drug binding to its intended target. This model allows researchers to predict how making small changes to the drug can improve its potency and selectivity.

But the targets, usually proteins, are not static. Each atom in the system is in constant motion, and these motions occur on a broad range of timescales ranging from femtoseconds (bond vibrations) to hours (drug residence times). Bridging this gap with molecular simulation cannot be done in a straightforward way. It requires new computational approaches like the ones being developed in the research group led by Alex Dickson, Assistant Professor in the Michigan State University Department of Computational Mathematics, Science and Engineering; Department of Biochemistry & Molecular Biology.

GPUs have changed the game by significantly lowering the cost of molecular simulation. Microsecond timescale simulations are now routine, and moderately sized GPU clusters can collect simulation data over hundreds of microseconds by running in parallel, which can be synthesized into a single statistical model that describes drug binding and release. With the sampling methods they are developing, Dickson’s group aims to push this even further, to simulate unbinding transitions that occur on extremely rare time scales on the molecular level but are vitally relevant to drug design.

Dickson expects that as methods continue to develop, simulations of drug binding will become widely used both in academia and the pharmaceutical industry, contributing to the rational design of new drug therapies. These simulations also open up a whole new world of dynamic structural information that can be used to identify new binding sites that can be targeted with new molecules. According to Processor Dickson, “We are only beginning to understand the physics behind the binding of proteins to small molecules, which has roots that are billions of years old, and has a profound impact on the future of our species.”

]]> 0
NVIDIA Docker: GPU Server Application Deployment Made Easy Tue, 28 Jun 2016 02:35:23 +0000 Over the last few years there has been a dramatic rise in the use of containers for deploying data center applications at scale. The reason for this is simple: containers encapsulate an application’s dependencies to provide reproducible and reliable execution of applications and services without the overhead of a full virtual machine. If you have ever spent a day provisioning a server with a multitude of packages for a scientific or deep learning application, or have put in weeks of effort to ensure your application can be built and deployed in multiple linux environments, then Docker containers are well worth your time.

Figure 1: Docker containers encapsulate applications’ dependencies to provide reproducible and reliable execution. The NVIDIA Docker plugin enables deployment of GPU-accelerated applications across any Linux GPU server with NVIDIA Docker support.

At NVIDIA, we use containers in a variety of ways including development, testing, benchmarking, and of course in production as the mechanism for deploying deep learning frameworks through the NVIDIA DGX-1’s Cloud Managed Software. Docker has been game-changing in how we manage our workflow. Using Docker, we can develop and prototype GPU applications on a workstation, and then ship and run those applications anywhere that supports GPU containers.

In this article, we will introduce Docker containers; explain the benefits of the NVIDIA Docker plugin; walk through an example of building and deploying a simple CUDA application; and finish by demonstrating how you can use NVIDIA Docker to run today’s most popular deep learning applications and frameworks including DIGITS, Caffe, and TensorFlow.

Last week at DockerCon 2016, Felix and Jonathan presented a talk “Using Docker for GPU-Accelerated Applications”. Here are the slides.

Introduction to Docker

A Docker container is a mechanism for bundling a Linux application with all of its libraries, data files, and environment variables so that the execution environment is always the same, on whatever Linux system it runs and between instances on the same host. Docker containers are user-mode only, so all kernel calls from the container are handled by the host system kernel. On its web site, Docker describes containers this way:


It’s important to distinguish containers and hypervisor-based virtual machines (VMs). VMs allow multiple copies of the operating system, or even multiple different operating systems, to share a machine. Each VM can host and run multiple applications. By comparison, a container is designed to virtualize a single application, and all containers deployed on a host share a single OS kernel, as Figure 2 shows. Typically, containers spin up faster, run the application with bare-metal performance, and are simpler to manage since there is no additional overhead to making an OS kernel call.

Figure 2: While VMs encapsulate the entire OS and any applications, containers encapsulate individual applications and their dependencies for portable deployment, but share the same host OS between containers.

Docker provides both hardware and software encapsulation by allowing multiple containers to run on the same system at the same time each with their own set of resources (CPU, memory, etc) and their own dedicated set of dependencies (library version, environment variables, etc.). Docker also provides portable Linux deployment: Docker containers can be run on any Linux system with kernel is 3.10 or later.  All major Linux distros have supported Docker since 2014. Encapsulation and portable deployment are valuable to both the developers creating and testing applications and the operations staff who run them in data centers.

Docker provides many more important features.

  • Docker’s powerful command-line tool, `docker build`, creates Docker images from source code and binaries, using the description provided in a “Dockerfile”.
  • Docker’s component architecture allows one container image to be used as a base for other containers.
  • Docker provides automatic versioning and labeling of containers, with optimized assembly and deployment. Docker images are assembled from versioned layers so that only the layers missing on a server need to be downloaded.
  • Docker Hub is a service that makes it easy to share docker images publicly or privately.
  • Containers can be constrained to a limited set of resources on a system (e.g one CPU core and 1GB of memory).

Docker provides a layered file system that conserves disk space and forms the basis for extensible containers.

Why NVIDIA Docker?

Docker containers are platform-agnostic, but also hardware-agnostic. This presents a problem when using specialized hardware such as NVIDIA GPUs which require kernel modules and user-level libraries to operate.  As a result, Docker does not natively support NVIDIA GPUs within containers.

One of the early work-arounds to this problem was to fully install the NVIDIA drivers inside the container and map in the character devices corresponding to the NVIDIA GPUs (e.g. /dev/nvidia0) on launch. This solution is brittle because the version of the host driver must exactly match the version of the driver installed in the container. This requirement drastically reduced the portability of these early containers, undermining one of Docker’s more important features.

To enable portability in Docker images that leverage NVIDIA GPUs, we developed nvidia-docker, an open-source project hosted on Github that provides the two critical components needed for portable GPU-based containers:

  1. driver-agnostic CUDA images; and
  2. a Docker command line wrapper that mounts the user mode components of the driver and the GPUs (character devices) into the container at launch.

nvidia-docker is essentially a wrapper around the docker command that transparently provisions a container with the necessary components to execute code on the GPU. It is only absolutely necessary when using nvidia-docker run to execute a container that uses GPUs. But for simplicity in this post we use it for all Docker commands.

Installing Docker and NVIDIA Docker

Before we can get started building containerized GPU applications, let’s make sure you have the prerequisite software installed and working properly. You will need:

To help with the installation, we have created an ansible role to perform the docker and nvidia-docker install. Ansible is a tool for automating configuration management of machines and deployment of applications.

To test that you are ready to go, run the following command. You should see output similar to what is shown.

ryan@titanx:~$ nvidia-docker run --rm hello-world
Using default tag: latest
latest: Pulling from library/hello-world
a9d36faac0fe: Pull complete 
Digest: sha256:e52be8ffeeb1f374f440893189cd32f44cb166650e7ab185fa7735b7dc48d619
Status: Downloaded newer image for hello-world:latest

Hello from Docker.
This message shows that your installation appears to be working correctly.

[... simplified output ...]

Now that everything is working, let’s get started with building a simple GPU application in a container.

Build Containerized GPU Applications

To highlight the features of Docker and our plugin, I will build the deviceQuery application from the CUDA Toolkit samples in a container. This three-step method can be applied to any of the CUDA samples or to your favorite application with minor changes.

  1. Set up and explore the development environment inside a container.
  2. Build the application in the container.
  3. Deploy the container in multiple environments.

Development Environment

The Docker equivalent of installing the CUDA development libraries is the following command:

nvidia-docker pull nvidia/cuda

This command pulls the latest version of the nvidia/cuda image from Docker Hub, which is a cloud storage service for container images. Commands can be executed in this container using docker run. The following is an invocation of nvcc --version in the container we just pulled.

ryan@titanx:~$ nvidia-docker run --rm -ti nvidia/cuda nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2015 NVIDIA Corporation
Built on Tue_Aug_11_14:27:32_CDT_2015
Cuda compilation tools, release 7.5, V7.5.17

If you need CUDA 6.5 or 7.0, you can specify a tag for the image. A list of available CUDA images for Ubuntu and CentOS can be found on the nvidia-docker wiki. Here is a similar example using CUDA 7.0.

ryan@titanx:~$ nvidia-docker run --rm -ti nvidia/cuda:7.0 nvcc --version
Unable to find image 'nvidia/cuda:7.0' locally
7.0: Pulling from nvidia/cuda
6c953ac5d795: Already exists 
[ … simplified layers -- ubuntu base image … ]
68bad08eb200: Pull complete 
[ … simplified layers -- cuda 7.0 toolkit … ] 

Status: Downloaded newer image for nvidia/cuda:7.0
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2015 NVIDIA Corporation
Built on Mon_Feb_16_22:59:02_CST_2015
Cuda compilation tools, release 7.0, V7.0.27

Since only the nvidia/cuda:7.5 image is locally present on the system, the docker run command above folds the pull and run operations together. You’ll notice that the pull of the nvidia/cuda:7.0 image is faster than the pull of the 7.5 image. This is because both container images share the same base Ubuntu 14.04 image, which is already present on the host. Docker caches and reuses image layers, so only the new layers needed for the CUDA 7.0 toolkit are downloaded. This is an important point: Docker images are built layer-by-layer and the layers can be shared with multiple images to save disk space on the host (as well as deployment time). This extensibility of Docker images is a powerful feature that we will explore later in the article. (Read more about Docker’s layered filesystem).  

Finally, you can explore the development image by running a container that executes a bash shell.

nvidia-docker run --rm -ti nvidia/cuda:7.5 bash

Feel free to play around in this sandboxed environment. You can install packages via apt or examine the CUDA libraries in /usr/local/cuda, but note that any changes you make to the running container will be lost when you exit (ctrl + d or the exit command).

This non-persistence is in fact a feature of Docker.  Every instance of a Docker container starts from the same initial state defined by the image.  In the next section we will look at a mechanism for extending and adding new content to an image.

Building an Application

Let’s build the deviceQuery application in a container by extending the CUDA image with layers for the application. One mechanism used for defining layers is the Dockerfile. The Dockerfile acts as a blueprint in which each instruction in the file adds a new layer to the image. Let’s look at the Dockerfile for the deviceQuery application.

# FROM defines the base image
FROM nvidia/cuda:7.5

# RUN executes a shell command
# You can chain multiple commands together with && 
# A \ is used to split long lines to help with readability
# This particular instruction installs the source files 
# for deviceQuery by installing the CUDA samples via apt
        cuda-samples-$CUDA_PKG_VERSION && \

# set the working directory 
WORKDIR /usr/local/cuda/samples/1_Utilities/deviceQuery

RUN make

# CMD defines the default command to be run in the container 
# CMD is overridden by supplying a command + arguments to 
# `docker run`, e.g. `nvcc --version` or `bash`
CMD ./deviceQuery

For a full list of available Dockerfile instructions see the Dockerfile reference page.

Run the following docker build command in a folder with only the Dockerfile to build the application’s container image from the Dockerfile blueprint.

nvidia-docker build -t device-query .

This command generates a docker container image named device-query which inherits all the layers from the CUDA 7.5 development image.

Running a Containerized CUDA Application

We are now ready to execute the device-query container on the GPU.  By default, nvidia-docker maps all of the GPUs on the host into the container.

ryan@titanx:~$ nvidia-docker run --rm -ti device-query
./deviceQuery Starting...


Detected 2 CUDA Capable device(s)

Device 0: "GeForce GTX TITAN X"
  CUDA Driver Version / Runtime Version          8.0 / 7.5
  CUDA Capability Major/Minor version number:    5.2
  Total amount of global memory:                 12287 MBytes (12883345408 bytes)
  (24) Multiprocessors, (128) CUDA Cores/MP:     3072 CUDA Cores
[... simplified output …]

Device 1: "GeForce GTX TITAN X"
  CUDA Driver Version / Runtime Version          8.0 / 7.5
  CUDA Capability Major/Minor version number:    5.2
  Total amount of global memory:                 12288 MBytes (12884705280 bytes)
[... simplified output …]

> Peer access from GeForce GTX TITAN X (GPU0) -> GeForce GTX TITAN X (GPU1) : Yes
> Peer access from GeForce GTX TITAN X (GPU1) -> GeForce GTX TITAN X (GPU0) : Yes

CUDA Runtime Version = 7.5, NumDevs = 2, 
Device0 = GeForce GTX TITAN X, 
Device1 = GeForce GTX TITAN X
Result = PASS

nvidia-docker also provides resource isolation capabilities through the NV_GPU environment variable. The following example runs the device-query container on GPU 1.

ryan@titanx:~$ NV_GPU=1 nvidia-docker run --rm -ti device-query
./deviceQuery Starting...


Detected 1 CUDA Capable device(s)

Device 0: "GeForce GTX TITAN X"
  CUDA Driver Version / Runtime Version          8.0 / 7.5
  CUDA Capability Major/Minor version number:    5.2
  Total amount of global memory:                 12288 MBytes (12884705280 bytes)

[... simplified output …]

CUDA Runtime Version = 7.5, NumDevs = 1, 
Device0 = GeForce GTX TITAN X
Result = PASS

Because we only mapped GPU 1 into the container, the deviceQuery application can only see and report on one GPU. Resource isolation allows you to specify which GPUs your containerized application is allowed to use. For best performance, make sure to account for the topology of the PCI tree when selecting a subset of GPUs. (To better understand the importance of the PCIe topology for multi-GPU performance see this blog post about efficient collective communication with the NCCL library.)

Deploying a GPU Container

In our example so far, we built and ran device-query on a workstation with two NVIDIA Titan X GPUs. Now, let’s deploy that container on a DGX-1 server. Deploying a container is the process of moving the container images from where they are built to where they will be run. One of a variety of ways to deploy a container is Docker Hub, a cloud service used to host container images similar to how Github hosts git repositories.  The nvidia/cuda images we used to build device-query are hosted by Docker Hub.

Before we can push device-query to Docker Hub, we need to tag it with the a Docker Hub username/account.

ryan@titanx:~$ nvidia-docker tag device-query ryanolson/device-query

Now pushing the tagged images to Docker Hub is simple.

ryan@titanx:~/docker/blog-post$ nvidia-docker push ryanolson/device-query
The push refers to a repository []
35f4e4c58e9f: Pushed 
170db279f3a6: Pushed
3f6d94182a7e: Mounted from nvidia/cuda 
[ … simplifed multiple ‘from nvidia/cuda’ layers … ]

For alternatives to using Docker Hub or other container hosting repositories, check out the docker save and docker load commands.

The device-query image we pushed to Docker Hub is now available to any docker-enabled server in the world. To deploy device-query on the DGX-1 in my lab, I simply pull and run the image. This example uses the CUDA 8.0 release candidate, which you can access now by joining the NVIDIA Accelerated Computing Developer Program.

lab@dgx-1-07:~$ nvidia-docker run --rm -ti ryanolson/device-query
Using default tag: latest
latest: Pulling from ryanolson/device-query
6c953ac5d795: Already exists 
[ … simplified multiple ‘already exists’ layers … ]
1cc994928295: Pull complete 
a9e6b6393938: Pull complete 
Digest: sha256:0b42703a0785cf5243151b8cba7cc181a6c20a3c945a718ce712476541fe4d70
Status: Downloaded newer image for ryanolson/device-query:latest

./deviceQuery Starting...


Detected 8 CUDA Capable device(s)

Device 0: "Tesla P100-SXM2-16GB"
  CUDA Driver Version / Runtime Version          8.0 / 7.5
  CUDA Capability Major/Minor version number:    6.0
  Total amount of global memory:                 16281 MBytes (17071669248 bytes)
  (56) Multiprocessors, (64) CUDA Cores/MP:     3584 CUDA Cores

[... simplified output … ]

Device 7: "Tesla P100-SXM2-16GB"
  CUDA Driver Version / Runtime Version          8.0 / 7.5
  CUDA Capability Major/Minor version number:    6.0
  Total amount of global memory:                 16281 MBytes (17071669248 bytes)

[ … simplified output … ]

CUDA Runtime Version = 7.5, NumDevs = 8, 
Device0 = Tesla P100-SXM2-16GB, 
Device1 = Tesla P100-SXM2-16GB, 
Device2 = Tesla P100-SXM2-16GB, 
Device3 = Tesla P100-SXM2-16GB, 
Device4 = Tesla P100-SXM2-16GB, 
Device5 = Tesla P100-SXM2-16GB, 
Device6 = Tesla P100-SXM2-16GB, 
Device7 = Tesla P100-SXM2-16GB
Result = PASS

Docker for Deep Learning and HPC

Now that you have seen how easy it is to build and run a simple CUDA application, how hard is it to run a more complex application like DIGITS or Caffe? What if I told you it was as easy as a single command?

nvidia-docker run --name digits --rm -ti -p 8000:34448 nvidia/digits

This command launches the nvidia/digits container and maps port 34448 in the container which is running the DIGITS web service to port 8000 on the host.

If you run that command on the same machine on which you are reading this, then try opening http://localhost:8000 now. If you executed the nvidia-docker command on a remote machine, you will need to access port 8000 on that machine. Just replace localhost with the hostname or IP address of the remote server if you can access it directly.

With one command, you have NVIDIA’s DIGITS platform up and running and accessible from a browser, as Figure 3 shows.

Figure 3. This screenshot shows a local web browser connecting to DIGITS running in a container on a remote host, using ssh tunneling to forward localhost:8000 to port 8000 on the remote host running the container.
Figure 3. This screenshot shows a local web browser connecting to DIGITS running in a container on a remote host, using ssh tunneling to forward localhost:8000 to port 8000 on the remote host running the container. (Click for a larger screenshot.)

This is why we love Docker at NVIDIA. As developers of open-source software like DIGITS, we want users like you to be able to use our latest software with minimal effort.

And we’re not the only ones: Google and Microsoft use our CUDA images as the base images for TensorFlow and CNTK, respectively.  Google provides pre-built Docker images of TensorFlow through their public container repository, and Microsoft provides a Dockerfile for CNTK that you can build yourself.

Let’s see how easy it is to launch a more complex application like TensorFlow which requires Numpy, Bazel and myriad other dependencies. Yep, it’s just a single line! You don’t even need to download and build TensorFlow, you can use the image provided on Docker Hub directly.

nvidia-docker run --rm -ti tensorflow/tensorflow:r0.9-devel-gpu

After running this command, you can test TensorFlow by running its included MNIST training script:

root@ab5f46ba17d2:~# python -m tensorflow.models.image.mnist.convolutional
[tensorflow/stream_executor/] successfully opened CUDA library locally
[tensorflow/stream_executor/] successfully opened CUDA library locally
[tensorflow/stream_executor/] successfully opened CUDA library locally
[tensorflow/stream_executor/] successfully opened CUDA library locally
[tensorflow/stream_executor/] successfully opened CUDA library locally
Successfully downloaded train-images-idx3-ubyte.gz 9912422 bytes.
Successfully downloaded train-labels-idx1-ubyte.gz 28881 bytes.
Successfully downloaded t10k-images-idx3-ubyte.gz 1648877 bytes.
Successfully downloaded t10k-labels-idx1-ubyte.gz 4542 bytes.

Start using NVIDIA Docker Today

In this post we covered the basics of building a GPU application in a container by extending the nvidia/cuda images and deploying our new container on multiple different platforms. When you containerize your GPU application, get in touch with us using the comments below so we can add your project to the the list of projects using nvidia-docker.

To get started and learn more, check out the nvidia-docker github page and our 2016 Dockercon talk.

]]> 7
Production Deep Learning with NVIDIA GPU Inference Engine Mon, 20 Jun 2016 06:01:59 +0000 Figure 1. NVIDIA GPU Inference Engine (GIE) provides even higher efficiency and performance for neural network inference. Tests performed using GoogLenet. CPU-only: Single-socket Intel Xeon (Haswell) E5-2698 v3@2.3GHz with HT. GPU: NVIDIA Tesla M4 + cuDNN 5 RC. GPU + GIE: NVIDIA Tesla M4 + GIE.
GPU: NVIDIA Tesla M4 + cuDNN 5 RC.


NVIDIA GPU Inference Engine (GIE) is a high-performance deep learning inference solution for production environments. Power efficiency and speed of response are two key metrics for deployed deep learning applications, because they directly affect the user experience and the cost of the service provided. GIE automatically optimizes trained neural networks for run-time performance, delivering up to 16x higher performance per watt on a Tesla M4 GPU compared to the CPU-only systems commonly used for inference  today.

Figure 1 shows GIE inference performance per watt of the relatively complex GoogLeNet running on a Tesla M4. GIE can deliver 20 Images/s/Watt on the simpler AlexNet benchmark.

In this post, we will discuss how you can use GIE to get the best efficiency and performance out of your trained deep neural network on a GPU-based deployment platform.

Deep Learning Training and Deployment

Solving a supervised machine learning problem with deep neural networks involves a two-step process.

  1. The first step is to train a deep neural network on massive amounts of labeled data using GPUs. During this step, the neural network learns millions of weights or parameters that enable it to map input data examples to correct responses. Training requires iterative forward and backward passes through the network as the objective function is minimized with respect to the network weights. Often several models are trained and accuracy is validated against data not seen during training in order to estimate real-world performance.
  2. The next step–inference–uses the trained model to make predictions from new data. During this step, the best trained model is used in an application running in a production environment such as a data center, an automobile, or an embedded platform. For some applications, such as autonomous driving, inference is done in real time and therefore high throughput is critical.

To learn more about the differences between training and inference, see Michael Andersch’s post on inference with GPUs.

The target deployment environment introduces various challenges that are typically not present in the training environment. For example, if the target is an embedded device using the trained neural network to perceive its surroundings, then the forward inference pass through the model has a direct impact on the overall response time and the power consumed by the device.  The key metric to optimize is power efficiency: the inference performance per watt.

Performance per watt is also the critical metric in maximizing data center operational efficiency.  In this scenario, the need to minimize latency and energy used on large volumes of geographically and temporally disparate requests limits the ability to form large batches.

Introducing GPU Inference Engine

GIE is a high-performance inference engine designed to deliver maximum inference throughput and efficiency for common deep learning applications such as image classification, segmentation, and object detection. GIE optimizes your trained neural networks for run-time performance and delivers GPU-accelerated inference for web/mobile, embedded and automotive applications.

Figure 2. The GPU Inference Engine Workflow. Step 1: Optimize a trained neural network. Step 2: Perform real-time inference with the GIE runtime.
Figure 2. The GPU Inference Engine Workflow.

There are two phases in the use of GIE: build and deployment (See Figure 2). In the build phase, GIE performs optimizations on the network configuration and generates an optimized plan for computing the forward pass through the deep neural network. The plan is an optimized object code that can be serialized and stored in memory or on disk.

The deployment phase generally takes the form of a long running service or user application that accepts batches of input data, performs inference by executing the plan on the input data and returns batches of output data (classification, object detection, etc). With GIE you don’t need to install and run a deep learning framework on the deployment hardware. Discussion of the batching and pipeline of the inference service is a topic for another post; instead we will focus on how to use GIE for inference.

GIE Build Phase

The GIE runtime needs three files to deploy a classification neural network:

  1. a network architecture file (deploy.prototxt),
  2. trained weights (net.caffemodel), and
  3. a label file to provide a name for each output class.

In addition, you must define the batch size and the output layer. Code Listing 1 illustrates how to convert a Caffe model to a GIE object. The builder (lines 4-7) is responsible for reading the network information. Alternatively, you can use the builder to define the network information if you don’t provide a network architecture file (deploy.prototxt).

GIE supports the following layer types.

  • Convolution: 2D
  • Activation: ReLU, tanh and sigmoid
  • Pooling: max and average
  • ElementWise: sum, product or max of two tensors
  • LRN: cross-channel only
  • Fully-connected: with or without bias
  • SoftMax: cross-channel only
  • Deconvolution
IBuilder* builder = createInferBuilder(gLogger);

// parse the caffe model to populate the network, then set the outputs
INetworkDefinition* network = builder->createNetwork();

CaffeParser parser;
auto blob_name_to_tensor = parser.parse(“deploy.prototxt”,

// specify which tensors are outputs

// Build the engine
builder->setMaxWorkspaceSize(1 << 30); 
ICudaEngine* engine = builder->buildCudaEngine(*network);

You can also use the GIE C++ API to define the network without the Caffe parser, as Listing 2 shows. You can use the API to define any supported layer and its parameters. You can define any parameter that varies between networks, including convolution layer weight dimensions and outputs as well as the window size and stride for pooling layers.

ITensor* in = network->addInput(“input”, DataType::kFloat, Dims4{…});
IPoolingLayer* pool = network->addPooling(in, PoolingType::kMAX, …);

After defining or loading the network, you must specify the output tensors as line 13 of Listing 1 shows; in our example the output is “prob” (for probability). Next, define the batch size (line 16), which can vary depending on the deployment scenario. Listing 1 uses a batch size of 1 but you may choose larger batch sizes to fit your application needs and system configuration. Underneath, GIE performs layer optimizations to reduce inference time. While this is transparent to the API user, analyzing the network layers requires memory, so you must specify the maximum workspace size (line 17).

The last step is to call buildCudaEngine to perform layer optimization and build the engine with the optimized network based on your provided inputs and parameters. Once the model is converted to a GIE object it is deployable and can either be used on the host device or saved and used elsewhere.

GIE performs several important transformations and optimizations to the neural network graph. First, layers with unused output are eliminated to avoid unnecessary computation. Next, where possible convolution, bias, and ReLU layers are fused to form a single layer. Figure 4 shows the result of this vertical layer fusion on the original network from Figure 3 (fused layers are labeled CBR in Figure 4). Layer fusion improves the efficiency of running GIE-optimized networks on the GPU.

Figure 3. An example convolutional neural network with multiple convolutional and activation layers.
Figure 4. An example of vertical layer fusion on a convolutional neural network. Here, convolutional layers are combined with subsequent bias and activation (ReLU) layers.

Another transformation is horizontal layer fusion, or layer aggregation, along with the required division of aggregated layers to their respective outputs, as Figure 5 shows. Horizontal layer fusion improves performance by combining layers that take the same source tensor and apply the same operations with similar parameters, resulting in a single larger layer for higher computational efficiency. The example in Figure 5 shows the combination of 3 1×1 CBR layers from Figure 4 that take the same input into a single larger 1×1 CBR layer. Note that the output of this layer must be disaggregated to feed into the different subsequent layers from the original input graph.

Figure 5. An example of horizontal layer fusion on a convolutional neural network. Here, multiple 1x1 CBR layers from Figure 4 are fused “horizontally”, or across similar layers in the graph that share the same input.
Figure 5. An example of horizontal layer fusion on a convolutional neural network. Here, multiple 1×1 CBR layers from Figure 4 are fused “horizontally”, or across similar layers in the graph that share the same input.

GIE performs its transformations during the build phase transparently to the API user after the GIE parser reads in the trained network and configuration file, as Listing 1 shows.

GIE Deploy Phase

The inference builder (IBuilder) buildCudaEngine method returns a pointer to a new inference engine runtime object (ICudaEngine).  This runtime object is ready for immediate use; alternatively, its state can be serialized and saved to disk or to an object store for distribution. The serialized object code is called the Plan.

As mentioned earlier, the full scope of batching and streaming data to and from the runtime inference engine is beyond the scope of this article.  Listing 3 demonstrates the steps required to use the inference engine to process a batch of input data to generate a result.

// The execution context is responsible for launching the 
// compute kernels
IExecutionContext *context = engine->createExecutionContext();

// In order to bind the buffers, we need to know the names of the 
// input and output tensors.
int inputIndex = engine->getBindingIndex(INPUT_LAYER_NAME),
int outputIndex = engine->getBindingIndex(OUTPUT_LAYER_NAME);

// Allocate GPU memory for Input / Output data
void* buffers = malloc(engine->getNbBindings() * sizeof(void*));
cudaMalloc(&buffers[inputIndex], batchSize * size_of_single_input);
cudaMalloc(&buffers[outputIndex], batchSize * size_of_single_output);

// Use CUDA streams to manage the concurrency of copying and executing
cudaStream_t stream;

// Copy Input Data to the GPU
cudaMemcpyAsync(buffers[inputIndex], input, 
                batchSize * size_of_single_input, 
                cudaMemcpyHostToDevice, stream);

// Launch an instance of the GIE compute kernel
context.enqueue(batchSize, buffers, stream, nullptr);

// Copy Output Data to the Host
cudaMemcpyAsync(output, buffers[outputIndex], 
                batchSize * size_of_single_output, 
                cudaMemcpyDeviceToHost, stream));

// It is possible to have multiple instances of the code above
// in flight on the GPU in different streams.
// The host can then sync on a given stream and use the results

GIE Performance

At the end of the day, the success of GIE comes down to the performance it provides for inference. To measure the performance benefits we compared the per-layer timings of the GoogLeNet network using Caffe and GIE on NVIDIA Tesla M4 GPUs with a batch size of 1 averaged over 1000 iterations with GPU clocks fixed in the P0 state.

Figure 6. GIE + GPU vs. Caffe + GPU GoogLeNet layer execution time (lower is better).

The bar graph (Figure 6) is sorted to show the 10 most computationally expensive GoogLeNet layers (as run by Caffe) ordered from left to right as light green bars. The dark green bars represent the same layers run using GIE (lower is better). Since GIE can combine layers both vertically and horizontally into a single optimized kernel, the Caffe timing shown for each bar is the sum of Caffe kernels corresponding to each fused GIE kernel, while the GIE timing for each bar is for a single fused and optimized kernel.

Bars with two layers correspond to vertically fused layers, namely CBR (convolution + bias + activation/relu). Bars with four layers (or whose names contain “||” separating two CBR layer names) correspond to two CBRs that are horizontally fused, meaning two CBRs that share the same input tensor and thus gain the advantage of cache reuse on a single-pass of the input tensor vs. two separate kernel launches with the same input tensor. Unsurprisingly, the GIE kernels with four fused layers show some of the largest relative speedups and contribute to ~30% of the overall speedup.  The remainder of the speedup predominately comes from the two vertically fused CBR layers, which on average have a lower relative speedup, but comprise the bulk of the computation.

Maximize Performance and Efficiency with GIE

The NVIDIA GPU Inference Engine enables you to easily deploy neural networks to add deep learning based capabilities to your products with the highest performance and efficiency. GIE supports networks trained using popular neural network frameworks including Caffe, Theano, Torch and Tensorflow. During the build phase GIE identifies opportunities to optimize the network, and in the deployment phase GIE runs the optimized network in a way that minimizes latency and maximizes throughput.

If you are running web or mobile applications that are backed by data center servers, GIE’s low overhead means that you can deploy more varied and complex models to add intelligence to your product that will delight your users. If you are using deep learning to create the next generation of smart devices, GIE helps you deploy networks with high performance, high accuracy, and high energy efficiency.

Moreover, GIE enables you to leverage the power of GPUs to perform neural network inference using mixed-precision FP16 data. Performing neural network inference using FP16 can reduce memory usage by half and provide higher performance on Tesla P100 and Jetson TX1 GPUs.

GIE is currently being evaluated under an Early Access (EA) Program. To be notified when GIE is ready for public release or if you are interested in participating in the EA program, please visit the GIE product page to contact us today. To learn more about neural network inference on GPUs, see Michael Andersch’s post on inference with GPUs.

]]> 0
Inside Pascal: NVIDIA’s Newest Computing Platform Mon, 20 Jun 2016 06:00:44 +0000 At the 2016 GPU Technology Conference in San Jose, NVIDIA CEO Jen-Hsun Huang announced the new NVIDIA Tesla P100, the most advanced accelerator ever built. Based on the new NVIDIA Pascal GP100 GPU and powered by ground-breaking technologies, Tesla P100 delivers the highest absolute performance for HPC, technical computing, deep learning, and many computationally intensive datacenter workloads. Today at ISC 2016 in Frankfurt, NVIDIA expanded the Tesla P100 product family with the new Tesla P100 for PCIe-based servers.




Tesla P100: Extreme Performance and Features for GPU Computing


  • Unified Memory and Compute Preemption—significantly improved programming model;

Tesla P100 Available for NVLink and PCIe Servers

Tesla P100 PCIe GPU Accelerator Board
Tesla P100 PCIe GPU Accelerator Board

Tesla P100 accelerators will be available in two forms: A traditional GPU accelerator board for PCIe-based servers, and an SXM2 module for NVLink-optimized servers. P100 for PCIe-based servers allows HPC data centers to deploy the most advanced GPUs within PCIe-based nodes to support a mix of CPU and GPU workloads. P100 for NVLink-optimized servers provides the best performance and strong-scaling  for hyperscale and HPC data centers running applications that scale to multiple GPUs, such as deep learning. The table below provides the complete specifications of both Tesla P100 accelerators.

The Pascal GP100 Architecture: Faster in Every Way



Pascal GP100 Block Diagram
Pascal GP100 Block Diagram

Delivering higher performance and improving energy efficiency are two key goals for new GPU architectures. A number of changes to the SM in the Maxwell architecture improved its efficiency compared to Kepler. Pascal builds on this and incorporates additional improvements that increase performance per watt even further over Maxwell. While TSMC’s 16nm Fin-FET manufacturing process plays an important role, many GPU architectural modifications were also implemented to further reduce power consumption while maintaining high performance.

The following table provides a high-level comparison of Tesla P100 specifications compared to previous-generation Tesla GPU accelerators.

Tesla Products Tesla K40 Tesla M40 Tesla P100 (NVLink) Tesla P100 (PCIe)
GPU / Form Factor Kepler
GK110 / PCIe
GM200 / PCIe
GP100 / SXM2
GP100 / PCIe
SMs 15 24 56 56
TPCs 15 24 28 28
FP32 CUDA Cores / SM 192 128 64 64
FP32 CUDA Cores / GPU 2880 3072 3584 3584
FP64 CUDA Cores / SM 64 4 32 32
FP64 CUDA Cores / GPU 960 96 1792 1792
Base Clock 745 MHz 948 MHz 1328 MHz 1126 MHz
GPU Boost Clock 810/875 MHz  1114 MHz 1480 MHz 1303 MHz
FP32 GFLOPs[1] 5040 6844 10608 9340
FP64 GFLOPs[1] 1680 213 5304 4670
Texture Units 240 192 224 224
Memory Interface 384-bit GDDR5 384-bit GDDR5 4096-bit HBM2 3072-bit HBM2 (12GB)
4096-bit HBM2 (16GB)
Memory Bandwidth 288 GB/s 288 GB/s 720 GB/s 540 GB/s (12GB)
720 GB/s (16GB)
Memory Size Up to 12 GB Up to 24 GB 16 GB 12 GB or 16 GB
L2 Cache Size 1536 KB 3072 KB 4096 KB 4096 KB
Register File Size / GPU 3840 KB 6144 KB 14336 KB 14336 KB
TDP 235 Watts 250 Watts 300 Watts 250 Watts
Transistors 7.1 billion 8 billion 15.3 billion 15.3 billion
GPU Die Size 551 mm² 601 mm² 610 mm² 610 mm²
Manufacturing Process 28-nm 28-nm 16-nm 16-nm

[1] The GFLOPS in this chart are based on GPU Boost Clocks.

Pascal Streaming Multiprocessor

The Pascal GP100 Streaming Multiprocessor
The Pascal GP100 Streaming Multiprocessor

GP100’s SM incorporates 64 single-precision (FP32) CUDA Cores. In contrast, the Maxwell and Kepler SMs had 128 and 192 FP32 CUDA Cores, respectively. The GP100 SM is partitioned into two processing blocks, each having 32 single-precision CUDA Cores, an instruction buffer, a warp scheduler, and two dispatch units. While a GP100 SM has half the total number of CUDA Cores of a Maxwell SM, it maintains the same register file size and supports similar occupancy of warps and thread blocks.

GP100’s SM has the same number of registers as Maxwell GM200 and Kepler GK110 SMs, but the entire GP100 GPU has far more SMs, and thus many more registers overall. This means threads across the GPU have access to more registers, and GP100 supports more threads, warps, and thread blocks in flight compared to prior GPU generations.

Overall shared memory across the GP100 GPU is also increased due to the increased SM count, and aggregate shared memory bandwidth is effectively more than doubled. A higher ratio of shared memory, registers, and warps per SM in GP100 allows the SM to more efficiently execute code. There are more warps for the instruction scheduler to choose from, more loads to initiate, and more per-thread bandwidth to shared memory (per thread).

Compared to Kepler, Pascal’s SM features a simpler datapath organization that requires less die area and less power to manage data transfers within the SM. Pascal also provides superior scheduling and overlapped load/store instructions to increase floating point utilization. The new SM scheduler architecture in GP100 improves upon the advances of the Maxwell scheduler and is even more intelligent, providing increased performance and reduced power consumption. Each warp scheduler (one per processing block) is capable of dispatching two warp instructions per clock.

High-Performance Double Precision

Because of the importance of high-precision computation for technical computing and HPC codes, a key design goal for Tesla P100 is high double-precision performance. Each GP100 SM has 32 FP64 units, providing a 2:1 ratio of single- to double-precision throughput. Compared to the 3:1 ratio in Kepler GK110 GPUs, this allows Tesla P100 to process FP64 workloads more efficiently.

Like previous GPU architectures, GP100 supports full IEEE 754‐2008 compliant single- and double‐precision arithmetic, including support for the fused multiply‐add (FMA) operation and full speed support for denormalized values.

FP16 Arithmetic Support for Faster Deep Learning

The rapidly growing field of deep learning is one of the fastest growing fields of computing. Deep learning has proven to provide a high-level of accuracy and adaptability in applications spanning automatic image captioning, autonomous driving object recognition, natural language understanding and translation, and even computer-generated art. For an in-depth introduction to deep learning, check out the Deep Learning in a Nutshell series here on Parallel Forall.

Unlike other technical computing applications that require high-precision floating-point computation, deep neural network architectures have a natural resilience to errors due to the backpropagation algorithm used in their training. Storing FP16 data compared to higher precision FP32 or FP64 reduces memory usage of the neural network, allowing training and deployment of larger networks. Using FP16 computation improves performance up to 2x compared to FP32 arithmetic, and similarly FP16 data transfers take less time than FP32 or FP64 transfers.


Improved Atomics

Atomic memory operations are important in parallel programming, allowing concurrent threads to correctly perform read-modify-write operations on shared data. Kepler significantly increased the throughput of atomic operations to global memory compared to the earlier Fermi architecture; however, both Fermi and Kepler implemented shared memory atomics using an expensive lock/update/unlock pattern.

Maxwell improved this by implementing native hardware support for shared memory atomic operations for 32-bit integers, and native shared memory 32-bit and 64-bit compare-and-swap (CAS), which can be used to implement other atomic functions with reduced overhead (compared to the Fermi and Kepler methods which were implemented in software).

GP100 further improves atomics by providing an FP64 atomic add instruction for values in global memory. The `atomicAdd()“ function in CUDA now applies to 32 and 64-bit integer and floating-point data. Previously, FP64 atomic addition had to be implemented using a compare-and-swap loop, which is generally slower than a native instruction.

Compute Capability 6.0

GP100 supports the new Compute Capability 6.0. The following table compares parameters of different Compute Capabilities for NVIDIA GPU architectures.

GPU Kepler GK110 Maxwell GM200 Pascal GP100
Compute Capability 3.5 5.3 6.0
Threads / Warp 32 32 32
Max Warps / Multiprocessor 64 64 64
Max Threads / Multiprocessor 2048 2048 2048
Max Thread Blocks / Multiprocessor 16 32 32
Max 32-bit Registers / SM 65536 65536 65536
Max Registers / Block 65536 32768 65536
Max Registers / Thread 255 255 255
Max Thread Block Size 1024 1024 1024
CUDA Cores / SM 192 128 64
Shared Memory Size / SM Configurations (bytes) 16K/32K/48K 96K 64K

Dramatic Improvements in Memory

Tesla P100 Accelerator (Front)
Tesla P100 Accelerator (Front)

Many applications today are bottlenecked by memory bandwidth, especially in high-performance computing. Today, developers of high-performance software across all types of processors devote a lot of effort to optimizing code for efficient memory accesses, and to keep data in the parts of the memory hierarchy closest to the computational units. Some applications—for example in deep learning where many-layered neural networks are trained using massive data sets—are limited more by memory capacity. So memory poses two challenges to computing performance: bandwidth and capacity.

Tesla P100 Accelerator (Back)
Tesla P100 Accelerator (Back)

Tesla P100 tackles both memory challenges using stacked memory, a technology which enables multiple layers of DRAM components to be integrated vertically on the package along with the GPU. Tesla P100 is the first GPU accelerator to use High Bandwidth Memory 2 (HBM2). HBM2 memory provides much greater bandwidth, more than twice the capacity, and higher energy efficiency, compared to current off-package GDDR5.

Rather than requiring numerous discrete memory chips surrounding the GPU as in traditional GDDR5 GPU board designs, HBM2 includes one or more vertical stacks of multiple memory dies. The memory dies are linked using tiny wires that are called through-silicon vias and microbumps. One 8 Gb HBM2 die contains over 5,000 through-silicon via holes. A passive silicon interposer is then used to connect the memory stacks and the GPU die. The combination of HBM2 stacks, GPU die, and silicon interposer are packaged in a single 55mm x 55mm BGA package.

Tesla P100 accelerators have four 4-die HBM2 stacks, for a total of 16 GB of memory, and 720 GB/s peak bandwidth, which is 3 times higher than the Tesla M40 memory bandwidth.

ECC Memory

Another HBM2 benefit is native support for error correcting code (ECC) funtionality, which provides higher reliability for technical computing applications that are sensitive to data corruption, such as in large-scale clusters and supercomputers, where GPUs process large datasets with long application run times.

ECC technology detects and corrects single-bit soft errors before they affect the system. In comparison, GDDR5 does not provide internal ECC protection of the contents of memory and is limited to error detection of the GDDR5 bus only: Errors in the memory controller or the DRAM itself are not detected.

GK110 Kepler GPUs offered ECC protection for GDDR5 by allocating some of the available memory for explicit ECC storage. 6.25% of the overall GDDR5 is reserved for ECC bits. In the case of a 12 GB Tesla K40 (for example), 750 MB of its total memory is reserved for ECC operation, resulting in 11.25 GB (out of 12 GB) of available memory with ECC turned on for Tesla K40. Also, accessing ECC bits causes a small decrease in memory bandwidth compared to the non-ECC case. Since HBM2 supports ECC natively, Tesla P100 does not suffer from the capacity overhead, and ECC can be active at all times without a bandwidth penalty. Like the GK110 GPU, the GP100 GPU’s register files, shared memories, L1 cache, L2 cache, and the Tesla P100 accelerator’s HBM2 DRAM are protected by a Single‐Error Correct Double‐Error Detect (SECDED) ECC code.

NVLink High Speed Interconnect

NVLink is NVIDIA’s new high-speed interconnect technology for GPU-accelerated computing. Supported on SXM-2 based Tesla P100 accelerator boards, NVLink significantly increases performance for both GPU-to-GPU communications, and for GPU access to system memory.

Today, multiple GPUs are common in workstations as well as the nodes of HPC computing clusters and deep learning training systems. A powerful interconnect is extremely valuable in multiprocessing systems. Our vision for NVLink was to create an interconnect for GPUs that would offer much higher bandwidth than PCI Express Gen 3 (PCIe), and be compatible with the GPU ISA to support shared memory multiprocessing workloads.

Support for the GPU ISA means that programs running on NVLink-connected GPUs can execute directly on data in the memory of another GPU as well as on local memory. GPUs can also perform atomic memory operations on remote GPU memory addresses, enabling much tighter data sharing and improved application scaling.

NVLink uses NVIDIA’s new High-Speed Signaling interconnect (NVHS). NVHS transmits data over a differential pair running at up to 20 Gb/sec. Eight of these differential connections form a “Sub-Link” that sends data in one direction, and two sub-links—one for each direction—form a “Link” that connects two processors (GPU-to-GPU or GPU-to-CPU). A single Link supports up to 40 GB/sec of bidirectional bandwidth between the endpoints. Multiple Links can be combined to form “Gangs” for even higher-bandwidth connectivity between processors. The NVLink implementation in Tesla P100 supports up to four Links, allowing for a gang with an aggregate maximum theoretical bandwidth of 160 GB/sec bidirectional bandwidth.

The figure below shows an 8-GPU Hybrid Cube Mesh that includes two fully NVLink-connected quads of GPUs, with NVLink connections between the quads, and GPUs within each quad connected to their respective CPUs directly through PCIe.

Eight GPU hybrid cube mesh architecture with NVLink.
Eight GPU hybrid cube mesh architecture with NVLink.

While NVLink primarily focuses on connecting multiple NVIDIA Pascal GP100 GPUs together it can also connect Pascal GP100 GPUs with IBM Power CPUs with NVLink support. The following figure highlights an example of a four-GPU system with dual NVLink-capable CPUs connected with NVLink. In this configuration, each GPU has 120 combined GB/s bidirectional bandwidth to the other 3 GPUs in the system, and 40 GB/s bidirectional bandwidth to a CPU.

A 4-GPU system with NVLink-capable CPUs connected to the GPUs via NVLink.

Increasing Developer Productivity with Unified Memory

Unified Memory is an important feature of the CUDA programming model that greatly simplifies programming and porting of applications to GPUs by providing a single, unified virtual address space for accessing all CPU and GPU memory in the system. Pascal GP100 features provide a significant advancement for GPU computing by expanding the capabilities and improving the performance of Unified Memory.

CUDA 6 Unified Memory with a Kepler GPU.
CUDA 6 Unified Memory with a Kepler GPU.

CUDA 6 introduced Unified Memory, which creates a pool of managed memory that is shared between the CPU and GPU, bridging the CPU-GPU divide. Managed memory is accessible to both the CPU and GPU using a single pointer. The CUDA system software automatically migrates data allocated in Unified Memory between GPU and CPU, so that it looks like CPU memory to code running on the CPU, and like GPU memory to code running on the GPU. For details of how Unified Memory in CUDA 6 and later simplifies porting code to the GPU, see the post “Unified Memory in CUDA 6”.

CUDA 6 Unified Memory was limited by the features of the Kepler and Maxwell GPU architectures: all managed memory touched by the CPU had to be synchronized with the GPU before any kernel launch; the CPU and GPU could not simultaneously access a managed memory allocation; and the Unified Memory address space was limited to the size of the GPU physical memory.CUDA_6_unified_memory


Pascal GP100 Unified Memory

Expanding on the benefits of CUDA 6 Unified Memory, Pascal GP100 adds features to further simplify programming and sharing of memory between CPU and GPU, and allowing easier porting of CPU parallel compute applications to use GPUs for tremendous speedups. Two main hardware features enable these improvements: support for large address spaces and page faulting capability.

Pascal GP100 Unified Memory is not limited by the physical size of GPU memory.

GP100 extends GPU addressing capabilities to enable 49-bit virtual addressing. This is large enough to cover the 48-bit virtual address spaces of modern CPUs, as well as the GPU’s own memory. Therefore, GP100 Unified Memory allows programs to access the full address spaces of all CPUs and GPUs in the system as a single virtual address space, unlimited by the physical memory size of any one processor.

Memory page faulting support in GP100 is a crucial new feature that provides more seamless Unified Memory functionality. Combined with the system-wide virtual address space, page faulting provides several benefits. First, page faulting means that the CUDA system software doesn’t need to synchronize all managed memory allocations to the GPU before each kernel launch. If a kernel running on the GPU accesses a page that is not resident in its memory, it faults, allowing the page to be automatically migrated to the GPU memory on-demand. Alternatively, the page may be mapped into the GPU address space for access over the PCIe or NVLink interconnects (mapping on access can sometimes be faster than migration). Note that Unified Memory is system-wide: GPUs (and CPUs) can fault on and migrate memory pages either from CPU memory or from the memory of other GPUs in the system.

With the new page fault mechanism, global data coherency is guaranteed with Unified Memory. This means that with GP100, the CPUs and GPUs can access Unified Memory allocations simultaneously. This was illegal on Kepler and Maxwell GPUs, because coherence could not be guaranteed if the CPU accessed a Unified Memory allocation while a GPU kernel was active. Note, as with any parallel application, developers need to ensure correct synchronization to avoid data hazards between processors.

Finally, on supporting platforms, memory allocated with the default OS allocator (e.g. ‘malloc’ or ‘new’) can be accessed from both GPU code and CPU code using the same pointer (see code example below). On these systems, Unified Memory is the default: there is no need to use a special allocator or for the creation of a special managed memory pool. Moreover, GP100’s large virtual address space and page faulting capability enable applications to access the entire system virtual memory. This means that applications can oversubscribe the memory system: in other words they can allocate, access, and share arrays larger than the total physical capacity of the system, enabling out-of-core processing of very large datasets.

With operating system support, Pascal is capable of supporting unified memory with the default system allocator. Here, malloc is all that is needed to allocate memory accessible from any CPU or GPU in the system.

Certain operating system modifications are required to enable Unified Memory with the system allocator. NVIDIA is collaborating with Red Hat and working within the Linux community to enable this powerful functionality.

NVIDIA® DGX-1™ Deep Learning Supercomputer


Data scientists and artificial intelligence (AI) researchers require accuracy, simplicity, and speed from their Deep Learning systems. Faster training and iteration ultimately means faster innovation and faster time to market. The NVIDIA DGX-1 is the world’s first purpose-built server for Deep Learning, with fully integrated hardware and software that can be deployed quickly and easily. Its revolutionary performance of up to 170 FP16 TFLOP/s significantly accelerates training time, making the NVIDIA DGX-1 the first “AI supercomputer in a box”.

The NVIDIA DGX-1 server is the first server using Tesla P100 accelerators interconnected with NVLink. Available in an eight (8) Tesla P100 accelerator configuration, the DGX-1 system is built with high performance/high reliability components in a 3U rack-mountable chassis for standalone use or cluster integration.

The 8-GPU configuration features two NVLink fully-connected P100 GPU quads that are tied together by four additional NVLinks in a Hybrid Cube Mesh topology (See the 8-GPU NVLink diagram above). Every GPU in a quad is also directly connected via PCIe to a PCIe switch that connects to a CPU. The bottom line is that an NVIDIA DGX-1 server with eight Tesla P100 accelerators can deliver over 12x the Deep Learning performance compared to previous GPU-accelerated solutions.

Combining powerful hardware with software tailored to Deep Learning, the NVIDIA DGX-1 enables developers and researchers with a turnkey solution for high-performance GPU-accelerated Deep Learning application development, testing, and network training.

Pascal GP100 vs Xeon - Alexnet Training Time.
Pascal GP100 vs Xeon – Alexnet Training Time.

Tesla P100: Built for HPC and Deep Learning

As you can see, the new NVIDIA Tesla P100 accelerator is a performance powerhouse with revolutionary new features for technical computing and deep learning. Faster in every way than its predecessors, Tesla P100 provides massive leaps in computational throughput, memory bandwidth and capacity, interconnect performance, and programmability.

In this blog post I’ve really only scratched the surface—there are a many more new features in the Pascal GP100 architecture and in Tesla P100, including new instructions, powerful features such as compute preemption, and more that I couldn’t fit into this post. To learn more, read the Tesla P100 white paper and check out our talk “Inside Pascal” from GTC 2016.

Note, this post was first published Apr 5, 2016, and updated June 19, 2016 with information on Tesla P100 for PCIe.

]]> 45


Figure 1: Applications of Graph Partitioning
Figure 1: Applications of Graph Partitioning


Also, graph partitioning and clustering aims to find a splitting of a graph into subgraphs based on a specific metric. In particular, spectral graph partitioning and clustering relies on the spectrum—the eigenvalues and associated eigenvectors—of the Laplacian matrix corresponding to a given graph. Next, I will formally define this problem, show how it is related to the spectrum of the Laplacian matrix, and investigate its properties and tradeoffs.


Let a graph G=(V,E) be defined by its vertex set V and edge set E. The vertex set V = \{1,...,n\} represents n nodes in a graph, with each node identified by a unique integer number i \in V. The edge set E = \{(i_{1},j_{1}),...,(i_{e},j_{e})\} represents e edges in a graph, with each edge from node i to j identified by a pair (i,j) \in E.

Figure 2: Minimum balanced cut of a graph G=(V,E) induced by a set of vertices S={1,2}

Applications often need to find a splitting of the graph into subgraphs of similar size connected by as few edges as possible. This statement is often formulated as the problem of finding a set of vertices S \in V that induces a minimum balanced cut of a graph G=(V,E) in the sense of a cost function

\rho(S) = \min_{S} \frac{|\partial(S)|}{|S||\bar{S}|} ,

where \bar{S} denotes complement of set S with respect to V, |.| denotes the cardinality (number of elements) of a set, and \partial(S)= \{(i,j): i \in S \land j \notin S\} denotes the boundary of a set.

For example, a minimum balanced cut of a graph G=(V,E), with V=\{1,2,3,4,5\} and E=\{(1,2),(2,4),(3,4),(3,5),(4,5)\} induced by a set of vertices S=\{1,2\} with \bar{S}=\{3,4,5\} and boundary \partial(S)=\{(2,4)\} is shown in Figure 2.

It is important to point out that both partitioning and clustering aim to split the original graph into multiple sub-graphs. However, in partitioning the number of partitions and often their size is fixed, while in clustering the fact that there are no partitions can be a result in itself [1]. Also, the optimality of the splitting can be measured by different cost functions, including modularity, betweenness centrality, or flow.

Cost Function

I will focus on the popular ratio and normalized cut cost functions, which are variations of the minimum balanced cut of a graph G=(V,E). The ratio and normalized cut cost functions are defined as

\tilde{\rho}(S_{1},...,S_{p}) = \min_{S_{1},...,S_{p}} \sum_{k=1}^{p} \frac{|\partial(S_{k})|}{|S_{k}|}


\tilde{\eta}(S_{1},...,S_{p}) = \min_{S_{1},...,S_{p}} \sum_{k=1}^{p} \frac{|\partial(S_{k})|}{\text{vol}(S_{k})} ,

respectively, where set S_{k} \subseteq V for k=1,...,p, the intersection of sets S_{k} \cap S_{j} = \emptyset for k \neq j, and \text{vol}(S)=\sum_{i \in S} \text{d}(i) denotes the volume of a set, where d(i) is the degree (number of edges) of node i.

These cost functions are simpler than they look. Notice that the numerator refers to the number of edges cut between partitions and the denominator is related to the number of elements assigned to a particular partition. For example, in distributed sparse matrix-vector multiplication, the numerator relates to the number of elements that must be sent between partitions and the denominator to the work done per partition, measured in terms of the number of rows in the former or number of non-zero elements in the latter equation.

Also, it turns out that it is possible to express these cost functions as

\tilde{\rho}(S_{1},...,S_{p}) = \min_{S_{1},...,S_{p}} \sum_{k=1}^{p} \textbf{u}_{k}^{T} L \textbf{u}_{k} = \min_{U^{T}U=I} \text{Tr} (U^{T} L U)



where tall matrix U=[\textbf{u}_{1},...,\textbf{u}_{p}], I \in \mathbb{R}^{p \times p} is the identity matrix and \text{Tr}(.) is the trace of a matrix (the sum of its diagonal elements). Here vectors \textbf{u}_{k} take only discrete values, with non-zeroes corresponding to indices in the set S_{k}, while L is the Laplacian matrix, which will be defined next.

Laplacian Matrix

The Laplacian matrix is defined as L=D-A, where A \in \mathbb{R}^{n \times n} is the adjacency matrix of the graph G=(V,E) and the diagonal matrix D = \text{diag}(A\textbf{e}), where vector \textbf{e}=(1,...,1)^{T}.

For example, the Laplacian matrix for the graph G=(V,E) shown ing Figure 2 can be written

\begin{bmatrix} 1 & -1 & & \\ -1 & 2 & & -1 & \\ & & 2 & -1 & -1 \\ & -1 & -1 & 3 & -1 \\ & & -1 & -1 & 2 \end{bmatrix} = \begin{bmatrix} 1 & & & & \\ & 2 & & & \\ & & 2 & & \\ & & & 3 & \\ & & & & 2 \end{bmatrix} - \begin{bmatrix} & 1 & & & \\ 1 & & & 1 & \\ & & & 1 & 1 \\ & 1 & 1 & & 1 \\ & & 1 & 1 & \end{bmatrix}

The Laplacian matrix has very interesting properties. To illustrate it, let the set S=\{1,2\} as shown in Figure 2 and vector \textbf{x} be such that it has 1 at position i \in S and zero otherwise. Then, I can express the cardinality |S|, volume \text{vol}(S) as well as the cardinality of the boundary |\partial(S)|  of the set S using the vector \textbf{x} and the Laplacian matrix L in the following way.

|S|=\textbf{x}^{T}\textbf{x}=  \begin{bmatrix} 1 & 1 & & & &\end{bmatrix}\begin{bmatrix} 1 \\ 1 \\  \\  \\ \\  \end{bmatrix} = 2

\text{vol}(S) =\textbf{x}^{T}D\textbf{x} = \begin{bmatrix} 1 & 1 & & & &\end{bmatrix}\begin{bmatrix} 1 & & & & \\ & 2 & & & \\ & & 2 & & \\ & & & 3 & \\ & & & & 2 \end{bmatrix}\begin{bmatrix} 1 \\ 1 \\  \\  \\ \\  \end{bmatrix} = 3

|\partial(S)|=\textbf{x}^{T}L\textbf{x} = \begin{bmatrix} 1 & 1 & & & &\end{bmatrix}\begin{bmatrix} 1 & -1 & & \\ -1 & 2 & & -1 & \\ & & 2 & -1 & -1 \\ & -1 & -1 & 3 & -1 \\ & & -1 & -1 & 2 \end{bmatrix}\begin{bmatrix} 1 \\ 1 \\  \\  \\ \\  \end{bmatrix} = 1

This illustrates why in the previous section I could express all the terms in the ratio and normalized cost functions in terms of vector \textbf{x} and the Laplacian matrix L. A more detailed explanation is given in our technical report [2].

Key Idea of Spectral Scheme

Notice that obtaining the minimum of the cost function by finding the best non-zero discrete values for the vector \textbf{u} is no easier than finding the best indices for the set S. The two formulations of the cost functions are equivalent and both are NP-hard problems.

The key idea of spectral partitioning and clustering is not to look for the discrete solution directly, but instead to do it in two steps.

First, relax the discrete constraints and let vector \textbf{u} take real instead of discrete values. In this case, following the linear algebra Courant-Fischer theorem (sometimes referred to as the Min-Max theorem), the minimum of the cost function is obtained by the eigenvectors associated with the p smallest eigenvalues of the Laplacian matrix.

Second, map the obtained real values back to the discrete ones to find the solution of interest. This step can be done using simple heuristics, such as sorting and looking for a gap in the real values, or using more advanced multi-dimensional algorithms, such as the k-means algorithm. In the former case all real values in between gaps, while in the latter case all real values clustered around a particular centroid, will be assigned the same discrete value and therefore will belong to the same particular partition or cluster.

There is no guarantee that the two-step approach will find the best solution, but in practice it often finds a good enough approximation and works reasonably well.

Figure 3 provides a visual outline of the process and Algorithm 1 presents the algorithm in pseudocode.

Figure 3: Spectral Partitioning Outline
Figure 3: Spectral Partitioning Outline
Let G=(V,E) be an input graph
Let A be the adjacency matrix of G
Let diagonal matrix D = diag(Ae)
Set the Laplacian matrix L=D-A
Solve the eigenvalue problem L u = λu
Use heuristics to transform real into discrete values

Eigenvalue Problem

The solution of the eigenvalue problem is often the most time consuming part of spectral graph partitioning /clustering. There are many eigenvalue solvers that can be used to solve it, including Lanczos, Tracemin, Jacobi-Davidson and LOBPCG. In particular, Figures 3 and 4 show experimental results comparing the performance and quality, respectively, of the Lanczos and LOBPCG methods when looking for the 30 smallest eigenvectors of a few matrices from the DIMACS graph collection. Although Lanczos is often the fastest eigenvalue solver, when incomplete-LU with 0 fill-in (ILU0) is available, the preconditioned LOBPCG eigenvalue solver can be competitive and often computes a superior-quality solution.

Figure 3: Performance of Different Eigenvalue Solvers on the GPU
Figure 3: Performance of Various Eigenvalue Solvers on the GPU.
Figure 4: Quality of Various Eigenvalue Solvers on the GPU.
Figure 4: Quality of Various Eigenvalue Solvers on the GPU.


Now I’ll compare the spectral scheme on the GPU with the spectral scheme implemented on the CPU in the CHACO software package. The experiments are performed on a workstation with a 3.2 GHz Intel Core i7-3930K CPU and an NVIDIA Tesla K40c GPU.

The schemes are very similar, but not identical because CHACO has a slightly different implementation of the algorithms, and also attempts to provide a load-balanced cut within a fixed threshold \epsilon, so that for instance |S_{k}| \le (1+\epsilon)|V|/p. Therefore, CHACO’s cost function is similar to the ratio cut, but the clustering at the end is biased towards providing a load-balanced partitioning, while still minimizing the edge cuts. Also, CHACO implements spectral bisection, so when comparing to it I split the graph into only two partitions.

Figures 5 and 6 show the performance and quality of both spectral schemes, respectively. Notice that the GPU spectral scheme using Lanczos often obtains the solution faster, but with variable quality compared to the CPU spectral scheme in CHACO, which also uses a variation of the Lanczos method. On the other hand when using preconditioned LOBPCG the GPU implementation is usually faster, and most of the time obtains a higher quality solution as measured by the cost functions. The detailed results of these experiments can be found in our technical report [2].

Figure 5: Speedup of spectral graph partitioning on GPU vs. CPU (CHACO).
Figure 6: Quality of Spectral Graph Partitioning on GPU vs. CPU (CHACO)
Figure 6: Quality of Spectral Graph Partitioning on GPU vs. CPU (CHACO).

Finally, as mentioned earlier there exist many different partitioning and clustering strategies. In particular, some of the popular approaches for providing a balanced cut of a graph use multi-level schemes, implemented in software packages such as METIS. Both spectral and multi-level schemes are global methods that work on the entire graph, in contrast to local heuristics, such as the Kernighan-Lin algorithm.

It is interesting to compare the quality of spectral and multi-level schemes in terms of the the edge-cut and cost function obtained by them. The numerical experiments shown in Figures 7 and 8 plot the ratio of these quantities (the cost obtained by METIS divided by the cost obtained by the GPU spectral scheme), for 30 partitions. The result trends indicate that the behavior of the spectral and multi-level schemes is starkly different for two classes of problems: (i) meshes arising from discretization of PDEs and (ii) social network graphs that often have power-law-like distributions of edges per node. My conjecture is that the difference in quality between these schemes results from the fact that  multi-level schemes often rely on local information to construct a graph hierarchy that is used to partition the graph.

Notice that for PDEs the quality of the partitioning obtained by both schemes is essentially the same, while for networks with high degree nodes, such as social networks, spectral schemes can obtain significantly higher quality partitions. Even though in our experiments the time taken by the spectral schemes is often larger than that taken by the multi-level schemes, I think that spectral schemes can be a good choice in applications where quality is important. For example, in sparse linear algebra applications, even modest improvements in the quality of partitioning can lead to a significant impact on the overall application performance, so the extra partitioning cost of a spectral scheme may be worthwhile.

Figure 7: Spectral vs. Multi-level Schemes (METIS) for PDEs.
Figure 8: Spectral vs. Multi-level Schemes (METIS) for Networks.


I hope that after reading this blog post you have learned some of the intuition behind the spectral graph partitioning/clustering scheme and how it compares to other similar algorithms. A more formal treatment of the subject, with precise derivation of the theoretical results and detailed numerical experiments, can be found in our technical report [2].

The numerical experiments show that spectral partitioning on GPUs can outperform spectral partitioning on the CPU by up to 7x. Also, it is clear that multi-level schemes are a good choice for partitioning meshes arising from PDEs, while spectral schemes can achieve high quality partitioning and clustering on network graphs with high-degree nodes, such as graphs of social networks.

If you need to accelerate graph algorithms in your applications, check out the new GPU-accelerated nvGRAPH library. You can also read more about nvGRAPH in the post “CUDA 8 Features Revealed”. We are considering adding spectral partitioning to nvGRAPH in the future. Please let us know in the comments if you would find this useful.

A Note on Drawing Graphs

Finally, note that eigenvectors of the Laplacian matrix have many other applications. For example, they can be used for drawing graphs. In fact, the graph drawings in this blog were done with them. The interpretation of eigenvectors for this application has been studied in [3].


[1] M.E.J. Newman, “Modularity and Community Structure in Networks”, Proc. National Academy of Science, Vol. 103, pp.  8577–8582, 2006.

[2] M. Naumov and T. Moon, “Parallel Spectral Graph Partitioning”, NVIDIA Research Technical Report, NVR-2016-001, March, 2016.

[3] Y. Koren, “Drawing Graphs by Eigenvectors: Theory and Practice”, Computers & Mathematics with Applications, Vol. 49, pp. 1867–1888, 2005.

]]> 0

Brad: Can you talk a bit about your current research?



Suppose we obtained m users’ ratings on  items (say, movies). If user u rated item v, we use r_{uv} as the non-zero element of R at position (u, v). We want to minimize the following cost function JTo avoid overfitting, we use weighted-λ-regularization proposed in [1], where n_{x_u} and r_uv denote the number of total ratings on user u and item v, respectively.

J = \sum_{u,v} (r_{uv} - x^T_u \theta_v)^2 + \lambda (\sum_{u} n_{x_u} \|x_u\|^2 + \sum_{v} n_{theta_v} \| \theta_v \|^2)

Many optimization methods, including Alternating Least Squares (ALS) and Stochastic Gradient Descent (SGD) [1] have been applied to minimize J. We adopt the ALS approach, which first optimizes X while fixing \theta, and then optimizes theta while fixing X. That is, in one iteration we need to solve these two equations alternatively:

\sum_{r_{uv} \neq 0} (\theta_v \theta_v^T + \lambda I) \cdot x_u = \theta^T \cdot R_{u*}^T \text{    (1)}

\sum_{r_{uv} \neq 0} (x_u x_u^T + \lambda I)  \cdot \theta_v = X^T \cdot R_{*v} \text{      (2)}

B: What are some of the biggest challenges in this project?

W: Recently, the GPU has emerged as an accelerator for parallel algorithms. It has massive computing power (typically 10x higher floating-point operations per second (FLOP/s) versus a CPU) and memory bandwidth (typically 5x versus a CPU), but with limited amount of control logic and memory capacity. Particularly, the GPU’s success in deep learning inspired us to try GPUs for MF. In deep learning, the computation is mainly dense matrix multiplication which is compute bound. As a result, a GPU can train deep neural networks 10x as fast as a CPU by saturating its FLOP/s. However, unlike deep learning, a MF problem involves sparse matrix manipulation which is usually memory bound. Given this, we want to explore a MF algorithm and a system that can still leverage the GPU’s compute and memory capability.

Given the widespread use of MF, a scalable and speedy implementation is very important. There are two challenges:

  1. On a single GPU, MF deals with sparse matrices, which makes it difficult to utilize the GPU’s compute power.
  2. We need to scale to multiple GPUs to deal with large, industry-scale problems (hundreds of millions of non-zero ratings).

B: What NVIDIA technologies are you using to overcome the challenges?

W: We developed cuMF, a CUDA-based matrix factorization library that optimizes the alternating least square (ALS) method to solve very large-scale MF. cuMF uses CUDA 7.0+, cuBLAS and cuSPARSE, and has been tested on both Kepler (Tesla K40/K80) and Maxwell (Titan X) GPUs.

cuMF achieves excellent scalability and performance by innovatively applying the following techniques on GPUs:

(1) On a single GPU, we optimize memory access in ALS by various techniques including reducing discontiguous memory access, retaining hotspot variables in faster memory, and aggressively using registers. These techniques allow cuMF to get closer to the roofline performance of a single GPU. See Figure 2.

Figure 2
Figure 2. Memory optimization in cuMF. To solve X using Θ, cuMF needs to obtain the sum from Equation 1 for for different xus. CuMF reads θv through texture memory, so that θvs used by different users can be cached and reused. θv are loaded into shared memory, and the partial summation result of Equation 1 is tiled and stored in registers until the summation over all rated items is done.

(2) On multiple GPUs, we add data parallelism to ALS’s inherent model parallelism. Data parallelism needs a faster reduction operation among GPUs, leading to (3).

In distributed machine learning, model parallelism and data parallelism are two common schemes. Model parallelism partitions parameters among multiple learners with each one learns a subset of parameters. Data parallelism partitions the training data among multiple learners so that each one learns all parameters from its partial observation. These two schemes can be combined when both the number of model parameters and the training data are large.

Model parallelism. As seen from equation 1, the solution of each x_u is independent so model parallelism is straightforward.

Data parallelism. To tackle large-scale problems, on top of the existing model parallelism, we design a data-parallel approach. When  is big and cannot stay in one GPU, we re-write the summation matrix in eq. 1 to its data-parallel form as:

\sum_{r_{uv} \neq 0} (\theta_v \theta_v^T + \lambda I) = \sum_{i=1}^{p} \sum_{r_{uv} \neq 0}^{\text(GPU)_i} (\theta_v \theta_v^T + \lambda I)

That is, instead of transferring all \theta_vs to one GPU to calculate the sum, it calculates a local  sum on each GPU with only the local \theta_vs, and reduces (aggregates) many local \theta_vs later (See Figure 3).

Figure 3. Parallelism on multiple GPUs. ΘT is partitioned evenly and vertically, and stored on p GPUs (p=4 in the figure). X is partitioned evenly and horizontally, and solved in batches, achieving model parallelism. Each X batch is solved in parallel on p GPUs, each with ΘT 's partition on it, achieving data parallelism.

(3) We also developed an innovative topology-aware, parallel reduction method to fully leverage the bandwidth between GPUs. This way, cuMF ensures that multiple GPUs are efficiently utilized simultaneously.

The experimental results demonstrate that with up to four Titan X GPUs on one machine, cuMF is (1) competitive compared with multi-core methods, on medium-sized problems; (2) much faster than vanilla GPU implementations without memory optimization; (3) 6 to 10 times as fast, and 33 to 100 times as cost-efficient as distributed CPU systems, on large-scale problems; (4) more significantly, able to solve the largest matrix factorization problem ever reported.

Table 1. Speed and cost of cuMF on one machine with four GPUs, compared with three different CPU systems in the cloud. Note: experiment details are in our paper [2].
Figure 2. cuMF as Spark accelerator.
Figure 4. cuMF as Spark accelerator.

CuMF can be used standalone, or to accelerate the ALS implementation in Spark MLlib. We modified Spark’s ml/recommendation/als.scala (code) to detect GPUs and offload the ALS forming and solving to GPUs, while retain shuffling on Spark RDDs.

This approach has several advantages. First, existing Spark applications relying on mllib/ALS need no change. Second, we leverage the best of Spark (to scale-out to multiple nodes) and GPUs (to scale-up in one node).

B: What is the impact of your research to the larger data science community?

W: First, by using cuMF for matrix factorization, a large category of workloads that derive latent features from observations can be accelerated. These workloads include recommendation, topic modeling and word embedding.

Second, cuMF sheds light on accelerating machine learning algorithms involving sparse and graph data. We are glad to see that NVIDIA is also moving toward this direction by releasing nvGRAPH in the forthcoming CUDA 8.

B: When and why did you start looking at using NVIDIA GPUs?

W: We started looking at using NVIDIA GPUs to solve this problem in late 2014. We thought NVIDIA GPUs massive floating point throughput and large memory bandwidth would help in accelerating MF.

What is your GPU Computing experience and how have GPUs impacted your research?

It was a very smooth experience learning and using CUDA. I learned a lot from the two online parallel programming courses that both use CUDA — Intro to Parallel Programming and Heterogeneous Parallel Programming. There are a few excellent books that I can refer to, including Programming Massively Parallel Processors: A Hands-on Approach, by David Kirk and Wen-mei Hwu. In addition, I always get useful information and answers from the NVIDIA Forums.

GPUs have helped cuMF to be the state-of-the-art MF solution in terms of speed. Our research result is to be published in HPDC [2], a premier high performance computing conference.

B: What’s next for cuMF? 

W: We have open-sourced cuMF. In future work we plan to make cuMF more user-friendly by providing a Python interface. We also want to further enhance it by using the latest hardware and software features such as NVLink, FP16 and NCCL.

B: In terms of technology advances, what are you looking forward to in the next five years?

W: In my personal opinion, in the next five years HPC technology is going to impact many aspects of machine learning. HPC has enabled and will keep enabling much more sophisticated machine learning, dealing with even bigger data, in data center and on devices.

Learn More!

You can read the HPDC paper describing cuMF linked in reference [2], below, or watch the GTC 2016 talk. If you’d like to try out cuMF, check out the source code on Github.

If you are interested in data analytics applications on GPU, you may be interested in the new nvGraph library.


[1] Y. Koren, R. M. Bell, and C. Volinsky, Matrix factorization techniques for recommender systems. Computer, vol. 42, no. 8, pp. 30- 37, 2009.

[2] Wei Tan, Liangliang Cao, Liana Fong. Faster and Cheaper: Parallelizing Large-Scale Matrix Factorization on GPUs. HPDC, Kyoto, Japan, 2016

]]> 1
Train Your Reinforcement Learning Agents at the OpenAI Gym Wed, 27 Apr 2016 17:00:51 +0000 Today OpenAI, a non-profit artificial intelligence research company, launched OpenAI Gym, a toolkit for developing and comparing reinforcement learning algorithms. It supports teaching agents everything from walking to playing games like Pong or Go.

John Schulman is a researcher at OpenAI.
John Schulman is a researcher at OpenAI.


What is OpenAI? What is your mission?






While reinforcement learning is focused on making good decisions, supervised and unsupervised learning are mostly focused on making predictions. However, there are lots of connections between these areas, some of which are active topics of research. Besides its different focus, the sequential nature of reinforcement learning also sets it apart from most supervised learning problems. In reinforcement learning, the agent’s decisions affect what input data it receives, i.e., what situations it ends up in. That makes it harder to develop stable algorithms and necessitates exploration—the agent needs to actively seek out unknown territory where it might receive large rewards.

What is OpenAI Gym, and how will it help advance the development of AI?

OpenAI Gym is a toolkit for developing and comparing reinforcement learning algorithms. It includes a curated and diverse collection of environments, which currently include simulated robotics tasks, board games, algorithmic tasks such as addition of multi-digit numbers, and more. We expect this collection to grow over time, and we expect users to contribute their own environments. These environments all expose a common interface, making it possible to write generic algorithms that can be applied to many different environments.

OpenAI Gym also has a site where people can post their results on these environments and share their code. The goal is to make it easy for people to iterate on and improve RL algorithms, and get a sense for which algorithms really work.

Just to give you a sense of what the code looks like, here’s how you can create one of our environments (the classic cart-pole task, where the aim is to balance a vertical pole on a rolling cart), simulate some random actions, and then submit the results to the scoreboard. (In reality, you’ll only want to submit your results if you’ve implemented a learning algorithm.)

import gym
env = gym.make("CartPole-v0") # Creates the classic cart-pole environment
env.monitor.start("/tmp/gym-results", algorithm_id="random") 
# ^^^ starts logging statistics and recording videos
for _ in xrange(1000):
    action = env.action_space.sample() # random action
    observation, reward, done, info = env.step(action)
gym.upload("/tmp/gym-results", api_key="JMXoHnlRtm86Fif6FUw4Qop1DwDkYHy0")
# ^^^ upload stats to the website

This snippet does not include any learning or training—that would require much more code. Soon we are hoping to show clean implementations of various important algorithms on OpenAI Gym environments. If you’re interested, keep an eye on our website.

How are neural networks used in reinforcement learning?

To answer this, I’ll need to talk a little bit about what RL algorithms learn. Some reinforcement learning algorithms are focused on learning a policy, which is a function that takes in observations (e.g., camera images) and outputs actions (e.g., motor torques). Other algorithms are focused on learning value functions, which measure the goodness of states (i.e., the state of the world) and actions. Since we usually don’t have access to the full state of the world, we typically use one or more past observations instead. A Q-function (a type of value function) measures the goodness of a state-action pair (s, a), i.e., Q(s, a) tells you “how much reward am I going to receive, if I’m in state s, and I choose action a”. Given this Q-function, you can simply choose the action with highest expected reward. In other words, the Q-function defines a policy. The following video shows a Deep Q-Network (DQN) trained to play Breakout on OpenAI Gym.

Policy-based algorithms and Q-function-based algorithms are very similar at their core, and we can use neural networks to represent the policies and Q-functions. For example, when playing Atari games, the input to these networks is an image of the screen, and there is a discrete set of actions, e.g. {NOOP, LEFT, RIGHT, FIRE}. As a Q-function for this task, you can use a convolutional neural network that takes the screen image as input and outputs a real number for each of the 4 actions, specifying the goodness of that action (Mnih et al, 2013). As a policy for this task, you can use a convolutional neural network with similar architecture, which instead outputs the probability of each action. (Guo et al., 2014Schulman et al., 2015).

Architectures used in (Schulman et al., 2015), for Atari game playing and simulated robot control.
Architectures used in (Schulman et al., 2015), for simulated robot control (top) and Atari game playing (bottom).

What makes OpenAI Gym unique? Are there other similar open source environments currently available?

There are various open-source collections of environments, including but not limited to RL-Glue, RLPy, and Arcade Learning Environment. We’ve drawn inspiration and code from these libraries. OpenAI Gym also integrates work done recently by researchers at UC Berkeley on benchmarking deep reinforcement learning algorithms. A paper describing this benchmark is available on the ArXiv and will be presented at ICML this year.

OpenAI Gym goes beyond these previous collections by including a greater diversity of tasks and a greater range of difficulty (including simulated robot tasks that have only become plausibly solvable in the last year or so). Furthermore, OpenAI Gym uniquely includes online scoreboards for making comparisons and sharing code.

Who will use OpenAI Gym? How will AI researchers benefit from RL-Gym?

We hope to make OpenAI Gym accessible to people with different amounts of background. Users with no background in RL can download the codebase and get started experimenting in minutes. They can visit the scoreboards for different environments and download the code linked to solutions. Then they can verify the solutions (an important and useful service!) and tinker with them.

AI researchers will be able to use the included environments for their RL research. Each environment is semantically versioned, making it possible to report results in papers and have them remain meaningful. Researchers will be able to compare their algorithms’ performance with others on the scoreboard, and find code for the leading algorithms.

Do you have plans to accelerate OpenAI Gym with NVIDIA GPUs? How will GPUs benefit your work?

GPUs are becoming indispensable for learning problems that involve large neural networks. We will be using GPUs for training networks on the larger-scale tasks, and we expect many of our users to do the same.

Also, our current environments include robotics simulations, and we expect to include more complicated and realistic 3D environments in the future. GPUs will enable high-quality rendering and physics simulation in these more complicated 3D environments.

Would more realistically rendered environments be useful to transfer learning to the real-world?

Yes, I believe that photorealistic rendering is good enough to allow robots to be trained in simulation, so that the learned policies will transfer to the real world. Lots of exciting possibilities lie ahead.

What’s next for OpenAI and OpenAI Gym?

We’ll soon be publishing results from some of our ongoing research projects in unsupervised learning and reinforcement learning. We’re excited to see what users do with OpenAI Gym, and plan to continue building it into a tool that’s great for the research community as well as newcomers to the field.

To learn more see the OpenAI Gym announcement blog post. For an introduction to Deep Learning, check out the Deep Learning in a Nutshell series.

]]> 3


Figure 1: Illustration of the All-Reduce collective.
Figure 1: Illustration of the All-Reduce collective.

Collective Communication


Another common collective is all-gather, where each of K processors begins with an independent array of N values, and collects data from all other processors to form a result of dimension N*K, as Figure 2 shows.

Figure 2: Illustration of the All-Gather collective.
Figure 2: Illustration of the All-Gather collective.

Broadcast is a third example. Here an N-element buffer on one processor is copied to all other processors, as Figure 3 shows.

Figure 3: Illustration of the Broadcast collective.
Figure 3: Illustration of the Broadcast collective.

All of the above collectives can be performed either “out-of-place,” with separate input and output buffers, or “in-place,” with overlapping input and output.

Figure 4: A common PCIe topology for 4 GPUs attached to a single CPU. Red arrows represent PCIe x16 connections.
Red arrows represent PCIe x16 connections.

There are numerous approaches to implementing collectives efficiently. However, it is critical that our implementation takes the topology of interconnects between processors into account. For example, consider a broadcast of data from GPU0 to all other GPUs in the PCIe tree topology pictured below.

A two-step tree algorithm is a common choice in this situation. In the first step the data is sent from the GPU0 to a second GPU, and in the second step both of these send data to the remaining processors. However, we have a choice. We can either send data from GPU0 to GPU1 in the first step and then GPU0 to GPU2 and GPU1 to GPU3 in the second, or we can perform the initial copy form GPU0 to GPU2 and then GPU0 to GPU1 and GPU2 to GPU3 in the second step. Examining the topology, it is clear that the second option is preferred, since sending data simultaneously from GPU0 to GPU2 and GPU1 to GPU3 would cause contention on the upper PCIe links, halving the effective bandwidth for this step. In general, achieving good performance for collectives requires careful attention to the interconnect topology.

Figure 5: Ring order of GPUs in PCIe tree.
Figure 5: Ring order of GPUs in PCIe tree.

To optimize Broadcast bandwidth, an even better approach is to treat the PCIe topology above as a ring.

The broadcast is then performed by relaying small chunks of the input around the ring from GPU0 to GPU3. Interestingly, ring algorithms provide near optimal bandwidth for nearly all of the standard collective operations, even when applied to “tree-like” PCIe topologies. But note that selecting the correct ring order remains important.

GPU Collectives with NCCL

In order to provide maximum bandwidth, NCCL implements ring-style collectives. NCCL implicitly indexes the GPUs into the optimal ring order under the hood. This provides great performance for your application while freeing you from having to worry about specific hardware configurations.

Many collectives requires a buffer for intermediate results. In order to minimize memory overhead to a few MB on each GPU, NCCL splits large collectives into many small chunks. It would be highly inefficient to launch separate kernels and cudaMemcpy calls for every step and chunk of a collective algorithm. Instead, each collective is implemented by a monolithic CUDA kernel. NCCL makes extensive use of GPUDirect Peer-to-Peer direct access to push data between processors. Where peer-to-peer direct access is not available (e.g., when traversing a QPI interconnect), the pushed data is staged through a buffer in pinned system memory. Similarly, synchronization is performed by polling on volatile variables in device (or pinned system) memory.

Internally, NCCL implements each collective in terms of three primitives: Copy, Reduce, and ReduceAndCopy. Each of these is optimized to efficiently transfer fine-grained slices of data (4-16KB) between GPUs. The kernels have also been specifically optimized to achieve maximum bandwidth at low occupancy. As a result, NCCL can saturate a PCIe 3.0 x16 interconnect using a single block of CUDA threads. This leaves the bulk of the GPU free to execute compute tasks concurrently with the communication.

NCCL currently supports the all-gather, all-reduce, broadcast, reduce, and reduce-scatter collectives. Any number of GPUs can be used, as long as they reside in a single node.

Using NCCL

The open-source release of NCCL is available on Github. The library should build on any common Linux distribution and is compatible with CUDA 7.0 and later. CUDA supports direct access only for GPUs of the same model sharing a common PCIe root hub. GPUs not fitting these criteria are still supported by NCCL, though performance will be reduced since transfers are staged through pinned system memory.

The NCCL API closely follows MPI. Before performing collective operations, a communicator object must be created for each GPU. The communicators identify the set of GPUs that will communicate and maps the communication paths between them. We call the set of associated communicators a clique. There are two ways to initialize communicators in NCCL. The most general method is to call ncclCommInitRank() once for each GPU.

ncclResult_t ncclCommInitRank(ncclComm_t* comm, 
                              int nGPUs,
                              ncclUniqueId cliqueId,
                              int rank);

This function assumes that the GPU belonging to the specified rank has already been selected using cudaSetDevice(). nGPUs is the number of GPUs in the clique. cliqueId allows the ranks of the clique to find each other. The same cliqueId must be used by all ranks. To achieve this, call ncclGetUniqueId() in one rank and broadcast the resulting uniqueId to the other ranks of the clique using the communication framework of your choice (for example, MPI_Bcast).

The last argument to ncclCommInitRank() specifies the index of the current GPU within the clique. It must be unique for each rank in the clique and in the range [0, nGPUs). This index is used, for example, to identify the source GPU for a broadcast operation, or to order the contributions to an all-gather.

Upon successful initialization, ncclCommInitRank() returns ncclSuccess and *comm is set to the new NCCL communicator object.

Internally, ncclInitRank() performs a synchronization between all communicators in the clique. As a consequence, it must be called from a different host thread for each GPU, or from separate processes (e.g., MPI ranks). Specifically, calling ncclInitRank() in a single threaded loop over GPUs will result in deadlock.

The second way to initialize the communicators is to use ncclCommInitAll(). This is essentially a convenience routine that spares you the effort of spawning extra host threads to initialize NCCL in an otherwise single-threaded application.

ncclResult_t ncclCommInitAll(ncclComm_t* comms,
                             int nGPUs,
                             int* devList);

The comms argument now refers to an array of ncclComm_t objects, one for each of the nGPUs ranks in the clique. devList specifies which CUDA device gets associated with each rank. With the communicator object initialized, you can call the collectives through their host functions, such as ncclAllReduce().

ncclResult_t ncclAllReduce(void* sendoff,
                           void* recvbuff,
                           int count,
                           ncclDataType_t type,
                           ncclRedOp_t op,
                           ncclComm_t comm,
                           cudaStream_t stream);

Documentation for each collective is provided in nccl.h. Most arguments have analogues in the MPI collectives API. The notable exception is the stream argument. Similar to many CUDA “Async” routines, NCCL collectives schedule the operation in a stream but may return before the collective is complete. By queuing collectives and other compute kernels in separate streams, the GPU is able to overlap collective communication with more compute-intensive workloads. To maximize overlap, schedule NCCL collectives on high-priority streams to allow them to slip in among compute intensive grids.

ncclAllReduce() must be called once for each communicator in the clique, each call providing its own send and receive buffers, etc. NCCL collectives assume that the appropriate CUDA context is already current. Because NCCL is asynchronous, a simple loop can be used to initiate a collective from a single threaded application.

for (int gpu=0; gpu<nGPUs; ++gpu) {


The performance that NCCL collectives achieve depends on the exact topology of the computer. In the best case, all GPUs share peer access. This is most common for workstations equipped with a few GPUs. Larger servers usually sport dual CPUs, with some GPUs segregated on different IO hubs. Figure 6 shows NCCL bandwidth for various collectives measured on an NVIDIA Digits DevBox equipped with 4 GeForce GTX Titan X GPUs arranged as in Figure 4 above.

Figure 6: Link bandwidth achieved by various NCCL collectives.
Figure 6: Link bandwidth achieved by various NCCL collectives.

The red bar at 10.4 GB/s represents the bandwidth achieved by a large cudaMemcpy between two of the GPUs (specifically GPU 0 and 2 in Figure 4). NCCL sustains a high percentage of this peak bandwidth while performing communication among all four GPUs.

Future Directions & Final Remarks

The goal of NCCL is to deliver topology-aware collectives that can improve the scalability of your multi-GPU applications. By using NCCL you can get great performance without having to think about low-level hardware details.

Simple code examples for both single-process and MPI applications are distributed with NCCL. As a research project, we welcome your feedback as we continue to evolve the project. Please try it out and let us know what you think! For an update on the latest developments, come see my NCCL talk at GTC.

]]> 2
Optimizing Recurrent Neural Networks in cuDNN 5 Wed, 06 Apr 2016 21:15:57 +0000 Figure 1: cuDNN 5 + Torch speedup vs. Torch GPU implementation, M40, Intel® Xeon® Processor E5-2698 Network A: RNN size 2560, Wordvec size 2560, num layers 1, Seq length 200, max epochs 1 Network B: RNN size 256, num layers 3, max epochs 50, batch size 64 Network C: RNN size 256, Wordvec size 256, num layers 1, Seq length 1000, max epochs 1


  • Faster forward and backward convolutions using the Winograd convolution algorithm;
  • 3D FFT Tiling;
  • Spatial Transformer Networks;




Step 1: Optimizing a Single Iteration

The following equations govern the forward propagation of data through an LSTM unit. Figure 2 shows a diagram of an LSTM unit.

c_{t}=f_{t}\circ c_{t-1}+i_{t}\circ c'_{t}

Figure 2: A diagram of an LSTM unit
Figure 2: A diagram of an LSTM unit

From a computational perspective this boils down to eight matrix-matrix multiplications (GEMMs)—Four with input i, four with input h—and lots of point-wise operations.

The starting point for this case-study is an LSTM implemented operation-by-operation. For each iteration, for each layer, the implementation calls cuBLAS sgemm to perform each of the eight GEMMs, and hand-written CUDA kernels to call each of the point-wise operations. Pseuduocode for the method follows.

for layer in layers:
  for iteration in iterations:
     perform 4 SGEMMs on input from last layer
     perform 4 SGEMMs on input from last iteration
     perform point-wise operations

As a benchmark I measure run time per step, per layer on a Tesla M40 GPU. My benchmark LSTM has 512 hidden units and computes mini batches of size 64. The performance of this baseline implementation is fairly poor, achieving approximately 350 GFLOPS on the M40. The peak performance of this GPU is around 6000 GFLOPs, so there is a lot of room to improve. Let’s get started.

Optimization 1: Combining GEMM Operations

GPUs have very high peak floating-point throughput, but they need a lot of parallelism to approach this peak. The more parallel work you give them, the higher the performance they can achieve. Profiling this LSTM code shows that the GEMM operations use significantly fewer CUDA thread blocks than there are SMs on the GPU, indicating the GPU is massively underused.

GEMM is typically parallelized over the output matrix dimensions, with each thread computing many output elements for maximum efficiency. In this case each of the eight output matrices comprises 512×64 elements, which results in only four thread blocks. Ideally can run significantly more blocks than the GPU has SMs, and to maximize the theoretical occupancy for this kernel at least four blocks per SM (or 96 in total) are needed. (See the CUDA Best Practices guide for more on occupancy.) 

If n independent matrix multiplications share the same input, then they can be combined into one larger matrix multiplication with an output n times larger. The first optimization is therefore to combine the four weight matrices operating on the recurrent step into one weight matrix, and to combine the four weight matrices operating on the input into another. This gives us two matrix multiplications instead of eight, but each is four times the size and has four times the parallelism (16 blocks per GEMM). This optimization is fairly common in most framework implementations: it’s a very easy change that leads to a good speedup: the code runs roughly 2x faster.

Optimization 2: Streaming GEMMs

Even with the larger combined GEMMs, performance is still limited by lack of parallelism: there are 16 blocks instead of four, but the target is to have at least 96. The two remaining GEMMs are independent so they can be computed concurrently using CUDA streams. This doubles the number of possible concurrent blocks to 32.

Optimization 3: Fusing Point-wise Operations

Figure 3 shows that now a lot of time is spent in point-wise operations. There’s no need to do these in separate kernels; fusing them into a single kernel reduces data transfers to and from global memory and significantly reduces kernel launch overhead. 

Figure 3: NVIDIA Visual Profiler output showing the operations in a single cell. Pointwise operations were originally taking more than half of the runtime, after optimization point-wise operations are only a small fraction.

At this point I’m fairly happy with the performance of a single iteration: the majority of the computation is in the GEMMs, and they’re exposing as much parallelism as they can. This implementation is about 5x faster than the baseline implementation, but there are more improvements to come.

for layer in layers:
  for iteration in iterations:
     perform sgemm on input from last layer in stream A
     perform sgemm on input from last iteration in stream B
     wait for stream A and stream B
     perform point-wise operations in one kernel

Step 2: Optimizing Many Iterations

In an RNN the operations for a single iteration are repeated many times. This means it’s important to have those operations running efficiently, even if this comes at an upfront cost.

Optimization 4: Pre-Transposing the Weight Matrix

When performing a GEMM the standard BLAS API allows you to transpose either of the two input matrices. Some of the four combinations of transpose/not-transposed run slightly faster or slower than others. Depending on the way that the equations are mapped to the computation, a slower version of the GEMM may be used. By performing a transpose operation up-front on the weight matrix, each step can be made slightly faster. This comes at the cost of the transpose, but that is fairly cheap, so if the transposed matrix is to be used for more than a few iterations it is often worth it.

Optimization 5: Combining Input GEMMs

In many cases all of the inputs are available at the start of the RNN computation. This means that the matrix operations working on these inputs can be started immediately. It also means that they can be combined into larger GEMMs. While at first this may seem like a good thing (there’s more parallelism in the combined GEMMs), propagation of the recurrent GEMMs depends upon the completion of the input GEMMs. So there’s a tradeoff: combining input GEMMs gives more parallelism in that operation, but also prevents overlap with the recurrent GEMMs. The best strategy here depends a lot on the RNN hyperparameters. Combining two input GEMMs works best in this case.

for layer in layers:
  transpose weight matrices
  for iteration in iterations / combination size:
     perform sgemm on combined input from last layer in stream A
     for sub-iteration in combination size:
        perform sgemm on input from last iteration in stream B
        wait for stream A 
     wait for stream B
     for sub-iteration in combination size;
        perform pointwise operations in one kernel

Step 3: Optimization with Many Layers

The final step is to consider optimization between layers. Once again, there’s quite a lot of parallelism to be found here. Figure 4 shows the dependency graph for an RNN. As iteration n of a given layer only depends on iteration n-1 of that layer and iteration n of the previous layer it is possible to start on a layer before you’ve finished on the previous layer. This is really powerful: if there are two layers there is twice as much parallelism available.

Figure 4: As dependencies are resolved a wavefront of operations moves through the network.

Going up from one to four layers, throughput increases by roughly 1.7x: from 2.3 TFLOPs to 3.9 TFLOPs. At this point the gains from exposing more parallelism are starting to become more limited. Compared to the original implementation which only had four blocks running at any given time, this implementation can run up to 128 blocks concurrently. This is enough to make use of all of the M40’s resources, achieving nearly 70% of peak floating point performance and running more than 10x faster than the original implementation.

The following table shows the performance achieved after each of the optimizations I have described and the speedup vs. the baseline code.

Optimization GFLOPS Speedup
Baseline 349 (1.0x)
Combined GEMMs 724 2.1x
GEMM Streaming 994 2.8x
Fused point-wise operations 1942 5.5x
Matrix pre-transposition 2199 6.3x
Combining Inputs 2290 6.5x
Four layers 3898 11.1x


Propagating gradients through the backward pass is very similar to propagating values forward. Once the gradients have been propagated the weights can be updated in a single call spanning all iterations: there are no longer any recurrent dependencies. This results in a very large efficient matrix multiplication.


To get the best performance out of Recurrent Neural Networks you often have to expose much more parallelism than direct implementation of the equations provides. In cuDNN we’ve applied these optimizations to four common RNNs, so I strongly recommend that you use cuDNN 5 if you are using these RNNs in your sequence learning application.

For more information:

Other resources:

]]> 4