Parallel Computing with Mathematica

Parallel Computing Toolkit and gridMathematica

Roman E. Maeder
Wolfram Research

Topics

Parallel Computing Toolkit

gridMathematica

Examples of Applications

Support for Parallel Computing

What is needed to add parallel computing to a programming environment, such as Mathematica?

Parallel computing ( a system programmer's view)

Parallel programming requires support for:

In this context

Parallel algorithms are expressed directly in Mathematica's programming language, no compilation is necessary.

Such an expression may be a simple calculation or a larger program.

Parallel Computing Toolkit

The Parallel Computing Toolkit is a Mathematica application package that provides everything that is needed for developing and running parallel computations.

Features of Parallel Computing Toolkit

The Application Package

All required support code is contained in an application package, which can simply be imported into a running kernel, the master kernel, which controls the parallel computation.

Needs["Parallel`"]

Note that the master kernel is normally connected to a Mathematica frontend, through which the user interacts with the parallel computation.

The master kernel interacts with any number of remote slave kernels, which perform part of a parallel computation. Starting, initialization, monitoring, and stopping slave kernels is done through the master kernel.

Requirements

All scheduling happens in the master kernel. No extra software is required on the slave kernels. No support system for parallel computing, such as PVM or MPI is required.

Slave kernels are normally started once and are then used for any number of computations. The Parallel Computing Toolkit contains its own scheduler and load balancing tools.

Supported Configurations

Master and slave kernels can run on a multi-processor machine, a tightly coupled network of processors (a cluster), a network of workstations or any remote machines accessible through TCP/IP.

Presentation Setup (Sun X4600 at Wolfram Technology Conference 2006)

In this demonstration, the frontend runs on my laptop computer, the master kernel and slave kernels run on a 8 Opteron/16 core Sun X4600 machine.

Starting Remote Kernels

Any connection method supported by MathLink can be used to start remote kernels.

LaunchSlave["localhost"]

LaunchSlave["gridhead.wolfram.com","ssh -l maeder `1` '/usr/local/bin/math -mathlink'",
ConnectionType -> LinkLaunch]

LaunchSlave["xserve1","ssh -n `1` '/usr/local/bin/math51 -mathlink -linkmode Connect -linkname `2` -noinit >&/dev/null &'",LinkHost -> "192.168.1.1"]

Simple Calculations

RemoteEvaluate[] sends the same expression to all slave kernels.
Here, we ask them for their unique processor ID.

RemoteEvaluate[$ProcessorID]

After all computations are done, the slave kernels should be shut down (we will restart more kernels in a minute).

CloseSlaves[];

Setup

The normally available slave kernels can be set up as defaults.

$AvailableMachines

Here is the command template needed for connecting to kernels on remote machines.

$RemoteCommand

Additional kernels on the local multiprocessor machine are started like this.

$mathkernel

It is now easy to start them up.

LaunchSlaves[]

This command displays a number of properties of the remote kernels, such as the machine on which they are running, the Mathematica version, and CPU and memory usage so far.

TableForm[RemoteEvaluate[{$ProcessorID, $MachineName, $Version, TimeUsed[],  
              MaxMemoryUsed[]}],TableHeadings -> {None, {"ID", "host", "Mathematica \
Version",  "CPU Time", "Memory"}},TableDepth -> 2]

Parallel Commands

Parallel Versions of Mathematica Functions

These functions have been implemented for parallel evaluation in PCT.

?Parallel*

Rather than go on with more Parallel* special versions, a more general approach has been developed for Version 2. It works for all tasks that operate uniformly on a (long) expression, such as a list.

ParallelEvaluate

A general mechanism, ParallelEvaluate[], allows the distribution of a computation to available processors. It dispatches parts of the computation to remote processors and is flexible to handle a large number of different cases.

?ParallelEvaluate

It computes f[{(Subscript[e, 1]), …, (Subscript[e, n])}] as follows:

Examples of different kinds of parallel computations

ParallelEvaluate[{1 + 2, 2 + 3, 3 + 4, 4 + 5, 5 + 6}]

ParallelEvaluate[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}, 
      PrimeQ]

This is a parallel version of

Map[PrimeQ, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}]

ParallelEvaluate[1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10]

Less obvious example:

