Transcript Robust High Performance Optimization for

Robust High Performance
Optimization for Clustering,
Multi-Dimensional Scaling
and Mixture Models
CGB Indiana University Lunchtime Talk
January 22 2008
Geoffrey Fox discussing work with:
Seung-Hee Bae, Neil Devadasan, Rajarshi Guha, Marlon Pierce,
Xiaohong Qiu, David Wild, Huapeng Yuan
Community Grids Laboratory, Research Computing UITS,
School of informatics and POLIS Center Indiana University
George Chrysanthakopoulos, Henrik Frystyk Nielsen
Microsoft Research, Redmond WA
http://www.infomall.org/multicore
[email protected], http://www.infomall.org
1
Abstract








We first review the pros and cons of various approaches to non linear
optimization in the presence of local minima, ill conditioned matrices and
ambiguous choice of appropriate number of degrees of freedom (over and
under fitting).
We define constraints on approaches from need to run well in parallel on
systems of multicore CPU's.
We present a uniform approach to data clustering and Gaussian mixture model
ling that uses deterministic (not Monte Carlo) annealing to mitigate the local
minima problem and naturally relates the appropriate number of parameters
(clusters or mixture components) to the scale at which problem is examined.
New clusters (mixtures) are introduced at phase transitions as the annealing
temperature is lowered and second derivative matrix becomes singular.
We contrast three ways of visualizing this structure in low (2) dimensions with
Principal Component Analysis PCA, Generative Topographic Mapping GTM
and Multi-Dimensional Scaling MDS using annealing to regularize GTM and
MDS.
Currently we have implemented in preliminary fashion deterministic annealing
clustering and GTM in a fashion that runs well on multicore systems.
We have applied these techniques to Geographical Information Systems
(clustering demographic data in 2D) and Cheminformatics in 1052 and lower
dimensions.
We would like to understand other applications that can constrain and test
these techniques.
Too much Computing?



Historically both grids and parallel computing have tried to
increase computing capabilities by
• Optimizing performance of codes at cost of re-usability
• Exploiting all possible CPU’s such as Graphics coprocessors and “idle cycles” (across administrative
domains)
• Linking central computers together such as NSF/DoE/DoD
supercomputer networks without clear user requirements
Next Crisis in technology area will be the opposite problem –
commodity chips will be 32-128way parallel in 5 years time
and we currently have no idea how to use them on commodity
systems – especially on clients
• Only 2 releases of standard software (e.g. Office) in this
time span so need solutions that can be implemented in
next 3-5 years
Intel RMS analysis: Gaming and Generalized decision
support (data mining) are ways of using these cycles
Intel’s Projection
Too much Data to the Rescue?




Multicore servers have clear “universal parallelism” as many
users can access and use machines simultaneously
Maybe also need application parallelism (e.g. datamining) as
needed on client machines
Over next years, we will be submerged of course in data
deluge
• Scientific observations for e-Science
• Local (video, environmental) sensors
• Data fetched from Internet defining users interests
Maybe data-mining of this “too much data” will use up the
“too much computing” both for science and commodity PC’s
• PC will use this data(-mining) to be intelligent user
assistant?
• Must have highly parallel algorithms
Intel’s Application Stack
Total
Asian
Hispanic
Renters
30 Clusters
30 Clusters
10 Clusters
Clustering algorithm annealing by decreasing distance scale and gradually finds more cluster
resolution improved
Here we see 10 clusters increasing to 30 as algorithm progresses
Parallel Multicore
Deterministic Annealing Clustering
Parallel Overhead
on 8 Threads on Intel 8 core
0.45
10 Clusters
0.4
Speedup = 8/(1+Overhead)
= 7.6 if Overhead is 0.05
0.35
Overhead = Constant1 + Constant2/n
0.3
Constant1 = 0.05 to 0.1 (Windows)
due to thread
runtime fluctuations
0.25
0.2
20 Clusters
0.15
0.1
0.05
10000/(Grain Size n = points per core)
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Basic Nonlinear Optimization I




Nonlinear optimization (in either 2 or maximum likilihood) is
subject to well known problems including diverging or converging
but to non optimal solution (local minima)
If we have n parameters i and need to minimize
F = 2 or – ln L
Then naively one expands
F = F(0) + k=1n k F(0) /i +  k,l=1n k l 2F(0) /kl
• with  =  - 0



