Data Analytics and Simulation in Parallel with MATLAB*

1© 2016 The MathWorks, Inc.
Jos Martin
jos.martin@mathworks.com
Senior Engineering Manager – Parallel Computing
Parallel Computing with MATLAB and Simulink

2
Introduction
Ø  Why Parallel Computing
•  Need faster insight to bring competitive products to market quickly
•  Computing infrastructure is broadly available (Multicore Desktops, GPUs, Clusters)
Ø  With MathWorks Parallel Computing Tools
•  Leverage computational power of available hardware
•  Accelerate workflows with minimal to no code changes to your original code
•  Seamlessly scale from your desktop to clusters or on the cloud
•  Save engineering and research time and focus on insight

3
Agenda
Ø  Parallel Computing Paradigm
Ø  Task Parallelism
Ø  Data Parallelism
Ø  Summary

4
Agenda
Ø  Summary

5
Parallel Computing Paradigm
Multicore Desktops
Multicore Desktop
Core 5
Core 1 Core 2
Core 6MATLAB Desktop
(client)
Worker Worker
Worker Worker

6
Cluster Hardware
Cluster of computers
Core 5
Core 1 Core 2
Core 6
MATLAB Desktop
(client)
Core 5
Core 1 Core 2
Core 6
Core 5
Core 1 Core 2
Core 6 Core 5
Core 1 Core 2
Core 6
Worker Worker
Worker Worker
Worker Worker
Worker Worker
Worker Worker Worker Worker
Worker Worker Worker Worker

7
NVIDIA GPUs
Using NVIDIA GPUs
MATLAB Desktop
(client)
GPU cores
Device Memory

8
Cluster Computing Paradigm
§  Prototype on the desktop
§  Integrate with existing
infrastructure
§  Access directly through
MATLAB
User Desktop Headnode
Compute
Nodes
Parallel Computing Toolbox
MATLAB
MATLAB Distributed Computing Server

9
Parallel-enabled Toolboxes (MATLAB® Product Family)
Enable parallel computing support by setting a flag or preference
Optimization
Parallel estimation of
gradients
Statistics and Machine Learning
Resampling Methods, k-Means
clustering, GPU-enabled functions
Neural Networks
Deep Learning, Neural Network
training and simulation
Image Processing
Batch Image Processor, Block
Processing, GPU-enabled functions
Computer Vision
Parallel-enabled functions
in bag-of-words workflow
Signal Processing and
Communications
GPU-enabled FFT filtering,
cross correlation, BER
simulations
Other Parallel-enabled Toolboxes

10
Parallel-enabled Toolboxes (Simulink® Product Family)
Enable parallel computing support by setting a flag or preference
Simulink Control Design
Frequency response estimation
Simulink/Embedded Coder
Generating and building code
Simulink Design Optimization
Response optimization, sensitivity
analysis, parameter estimation
Communication Systems Toolbox
GPU-based System objects for
Simulation Acceleration
Other Parallel-enabled Toolboxes

11
Agenda
Ø  Summary

12
parfor
Definition
Code in a parfor loop is guaranteed by the programmer to be execution order
independent
Why is that important?
We can execute the iterates of the loop in any order, potentially at the same time on
many different workers.

13
parfor – how it works
§  Static analysis to deduce variable classification
–  Works out what data will need to be passed to which iterates
§  A loop from 1:N has N iterates which we partition into a number of intervals
–  Each interval will likely have a different number of iterates
§  Start allocating the intervals to execute on the workers
§  Stitch the results back together
–  Using functions from the static analysis

14
Variable Classification
reduce = 0; bcast = …; in = …;
parfor i = 1:N
temp = foo1(bcast, i);
out(i) = foo2(in(i), temp);
reduce = reduce + foo3(temp);
end

15
Loop variable
parfor i = 1:N
end

16
Making extra parallelism
§  No one loop appears to have enough iterations to go parallel effectively
–  Imagine if Na, Nb, and Nc are all reasonably small
for ii = 1:Na
for jj = 1:Nb
for kk = Nc
end
end
end
Na * Nb * Nc == quiteBigNumber
mergeLoopsDemo

17
Making extra parallelism
[N, fun] = mergeLoopRanges([Na Nb Nc]);
parfor xx = 1:N
[ii,jj,kk] = fun(xx);
doOriginalLoopCode
end

18
Sliced Variable
parfor i = 1:N
end

19
Broadcast variable
parfor i = 1:N
end

20
Reusing data
D = makeSomeBigData;
for ii = 1:N
parfor jj = 1:M
a(jj) = func(D, jj);
end
end

