# Parallel Framework for Dimensionality Reduction of Large-Scale Datasets

Authors: Sai Kiranmayee Samudrala; Jaroslaw Zola; Srinivas Aluru; Baskar Ganapathysubramanian
Journal: Scientific Programming
 (2015)

---

## Abstract

Dimensionality reduction refers to a set of mathematical techniques used to
reduce complexity of the original high-dimensional data, while preserving its
selected properties. Improvements in simulation strategies and experimental
data collection methods are resulting in a deluge of heterogeneous and
high-dimensional data, which often makes dimensionality reduction the only
viable way to gain qualitative and quantitative understanding of the
data. However, existing dimensionality reduction software often does not scale to
datasets arising in real-life applications, which may consist of thousands of
points with millions of dimensions. In this paper, we propose a parallel
framework for dimensionality reduction of large-scale data. We identify key
components underlying the spectral dimensionality reduction techniques, and
propose their efficient parallel implementation. We show that the resulting
framework can be used to process datasets consisting of millions of points when
executed on a 16,000-core cluster, which is beyond the reach of currently
available methods. To further demonstrate applicability of our framework we
perform dimensionality reduction of 75,000 images representing morphology
evolution during manufacturing of organic solar cells in order to identify how
processing parameters affect morphology evolution.


---

## Body

## 1. Introduction

Computational analysis of high-dimensional data continues to be a challenging problem, spurring the development of numerous computational techniques. An important and emerging class of methods for dealing with such data is dimensionality reduction. In many applications, features of interest can be preserved while mapping the high dimensionality data to a small number of dimensions. These mappings include popular techniques such as principle component analysis (PCA) [1] and complex nonlinear maps such as Isomap [2] and kernel PCA [3].Linear manifold learning techniques, for example, PCA or multidimensional scaling [4–7], existed as orthogonalization methods for several decades. Nonlinear methods like Isomap, LLE (locally linear embedding) [8], and Hessian LLE [9] were discovered recently. Another class of methods that emerged in the past few years are the unsupervised learning techniques, including artificial neural networks for Sammon’s nonlinear mapping [10], Kohenen’s or self organizing maps (SOM) [11], and curvilinear component analysis [12]. Modifications to the existing algorithms of manifold learning, to improve either their efficiency or performance, were another area where efforts were focused [13–16]. For example, Landmark Isomap [17] is a modification to the original Isomap method to extend its usage to larger datasets by picking a few representative points and applying Isomap technique to them. Along with the emergence of new manifold learning techniques, different sequential implementations of these techniques, targeting various hardware platforms and various programming languages, have been developed [18, 19].Dimensionality reduction techniques are often compute intensive and do not easily scale to large datasets. Recent advances in high-throughput measurements using physical entities, such as sensors, or results of complex numerical simulations are generating data of extremely high dimensionality. It is becoming increasingly difficult to process such data sequentially.In this paper, we propose a parallel framework for dimensionality reduction. Rather than focusing on a particular method, we consider the class of spectral dimensionality reduction methods. Till date, few efforts have been made in developing parallel implementations of these methods, other than development of a parallel version of PCA [20, 21], parallelization of multidimensional scaling (MDS) for genomic data [22], and algorithm development for GPU platforms [23, 24].We perform a systematic analysis of spectral dimensionality reduction techniques and provide their unified view that can be exploited by dimensionality reduction algorithm designers. We identify common computational building blocks required for implementing spectral dimensionality reduction methods and use these abstractions to derive a common parallel framework. We implement such a framework and show that it can handle large datasets and it scales to thousands of processors. We demonstrate advantages of our software by analyzing 75,000 images of morphology evolution during manufacturing of organic solar cells, which enables us to visually inspect and correlate fabrication parameters with morphology.The remainder of this paper is organized as follows. In Section2 we introduce the dimensionality reduction problem and describe basic spectral dimensionality reduction techniques, highlighting their computational kernels. In Section 3 we provide a detailed description of our parallel framework including algorithmic solutions. Finally, in Section 4 we present experimental results, and we conclude the paper in Section 5.


## 2. Definitions and Methods Overview

The problem of dimensionality reduction can be formulated as follows: Consider a setX=x0,x1,…,xn-1 of n points, where xi∈RD, and D≫1. We are interested in finding a set Y=y0,y1,…,yn-1, such that yi∈Rd, d≪D, and ∀i,j|xi-xj|h=|yi-yj|h. Here, |a-b|h denotes a specific norm that captures properties we want to preserve during dimensionality reduction [25]. For instance, by defining h as Euclidean norm we preserve Euclidean distance, thus obtaining a reduction equivalent to the standard technique of principal component analysis (PCA) [1]. Similarly, defining h to be the angular distance (or conformal distance [26]) results in locally linear embedding (LLE) [8] that preserves local angles between points. In a typical application [27, 28], xi represents a state of the analyzed system, for example, temperature field and concentration distribution. Such state description can be derived from experimental sensor data or can be the result of a numerical simulation. However, irrespective of the source, it is characterized by high dimensionality, that is, D is typically of the order of 106 [29]. While xi represents just a single state of the system, common data acquisition setups deliver large collections of such observations, which correspond to the temporal or parametric evolution of the system [27]. Thus, the cardinality n of the resulting set X is usually large (n~104-105). Intuitively, information obfuscation increases with data dimensionality. Therefore, in the process of dimensionality reduction (DR) we seek as small a dimension d as possible, given constraints induced by the norm |a-b|h [25]. Routinely, d<4 as it permits, for instance, visualization of the set Y.DR techniques have been extensively researched over the last decade [25]. In particular, methods based on spectral data decomposition have been very successful [1, 2, 9] and have been widely adopted. Early approaches in this category exploited simple linear structure of the data, for example, PCA or multidimensional scaling (MDS) [30]. More recently, techniques that can unravel complex nonlinear structures in the data, for example, Isomap [2], LLE, and kernel PCA [3], have been developed. While all these methods have been proposed taking into account specific applications [19, 25], their underlying formulations share similar algorithmic mechanisms. In what follows we provide a more detailed overview of spectral DR techniques and we identify their common computational kernels that form the basis for our parallel framework.


### 2.1. Spectral Dimensionality Reduction