As long as number of observables N is large in sum
2 = i=1N qi2() or – ln L = i=1N ln qi()
One can approximate
2F(0) /kl  i=1N (1/ qi)qi(0) /k qi(0) /l
dropping term proportional to 2qi(0) /k l
This shows that second derivative matrix is positive (semi-)
definite
Basic Nonlinear Optimization II

Taylor expansion of F (Newton’s method) is just
F ( )  F ( 0 )  θ b  0.5 Δθ A Δθ
T
T
where b is first derivative and A second derivative
with minima given by
1
Δ   A b


This formula can rarely be used as in practice matrix A is
poorly conditioned corresponding to linear combination of 
being (ill/un) determined
This easy to understand by transforming basis to those of
eigenvectors of A which is a symmetric positive (semi) definite
matrix
Basic Nonlinear Optimization III
0
….
1
Suppose A  0 2
…….
0 0
0
0 with 1  2  ......  n
n
Then if b refered to eigenbasis , we get
 k  bk / k and when k small, then  k is unreliable as neglected terms
in expansion are larger tha n included second order term s.


In my experience, the lowest eigenvalues k of A are zero to
rounding error, so corresponding shifts k are unreliable and
incrementing  = 0 +  gives a nonsense point and F rapidly
diverges as one iterates
There are several approaches and they actually work and with
modification, Newton’s method is very reliable and one can
expect to get good results
• Unless n is very large (tens of thousands?) when it can be
computationally unattractive
Reliable Newton’s Method I

WNPF (Well Nigh Perfect Fit – code lost): Here one
diagonalizes A and modifies equation k = bk/k, by
zeroing shifts k if k < cut where is cut is dynamically
adjusted and I typically started with cut =  max with
initially  ~ 0.1
•  adjusted up or down depending on success of fit
• Only shift in well determined directions but do many
directions at once

Converges in 10-20 iterations
See my first paper
65 Citations
Reliable Newton’s Method II

Manxcat.f (Code still exists): This can be too expensive
so when I needed millions of such optimizations I used
Marquardt’s method which corresponds to
• k  k + Q accomplished by
• Replacing A  Q  Identity matrix + A
• so removes effect of 1/k for small k without expense of
determining eigenvalues/vectors
• Has nonzero but small shifts in ill determined directions
• Q is dynamically changed in a complicated fashion on a
logarithmic scale


Similar performance to WNPF but less flexible
Code ran million optimizations without aborting once!
Gradient Descent











Many optimizations do not use second order Taylor expansion
but rather look at a shift
 =  G(0) so estimate of minimum is
 = 0 +  G(0)
Where  may be varied and G has a more or less fancy way of determining it
However usually G is related to steepest descent direction
i.e. G = - F which is direction F decreases fastest along
One can often find a rigorous bound so that one can show F definitely
decreases
However essentially always a local minima as only explore directions of
largest eigenvalues (still could be a good minima)
Note unlikely to be a minimum in “other directions” so really a “partial” not
a “local” minima
Convergence tends to be slow as error lower order in F than for Newton’s
method but derivatives do not need to be calculated
Computation  n N for n parameters and N data points
Newton’s method computation  n2 N (compute A) plus n3 (Solve for )
Simulated Annealing





Here we consider function F() to be minimized as an
Energy E() and randomly choose points in a
neighborhood of 0
One accepts a new value of  with a probability P
depending on temperature T and energies E() and
E(0)
Typical choice is P=1 if E() < E(0), and
P = exp([E(0) − E()]/T) if E() > E(0)
Temperatures are decreased gradually
This method avoids local minima due to allowing
upward fluctuations at high temperature but execution
time is huge although each iteration has time
proportional to n N
CICC Chemical Informatics and Cyberinfrastructure
Collaboratory Web Service Infrastructure
Cheminformatics Services
Statistics Services
Database Services
Core functionality
Fingerprints
Similarity
Descriptors
2D diagrams
File format conversion
Computation functionality
Regression
Classification
Clustering
Sampling distributions
3D structures by
CID
SMARTS
3D Similarity
Docking scores/poses by
CID
SMARTS
Protein
Docking scores
Applications
Applications
Docking
Predictive models
Filtering
Feature selection
Druglikeness
2D plots
Toxicity predictions
Arbitrary R code (PkCell)
Mutagenecity predictions
PubChem related data by
Anti-cancer activity predictions
Need to make
Pharmacokinetic parameters
CID, SMARTS
all this parallel
OSCAR Document Analysis
InChI Generation/Search
Computational Chemistry (Gamess, Jaguar etc.)
Core Grid Services
Service Registry
Job Submission and Management
Local Clusters
IU Big Red, TeraGrid, Open Science Grid
Varuna.net
Quantum Chemistry
Portal Services
RSS Feeds
User Profiles
Collaboration as in Sakai
Deterministic Annealing for Data Mining