21
Reusing data
D = parallel.pool.Constant( @makeSomeBigData );
for ii = 1:N
parfor jj = 1:M
a(jj) = func(D.value, jj);
end
end
% Alternatively you can send data once and re-use it
oD = parallel.pool.Constant( someLargeData )

22
Common parallel program
set stuff going
while not all finished {
for next available result do something;
}

23
parfeval
§  Allows asynchronous programming
f = parfeval(@func, numOut, in1, in2, …)
§  The return f is a future which allows you to
–  Wait for the completion of calling func(in1, in2, …)
–  Get the result of that call
–  … do other useful parallel programming tasks …

24
Fetch Next
§  Fetch next available unread result from an array of futures.
[idx, out1, ...] = fetchNext(arrayOfFutures)
§  idx is the index of the future from which the result is fetched
§  Once a particular future has returned a result via fetchNext it will never do so again
–  That particular result is considered read, and will not be re-read

25
Common parallel program (MATLAB)
% Set stuff going
for ii = N:-1:1
fs(ii) = parfeval(@stuff, 1);
end
% While not all finished
for ii = 1:N
% For next available result
[whichOne, result] = fetchNext(fs);
doSomething(whichOne, result);
end
parfevalWaitbarDemo

26
parfevalOnAll
§  Frequently you want setup and teardown operations
–  which execute once on each worker in the pool, before and after the actual work
§  Execution order guarantee:
It is guaranteed that relative order of parfeval and parfevalOnAll as executed on the
client will be preserved on all the workers.

27
Agenda
Ø  Summary

28
Remote arrays in MATLAB
MATLAB provides array types for data that is not in “normal” memory
distributed array
(since R2006b)
Data lives in the combined memory of a cluster of
computers
gpuArray
(since R2010b)
Data lives in the memory of the GPU card
tall array
(since R2016b)
Data lives on disk, maybe spread across many disks
(distributed file-system)

29
Normal array – calculation happens in main memory:
Rule: take the calculation to where the data is
x = rand(...)
x_norm = (x – mean(x)) ./ std(x)

30
gpuArray – all calculation happens on the GPU:
x = gpuArray(...)
Rule: take the calculation to where the data is
distributed – calculation is spread across the cluster:
x = distributed(...)
tall – calculation is performed by stepping through files:
x = tall(...)

31
How big is big?
What does “Big Data” even mean?
“Any collection of data sets so large and complex that it becomes difficult to
process using … traditional data processing applications.”
(Wikipedia)
“Any collection of data sets so large that it becomes difficult to process using
traditional MATLAB functions, which assume all of the data is in memory.”
(MATLAB)

32
How big is big?
The Large Hadron Collider reached peak
performance on 29 June 2016
§  2076 bunches of 120 billion protons currently
circulating in each direction
§  ~1.6x1014 collisions per week, >30 petabytes of
data per year
§  too big to even store in one place
§  used to explore interesting science, but taking
researchers a long time to get through
In 1085 William 1st commissioned a survey
of England
§  ~2 million words and figures collected over two
years
§  too big to handle in one piece
§  collected and summarized in regional pieces
§  used to generate revenue (tax), but most of the
data then sat unused
Image courtesy of CERN.
Copyright 2011 CERN.

33
Tall arrays (new R2016b)
§  MATLAB data-type for data that doesn’t fit into memory
§  Ideal for lots of observations, few variables (hence “tall”)
§  Looks like a normal MATLAB array
–  Supports numeric types, tables, datetimes, categoricals, strings, etc…
–  Basic maths, stats, indexing, etc.
–  Statistics and Machine Learning Toolbox support (clustering, classification, etc.)

34
Cluster of
Machines
Memory
Single
Machine
Memory
§  Data is in one or more files
§  Typically tabular data
§  Files stacked vertically
§  Data doesn’t fit into memory
(even cluster memory)

35
Cluster of
Machines
Memory
Single
Machine
Memory
§  Use datastore to define files
§  Allows access to small pieces of
data that fit in memory.
Datastore

ds = datastore('*.csv')
while hasdata(ds)
piece = read(ds);
% Process piece
end

36
tall
array

Cluster of
Machines
Memory
Single
Machine
Memory
§  Create tall table from datastore
§  Operate on whole tall table just like
ordinary table
Datastore

ds = datastore('*.csv')
tt = tall(ds)
summary(tt)
max(tt.EndTime – tt.StartTime)
Single
Machine
MemoryProcess

37
tall
array

Cluster of
Machines
Memory
Single
Machine
Memory
§  With Parallel Computing Toolbox,
process several pieces at once
Datastore

Single
Machine
MemoryProcess

Single
Machine
MemoryProcess