The goal of DR is to identify a low-dimensional representationY of the original dataset X, that preserves certain predefined properties. The key idea underpinning spectral DR can be explained as follows. We encode desired information about X, that is, topology or distance, in its entirety by considering all pairs of points in X. This encoding is represented as a matrix An×n. Next, we subject matrix A to unitary transformation V, that is, transformation that preserves norm of A, to obtain its sparsest form Λ, where A=VΛVT. Here, Λn×n is a diagonal matrix with rapidly diminishing entries. As a result, it is sufficient to consider only d entries of Λ to capture all the information encoded in A. These d entries constitute the set Y. The above procedure hinges on the fact that unitary transformations preserve original properties of A [31]. Note also that it requires a method to construct matrix A in the first place. Indeed, what differentiates various spectral methods is the way information is encoded in A.We summarize the general idea of spectral DR in Algorithm1. In the first four steps we construct the matrix A. As indicated, this matrix encodes information about the property that we wish to preserve in the process of DR. To obtain A we first identify the k nearest neighbors (KNN) of each point xi∈X. Note that currently all studied methods use KNN defined in Euclidean space. This enables us to define a weighted graph G that encapsulates, both distance and topological, properties of the set X. Given graph G, we can construct a function FG:X×X→R to isolate the desired property. For instance, consider the Isomap algorithm in which the geodesic distance is maintained. In this case, FG returns the length of the shortest path between xi and xj in G. Note that for some methods FG is very simple; for example, for PCA it is equivalent to a distance measure ω, FG(xi,xj)=ω(xi,xj), while for other methods FG can be more involved. Differences between various DR methods and their corresponding function FG are outlined in Table 1. The property extracted by function FG is stored in an auxiliary matrix W, which is next normalized to obtain matrix A. This process of normalization is a simple algebraic transformation, which ensures that A is centered and hence that the final low-dimensional set of points Y contains the origin and is not an affine translation [31]. Subsequently, A is spectrally decomposed into its eigenvalues that constitute the sparsest representation of A. Resulting eigenvectors and eigenvalues are then postprocessed to extract the set Y of low-dimensional points.Table 1
Comparison of selected spectral dimensionality reduction methods.
PCA
Isomap
LLE
Parameterk in KNN
n
~d
~d
FunctionFG
ω
(
x
i
,
x
j
)
Length of the shortest path betweenxi and xj in G
α
i
j if (xi,xj)∈E, 0 otherwise, where xi=∑xl:(xi,xl)∈E‍αilxl
Normalization
W
i
j
*
=
W
i
j
2, A=HTW*H
W
i
j
*
=
W
i
j
2,  A=HTW*H
A
=
(
I
-
W
)
-
1
(
I
-
W
)
-
T
Note:I is the identity matrix, and H=I-(1/n)1n×n.Algorithm 1:Spectral dimensionality reduction.
Input: Set X=x0,x1,…,xn-1, xi∈RD, and the target dimension d.
Output: Set Y=y0,y1,…,yn-1, yi∈Rd.
(1)  For eachxi∈X find its k nearest neighbors.
(2) Define directed weighted graphG=(X,E,ω),
where(xi,xj)∈E iff xj is a neighbor of xi,
andω(xi,xj) is a distance measure,
usuallyω(xi,xj)=|xi-xj|2.
(3) LetWij=FG(xi,xj), where FG extracts specific property
from graphG.
(4) NormalizeW to obtain matrix A.
(5) Find eigenvectors ofA, A=VΛVT.
(6) Identify latent dimensionalityd.
(7)Y is represented by the first d rows of V.The abstract representation of spectral DR methods in Algorithm1 is based on a thorough analysis of existing techniques [1, 2, 8]. While this representation is compact, it offers sufficient flexibility to, for instance, design new dimensionality reduction procedures. At the same time it provides clear separation of individual computational steps and explicates data flow in any DR process. We exploit both these facts when designing our parallel framework.


### 2.2. Performance Analysis of Dimensionality Reduction Methods

We used the above presentation of DR methods to identify their basic computational kernels. To better understand how these kernels contribute to the overall performance of different DR methods, we performed a set of experiments using domain specific implementation in Matlab. Experiments were carried out for varyingn and a fixed D=1000 on a workstation with 8 GB of RAM and an Intel 3.2 GHz processor. Obtained results are presented in Table 2.Table 2
Run time (in seconds) of different DR components (d=3).
n
100
1000
2000
4000
KNN
0.08640
1.34998
5.66768
27.91930
W in PCA
—
—
—
—
W in Isomap
0.06470
14.9030
130.130
1153.30
W in LLE
0.08960
0.12601
0.24609
0.49253
Normalize
0.00195
0.11875
0.74934
5.56630
Eigensolve
0.02916
0.05536
0.23267
0.85211
ExtractY
0.00020
0.00014
0.00016
0.00022As can be seen, the run time of analyzed methods is dominated by two steps, namely,KNN and construction of the auxiliary matrix W. Together they account for 99.8% of the total execution time for n=4000. In our implementation the KNN procedure depends on all-versus-all distance calculations. This is justified taking into account that D is very large and thus efficient algorithmic strategies for KNN, for example, based on hierarchical space decomposition [32], are infeasible. Consequently, complexity of this step is O(Dn2). The cost of computing matrix W depends explicitly on the definition of function FG. Among existing DR techniques this function is the most complex for the Isomap method. Recall that in the process of DR we are interested in preserving either distance or local topology characteristics. Local topology properties can be directly obtained from KNN [8, 9, 33], inducing computationally efficient definition of FG. Conversely, distance characteristics must conform to global constraints and therefore have higher computational complexity [34]. In case of Isomap, pairwise geodesic distances can be efficiently derived from all-pairs shortest path distances using, for example, Floyd-Warshall algorithm, with O(n3) worst-case complexity.Another significant DR component is normalization. Although implementation of this step varies between different methods it is invariably dominated by matrix-matrix multiplication. Therefore, we assume overall normalization complexity to beO(n3). The last important component is the eigenvalue solver. In general, complexity of this kernel varies depending on the particular solver used. Commonly employed algorithms include Lanczos method [35], Krylov subspace methods [36], or deflation-based power methods [37, 38]. The choice of method is driven by the structure of the matrix and the number of required eigenvalues. Standard distance preserving DR methods operate on dense symmetric matrices, while topology preserving methods involve sparse symmetric matrices. Accordingly, complexity of these techniques is usually O(dn2), where d is the number of desired eigenvalues.A final key factor we have to consider is memory complexity of the described kernels. Here, the main contributing structures are matricesW and A. These matrices are most often dense and in the majority of cases require O(n2) storage. Because KNN directly depends on distances between all pairs, it utilizes an n×n matrix as well. Finally, input dataset X requires O(Dn) memory.One important caveat that affects the above analysis is the relationship betweenD and n. In many applications D is significantly greater than n. This is not surprising taking into account that acquiring high resolution data (hence high-dimensional) is resource intensive. Therefore one may expect that with increasing D there is usually decrease of n. In our applications [29, 39] it is not uncommon that D=O(n2) or even D=O(n4). Consequently, the KNN step in Algorithm 1 becomes the most compute intensive while memory requirement is dominated by the input data. Observe that this trend is reflected in our experimental data.


## 2.1. Spectral Dimensionality Reduction