We are looking at deterministic annealing algorithms because although
heuristic
• They have clear scalable parallelism (e.g. use parallel BLAS)
• They avoid (some) local minima and regularize ill defined problems
in an intuitively clear fashion
• They are fast (no Monte Carlo)
• They determine number of clusters automatically (Most important?)
• Support iterative visualization approaches as increment number of
clusters systematically
Developed first by Durbin as Elastic Net for TSP
Extended by Rose (my student then; now at UCSB)) and Gurewitz
(visitor to C3P) at Caltech for signal processing and applied later to
many optimization and supervised and unsupervised learning methods.
See K. Rose, "Deterministic Annealing for Clustering, Compression,
Classification, Regression, and Related Optimization Problems,"
Proceedings of the IEEE, vol. 80, pp. 2210-2239, November 1998 Cited
by 286
A deterministic annealing approach to clustering K Rose, E Gurewwitz,
G Fox - Pattern Recognition Letters, 1990 Cited by 159
Statistical mechanics and phase transitions in clustering K Rose, E
Gurewitz, GC Fox - Physical Review Letters, 1990 - Cited by 276
High Level Theory





Deterministic Annealing can be looked at from a Physics,
Statistics and/or Information theoretic point of view
In physics view, you get result by using “mean field
approximation” but simplest is information theoretic
Let E be function that one wishes to minimize
Then we wish to find minima that is most probable given density
of states.
One minimizes Free Energy F= E-TS where S is Entropy and S is
just sum over plnp for each degree of freedom with probability p
• At large T, Entropy dominates while at small T Energy dominates
• Annealing lowers temperature so solution tracks continuously


Key idea is to only apply this to “missing data” i.e. to hidden
degree of freedom k that labels cluster or mixture
So S is a sum over these discrete degrees of freedom and one can
minimize F to find the form of p for given energy function
Deterministic Annealing for Clustering I
N Points X i and K Cluster Centers Y k
Pr( X i  Ck )  exp(  E ( X i , Y k ) / T ) / Z ( X i ) where
Z ( X i )  k 1 exp(  E ( X i , Y k ) / T )
K
E ( X i , Y k )  ( X i  Y k )2
Free Energy F  T i1 ln Z ( X i ))
N
Compare Simple Gaussian Mixture (K centers) with
Z ( X i )  k 1 Pk exp(  E ( X i , Y k ) /( 2 k ))
K
2
Deterministic Annealing for Clustering II
with Pr( X i  Ck )  exp(  E ( X i , Y k ) / T ) / Z ( X i , Y k )
old
Calculate Y k
new

N
i 1
X i Pr( X i  Ck )
old

N
i 1






Pr( X i  Ck )
This is an extended K-means algorithm or a simplified Gaussian
mixture
Start with a single cluster giving as solution Y1 as centroid
For some annealing schedule for T, iterate above algorithm
testing correlation matrix in Xi about each cluster center to see if
“elongated”
Split cluster if elongation “long enough”; splitting is a phase
transition in physics view
You do not need to assume number of clusters but rather a final
resolution T or equivalent
At T=0, uninteresting solution is N clusters; one at each point xi
Deterministic
Annealing
F({y}, T)
Solve Linear
Equations for
each temperature
Nonlinearity
effects mitigated
by approximating
with solution at
previous higher
temperature


Configuration {y}
Minimum evolving as temperature decreases
Movement at fixed temperature going to local
minima if not initialized “correctly
Views from
Past on
Physical
Optimization
You need to Specify a minimum temperature
Parallelism


