CGB Indiana University Lunchtime Talk January 22 2008 Geoffrey Fox discussing work with:

1. Robust High Performance Optimization for Clustering, Multi-Dimensional Scaling and Mixture Models CGB Indiana University Lunchtime TalkJanuary 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 CenterIndiana University George Chrysanthakopoulos, Henrik Frystyk Nielsen Microsoft Research, Redmond WA http://www.infomall.org/multicore gcf@indiana.edu, http://www.infomall.org

2. 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.

3. 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 co-processors 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

4. Intel’s Projection

5. 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

6. Intel’s Application Stack

7. Total Asian Hispanic Renters 30 Clusters 10 Clusters 30 Clusters

8. Clustering algorithm annealing by decreasing distance scale and gradually finds more clusters as resolution improved Here we see 10 clusters increasing to 30 as algorithm progresses

9. Parallel MulticoreDeterministic Annealing Clustering Parallel Overheadon 8 Threads on Intel 8 core Speedup = 8/(1+Overhead) = 7.6 if Overhead is 0.05 10 Clusters Overhead = Constant1 + Constant2/n Constant1 = 0.05 to 0.1 (Windows) due to threadruntime fluctuations 20 Clusters 10000/(Grain Size n = points per core)

10. 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 lnqi() • One can approximate 2F(0) /kl  i=1N (1/ qi)qi(0) /k qi(0) /ldropping term proportional to2qi(0) /k l • This shows that second derivative matrix is positive (semi-) definite

11. Basic Nonlinear Optimization II • Taylor expansion of F (Newton’s method) is just • 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

12. …. ……. Basic Nonlinear Optimization III • In my experience, the lowest eigenvalues kof 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

13. 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< cutwhere is cutis dynamically adjusted and I typically started with cut =  maxwith 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

14. 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 + Qaccomplished 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!

15. 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 = - Fwhich 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 )

16. 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

17. Need to make all this parallel OSCAR Document Analysis InChI Generation/Search Computational Chemistry (Gamess, Jaguar etc.) Varuna.net Quantum Chemistry CICC Chemical Informatics and Cyberinfrastructure Collaboratory Web Service Infrastructure Portal Services RSS Feeds User Profiles Collaboration as in Sakai Core Grid Services Service Registry Job Submission and Management Local Clusters IU Big Red, TeraGrid, Open Science Grid

18. 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

19. 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

20. Deterministic Annealing for Clustering I

21. Deterministic Annealing for Clustering II • 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

22. DeterministicAnnealing • Minimum evolving as temperature decreases • Movement at fixed temperature going to local minima if not initialized “correctly F({y}, T) Solve Linear Equations for each temperature Nonlinearity effects mitigated by approximating with solution at previous higher temperature Configuration {y}

23. Views from Past on Physical Optimization

24. You need to Specify a minimum temperature

26. Parallelism • All the methods – Deterministic Annealing, Mixture Models, GTM, PCA and Newton’s method involve • Sums over data pointsx=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

27. Where are we for Clustering? • We have deterministically annealed clustering running well on 8-core (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

28. 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 6-20 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.

29. In detail, different groups have different cluster centers

30. Age < 5 v. Age > 75 Age < 5 v. Age 25-34

31. Parallel Multicore Deterministic Annealing Clustering Parallel Overhead for large (2M points) Indiana Census clustering on 8 Threads Intel 8bThis fluctuating overhead due to 5-10% runtime fluctuations between threads “Constant1” Increasing number of clusters decreases communication/memory bandwidth overheads

32. Parallel Multicore Deterministic Annealing Clustering Parallel Overhead for subset of PubChem clustering on 8 Threads (Intel 8b) The fluctuating overhead is reduced to 2% (as bits not doubles)40,000 points with 1052 binary properties (Census is 2 real valued properties) “Constant1” Increasing number of clusters decreases communication/memory bandwidth overheads

33. 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

34. 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=1nweight(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/ijwhich emphasizes small distances ij

35. 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 Ykin original space and a “mean”Mabout which we do the expansion • For cluster splitting the mean Mis 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

36. 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

37. Classical MDS: SMACOFScaling 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 • Then one can establish a rigorous bound to show that one always decreases Stress (X) by iteration labeled by t if • This presumably will always find local minima?

38. 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)

39. 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=1MWm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

40. 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

41. Basic Gaussian Mixture Models • Mixture models write the probability of a point X(x) as • And set the Liklihood L = x=1N p(x,) where parameters  are centers Y(k), Probability Pkthat point generated by center k and covariance matrix (k)2 • Expectation Maximization is (from 1977) iterative solution to minimizing – lnL • It starts by estimating • And using this to calculate allthe  using formulae like: • 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 mixturesM Kloppenburg, P Tavan - Physical Review E, 1997 Cited by 18

42. DA / EM / DAEM / GTM I • The General Formula for “F = - lnL” or “Free Energy” F = E-TS is: • 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

43. DA / EM / DAEM / GTM II • The General Formula for “F = - lnL-TS” is: • Note x is a label and E(x) is data • For Deterministic AnnealingGaussian 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

44. DA / EM / DAEM / GTM III • The General Formula for “F = - lnL” is: • 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=1MWmm(X(k)) • Fix m(X) = exp( - 0.5 (X-m)2/2 ) but other choices possible • Vary Wm andbut fix values of M and K a priori • Y(k) E(x) Wm are vectors in original high D dimensional space • X(k) and m are vectors in 2 dimensional mapped space There is presumably a version of GTM using deterministic annealing using either a pure cluster basis – anneal in  orusing an annealed Gaussian mixture with temperature 1/ plus intrinsic widths

45. 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

46. 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.

47. 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

48. 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 CCRhttp://msdn.microsoft.com/robotics/ as it supports both MPI and dynamic threading style of parallelism

49. 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

50. 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