The goal of DR is to identify a low-dimensional representationY of the original dataset X, that preserves certain predefined properties. The key idea underpinning spectral DR can be explained as follows. We encode desired information about X, that is, topology or distance, in its entirety by considering all pairs of points in X. This encoding is represented as a matrix An×n. Next, we subject matrix A to unitary transformation V, that is, transformation that preserves norm of A, to obtain its sparsest form Λ, where A=VΛVT. Here, Λn×n is a diagonal matrix with rapidly diminishing entries. As a result, it is sufficient to consider only d entries of Λ to capture all the information encoded in A. These d entries constitute the set Y. The above procedure hinges on the fact that unitary transformations preserve original properties of A [31]. Note also that it requires a method to construct matrix A in the first place. Indeed, what differentiates various spectral methods is the way information is encoded in A.We summarize the general idea of spectral DR in Algorithm1. In the first four steps we construct the matrix A. As indicated, this matrix encodes information about the property that we wish to preserve in the process of DR. To obtain A we first identify the k nearest neighbors (KNN) of each point xi∈X. Note that currently all studied methods use KNN defined in Euclidean space. This enables us to define a weighted graph G that encapsulates, both distance and topological, properties of the set X. Given graph G, we can construct a function FG:X×X→R to isolate the desired property. For instance, consider the Isomap algorithm in which the geodesic distance is maintained. In this case, FG returns the length of the shortest path between xi and xj in G. Note that for some methods FG is very simple; for example, for PCA it is equivalent to a distance measure ω, FG(xi,xj)=ω(xi,xj), while for other methods FG can be more involved. Differences between various DR methods and their corresponding function FG are outlined in Table 1. The property extracted by function FG is stored in an auxiliary matrix W, which is next normalized to obtain matrix A. This process of normalization is a simple algebraic transformation, which ensures that A is centered and hence that the final low-dimensional set of points Y contains the origin and is not an affine translation [31]. Subsequently, A is spectrally decomposed into its eigenvalues that constitute the sparsest representation of A. Resulting eigenvectors and eigenvalues are then postprocessed to extract the set Y of low-dimensional points.Table 1
Comparison of selected spectral dimensionality reduction methods.
PCA
Isomap
LLE
Parameterk in KNN
n
~d
~d
FunctionFG
ω
(
x
i
,
x
j
)
Length of the shortest path betweenxi and xj in G
α
i
j if (xi,xj)∈E, 0 otherwise, where xi=∑xl:(xi,xl)∈E‍αilxl
Normalization
W
i
j
*
=
W
i
j
2, A=HTW*H
W
i
j
*
=
W
i
j
2,  A=HTW*H
A
=
(
I
-
W
)
-
1
(
I
-
W
)
-
T
Note:I is the identity matrix, and H=I-(1/n)1n×n.Algorithm 1:Spectral dimensionality reduction.
Input: Set X=x0,x1,…,xn-1, xi∈RD, and the target dimension d.
Output: Set Y=y0,y1,…,yn-1, yi∈Rd.
(1)  For eachxi∈X find its k nearest neighbors.
(2) Define directed weighted graphG=(X,E,ω),
where(xi,xj)∈E iff xj is a neighbor of xi,
andω(xi,xj) is a distance measure,
usuallyω(xi,xj)=|xi-xj|2.
(3) LetWij=FG(xi,xj), where FG extracts specific property
from graphG.
(4) NormalizeW to obtain matrix A.
(5) Find eigenvectors ofA, A=VΛVT.
(6) Identify latent dimensionalityd.
(7)Y is represented by the first d rows of V.The abstract representation of spectral DR methods in Algorithm1 is based on a thorough analysis of existing techniques [1, 2, 8]. While this representation is compact, it offers sufficient flexibility to, for instance, design new dimensionality reduction procedures. At the same time it provides clear separation of individual computational steps and explicates data flow in any DR process. We exploit both these facts when designing our parallel framework.


## 2.2. Performance Analysis of Dimensionality Reduction Methods

We used the above presentation of DR methods to identify their basic computational kernels. To better understand how these kernels contribute to the overall performance of different DR methods, we performed a set of experiments using domain specific implementation in Matlab. Experiments were carried out for varyingn and a fixed D=1000 on a workstation with 8 GB of RAM and an Intel 3.2 GHz processor. Obtained results are presented in Table 2.Table 2
Run time (in seconds) of different DR components (d=3).
n
100
1000
2000
4000
KNN
0.08640
1.34998
5.66768
27.91930
W in PCA
—
—
—
—
W in Isomap
0.06470
14.9030
130.130
1153.30
W in LLE
0.08960
0.12601
0.24609
0.49253
Normalize
0.00195
0.11875
0.74934
5.56630
Eigensolve
0.02916
0.05536
0.23267
0.85211
ExtractY
0.00020
0.00014
0.00016
0.00022As can be seen, the run time of analyzed methods is dominated by two steps, namely,KNN and construction of the auxiliary matrix W. Together they account for 99.8% of the total execution time for n=4000. In our implementation the KNN procedure depends on all-versus-all distance calculations. This is justified taking into account that D is very large and thus efficient algorithmic strategies for KNN, for example, based on hierarchical space decomposition [32], are infeasible. Consequently, complexity of this step is O(Dn2). The cost of computing matrix W depends explicitly on the definition of function FG. Among existing DR techniques this function is the most complex for the Isomap method. Recall that in the process of DR we are interested in preserving either distance or local topology characteristics. Local topology properties can be directly obtained from KNN [8, 9, 33], inducing computationally efficient definition of FG. Conversely, distance characteristics must conform to global constraints and therefore have higher computational complexity [34]. In case of Isomap, pairwise geodesic distances can be efficiently derived from all-pairs shortest path distances using, for example, Floyd-Warshall algorithm, with O(n3) worst-case complexity.Another significant DR component is normalization. Although implementation of this step varies between different methods it is invariably dominated by matrix-matrix multiplication. Therefore, we assume overall normalization complexity to beO(n3). The last important component is the eigenvalue solver. In general, complexity of this kernel varies depending on the particular solver used. Commonly employed algorithms include Lanczos method [35], Krylov subspace methods [36], or deflation-based power methods [37, 38]. The choice of method is driven by the structure of the matrix and the number of required eigenvalues. Standard distance preserving DR methods operate on dense symmetric matrices, while topology preserving methods involve sparse symmetric matrices. Accordingly, complexity of these techniques is usually O(dn2), where d is the number of desired eigenvalues.A final key factor we have to consider is memory complexity of the described kernels. Here, the main contributing structures are matricesW and A. These matrices are most often dense and in the majority of cases require O(n2) storage. Because KNN directly depends on distances between all pairs, it utilizes an n×n matrix as well. Finally, input dataset X requires O(Dn) memory.One important caveat that affects the above analysis is the relationship betweenD and n. In many applications D is significantly greater than n. This is not surprising taking into account that acquiring high resolution data (hence high-dimensional) is resource intensive. Therefore one may expect that with increasing D there is usually decrease of n. In our applications [29, 39] it is not uncommon that D=O(n2) or even D=O(n4). Consequently, the KNN step in Algorithm 1 becomes the most compute intensive while memory requirement is dominated by the input data. Observe that this trend is reflected in our experimental data.


## 3. Parallel Framework for Dimensionality Reduction

Dimensionality reduction very quickly becomes both memory and compute intensive, irrespective of the particular method. Memory consumption arises from the size of input data and the auxiliary matrices created in the process. The computational cost is dominated by pairwise computations and weight matrix construction. The goal of our framework is to scale DR methods to very large datasets that could be analyzed on large parallel machines.We designed our parallel DR package following the general outline presented in Algorithm1. Taking into account significant memory and computational complexity, we focused on distributed memory machines with MPI. To ensure modularity of the framework without sacrificing performance and scalability, we decided to employ a scheme in which processors are organized into a logical 2D mesh. In what follows, we assume a simple point-to-point communication model with latency ts and bandwidth 1/tm.


### 3.1. Constructing GraphG