All the methods – Deterministic Annealing, Mixture
Models, GTM, PCA and Newton’s method involve
Sums over data points x=1..N to calculate values and
derivatives
• One can divide points up between cores and efficiently
parallelize
• A little tricky to add results from separate data sets running
in parallel (problem is Cache coherence)

Matrix algebra such as finding eigenvalues
• Well known how to parallelize but need low latency i.e.
parallel not distributed system
Where are we for Clustering?


We have deterministically annealed clustering running well on 8core (2-processor quad core) Intel systems using C# and
Microsoft Robotics Studio CCR/DSS
Could also run on multicore-based parallel machines but didn’t
do this (is there a large Windows quad core cluster on
TeraGrid?)
• This would also be efficient on large problems

Applied to Geographical Information Systems (GIS) and census
data
• Could be an interesting application on future broadly deployed PC’s
• Visualize nicely on Google Maps (and presumably Microsoft Virtual Earth)


Applied to several Cheminformatics problems and have parallel
efficiency but visualization harder as in 150-1024 (or more)
dimensions
Will develop a family of such parallel (annealing) data-mining
tools where basic approach known for
• Clustering and Gaussian Mixtures (Expectation Maximization)
• Mapping High Dimensional Spaces – MDS and GTM
• Hidden Markov Methods
Clustering Data



Cheminformatics was tested successfully with small datasets and
compared to commercial tools
Cluster on properties of chemicals from high throughput
screening results to chemical properties (structure, molecular
weight etc.)
Applying to PubChem (and commercial databases) that have 620 million compounds
• Comparing traditional fingerprint (binary properties) with real-valued
properties

GIS uses publicly available Census data; in particular the 2000
Census aggregated in 200,000 Census Blocks covering Indiana
• 100MB of data

Initial clustering done on simple attributes given in this data
• Total population and number of Asian, Hispanic and Renters

Working with POLIS Center at Indianapolis on clustering of
SAVI (Social Assets and Vulnerabilities Indicators) attributes at
http://www.savi.org) for community and decision makers
• Economy, Loans, Crime, Religion etc.
In detail, different groups have
different cluster centers
Age < 5 v. Age > 75
Age < 5 v. Age 25-34
Parallel Multicore
Deterministic Annealing Clustering
Parallel Overhead for large (2M points) Indiana Census clustering
on 8 Threads Intel 8b
This fluctuating overhead due to 5-10% runtime fluctuations between threads
0.250
0.200
overhead
“Constant1”
0.150
0.100
0.050
Increasing number of clusters decreases
communication/memory bandwidth overheads
0.000
0
5
10
15
20
#cluster
25
30
35
Parallel Multicore
Deterministic Annealing Clustering
0.200
Parallel Overhead for subset of PubChem clustering on 8 Threads
(Intel 8b)
0.180
“Constant1”
The fluctuating overhead
is reduced to 2% (as bits not doubles)
40,000 points with 1052 binary properties
(Census is 2 real valued properties)
0.160
overhead
0.140
0.120
0.100
0.080
0.060
0.040
Increasing number of clusters decreases
communication/memory bandwidth overheads
0.020
0.000
0
2
4
6
8
10
#cluster
12
14
16
18
Visualizing High Dimensional Spaces




For GIS we have a 2D or 3D underlying space and visualization
of data correlations is clear
For cheminformatics we have hundreds (continuous) to
thousands (binary) degrees of freedom and difficult to evaluate
any data mining – in particular clustering
Need to map high dimensional space into lower (here two)
dimensional space with some constraints
This field has been well studied and there are at least two
approaches
• SOM (Self Organizing Maps) and GTM (Generative Topographic
Mapping) which are really clustering algorithms with a built in 2D
organization for clusters
• MDS (Multi Dimensional Scaling) which “just views this as an
optimization problem”.
• Principal Component Analysis (PCA) is here viewed as a special case of
MDS
MDS Multi Dimensional Scaling






This is rather straightforward (and perhaps good for that
reason).
Consider n points Y which are vectors in a high dimensional
space that we want to map into n points X in a low dimensional
space (bad notation)
Let us minimize:
Stress (X) = i<j=1n weight(i,j) (ij - d(Xi , Xj))2
ij is distance between original vectors Y in high dimensional
space but it can be any scoring of discrepancy between points i
and j.
d(Xi , Xj) is distance between mapped vectors
Here simplest choice is weight(i,j) = 1 but one can also look at
choices like Sammon’s mapping
weight(i,j) = 1/ij which emphasizes small distances ij
PCA Principal Component Analysis I