38
Example
New York taxi fares (150,000,000 rows (~25GB) per year)
>>
dataLocation
=
'hdfs://hadoop01glnxa64:54310/datasets/nyctaxi/';

>>
ds
=
datastore(
fullfile(dataLocation,
'yellow_tripdata_2015-‐*.csv')
);

>>
tt
=
tall(ds)

tt
=

M×19
tall
table

VendorID

tpep_pickup_datetime

tpep_dropoff_datetime

passenger_count

trip_distance

pickup_longitude

pic

________

_____________________

_____________________

_______________

_____________

________________

___

2

'2015-‐01-‐15
19:05:39'

'2015-‐01-‐15
19:23:42'

1

1.59

-‐73.994

40

1

'2015-‐01-‐10
20:33:38'

'2015-‐01-‐10
20:53:28'

1

3.3

-‐74.002

40.

1

'2015-‐01-‐10
20:33:38'

'2015-‐01-‐10
20:43:41'

1

1.8

-‐73.963

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:35:31'

1

0.5

-‐74.009

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:52:58'

1

3

-‐73.971

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:53:52'

1

9

-‐73.874

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:58:31'

1

2.2

-‐73.983

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:42:20'

3

0.8

-‐74.003

40.

:

:

:

:

:

:

:

:

:

:

:

:

:

:

Input data is tabular –
result is a tall table

39
>>
dataLocation
=
'hdfs://hadoop01glnxa64:54310/datasets/nyctaxi/';

>>
ds
=
datastore(
fullfile(dataLocation,
'yellow_tripdata_2015-‐*.csv')
);

>>
tt
=
tall(ds)

tt
=

M×19
tall
table

VendorID

tpep_pickup_datetime

tpep_dropoff_datetime

passenger_count

trip_distance

pickup_longitude

pic

________

_____________________

_____________________

_______________

_____________

________________

___

2

'2015-‐01-‐15
19:05:39'

'2015-‐01-‐15
19:23:42'

1

1.59

-‐73.994

40

1

'2015-‐01-‐10
20:33:38'

'2015-‐01-‐10
20:53:28'

1

3.3

-‐74.002

40.

1

'2015-‐01-‐10
20:33:38'

'2015-‐01-‐10
20:43:41'

1

1.8

-‐73.963

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:35:31'

1

0.5

-‐74.009

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:52:58'

1

3

-‐73.971

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:53:52'

1

9

-‐73.874

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:58:31'

1

2.2

-‐73.983

40.

1

'2015-‐01-‐10
20:33:39'

'2015-‐01-‐10
20:42:20'

3

0.8

-‐74.003

40.

:

:

:

:

:

:

:

:

:

:

:

:

:

:

Example
New York taxi fares (150,000,000 rows (~25GB) per year)
Number of rows is
unknown until all the
data has been read
Only the first few
rows are displayed

40
Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
0%
complete

Evaluation
0%
complete

Example
Once created, can process much like an ordinary table
%
Remove
some
bad
data

tt.trip_minutes
=
minutes(tt.tpep_dropoff_datetime
-‐
tt.tpep_pickup_datetime);

tt.speed_mph
=
tt.trip_distance
./
(tt.trip_minutes
./
60);

ignore
=
tt.trip_minutes
<=
1
|
...

%
really
short

tt.trip_minutes
>=
60
*
12
|
...

%
unfeasibly
long

tt.trip_distance
<=
1
|
...

%
really
short

tt.trip_distance
>=
12
*
55
|
...

%
unfeasibly
far

tt.speed_mph
>
55
|
...

%
unfeasibly
fast

tt.total_amount
<
0
|
...

%
negative
fares?!

tt.total_amount
>
10000;

%
unfeasibly
large
fares

tt(ignore,
:)
=
[];

%
Credit
card
payments
have
the
most
accurate
tip
data