The graph construction procedure is based on identifyingk nearest neighbors of each input point. Because of the high dimensionality of the input data it is advantageous to implement KNN in two steps, where we first compute all pairwise distances and then we identify neighbors in a simple scan. Note that these pairwise distances actually represent the distance measure ω (see Algorithm 1). Therefore we will consider ω to be an n×n distance matrix. Parallel pairwise computation is a well studied problem [40]. Here, we benefit from our earlier experience with accelerating pairwise computations on heterogeneous parallel processors [41].Letp=q2 denote the number of processors conceptualized as organized into a q×q virtual mesh. We decompose ω into blocks of (n/q)×(n/q) elements. Processor with coordinates (i,j) is responsible for computing elements of ω within block (i,j). This scheme requires that each processor stores two blocks of n/q points of the input dataset X that correspond to row-vectors and column-vectors used to compute respective part of the matrix ω. In our implementation, the distribution of the input dataset is performed by parallel I/O with initially preprocessed X. Note that to obtain a single element of ω we perform |a-b|2 norm computations, which are particularly well suited for vectorization. Therefore, we hand-tuned our code to benefit from SIMD extensions of modern processors.Given pairwise distances, the second step is to identify neighbors of individual points (i.e., vertices ofG). This step is executed only for methods where k<n, which virtually involves all methods other than PCA (see Table 1). As in the case of pairwise computations it can be efficiently parallelized using the following scheme. Initially, each processor creates a set of k candidate neighbors with respect to the block of matrix ω it stores. Specifically, processor with coordinates (i,j) searches for neighbors of the set of points xin/q,…,x(i+1)n/q by analyzing rows of its local block of ω. Because k is very small this operation can be performed using a simple linear scan. Next, all processors within the same row perform all-to-all communication to aggregate candidate neighbors. Here, the candidate neighborhood of point xl is assembled on the processor with coordinates lq/n,lp/nmod⁡q. This processor merges candidate neighborhood lists into the final set of k nearest neighbors. Observe that this completes the graph construction phase—graph G is now stored in the form of adjacency list distributed over p processors.The computational complexity of the entire procedure can be decomposed intoODn2/p cost of computing distance matrix, O(n/p)log⁡k cost of performing local KNN search, and O(ts+tm(nk/p))p cost of distributing local KNN information to all processors sharing a row. Hence, the total complexity is bounded by ODn2/p+(ts+tm(nk/p))p, which is optimal for p<n2.


### 3.2. Building Auxiliary MatrixW

Given graphG we proceed to the next step, which involves constructing the auxiliary matrix W from the information encapsulated in G. As is the case of ω we distribute W over q×q mesh of processors.Recall that the formulation of DR methods proposed in Algorithm1 ensures that the only step that is method dependent is the construction of matrix W. Consequently, any parallel implementation of this step will vary but it will reflect limitations inherent to the sequential counterpart. Specifically, topology preserving methods, such as, LLE, will involve only local data and hence will be amenable to embarrassing parallelism with limited or no communication. Conversely, distance preserving methods will inevitably require a global data view and thus potentially more sophisticated parallelization strategies. Following our previous claim regarding complexity of Isomap we focus our presentation on the parallel implementation of this particular method.The functionFG used in Isomap is based on the geodesic distance, which has been mathematically shown to be asymptotically equivalent to a graph distance in G (i.e., shortest path distance) [42]. However, the geodesic distance is a metric, while all-pairs shortest paths in directed graph G do not have to satisfy the symmetry condition. Therefore, to obtain W, special attention must be paid to how shortest path distances are used. More precisely, graph G must be transformed to ensure that it is symmetric. Note that after such transformation the graph is no longer regular; that is, certain nodes may have more than k neighbors (see Figure 1).Figure 1
GraphG before (a) and after (b) symmetrization for an example set of seven points and k=2.
(a)
(b)Taking into account the above requirements we obtain the following procedure of constructingW in parallel. First, all processors within the same row perform all-to-all communication to replicate graph G. As a result, each column of processors stores a copy of the entire graph G, that is, row-wise distributed between q processors in that column. Thanks to this, each processor can initialize its local part of W without further communication. After initialization W represents the distributed adjacency matrix of G, where Wij=ω(xi,xj) if xj is a neighbor of xi, and +∞ otherwise. In the next step symmetrization procedure is executed. Processors with coordinates (i,j) and (j,i), where i≠j, exchange respective blocks of W and select element-wise minimum value between blocks. Similar operation is performed locally by processors on the diagonal, that is, processors for which i=j. At this stage W can be used to identify all-pairs shortest paths. Several parallel algorithms have been proposed to address this problem, targeting the PRAM model [43–45], shared memory architectures [46], and multi/many cores [47–49], as well as distributed memory machines [45, 50]. Among the existing parallel strategies we decided to adopt the checkerboard version of the parallel Floyd algorithm [46]. Briefly, the method proceeds in n iterations, where in each iteration every processor performs On2/p operations to update its local block of W. All processors are synchronized in each iteration, owing to the fact that in iteration l, lth row and lth column of W have to be broadcasted. The broadcast is performed between processors that share the same row/column. After n iterations, matrix W stores all-pairs shortest path, which concludes the entire procedure.Complexity of this phase isOn3/p+nts+tm(n/p)log⁡(p) and is dominated by the parallel Floyd’s algorithm. While replication and symmetrization of G can be executed efficiently in On2/p+ts+tm(n2/p) time, all-pairs path searching involves extensive communication that slightly hinders scalability. Nevertheless, the algorithm remains scalable as long as p<n2/log2(n), which is true in the majority of real-life cases.


### 3.3. Matrix Normalization

The goal of normalization is to transform matrixW such that the resulting matrix A is both row and column centered; that is, ∑i‍Aij=0 and ∑j‍Aij=0. The normalization stage in all cases consists of matrix-matrix multiplication (see Table 1). However, in certain situations, especially in distance-preserving methods, explicit matrix multiplication can be avoided by taking advantage of structural properties of one of the matrices (e.g., the matrix H in Table 1). This is the case for, for example, PCA and Isomap, where we exploit the fact that matrices H and HT are given analytically and thus can be generated in-place on each processor that requires them to perform multiplication. Consequently, the communication pattern inherent to the standard parallel matrix-matrix multiplication algorithms is simplified to one parallel reduction in the final dot-product operation. The complexity of this approach is On3/p+(ts+tm(n2/p))p.


### 3.4. Finding Eigenvalues