This is closely related to analysis used in WNPF approach to
nonlinear fitting and DA method for determining whether to
split clusters. We form DD matrix PCA which depends on
vectors Yk in original space and a “mean” M about which we do
the expansion
n
PCA   K (Y k , M )(Y k  M )(Y k  M )
T
k 1


For cluster splitting the mean M is a cluster center and there is
a kernel which is probability that Yk belongs to a given cluster
For MDS, the kernel is often absent and the mean M is the
centroid of points Yk in original space
PCA Principal Component Analysis II


The PCA approach finds the eigenvalues and
eigenvectors of the DD matrix PCA and maps point
point Yk into its expansion in the first 2 (or choose here
lower dimensional space) PCA eigenvectors with
largest eigenvalues
This is optimal solution to minimizing stress under two
conditions
• One only looks at mappings Yk  Xk that are orthogonal
projections
• Choice of stress as i<j=1n (ij2- d(Xi , Xj)2)


For example there are “better” linear transformations
that add scaling of PCA vectors
Further PCA clearly weights large discrepancies most
Classical MDS: SMACOF
Scaling by minimizing a complicated function

This is a variant of steepest descent approach to minimizing
Stress. One sets two matrices V and B in nL dimensional space
(L=2 is dimension of target space of mapping).
• n is number of vectors Yk to be mapped
Vi , j  weight (i, j ) if i  j and Vi ,i  i j Vi , j
Bi , j ( Z )  weight (i, j )  δij /d( Z i ,Z j ) if i  j
and Bi ,i  i j Bi , j

Then one can establish a rigorous bound to show that one always
decreases Stress (X) by iteration labeled by t if
One iterates the Guttman Transform
X

[t ]
1
 V B( X
[ t 1]
)X
[ t 1]
This presumably will always find local minima?
More obvious MDS for Clustering or
Mixture Model Mapping



The number of parameters in MDS is typically modest as in
simple cases just 2 (dimension of target space)  Number of
Clusters (or Number of points to be mapped)
Thus Newton’s method should work well and it is trivial to
derive first and second order derivatives for (X)
As well as “general” approaches to regularization, one can
enhance convergence by
• Initializing optimization with SMACOF iterations
• Noting that DA clustering algorithm gives us results with
number of clusters increasing by one each time. Thus we
can initialize optimization for n clusters with results
from MDS with n-1 clusters (DA also tells you which
cluster center split)
GTM Generative Topographic Mapping




This addresses a slightly different problem of mapping ALL
points in a high dimensional space. It approaches this by
simultaneously clustering in high dimensional space and
mapping cluster (centers) to low dimension space
• We ignore clustering and view as a way to getting mapping by
looking at all points in space not just cluster centers
The key idea is a nonlinear mapping
• Y(k) = m=1M Wmm(X(k))
• m(X) = exp( - 0.5 (X-m)2/2 ) for fixed m and 
Y(k) are cluster centers in high dimensional space which map
from X(k) in lower dimensional space
If  largish compared to separation between basis vectors X(k)
then topology of Y(k) and X(k) match but distances will not
GTM for visualizing 2 Clusters in 155 Dimensions








Here GTM just
used to visualize.
Clustering done
separately
Deliberately “easy”
problem!
Note although
formal clustering
gives 2 clusters
GTM visualization
uses
K=225 clusters
M=64 basis
functions
N=335 points
Basic Gaussian Mixture Models

Mixture models write the probability of a point X(x) as
K
p( x, )   Pk exp( 0.5( X ( x)  Y (k )) 2 /  (k ) 2 ) /( 2 (k ) 2 ) D / 2
k 1



And set the Liklihood L = x=1N p(x,) where parameters  are centers Y(k),
Probability Pk that point generated by center k and covariance matrix (k)2
Expectation Maximization is (from 1977) iterative solution to minimizing –
lnL
It starts
Prob(k | X ( x), ) 
by estimating
2
2
2 D/2
Pk exp( 0.5( X ( x)  Y (k )) /  (k ) ) /( 2 (k ) )





/ p ( x,  )
And using this to calculate all
N
Y
(
k
)