keep
=
tt.payment_type
==
{'Credit
card'};

tt
=
tt(keep,:);

%
Show
tip
distribution

histogram(
tt.tip_amount,
0:25
)

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
16%
complete

Evaluation
14%
complete

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
33%
complete

Evaluation
30%
complete

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
50%
complete

Evaluation
45%
complete

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
66%
complete

Evaluation
63%
complete

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
83%
complete

Evaluation
75%
complete

Evaluating
tall
expression
using
the
Parallel
Pool
'local':

-‐
Pass
1
of
1:
Completed
in
4.9667
min

Evaluation
completed
in
5
min

Data only read once,
despite 21 operations

41
Scaling up
If you just have MATLAB:
§  Run through each ‘chunk’ of data one by one
If you also have Parallel Computing Toolbox:
§  Use all local cores to process several ‘chunks’ at once
If you also have a cluster with MATLAB Distributed
Computing Server (MDCS):
§  Use the whole cluster to process many ‘chunks’ at once

42
Scaling up
Working with clusters from MATLAB desktop:
§  General purpose MATLAB cluster
–  Can co-exist with other MATLAB workloads (parfor, parfeval,
spmd, jobs and tasks, distributed arrays, …)
–  Uses local memory and file caches on workers for efficiency
§  Spark-enabled Hadoop clusters
–  Data in HDFS
–  Calculation is scheduled to be near data
–  Uses Spark’s built-in memory and disk caching

43
Distributed Arrays
§  Unlike tall these need to fit in cluster memory
§  Full range of numerical algorithms
–  Linear Algebra, eigenvalues, etc
§  Sparse support
§  Want to write your own function?
–  Full access to local part of the data
–  Ability to rebuild distributed array from local parts

44
Standard Benchmarks
HPL
spmd
A = rand(10000, codistributor2dbc);
b = rand(10000, 1, codistributor2dbc);
x = Ab;
end
FFT
D = rand(1e8, 1, 'distributed');
F = fft(D);
STREAM Triad
B = rand(1e8, 1, 'distributed');
C = rand(1e8, 1, 'distributed');
q = rand;
A = B + q*C;
Create 2D Block Cyclic
distribution appropriate for
matrix solve

45
Single Program, Multiple Data (spmd)
§  Everyone executes the same program
–  Just with different data
–  Inter-lab communication library enabled (MPI)
–  Call labindex and numlabs to distinguish labs
§  Example
x = 1
spmd
y = x + labindex;
end

46
Variable types can change across spmd
x = 1;
assert( isa(x, 'double') )
spmd
assert( isa(x, 'double') )
y = labindex;
end
assert( isa(y 'Composite') )
y{1} == 1; % TRUE
spmd
assert( isa(y, 'double') )
end
Ordinary types are broadcast to all labs
Returned ordinary type is referenced
by a Composite
Composite can be dereferenced
on the client
Composite becomes the contained
ordinary type on a lab

47
Want to write an MPI Program?
Send / receive data between labs
labSend
labReceive
labSendReceive
labProbe
Who am I? Where am I?
labindex
numlabs
Global Operations across labs
gplus
gcat
gop
labBarrier
labBroadcast

48
RandomAccess – an interesting MPI program
spmd
% Initialize random number stream on each worker
randRA( (labindex-1) * m * 4 / numlabs, 'StreamOffset' );
t1 = tic;
for k = 1:nloops
% Make the local chunk of random data
list = randRA( b );
% Loop over the hyper-cube dimensions
for d = 0:logNumlabs-1
% Choose my partner for this dimension of the hypercube
partner = 1 + bitxor( (labindex-1), 2.^d );
% Choose my mask for this dimension of the hypercube
dim_mask = uint64( 2.^( d + logLocalSize ) );
% Choose which data to send and receive for this dimension
dataToSend = logical( bitand( list, dim_mask ) );
if partner <= labindex
dataToSend = ~dataToSend;
end
% Setup a list of data that will be sent, and list I will keep
send_list = list( dataToSend );
keep_list = list( ~dataToSend );
% Use send/receive to get some data that we should use next round
recv_list = labSendReceive( partner, partner, send_list );
% Our new list is the old list and what we've received
list = [keep_list, recv_list];
end
% Finally, after all rounds of communication, perform the table updates.
idx = 1 + double( bitand( localMask, list ) );
T(idx) = bitxor( T(idx), list );
end
% Calculate max time
t = gop( @max, toc( t1 ) );
end

49
RandomAccess – an interesting MPI program
spmd
t1 = tic;
for k = 1:nloops
% Make the local chunk of random data
list = randRA( b );
% Loop over the hyper-cube dimensions
for d = myDims(0:logNumlabs-1, labindex)
[send_list, keep_list] = partitionData(list, d);
% Use send/receive to transfer data
recv_list = labSendReceive(to(d), from(d), send_list);
% Our new list is the old list and what we've received
list = [keep_list, recv_list];
end
% Finally, perform the table updates.
T = updateTable(list);
end
% Calculate max time
t = gop( @max, toc( t1 ) );
end

50
Summary
§  Simple array types for data analytics and PGAS programming
§  Extensive parallel language suited to all types of problem
§  All runs in both interactive and batch (off-line) mode

Data Analytics and Simulation in Parallel with MATLAB*

Recommended

More Related Content

What's hot (20)

Similar to Data Analytics and Simulation in Parallel with MATLAB* (20)

More from Intel® Software (20)

Recently uploaded (20)

Data Analytics and Simulation in Parallel with MATLAB*