Computing eigenvalues is the final step in the dimensionality reduction process. Although parallel eigensolvers are readily available, they are usually designed for shared memory and multi/many core architectures [51–54]. This unfortunately makes them impractical for our purposes. At the same time, existing distributed memory solutions are not scalable and cannot handle large and dense data. For instance, one of the more popular packages,SLEPc [55], uses a simple 1D decomposition and in our tests did not scale to more than 4096 processors. A more recent library,elemental [56], which is still under development, offers 2D block-cyclic decomposition, but relies on a fixed block size (private communication). Finally,ELPA [57], which is the most promising library in this category, is also under development.For these reasons we decided to implement a custom eigenvalue solver that exploits special properties of matrixA (symmetric, positive semidefinite) and computes only the first d eigenvalues. Our solver is based on the power method [31] and matrix deflation and is outlined in Algorithm 2. Note that power methods are considered easy but not efficient to parallelize. At the same time, however, they are at heart of several important real-life systems, for instance, Google’s PageRank [58].Algorithm 2:2D-Block Parallel Power Method.
Input: Matrix An×n and the required numbed of eigenvalues d.
2D mesh of p=q×q processors.
Output: Set of eigenvalues and eigenvectors of A,
Λ0,…,d-1  = λ0,λ1,…,λd-1 and V0,…,d-1=v0,v1,…,vd-1.
(1) Let x be a column-wise distributed vector in Rn.
(2)for  i←0:d-1  do
(3) Initialize x randomly. Processors within the same column use the same seed.
(4) column←true
(5) while  not converged  do
(6)  Compute z=Ax locally.
(7)  if  column=true  then
(8)   Perform column-wise all-reduce to obtain z.
(9)  else
(10)     Perform row-wise all-reduce to obtainz.
(11)end if
(12)column←¬column
(13)x←z
(14)end while
(15)   Computeu=Az as in steps  (6)–(11).
(16)   Replicate entire vectoru and z on each processor.
(17)λi←z·uz·z
(18)vi←zz2
(19)   Deflate local block ofA: A←A-λiviviT.
(20)end forIn general, our approach follows the standard scheme of the power method (lines 3–18), repeatedd times to identify first d largest eigenvalues. After identification of an eigenvalue and its associated eigenvector, the matrix A is deflated—the contribution of the vector is removed from A (line 19). Observe that power method involves nested matrix-vector product (lines 6–11) that under normal circumstances would require parallel vector transposition. However, our parallel implementation benefits from the fact that A is symmetric, hence eliminating need for vector transposition. Indeed, the entire procedure consists of local matrix-vector product followed by all-reduce operation. Here, the reduction operation alternates between columns and rows as required to ensure that vector x is stored properly. Note that the power method is bounded by convergence criteria (line 5). In our case we use one of several popular conditions, which involves checking relative error between the current and previous estimate of the eigenvalue that can be performed every several iterations. We also note that convergence is significantly improved by using a matrix shifting strategy in the form A=A-δI, where δ is a positive number [59].Extracting eigenvalue and eigenvector in iterationi (lines 11–18) depends on vectors u and z, while deflation step involves λi. Therefore, it is advantageous to replicate both u and z in their entirety on each processor. We achieve this with all-to-all communication executed by processors within the same row. This allows us to execute the deflation step in parallel, with each processor updating its local block of matrix A. Thus, the complexity of a single iteration of the power method is On3/p+ts+tm(n/p)log⁡(p), while the deflation step is On2/p+tsp+tmn.To conclude this section we would like to emphasize that our solver operates under the same assumptions as any power method. It requires that the firstd eigenvectors of A are linearly independent, the initial vector x generated in ith iteration is not orthogonal to the eigenvector vi, and finally, the first d eigenvalues are nondegenerate [31]. Note that these conditions are not restrictive and are easily satisfied in the context of dimensionality reduction.


## 3.1. Constructing GraphG

The graph construction procedure is based on identifyingk nearest neighbors of each input point. Because of the high dimensionality of the input data it is advantageous to implement KNN in two steps, where we first compute all pairwise distances and then we identify neighbors in a simple scan. Note that these pairwise distances actually represent the distance measure ω (see Algorithm 1). Therefore we will consider ω to be an n×n distance matrix. Parallel pairwise computation is a well studied problem [40]. Here, we benefit from our earlier experience with accelerating pairwise computations on heterogeneous parallel processors [41].Letp=q2 denote the number of processors conceptualized as organized into a q×q virtual mesh. We decompose ω into blocks of (n/q)×(n/q) elements. Processor with coordinates (i,j) is responsible for computing elements of ω within block (i,j). This scheme requires that each processor stores two blocks of n/q points of the input dataset X that correspond to row-vectors and column-vectors used to compute respective part of the matrix ω. In our implementation, the distribution of the input dataset is performed by parallel I/O with initially preprocessed X. Note that to obtain a single element of ω we perform |a-b|2 norm computations, which are particularly well suited for vectorization. Therefore, we hand-tuned our code to benefit from SIMD extensions of modern processors.Given pairwise distances, the second step is to identify neighbors of individual points (i.e., vertices ofG). This step is executed only for methods where k<n, which virtually involves all methods other than PCA (see Table 1). As in the case of pairwise computations it can be efficiently parallelized using the following scheme. Initially, each processor creates a set of k candidate neighbors with respect to the block of matrix ω it stores. Specifically, processor with coordinates (i,j) searches for neighbors of the set of points xin/q,…,x(i+1)n/q by analyzing rows of its local block of ω. Because k is very small this operation can be performed using a simple linear scan. Next, all processors within the same row perform all-to-all communication to aggregate candidate neighbors. Here, the candidate neighborhood of point xl is assembled on the processor with coordinates lq/n,lp/nmod⁡q. This processor merges candidate neighborhood lists into the final set of k nearest neighbors. Observe that this completes the graph construction phase—graph G is now stored in the form of adjacency list distributed over p processors.The computational complexity of the entire procedure can be decomposed intoODn2/p cost of computing distance matrix, O(n/p)log⁡k cost of performing local KNN search, and O(ts+tm(nk/p))p cost of distributing local KNN information to all processors sharing a row. Hence, the total complexity is bounded by ODn2/p+(ts+tm(nk/p))p, which is optimal for p<n2.


## 3.2. Building Auxiliary MatrixW

Given graphG we proceed to the next step, which involves constructing the auxiliary matrix W from the information encapsulated in G. As is the case of ω we distribute W over q×q mesh of processors.Recall that the formulation of DR methods proposed in Algorithm1 ensures that the only step that is method dependent is the construction of matrix W. Consequently, any parallel implementation of this step will vary but it will reflect limitations inherent to the sequential counterpart. Specifically, topology preserving methods, such as, LLE, will involve only local data and hence will be amenable to embarrassing parallelism with limited or no communication. Conversely, distance preserving methods will inevitably require a global data view and thus potentially more sophisticated parallelization strategies. Following our previous claim regarding complexity of Isomap we focus our presentation on the parallel implementation of this particular method.The functionFG used in Isomap is based on the geodesic distance, which has been mathematically shown to be asymptotically equivalent to a graph distance in G (i.e., shortest path distance) [42]. However, the geodesic distance is a metric, while all-pairs shortest paths in directed graph G do not have to satisfy the symmetry condition. Therefore, to obtain W, special attention must be paid to how shortest path distances are used. More precisely, graph G must be transformed to ensure that it is symmetric. Note that after such transformation the graph is no longer regular; that is, certain nodes may have more than k neighbors (see Figure 1).Figure 1
GraphG before (a) and after (b) symmetrization for an example set of seven points and k=2.
(a)
(b)Taking into account the above requirements we obtain the following procedure of constructingW in parallel. First, all processors within the same row perform all-to-all communication to replicate graph G. As a result, each column of processors stores a copy of the entire graph G, that is, row-wise distributed between q processors in that column. Thanks to this, each processor can initialize its local part of W without further communication. After initialization W represents the distributed adjacency matrix of G, where Wij=ω(xi,xj) if xj is a neighbor of xi, and +∞ otherwise. In the next step symmetrization procedure is executed. Processors with coordinates (i,j) and (j,i), where i≠j, exchange respective blocks of W and select element-wise minimum value between blocks. Similar operation is performed locally by processors on the diagonal, that is, processors for which i=j. At this stage W can be used to identify all-pairs shortest paths. Several parallel algorithms have been proposed to address this problem, targeting the PRAM model [43–45], shared memory architectures [46], and multi/many cores [47–49], as well as distributed memory machines [45, 50]. Among the existing parallel strategies we decided to adopt the checkerboard version of the parallel Floyd algorithm [46]. Briefly, the method proceeds in n iterations, where in each iteration every processor performs On2/p operations to update its local block of W. All processors are synchronized in each iteration, owing to the fact that in iteration l, lth row and lth column of W have to be broadcasted. The broadcast is performed between processors that share the same row/column. After n iterations, matrix W stores all-pairs shortest path, which concludes the entire procedure.Complexity of this phase isOn3/p+nts+tm(n/p)log⁡(p) and is dominated by the parallel Floyd’s algorithm. While replication and symmetrization of G can be executed efficiently in On2/p+ts+tm(n2/p) time, all-pairs path searching involves extensive communication that slightly hinders scalability. Nevertheless, the algorithm remains scalable as long as p<n2/log2(n), which is true in the majority of real-life cases.