X ( x) Prob(k | X ( x), )
the  using formulae like:
x1
Note similarity to deterministic annealing cluster formulae

Deterministic annealing EM algorithm - N Ueda, R Nakano - Neural Networks, 1998
Cited by 124
A different approach which does not seem general is in
Deterministic annealing for density estimation by multivariate normal mixtures
M Kloppenburg, P Tavan - Physical Review E, 1997 Cited by 18

DA
/
EM
/
DAEM
/
GTM
I
The General Formula for “F = - lnL” or “Free Energy” F = E-TS is:
N
F  T  a ( x) ln{ k 1 g (k ) exp[ 0.5( E ( x)  Y (k )) 2 / (Ts(k ))]
K
x 1



Note x is a label and E(x) is data
For Deterministic Annealing
• a(x) = 1/N or generally p(x) with  p(x) =1
• g(k)=1
• s(k)=0.5
• T is annealing temperature varied down from  with final value of 1
• Vary Y(k) but can calculate Pk and (k) (even for matrix (k)) using IDENTICAL
formulae for Gaussian mixtures; iteration formula identical to Expectation
Maximization EM which is steepest descent
• K starts at 1 and is incremented by algorithm
For traditional Gaussian mixture models simplified to spherical distributions ((k) is
really a k  k symmetric correlation matrix)
• a(x) = 1
• g(k) = Pk/(2(k)2)D/2 where space D dimensional
• s(k) = (k)2
• T=1
• Vary Y(k) Pk and (k) but fix value of K a priori
DA / EM / DAEM / GTM II

The General Formula for “F = - lnL-TS” is:
N
F  T  a( x) ln{ k 1 g (k ) exp[ 0.5( E ( x)  Y (k )) 2 /(Ts (k ))]
K
x 1


Note x is a label and E(x) is data
For Deterministic Annealing Gaussian mixture models simplified
to spherical distributions ((k) is really a DD symmetric
correlation matrix)
• a(x) = 1
• g(k)={Pk/(2(k)2)D/2}1/T
• s(k)= (k)2
• T is annealing temperature varied down from  with final
value of 1
• Vary Y(k) Pk and (k)
• K starts at 1 and is incremented by algorithm
DA / EM / DAEM / GTM III

The General Formula for “F = - lnL” is:
N
F  T  a( x) ln{ k 1 g (k ) exp[ 0.5( X ( x)  Y (k )) 2 /(Ts (k ))]
K
x 1
Note x is a label and E(x) is data
 For GTM Generative Topographic Mapping
• a(x) = 1
• g(k) = (1/K)( /2)D/2 where space D dimensional
• s(k) = 1/ 
• T=1
• Y(k) = m=1M Wmm(X(k))
• Fix m(X) = exp( - 0.5 (X-m)2/2 ) but other choices possible
There is presumably a version of GTM using deterministic annealing
• Vary Wm and  but fix values of M and K a priori
using either a pure
cluster basis – anneal in  or
 Y(k) E(x) Wm are vectors in original high D dimensional space
using
an annealed Gaussian mixture with temperature 1/ plus intrinsic
 X(k) and m are vectors in 2 dimensional mapped space
widths

Some links for Parallel Computing

See http://www.connotea.org/user/crmc for references - select tag oldies for venerable links; tags like MPI
Applications Compiler have obvious significance
http://www.infomall.org/salsa for recent work
including publications
My tutorial on parallel computing

http://grids.ucs.indiana.edu/ptliupages/presentations/PC2007/index.html


Parallel Programming Model







If multicore technology is to succeed, mere mortals must be able to
build effective parallel programs on commodity machines
There are interesting new developments – especially the new Darpa
HPCS Languages X10, Chapel and Fortress
However if mortals are to program the 64-256 core chips expected in 5-7
years, then we must use near term technology and we must make it easy
• This rules out radical new approaches such as new languages
Remember that the important applications are not scientific computing
but most of the algorithms needed are similar to those explored in
scientific parallel computing
We can divide problem into two parts:
• “Micro-parallelism”: High Performance scalable (in number of
cores) parallel kernels or libraries
• Macro-parallelism: Composition of kernels into complete
applications
We currently assume that the kernels of the scalable parallel
algorithms/applications/libraries will be built by experts with a
Broader group of programmers (mere mortals) composing library
members into complete applications.
Multicore SALSA at CGL




Service Aggregated Linked Sequential Activities
Aims to link parallel and distributed (Grid) computing by
developing parallel applications as services and not as
programs or libraries
• Improve traditionally poor parallel programming
development environments
Developing set of services (library) of multicore parallel data
mining algorithms
Looking at Intel list of algorithms (and all previous experience),
we find there are two styles of “micro-parallelism”
• Dynamic search as in integer programming, Hidden Markov Methods
(and computer chess); irregular synchronization with dynamic threads
• “MPI Style” i.e. several threads running typically in SPMD (Single
Program Multiple Data); collective synchronization of all threads together

Most Intel RMS are “MPI Style” and very close to scientific
algorithms even if applications are not science
Scalable Parallel Components






How do we implement micro-parallelism?
There are no agreed high-level programming environments for
building library members that are broadly applicable.
However lower level approaches where experts define
parallelism explicitly are available and have clear performance
models.
These include MPI for messaging or just locks within a single
shared memory.
There are several patterns to support here including the
collective synchronization of MPI, dynamic irregular thread
parallelism needed in search algorithms, and more specialized
cases like discrete event simulation.
We use Microsoft CCR
http://msdn.microsoft.com/robotics/ as it supports both MPI
and dynamic threading style of parallelism
There is MPI style messaging and ..

OpenMP annotation or Automatic Parallelism of existing
software is practical way to use those pesky cores with existing
code
• As parallelism is typically not expressed precisely, one needs luck to get
good performance
• Remember writing in Fortran, C, C#, Java … throws away information
about parallelism

HPCS Languages should be able to properly express parallelism
but we do not know how efficient and reliable compilers will be
• High Performance Fortran failed as language expressed a subset of
parallelism and compilers did not give predictable performance

PGAS (Partitioned Global Address Space) like UPC, Co-array
Fortran, Titanium, HPJava
• One decomposes application into parts and writes the code for each
component but use some form of global index
• Compiler generates synchronization and messaging
• PGAS approach should work but has never been widely used – presumably
because compilers not mature
Summary of micro-parallelism


On new applications, use MPI/locks with explicit
user decomposition
A subset of applications can use “data parallel”
compilers which follow in HPF footsteps
• Graphics Chips and Cell processor motivate such
special compilers but not clear how many
applications can be done this way

OpenMP and/or Compiler-based Automatic
Parallelism for existing codes in conventional
languages
Composition of Parallel Components










The composition (macro-parallelism) step has many excellent solutions
as this does not have the same drastic synchronization and correctness
constraints as one has for scalable kernels
• Unlike micro-parallelism step which has no very good solutions
Task parallelism in languages such as C++, C#, Java and Fortran90;
General scripting languages like PHP Perl Python
Domain specific environments like Matlab and Mathematica
Functional Languages like MapReduce, F#
HeNCE, AVS and Khoros from the past and CCA from DoE
Web Service/Grid Workflow like Taverna, Kepler, InforSense KDE,
Pipeline Pilot (from SciTegic) and the LEAD environment built at
Indiana University.
Web solutions like Mash-ups and DSS
Many scientific applications use MPI for the coarse grain composition
as well as fine grain parallelism but this doesn’t seem elegant
The new languages from Darpa’s HPCS program support task
parallelism (composition of parallel components) decoupling
composition and scalable parallelism will remain popular and must be
supported.
Integration of Services and “MPI”/Threads






Kernels and Composition must be supported both inside chips (the multicore
problem) and between machines in clusters (the traditional parallel computing
problem) or Grids.
The scalable parallelism (kernel) problem is typically only interesting on true
parallel computers (rather than grids) as the algorithms require low
communication latency.
However composition is similar in both parallel and distributed scenarios and it
seems useful to allow the use of Grid and Web composition tools for the parallel
problem.
• This should allow parallel computing to exploit large investment in service
programming environments
Thus in SALSA we express parallel kernels not as traditional libraries but as
(some variant of) services so they can be used by non expert programmers
Bottom Line: We need a runtime that supports inter-service linkage and microparallelism linkage
CCR and DSS have this property
• Does it work and what are performance costs of the universality of
runtime?
• Messaging need not be explicit for large data sets inside multicore node.
However still use small messages to synchronize