ParallelEvaluate[{1, 2, 3, 4, 5, 6, 7},MemberQ[#, 5]&,Or]

ParallelEvaluate[{1, 2, 3, 4, 5, 6, 7, 8},Cases[#, _?OddQ]&]

Semi-Automatic Parallelization

As we saw, to parallelize

Count[{1, 2, 3, 4, 5, 6, 7},_?PrimeQ]

you have to say

ParallelEvaluate[{1, 2, 3, 4, 5, 6, 7},Count[#, _?PrimeQ]&,Plus]

The idea is always the same. Wouldn't this be nice:

ParallelEvaluate[Count[{1, 2, 3, 4, 5, 6, 7}, _?PrimeQ]]

It is quite simple:

Count/:ParallelEvaluate[
        Count[list_, params___]]:=
ParallelEvaluate[list, Count[#, params] &, 
        Plus]

All of these have been overloaded:

$ParallelCommands

ParallelEvaluate[] takes relative performance of remote processors into account when distributing a calculation. You can use any Mathematica calculation to measure the actual performance and calibrate the processor speed estimates.

Concurrency

The building blocks of concurrency:

Note: a process is simply an expression to be evaluated.

Starting Processes, Waiting for Results

Enter the given expression into the process queue, return a process id.

Queue[cmd]

Waits for one of the given processes to finish.

WaitOne[pid_list]

Wait for all processes or a single process, return list of results or single result, respectively.

Wait[pid_list], Wait[pid]

Examples

Queue a process and wait for it to finish.

j1 = Queue[1 + 1]

Wait[j1]

Queue several processes to test (2^i) - 1 for primality.

pids = Queue[PrimeQ[2^# - 1]] & /@ {1, 2, 3, 4, 5}

Nondeterministic waiting for a job and a convenient use of the return value to drain the queue. This command can be repeated as long as the list of remaining pids is nonempty. Each time, one result is returned.

{res, pid, pids} = WaitOne[pids]

Finally, wait for all remaining processes.

Wait[pids]

Nondeterminism

The WaitOne[] primitive implements nondeterminism. You can use this, for example, to try out different methods to solve a problem, or to search for data with certain properties.

Scheduling

Users can implement their own scheduling policy. PCT contains sample implementations of FIFO and priority queues.

Automatic Parallel Code Generation

Wait[Queue[expr]] is just the evaluation function (or identity).

Wait[Queue[1 + 2]]

To evaluate a function remotely, you can use

Wait[Composition[Queue,f][x]]

This observation gives us a powerful recipe for parallelizing functional operations. This example develops a parallel version of the Mathematica function Map[].

Map[f, {a, b, c}]

We compose f with Queue[]. The result is a list of process IDs that represent the individual function applications.

Map[Composition[Queue, f], {a, b, c}]

Wait[] runs the scheduler and returns with the results, which were computed in parallel on any available slave kernels.

Wait[%]

The powerful functional language allows us to derive parallel code automatically.

Shared Memory

How can processes communicate with each other?

Mathematica kernels do not share memory, even if they run on the same machine. Our computation model is always a distributed memory one.

Virtual Shared Memory: implementing shared memory on distributed memory machines with a software layer that handles all variable accesses, using message passing.

When a process needs the value of a global (shared) variable or wants to set its value, it sends a request back to the scheduler. The scheduler maintains the variables, udpates their values, and sends back the requested values.

Basic Examples

We will work with two slave kernels.

s1 = $Slaves[[1]]; s2 = $Slaves[[2]];

Declare shared variables.

SharedVariables[{x}]

Basic tests. One slave kernel reads the variable, another one increments its value.

x = 0;

RemoteEvaluate[x = x + 1,s1]

RemoteEvaluate[x,s2]

Increment the value in one indivisible step. This operation is important for synchronization.

RemoteEvaluate[x++,s2]

RemoteEvaluate[x,s1]

No special syntax is required to use shared variables. Just do things the normal way. If a variable is shared, all accesses are transparently sent to the master process.

ClearShared[x];

Synchronization

Exclusive access to a variable or other resource can be controlled by a lock. An update of a shared variable implements a lock and allows semaphores and critical sections.

Here is the typical code that implements a critical section with a lock. Each of the five parallel processes accesses the same shared variable y, first reading it, then writing a new value. Here is the code without locking.

SharedVariables[{y}]

y = 0;
ParallelMap[Module[{a},
		Pause[Random[]];
		a = y;
		Pause[0.5  Random[]];
		y = a + 1;
		a + 1]&,
	Range[5]]

As a consequence, the value of y is not the same it would be with sequential execution of the same program.

Locking Primitives

Exclusive access to variables can be achieved with one new primitive operation, a conditional atomic write. In machine code, it is often called test and set or read-modify-write.

TestAndSet[s, e]

Set s to e (just like s=e), but only if s is currently unset or its value is Null. Return the (old or new) value of s. The read-modify-write is atomic.

Example with Synchronization

Before the critical section, the process acquires a lock, which is released after the critical section. In this way no conflicting updates happen, as can be seen by the output, in which all values from 1 to 5 occur.

SharedVariables[{y, lock}]

y = 0;
ParallelMap[Module[{a},
	Pause[Random[]];
	While[TestAndSet[lock, #] =!= #, 
                    Pause[0.1  Random[]]];
	a = y;
	Pause[0.5  Random[]];
	y = a + 1;
	lock = Null;
	a + 1]&,
  Range[5]]

ClearSlaves[];

Advanced Use: Remote Procedure Calls

A shared downvalue can be thought of as a remote procedure call. Evaluating such a function call in a slave kernel carries the expression back to the master kernel, where it is evaluated instead. Here is an example for collecting the results of a search in a large search space.

Parallel Sow and Reap

Parallelizing Sow/Reap: sow on the slave kernels, reap all results in the master kernel. But the calls of Sow[] need to be executed on the master, not the slaves, where they are generated. We define an auxiliary shared downvalue sow[e].

sow[e_] := (Sow[e]; Null)
Protect[sow];
SharedDownValues[sow]

The extra Null prevents the value being needlessly sent back to the remote kerner where it came from.

Reap[
RemoteEvaluate[
        sow[$ProcessorID]]
]

Searching large implicit spaces

To search a large implicit space it is important not the generate all its elements beforehand and then traverse the space; rather, the elements should be generated iteratively and dealt with immediately (visitor paradigm).

This method does not scale well:

space=Table[i, {i, 10}];

Reap[
Scan[Function[i, If[PrimeQ[i], Sow[i]]], 
        space]
]

Instead use:

Reap[
Do[If[PrimeQ[i], Sow[i]], {i, 10}]
]

This example made me reconsider my decision not to implement ParallelScan and ParallelDo.

Parallel Search

ParallelTable in place of the nonexistant ParallelDo would return the whole structure, which defeats the purpose. We can cheat by using an iterator that returns a single value, such as ParallelSum.

Reap[
ParallelSum[
        If[PrimeQ[i], sow[i]; 1, 0], {i, 
          1000}]
]

Note that the result is not sorted!

gridMathematica

gridMathematica provides a quick way to set up and run large calculations by bundling Mathematica kernel licenses and Parallel Computing Toolkit in a way that is easy to install on parallel computers.

A typical installation of gridMathematica involves a master kernel, a license manager, and one slave kernel per available node. The master kernel handles all input, output, and scheduling. It can be controlled from any Mathematica front end or via batch files, either locally or via a remote connection. Once the remote kernels are launched, they are ready to receive commands from the master machine.

gridMathematica is platform independent and can be used on dedicated multiprocessor machines as well as on homogeneous and heterogeneous clusters. The only technical requirement, apart from the ability to run Mathematica, is a TCP/IP connection between the individual computing nodes. This connection allows customers to run the same code on any available machines without any porting work. It also makes it easy to build ad-hoc clusters out of underutilized computers or to take advantage of low-use periods.

Models of Cluster/Grid computations

The connections between the Mathematica frontend (if any), the master kernel, and the slave kernels are all MathLink connections. All of these can be local (between endpoints on the same machine), or networked. Typical configurations include:

Grid Management

Support for popular grid management and job control software is being developed by Wolfram Research and made part of gridMathematica releases as it becomes available.

See the presentation Grid Management Technologies and Mathematica today at 3:30pm for more information on this topic.

Comparison with other Environments

gridMathematica does not need any load-balancing or resource allocation software or parallel libraries, such as PVM or MPI.

All parallel code is written in Mathematica. No compilation is necessary and no common file system is required. Parallel programs can be developed and debugged interactively, and run interactively or in batch mode.

The current version of PCT does not support symmetric multiprocessing.

Example 1: Parallel Search and Nondeterminism

We want to find an irreducible polynomial modulo a prime p. Such polynomials are important in signal processing, data storage, and communication.

Sequential Code Development

This generates a random polynomial of a given degree mod p.

randPoly[n_Integer?Positive, x_, p_] := 
      Random[Integer, {1, p - 1}]*x^n + 
        Sum[Random[Integer, {0, p - 1}]*x^i, {i, 1, n - 1}] + 
        Random[Integer, {1, p - 1}]

It is irreducible if it does not have any factors other than a numerical content.

irreducibleQ[poly_, p_] := 
      With[{fl = FactorList[poly, Modulus -> p]}, 
        Length[fl] < 3 && fl[[-1, 2]] == 1]

Here is an example polynomial modulo 2.

poly = randPoly[10, x, 2]

We test whether it is irreducible.

irreducibleQ[poly, 2]

This loop generates new polynomials until it finds an irreducible one.

findOne[n_, p_, x_:  x] := 
      Module[{poly}, 
While[! irreducibleQ[
              poly = randPoly[n, x, p], p], 
          Null]; 
poly]

Here is an irreducible polynomial of degree 100 mod 2.

findOne[100, 2]

Parallelization

We run such a loop on each slave kernel and wait until one of them returns a result. The remaining computations are then aborted gracefully.

Here is the number of available remote kernels.

Length[$Slaves]

ExportEnvironment[randPoly, irreducibleQ, findOne]

With[{n=512,p=2},
Module[{id, ids},
ids=Table[Queue[findOne[n, p]],{Length[$Slaves]}];
res = $Failed;
While[res === $Failed  &&  Length[ids] > 0,  
{res, id, ids}=WaitOne[ids]
];
ResetSlaves[];
]]

Here is the result.

res

Clear[res, poly]

Example 2: Monte-Carlo Simulation

Components

Monte-Carlo simulation is used to approximate values for which no closed-form solution exists. It works by generating a large number of samples (or scenarios) and averaging the results.

In this example we want to find the theoretical value of a path-dependent financial derivative instrument. A scenario is a possible stock-price movement over a number of days.

Stock Price Movement

A common assumption is that the returns of a stock are log-normally distributed. A scenario consists of the daily stock prices under this assumption.

Needs["Statistics`ContinuousDistributions`"]

This function generates a list of n log-normally distributed values with given mean and standard deviation. These are the simulated returns of the stock.

RandomLogNormal[μ_, σ_, n_] := 
      RandomArray[LogNormalDistribution[μ, σ], n]

RandomLogNormal[0.01, 0.1, 10]

The stock prices themselves are then obtained by successive multiplication.

FoldList[Times, 1, %]

This function generates such a sequence of stock prices.

\(Scenario[μ_, σ_, n_] := 
	
      Rest[FoldList[Times, 1, 
          RandomLogNormal[μ\/n, σ\/\@n, n]]]\)

sample = Scenario[0.01, 0.1, 10]

Valuation

The payoff of an arithmetic-mean call option depends on the mean of the stock price over the lifetime of the option. Here, x represents the strike price.

avgPayoff[prices_, x_] := Max[Mean[prices] - x, 0]

avgPayoff[sample, 1.0]

Averaging

To find the value of the derivative instrument, we perform many simulations and take the average of the results.

MonteCarloValuation[payoff_, σ_, r_, n_, trials_] := 
      Sum[payoff[Scenario[r - (σ^2)/2, σ, n]], {trials}]/
        trials

MonteCarloValuation[avgPayoff[#, 1.0] &, 0.2, Log[1.1], 360, 
      1000]

Each time we evaluate the valuation function we get a different result.

MonteCarloValuation[avgPayoff[#, 1.0] &, 0.2, Log[1.1], 360, 
      1000]

Parallelization

Averaging is a linear operation, so we can simply perform parts of it in parallel and combine the results.

We can export all the code and then use MonteCarloValuation on the slave kernels.

RemoteNeeds["Statistics`ContinuousDistributions`"];

ExportEnvironment[{MonteCarloValuation, avgPayoff, Scenario, 
          RandomLogNormal}];

Running the Simulation

This is the number of kernels available.

Length[$Slaves]

To perform the same calculation on all slave kernels, we can use RemoteEvaluate[].

res = 
      RemoteEvaluate[
        MonteCarloValuation[avgPayoff[#, 1.0] &, 0.2, Log[1.1], 360, 
          2000]]

Because all slave kernels used the same number of simulations, we can simply average the results.

Mean[res]

Example 3: Global Interval Minimization

Introduction

Minimization finds the minimum value (and where it occurs) of a function of several variables. It is used in optimization problems.

Local minima

Examples: FindMinimum[]

Global techniques

Examples:
NMiminimize[]
Global Optimization by Loehle Enterprises
LGO by János Pintér

Interval methods

Algorithm

The domain rectangle is subdivided into (2^d) intervals (d is the dimension of the problem) until either

Exclude a rectangle from further consideration if it can be established that it cannot contain a global minimum. This is the case if the value of the objective function over the rectangle is larger than a known function value at a point lying in some other interval (a feasible point, shown in red). These rectangles are colored blue in this illustration.

{x0, x1} = {0.05, 0.2};
f[x_] := Sin[1/x];
fp = 
        Plot[f[x], {x, x0, x1}, DisplayFunction -> Identity];
fi = 
        IntervalPlot[f[x], {x, x0, x1}, Tolerance -> 0.1, 
          DisplayFunction -> Identity];
ri = fi[[1, 2]];
xx = 0.1016; fx = f[xx];
r0 = Cases[ri, Rectangle[{_, l_}, _] /; l > fx];
r1 = Complement[ri, r0];
Show[
        Graphics[{RGBColor[0.4, 0.8, 1], r0, RGBColor[0.2, 1, 0.6], r1}], fp, 
        Graphics[{{RGBColor[1, 0, 0], PointSize[0.03], 
              Point[{xx, fx}]}, {Dashing[{0.02, 0.02}], 
              Line[{{x0, fx}, {x1, fx}}]}}]];

Parallelization

The list of remaining rectangles and the best minimum so far are global values that need to be shared among all processors working on parts of the problem. The communication and synchronization overhead can be decreased by keeping local queues and synchronizing only occasionally.

The resulting code is quite sophisticated and is several pages long. Here is an example of its performance.

An example proposed by János Pintér

Here is the function of two variables whose global minimum we want to find.

f = 
        Function[{x, y}, 
          x^2 + y^2 + Log[1 + 20  Abs[Sin[x^2 + y^2]]] + 
            10  Abs[Sin[5  x^2 + 15  y^2]]];

Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, 
        PlotPoints -> 80];

This is the domain in which we search for the minimum.

domain = N[Interval /@ {{-4, 10}, {-4, 10}}]

Here is the precision with which we want to find the minimum.

tol = 1.0*^-6;

This command initializes the global shared variables and starts the parallel computation. When it finishes it reports performance statistics.

pq = priorityQueue[];
enQ[{#, f @@ #}] & /@ initials;
minmax = Infinity;  nbusy = Length[$Slaves];
With[{f = f, tol = tol, 
          ming = 2*Length[$Slaves] + 1}, 
fin = 
          RemoteEvaluate[
            driver[f, tol, ming]]
];

Here is the value of the minimum.

minmax

This is the collection of intervals in which the minimum is guaranteed to lie.

result

Show[rg, AspectRatio -> Automatic, 
        Frame -> True];

Example 4: Search for Reducible Polynomials

Introduction

If you want to verify that the polynomials (∑_(i = 1)^(n + 1))i (x^(i - 1)) for n = 1, 2, … are all irreducible (an open conjecture), factor the polynomials and then check that the length of the list of factors is 2, one factor being the overall numerical factor.

This computation will continue forever. To stop it, abort the local evaluation by pressing ⌘[.] or choosing Kernel ⊳ Abort Evaluation. After the abort, collect any waiting results, as follows.

Here is the definition of the polynomial in x with degree n.

poly[n_, x_] := (∑_(i = 1)^(n + 1))i (x^(i - 1))

Parallel Search

Here you make the definition known to all remote kernels.

ExportEnvironment[{poly}];

Now you can start the computation, requiring that it print each result as it continues. To stop the computation, abort it by choosing Kernel ⊳ Abort Evaluation or pressing ⌘[.].

d = 200;  ids = {}; results = {};
While[True,
While[$QueueLength == 
                0, 
AppendTo[ids, 
                Queue[{d}, {d, 
                    Length[FactorList[
                        poly[d, 
                          x]]]}]]; 
d++; 
If[! QueueRun[],  Break[]];
];
If[
              Length[ids] > 0, 
{res, id, ids} = 
                WaitOne[ids]; 
Print[
                res]; 
AppendTo[results, 
                res];
];
]

Cleanup

After the abort, you can either wait for the currently running calculations,

res0 = Wait[ids]

results = Join[results, res0];

or you can abort the remote kernels.

ResetSlaves[];

Unattended Calculation

The code can be extended in several areas to allow it to run unattended for longer times.

Results

A computation has been ongoing since December 2004 and an Orion DT-12. Polynomials up to degree 21 000 were tested, all unsuccessful, so far. Here is a graph of the factorization time for each degree.

Further Information

Resources on the Web

My own material on PCT: http://www.mathconsult.ch/math/pct/

Parallel Computing Toolkit Homepage: http://www.wolfram.com/products/applications/parallel/

Support Portal at Wolfram Research: http://support.wolfram.com/applicationpacks/parallel/

Conclusions

CloseSlaves[];

Mathematica and the Parallel Computing Toolkit offer a unique environment for the development, testing, and use of parallel programs. It brings parallel computation within reach of anybody with access to more than one processor.