## 3.3. Matrix Normalization

The goal of normalization is to transform matrixW such that the resulting matrix A is both row and column centered; that is, ∑i‍Aij=0 and ∑j‍Aij=0. The normalization stage in all cases consists of matrix-matrix multiplication (see Table 1). However, in certain situations, especially in distance-preserving methods, explicit matrix multiplication can be avoided by taking advantage of structural properties of one of the matrices (e.g., the matrix H in Table 1). This is the case for, for example, PCA and Isomap, where we exploit the fact that matrices H and HT are given analytically and thus can be generated in-place on each processor that requires them to perform multiplication. Consequently, the communication pattern inherent to the standard parallel matrix-matrix multiplication algorithms is simplified to one parallel reduction in the final dot-product operation. The complexity of this approach is On3/p+(ts+tm(n2/p))p.


## 3.4. Finding Eigenvalues

Computing eigenvalues is the final step in the dimensionality reduction process. Although parallel eigensolvers are readily available, they are usually designed for shared memory and multi/many core architectures [51–54]. This unfortunately makes them impractical for our purposes. At the same time, existing distributed memory solutions are not scalable and cannot handle large and dense data. For instance, one of the more popular packages,SLEPc [55], uses a simple 1D decomposition and in our tests did not scale to more than 4096 processors. A more recent library,elemental [56], which is still under development, offers 2D block-cyclic decomposition, but relies on a fixed block size (private communication). Finally,ELPA [57], which is the most promising library in this category, is also under development.For these reasons we decided to implement a custom eigenvalue solver that exploits special properties of matrixA (symmetric, positive semidefinite) and computes only the first d eigenvalues. Our solver is based on the power method [31] and matrix deflation and is outlined in Algorithm 2. Note that power methods are considered easy but not efficient to parallelize. At the same time, however, they are at heart of several important real-life systems, for instance, Google’s PageRank [58].Algorithm 2:2D-Block Parallel Power Method.
Input: Matrix An×n and the required numbed of eigenvalues d.
2D mesh of p=q×q processors.
Output: Set of eigenvalues and eigenvectors of A,
Λ0,…,d-1  = λ0,λ1,…,λd-1 and V0,…,d-1=v0,v1,…,vd-1.
(1) Let x be a column-wise distributed vector in Rn.
(2)for  i←0:d-1  do
(3) Initialize x randomly. Processors within the same column use the same seed.
(4) column←true
(5) while  not converged  do
(6)  Compute z=Ax locally.
(7)  if  column=true  then
(8)   Perform column-wise all-reduce to obtain z.
(9)  else
(10)     Perform row-wise all-reduce to obtainz.
(11)end if
(12)column←¬column
(13)x←z
(14)end while
(15)   Computeu=Az as in steps  (6)–(11).
(16)   Replicate entire vectoru and z on each processor.
(17)λi←z·uz·z
(18)vi←zz2
(19)   Deflate local block ofA: A←A-λiviviT.
(20)end forIn general, our approach follows the standard scheme of the power method (lines 3–18), repeatedd times to identify first d largest eigenvalues. After identification of an eigenvalue and its associated eigenvector, the matrix A is deflated—the contribution of the vector is removed from A (line 19). Observe that power method involves nested matrix-vector product (lines 6–11) that under normal circumstances would require parallel vector transposition. However, our parallel implementation benefits from the fact that A is symmetric, hence eliminating need for vector transposition. Indeed, the entire procedure consists of local matrix-vector product followed by all-reduce operation. Here, the reduction operation alternates between columns and rows as required to ensure that vector x is stored properly. Note that the power method is bounded by convergence criteria (line 5). In our case we use one of several popular conditions, which involves checking relative error between the current and previous estimate of the eigenvalue that can be performed every several iterations. We also note that convergence is significantly improved by using a matrix shifting strategy in the form A=A-δI, where δ is a positive number [59].Extracting eigenvalue and eigenvector in iterationi (lines 11–18) depends on vectors u and z, while deflation step involves λi. Therefore, it is advantageous to replicate both u and z in their entirety on each processor. We achieve this with all-to-all communication executed by processors within the same row. This allows us to execute the deflation step in parallel, with each processor updating its local block of matrix A. Thus, the complexity of a single iteration of the power method is On3/p+ts+tm(n/p)log⁡(p), while the deflation step is On2/p+tsp+tmn.To conclude this section we would like to emphasize that our solver operates under the same assumptions as any power method. It requires that the firstd eigenvectors of A are linearly independent, the initial vector x generated in ith iteration is not orthogonal to the eigenvector vi, and finally, the first d eigenvalues are nondegenerate [31]. Note that these conditions are not restrictive and are easily satisfied in the context of dimensionality reduction.


## 4. Experimental Results

To assess scalability of our framework and test its performance in real-life applications, we performed a set of experiments using the TACC Ranger cluster [60]. A single node of this machine is based on AMD processors working at 2.3 GHz and provides 16 cores with 32 GB of DDR2 RAM and 512 KB of L2 cache per core. Nodes are connected by a multistage Infiniband network that offers 1 GB/s bandwidth. To compile all test programs and the framework, we used the Intel C++ 10.1 compiler with standard optimization flags and MVAPICH 1.0.1 MPI implementation. In every test we ran one MPI process per CPU core, which we refer to as processor.


### 4.1. Scalability Tests

In the first set of experiments we measured how problem size influences performance of our solution. We created a collection of synthetic datasets consisting ofn={4096,8192,16384,32768} points with D=10000. Next, we performed Isomap dimensionality reduction using different numbers of processors. Obtained results are summarized in Table 3 and Figure 2.Table 3
Run time in seconds for differentp and varying problem sizes n. Due to memory limitations problem instance with n=32768 could not be solved on less than p=256 processors.
p
n
4096
8192
16384
32768
16
404.63
3492.28
45288.93
—
64
101.72
761.75
6906.64
—
256
33.99
263.24
1655.39
14613.33
1024
39.06
124.19
682.91
3964.65Figure 2
Relative speedup for different problem sizesn.The results show that our framework provides very good scalability for large problem sizes in the entire range of tested processor configurations. The super-linear speedup observed forn=16384 is naturally explained by cache performance. Observe that the dominating computational factors in our framework are operations like matrix-matrix and matrix-vector products, which are well suited to exploit memory hierarchy. In our current implementation, we use direct SIMD-based implementation of these routines. A slightly weaker performance for small problem sizes and large number of processors can be attributed to network latency that offsets computational gains.To further understand how different components of the framework perform, we measured their run time obtained for changing problem sizes. Table4 shows that all modules scale as we would expect based on their theoretical complexity. The most time consuming stages are construction of the auxiliary matrix W for Isomap and normalization. This is not surprising taking into account that both components scale as O(n3) and the parallel Floyd’s algorithm involves n rounds of communication. The abrupt performance decrease in the normalization stage, which can be observed for n=16384, can be attributed to cache performance. Recall that normalization depends on matrix-matrix multiplication and hence is inherently sensitive to data locality. The final remark concerns k nearest neighbors module and eigenvalue solver. The KNN scales linearly with the data dimension D (see Table 5) and both modules can be used as standalone replacements whenever KNN or d largest eigenvalues problem has to be solved.Table 4
Component-wise run time in seconds for varying problem sizes, andp=1024 and D=10000.
n
2048
4096
8192
16384
32768
KNN
0.623
1.389
5.721
22.254
86.706
W in Isomap
9.132
56.517
128.225
457.306
1697.124
Normalize
0.160
0.905
6.526
223.240
2546.11
Eigensolve
0.050
0.155
0.188
0.699
2.838Table 5
Run time in seconds ofKNN component for n=1024 and different number of processors and varying D.
p
D
100
1000
10000
100000
16
0.053
0.530
5.466
92.014
64
0.015
0.115
1.373
22.984
256
0.005
0.027
0.349
5.860
1024
0.002
0.007
0.682
1.875In the final test we compared our parallel eigensolver withSLEPc [55], one of the most popular and widely used libraries providing eigensolvers.SLEPc is an efficient and portable framework that offers an intuitive user interface. In many cases it is the first choice for solving large-scale sparse eigenvalue problems.Table6 shows that our implementation systematically outperformsSLEPc. This can be explained by two main factors: first, unlikeSLEPc, our implementation follows 2D data decomposition scheme, which offers better scalability, and second, we are seeking only the d largest eigenvalues.Table 6
Comparison of our eigensolver withSLEPc. For p=1024SLEPc failed to execute.
p
n
=
1024
n
=
4096
Our solver
SLEPc
Our solver
SLEPc
16
0.0444
4.8159
2.5315
0.7049
64
0.0088
2.1666
0.6056
0.8134
256
0.0705
8.5538
0.1251
2.4143
1024
0.0742
*
0.1320
*
4096
0.0411
N/A
0.2024
10.9992


### 4.2. Using Dimensionality Reduction to Explore Manufacturing Pathways

Solar cells, or plastic solar cells, manufactured from organic blends, that is, a blend of two polymers, represent a promising low-cost, rapidly deployable strategy for harnessing solar energy. While highly cost-effective and flexible, their low power conversion efficiency makes them less competitive on a commercial scale in comparison with conventional inorganic solar cells. A key aspect determining the power conversion efficiency of organic solar cells is the morphological distribution of the two polymer regions in the device. Recent studies reveal that significant improvement in power conversion efficiency is possible through better morphology control of the organic thin film layer during the manufacturing process [29, 61–66]. High-throughput exploration of the various manufacturing parameters, for example, evaporation rate, blend ratio, substrate patterning frequency, substrate patterning intensity, and solvent, can potentially unravel process-morphology relationships that can help tailor processing pathways to obtain enhanced morphologies. Note that such high-throughput analysis results in datasets that are too large to visually inspect for trends and relationships. A promising approach towards unraveling process-morphology relationships in this high-throughput data is via dimensionality reduction. Here, we showcase our parallel framework on this pressing scientific problem. In particular, we focus on using dimensionality reduction to understand the effects of substrate patterning [67, 68], described by patterning frequency and intensity, on morphology evolution. We note that nanotip photolithography patterning of the substrate has shown significant potential to guide morphology evolution [69].The input dataset consists ofn=75150 morphologies. Each morphology is a 2-dimensional snapshot which is vectorized to have dimensionality D=8326. Figure 3 shows several representative final morphologies obtained by varying the patterning frequency, lp, from 0.5 to 1.50, and the intensity of the attraction/repulsion, μ, from 1+1e-6 to 1+8e-4 (in the remaining we omit leading 1 for clarity). In all these cases, the lower surface of the domain is patterned to attract and/or repel specific classes of polymers, thus affecting the morphology. We performed dimensionality reduction on this dataset using p=16384 processors on TACC Ranger. The total run time was 1058 seconds.Figure 3
Snapshots of microstructures representing final morphologies from 50 different manufacturing processes.Figure4 plots the first 10 eigenvalues of the data. Note that the first three eigenvalues, and hence the first three principle components of the data, represent ~70% of the information content of the entire data. We therefore characterize each morphology in terms of this three-dimensional representation.Figure 4
Scree plot of the ten largest eigenvalues (a) and their proportional energy covered (b).
(a)
(b)Figure5 representsall the morphologies in such reduced space. In the plot to the left, the points are color coded according to the patterning frequency used, while in the plot to the right, the points are color coded according to the patterning intensity. This plot provides significant visual insight into the effects of patterning frequency and intensity. There exists a central plane of patterning frequency, where the morphology evolution is highly regulated irrespective of the patterning intensity (lp≤1). This is particularly valuable information as the patterning frequency is much easier to control than patterning intensity from a manufacturing perspective. For patterning frequencies above lp=1, the morphologies are highly sensitive to slight variations in both frequency and intensity. This is also clearly seen in Figure 6, where slight variations in the intensity give dramatically different final morphologies. Notice also that higher intensity does not necessarily give different morphologies. This is a very important insight that allows us to preclude further, potentially expensive, exploration of the phase space of increasing patterning intensity.Figure 5
Morphology evolution as captured by the first three principal components of the original data. Different colors represent different patterning frequencylp (left) and different patterning intensity μ (right).
(a)
(b)Figure 6
Morphology evolution forlp=1.50, and categorization of parametric space.Finally, the low-dimensional plots illustrate the ability to achieve the same morphology using different processing conditions. For instance, in Figure7, we isolate the morphology evolution under two processing conditions that result in identical morphologies. Such correlations, most sensitive regions versus least sensitive regions (Figure 6) and configurations resulting in identical morphologies, are enormously useful as we tailor processing pathways to achieve designer morphologies. This analysis illustrates the power of parallel dimensionality reduction methods to achieve this goal. We defer a comprehensive physics-based analysis of this dataset to a subsequent publication.Figure 7
Multiple pathways of morphology evolution.


## 4.1. Scalability Tests

In the first set of experiments we measured how problem size influences performance of our solution. We created a collection of synthetic datasets consisting ofn={4096,8192,16384,32768} points with D=10000. Next, we performed Isomap dimensionality reduction using different numbers of processors. Obtained results are summarized in Table 3 and Figure 2.Table 3
Run time in seconds for differentp and varying problem sizes n. Due to memory limitations problem instance with n=32768 could not be solved on less than p=256 processors.
p
n
4096
8192
16384
32768
16
404.63
3492.28
45288.93
—
64
101.72
761.75
6906.64
—
256
33.99
263.24
1655.39
14613.33
1024
39.06
124.19
682.91
3964.65Figure 2
Relative speedup for different problem sizesn.The results show that our framework provides very good scalability for large problem sizes in the entire range of tested processor configurations. The super-linear speedup observed forn=16384 is naturally explained by cache performance. Observe that the dominating computational factors in our framework are operations like matrix-matrix and matrix-vector products, which are well suited to exploit memory hierarchy. In our current implementation, we use direct SIMD-based implementation of these routines. A slightly weaker performance for small problem sizes and large number of processors can be attributed to network latency that offsets computational gains.To further understand how different components of the framework perform, we measured their run time obtained for changing problem sizes. Table4 shows that all modules scale as we would expect based on their theoretical complexity. The most time consuming stages are construction of the auxiliary matrix W for Isomap and normalization. This is not surprising taking into account that both components scale as O(n3) and the parallel Floyd’s algorithm involves n rounds of communication. The abrupt performance decrease in the normalization stage, which can be observed for n=16384, can be attributed to cache performance. Recall that normalization depends on matrix-matrix multiplication and hence is inherently sensitive to data locality. The final remark concerns k nearest neighbors module and eigenvalue solver. The KNN scales linearly with the data dimension D (see Table 5) and both modules can be used as standalone replacements whenever KNN or d largest eigenvalues problem has to be solved.Table 4
Component-wise run time in seconds for varying problem sizes, andp=1024 and D=10000.
n
2048
4096
8192
16384
32768
KNN
0.623
1.389
5.721
22.254
86.706
W in Isomap
9.132
56.517
128.225
457.306
1697.124
Normalize
0.160
0.905
6.526
223.240
2546.11
Eigensolve
0.050
0.155
0.188
0.699
2.838Table 5
Run time in seconds ofKNN component for n=1024 and different number of processors and varying D.
p
D
100
1000
10000
100000
16
0.053
0.530
5.466
92.014
64
0.015
0.115
1.373
22.984
256
0.005
0.027
0.349
5.860
1024
0.002
0.007
0.682
1.875In the final test we compared our parallel eigensolver withSLEPc [55], one of the most popular and widely used libraries providing eigensolvers.SLEPc is an efficient and portable framework that offers an intuitive user interface. In many cases it is the first choice for solving large-scale sparse eigenvalue problems.Table6 shows that our implementation systematically outperformsSLEPc. This can be explained by two main factors: first, unlikeSLEPc, our implementation follows 2D data decomposition scheme, which offers better scalability, and second, we are seeking only the d largest eigenvalues.Table 6
Comparison of our eigensolver withSLEPc. For p=1024SLEPc failed to execute.
p
n
=
1024
n
=
4096
Our solver
SLEPc
Our solver
SLEPc
16
0.0444
4.8159
2.5315
0.7049
64
0.0088
2.1666
0.6056
0.8134
256
0.0705
8.5538
0.1251
2.4143
1024
0.0742
*
0.1320
*
4096
0.0411
N/A
0.2024
10.9992


## 4.2. Using Dimensionality Reduction to Explore Manufacturing Pathways

Solar cells, or plastic solar cells, manufactured from organic blends, that is, a blend of two polymers, represent a promising low-cost, rapidly deployable strategy for harnessing solar energy. While highly cost-effective and flexible, their low power conversion efficiency makes them less competitive on a commercial scale in comparison with conventional inorganic solar cells. A key aspect determining the power conversion efficiency of organic solar cells is the morphological distribution of the two polymer regions in the device. Recent studies reveal that significant improvement in power conversion efficiency is possible through better morphology control of the organic thin film layer during the manufacturing process [29, 61–66]. High-throughput exploration of the various manufacturing parameters, for example, evaporation rate, blend ratio, substrate patterning frequency, substrate patterning intensity, and solvent, can potentially unravel process-morphology relationships that can help tailor processing pathways to obtain enhanced morphologies. Note that such high-throughput analysis results in datasets that are too large to visually inspect for trends and relationships. A promising approach towards unraveling process-morphology relationships in this high-throughput data is via dimensionality reduction. Here, we showcase our parallel framework on this pressing scientific problem. In particular, we focus on using dimensionality reduction to understand the effects of substrate patterning [67, 68], described by patterning frequency and intensity, on morphology evolution. We note that nanotip photolithography patterning of the substrate has shown significant potential to guide morphology evolution [69].The input dataset consists ofn=75150 morphologies. Each morphology is a 2-dimensional snapshot which is vectorized to have dimensionality D=8326. Figure 3 shows several representative final morphologies obtained by varying the patterning frequency, lp, from 0.5 to 1.50, and the intensity of the attraction/repulsion, μ, from 1+1e-6 to 1+8e-4 (in the remaining we omit leading 1 for clarity). In all these cases, the lower surface of the domain is patterned to attract and/or repel specific classes of polymers, thus affecting the morphology. We performed dimensionality reduction on this dataset using p=16384 processors on TACC Ranger. The total run time was 1058 seconds.Figure 3
Snapshots of microstructures representing final morphologies from 50 different manufacturing processes.Figure4 plots the first 10 eigenvalues of the data. Note that the first three eigenvalues, and hence the first three principle components of the data, represent ~70% of the information content of the entire data. We therefore characterize each morphology in terms of this three-dimensional representation.Figure 4
Scree plot of the ten largest eigenvalues (a) and their proportional energy covered (b).
(a)
(b)Figure5 representsall the morphologies in such reduced space. In the plot to the left, the points are color coded according to the patterning frequency used, while in the plot to the right, the points are color coded according to the patterning intensity. This plot provides significant visual insight into the effects of patterning frequency and intensity. There exists a central plane of patterning frequency, where the morphology evolution is highly regulated irrespective of the patterning intensity (lp≤1). This is particularly valuable information as the patterning frequency is much easier to control than patterning intensity from a manufacturing perspective. For patterning frequencies above lp=1, the morphologies are highly sensitive to slight variations in both frequency and intensity. This is also clearly seen in Figure 6, where slight variations in the intensity give dramatically different final morphologies. Notice also that higher intensity does not necessarily give different morphologies. This is a very important insight that allows us to preclude further, potentially expensive, exploration of the phase space of increasing patterning intensity.Figure 5
Morphology evolution as captured by the first three principal components of the original data. Different colors represent different patterning frequencylp (left) and different patterning intensity μ (right).
(a)
(b)Figure 6
Morphology evolution forlp=1.50, and categorization of parametric space.Finally, the low-dimensional plots illustrate the ability to achieve the same morphology using different processing conditions. For instance, in Figure7, we isolate the morphology evolution under two processing conditions that result in identical morphologies. Such correlations, most sensitive regions versus least sensitive regions (Figure 6) and configurations resulting in identical morphologies, are enormously useful as we tailor processing pathways to achieve designer morphologies. This analysis illustrates the power of parallel dimensionality reduction methods to achieve this goal. We defer a comprehensive physics-based analysis of this dataset to a subsequent publication.Figure 7
Multiple pathways of morphology evolution.


## 5. Conclusions

In this work we illustrate a systematic analysis of dimensionality reduction techniques and recast them into a unified view that can be exploited by dimensionality reduction algorithm designers. We subsequently identify the common computational building blocks required to implement a spectral dimensionality reduction method. We use this insight to design and implement a parallel framework for dimensionality reduction that can handle large datasets and scales to thousands of processors. We demonstrate the capability and scalability of this framework on several test datasets. We finally showcase the applicability and potential of the framework towards unraveling complex process-morphology relationships in the manufacturing of plastic solar cells.



---

*Source: 180214-2015-03-10.xml*