Data analytics using canonical correlation analysis and Monte Carlo simulation (2024)

Canonical correlation analysis

We begin with an overview of the CCA methodology. Suppose that one has a set of input variables \(\left\{ {{x_1},{x_2} \cdots ,{x_{{N_i}}}} \right\}\) and corresponding output variables \(\left\{ {{y_1},{y_2} \cdots ,{y_{{N_{\rm{o}}}}}} \right\}\), where N _i (N _o) are the number of input (output) variables. One seeks linear combinations (known as canonical variates) \(V = \mathop {\sum}\nolimits_i {{\alpha _i}{x_i}} \) and \(W = \mathop {\sum}\nolimits_i {{\beta _i}{y_i}} \), where α _i and β _i are known as canonical weights, such that these combinations are maximally correlated. If one regards these input and output variables as the components of two vectors, \(\vec x\) and \(\vec y\) _, respectively, then one can define a covariance matrix whose (i, j)-th element is cov(x _i , y _j ). CCA begins with the construction of this matrix (or the associated correlation matrix) Σ.^{14, 15} The matrix elements are estimated from a sample covariance matrix that is calculated from the data from the M experiments.¹⁴ Then, Σ may be written in block form as

The aim of CCA is to determine {α _i } and {β _i } such that the correlation between W and V is maximized. Since the covariance between W and V is cov (W, V) = \({\vec \beta ^T}\mathop {\sum}\nolimits_{\vec y\vec x} {\vec \alpha } \), one seeks to maximize the correlation

where the denominator is the product of the variances of W and V, respectively. This may be accomplished by finding the solution of a system of hom*ogeneous equations or, equivalently, by finding the eigenvectors and corresponding eigenvalues of two operators σ ₁ and σ ₂. The two operators are constructed from products of matrix blocks of Σ and given by

The eigenvalues of each operator are the same, and the square roots of these eigenvalues are the canonical correlations. Moreover, the corresponding eigenvectors may be used to determine the α _i and β _i and, thereby, the canonical variates. More specifically, if \({\vec \alpha ^*}\) and \({\vec \beta ^*}\) are eigenvectors of σ ₁ corresponding to its maximum eigenvalue, then V ^* and W ^* are the desired canonical variates that maximize the correlation.¹⁴

CCA with Monte Carlo simulation (CCAMC)

While the CCA is extremely useful in highlighting linear relationships among input and output variables, it may be that a non-linear model of the data is more appropriate in some circ*mstances. More formally, one wishes to find some function of the subset of input variables, {x ₁, x ₂ …, x _n }, where n < N _i , such that the largest eigenvalue λ _max(x ₁, x ₂, x ₃, x ₄, …, x _n ) in the spectrum of σ ₁ (or σ ₂) is maximized. In practice, this is often tedious, especially for n large. With this in mind, we outline here a strategy using Monte Carlo simulation¹⁶ to identify appropriate non-linear combinations of variables for a given problem, and then validate this approach for two test cases below. The CCAMC strategy has the virtue that it can be straightforwardly implemented in most cases with relatively little computational cost. In addition, parallel computations can be performed to explore relatively wide regions of parameter space, as described below.

The procedure is as follows. One identifies, perhaps from a CCA analysis of data, the aforementioned input variables {x ₁, x ₂ …, x _n }, where n < N _i , that contribute significantly to the input canonical variate. It may be suspected that some (perhaps non-linear) function of these variables is, in fact, strongly correlated with the output, but the functional form may not be readily apparent. In this case, it is sensible to consider families of trial functions that ideally form a basis in some infinite-dimensional vector space. For example, for the vector space of multivariate polynomial functions one may parametrize the unknown function f(\(\vec x\)) in terms of Kolmogorov-Gabor polynomials¹⁷

where c ₀, c _i , etc., are unknown expansion coefficients. This expansion is used extensively in, for example, polynomial neural networks to identify non-linear relationships between input and output variables.¹² Alternatively, if variable ratios are more appropriate, one may employ a multivariate generalization of Padé approximants.^18,19,20 For example, in the trivariate case, one can write

where m is an expansion parameter and the c _ijk and d _ijk are unknown expansion coefficients. For more than three variables, Eq.5 can, of course, be generalized.

Upon replacing an input variable by one of the above representations, one then seeks the set of variables that maximize the largest eigenvalue, λ _max, of either σ ₁ or σ ₂. This may be accomplished by employing Monte Carlo simulation using λ _max as an objective function. In brief, one chooses values for the variable set and then calculates λ _max from a CCA analysis. These starting values are most conveniently chosen from the representations given in either Eq.4 or 5 with a set of arbitrary coefficients. This selection and a subsequent CCA analysis leads to a corresponding starting value for λ _max. As in a conventional Monte Carlo simulation, one constructs a Markov chain of states. One begins by defining a ‘‘box” in coefficient space for trial variations of the coefficients. For convenience, the size of the box can be chosen so that approximately 50% of the variations are accepted. Such trial variations lead to concomitant changes in the CCA operators given in Eq.3, and thereby to changes in their spectra. If a change in the given coefficient set leads to an increase (decrease) in λ _max, then this set is retained (rejected). This procedure is repeated for many iterations, usually from a different initial state, and may be performed in parallel. As is intuitively reasonable, the selection of an initial state determines the convergence to the maximum eigenvalue, and the use of many different initial states mitigates against trapping in configuration space. In practice, we have found that convergence is achieved after a few thousand iterations for the examples considered here. In effect, we are performing the analog of a zero-temperature Monte Carlo simulation. It is also possible to define an effective temperature and to perform simulated annealing, as indicated below, if one is still concerned about trapping.

To determine whether a given pair of variates is indeed correlated, one ascertains the significance of the results from a hypothesis test in which one tests the null hypothesis that a given pair of variates is uncorrelated. The results of this test are then quantified in terms of a statistic, such as Wilks lambda,²¹ and an associated p-value. For relatively small p-values, there is contrary evidence of a correlation and the corresponding variate pair is acceptable.

In practice, there are several issues to be addressed when implementing this methodology. First, in calculating σ ₁ and σ ₂ using Eq.3 one must calculate inverses of the block matrices \({\Sigma _{\vec X\vec X}}\) and \({\Sigma _{\vec Y\vec Y}}\). In some cases, however, these matrices may be ill-conditioned at some point during the Monte Carlo simulation. To remedy this situation, one can either regularize the block matrices²² or, alternatively, replace the inverse by the Moore-Penrose pseudoinverse.²³ For a given matrix Ξ and appropriate identity matrix I, the pseudoinverse Ξ ⁺ satisfies the relation that the sum of the squares of the matrix elements of ΞΞ ⁺ − I is a minimum. We have chosen to employ this latter approach and found it to be a satisfactory remedy. Second, it is possible to become trapped near local extrema during the simulation. It is, therefore, advisable to perform many simulations, each comprising approximately 15,000–20,000 iterations, possibly starting from different initial conditions to explore a wide range of parameter space. Such simulations may, of course, be performed in parallel. In addition, simulated annealing may be employed to mitigate this problem.¹⁶ This approach would require the introduction of a fictitious temperature that would permit fluctuations. These fluctuations imply that states that minimize λ _max would sometimes be accepted and the system would escape from a local trap. The subsequent extraction of thermal energy would then allow the system to relax to (perhaps) another extremum. We have not employed annealing here as it has not been necessary in the various test cases that we have considered. Third, one may wish to explore the robustness of this methodology. To verify robustness, we have systematically eliminated less impactful input variables from the analysis, and there has been little to no change in the correlation coefficient, etc. Finally, the selection “a priori” of a best functional form (e.g., Eqs.4 and 5) is not crucial here. Given that the CCAMC methodology may be implemented in parallel, it is straightforward to examine many different functional forms (including, of course, others not given here) at the same time to identify which is best. This flexibility is another strength of this methodology.

Applications

To validate the CCAMC methodology outlined above, we consider here two experimental applications of relevance in applied physics and materials science and engineering, namely: (1) quantifying abnormal grain growth (AGG) in ceramic oxides, and (2) relating the detailed microstructural features in CuInSe₂ thin films, the so-called grain-boundary character distribution (GBCD), to the electrical and optoelectronic properties of these films. In each case, one seeks to identify a combination of input variables that is strongly correlated with an observable output.

Application 1: Abnormal grain growth

AGG occurs in a polycrystal when a small group of grains, typically having anisotropic boundary energies or mobilities, becomes relatively large by growing into the surrounding matrix of ‘‘normal” grains.^{24, 25} It is a ubiquitous phenomenon in many systems, especially thin films, and occurs, in some cases, due to the presence of impurity excesses (e.g., Ca or Si in alumina).²⁶ AGG is found in both metallic and ceramic material microstructures and may have a significant impact on their electromechanical properties. Figure1 shows a prototypical microstructure, as obtained from electron backscatter diffraction (EBSD), highlighting abnormal grains in a background of ‘‘normal grains” for an alumina sample.

In a recent paper, Lawrence et al.²⁷ introduced a series of metrics that are useful in quantifying various microstructural aspects of AGG, and then employed a CCA to highlight the influence of processing variables on these metrics in doped aluminas.²⁷ These metrics highlight specific features of the joint probability density function (pdf) p(G, a) of grain size, G, and aspect ratio, a, for a given microstructure, and, therefore, necessarily focus on grains in the tail of the distribution. More specifically, the metrics reflect the extreme events (i.e., large grains) that determine the shape of the tail in terms of conditional expectation values, and will be denoted here by {ϕ ₁, ϕ ₂, …, ϕ ₅}. The object of the study by Lawrence et al. was to understand how the processing parameters for the doped aluminas, including compositions c (for such dopants as MgO, CaO, Na₂O, and SiO₂), temperature T and time t influence the aforementioned metrics. It should be noted that impurities often segregate to grain boundaries or other defects and may thereby affect thermomechanical and kinetic properties.^28,29,30

In this study, we wish to identify which (possibly non-linear) combination(s) of processing variables is (are) most correlated with microstructural abnormality. To accomplish this aim, after the spark-plasma sintering of specialty alumina powders, we compiled statistics for microstructures associated with grain growth at some annealing temperature. The original work by Lawrence et al. was based on 33 samples; in this study we have augmented the original data set to obtain 68 independent samples. Each of the 68 independent specialty powders had a unique combination of processing variables, and the final microstructures had grain populations in the range of approximately 1000–9000 grains. For this analysis, it is convenient to define a set of input (i.e., processing) variables, \(\vec x\), and a set of output (i.e., metric) variables, \(\vec y\), as

where the composition subscripts denote dopant type.

A CCA was performed using the data set described above. From this analysis, one finds a single pair of canonical variates having a (maximum) correlation coefficient of 0.59 and a p-value of 0.10 for the hypothesis test using Rao’s F-test.¹⁴ In Fig.2a the corresponding canonical variates, V and W, are plotted for each of the data points, along with the associated regression line.

To find variates having a stronger correlation, one can either use physical intuition to identify a new combination of variables or employ, for example, the CCAMC method. Based on the earlier work of Lawrence et al., one such variable that reflects perceived correlations, and identified via physical intuition, is the ratio \(r = {c_{{\rm{MgO}}}}{\rm{/}}\left( {{c_{{\rm{CaO}}}} + {c_{{\rm{Si}}{{\rm{O}}_2}}}} \right)\). With this guidance and after augmenting the input variable set with r, a CCA analysis yields a pair of canonical variates having a (maximum) correlation coefficient of 0.69 with a p-value of 0.02. To determine whether the CCAMC procedure can yield further improvement, consider the Padé functional form given in Eq.5 with the variable set \(\left\{ {{c_{{\rm{MgO}}}},{c_{{\rm{CaO}}}},{c_{{\rm{Si}}{{\rm{O}}_{\rm{2}}}}}} \right\}\) and expansion exponent m = 1. These approximants replace r. These input variables were identified from the CCA as making significant contributions to the input canonical variate. Moreover, the inverse relationship between some pairs of variables suggests the Padé functional form. After approximately 100 simulations, each comprising 15,000 iterations, expansion coefficients were found yielding a maximum observed correlation coefficient of 0.86 with an associated p-value of 4 × 10⁻⁸. In Fig.2b, the corresponding canonical variates are plotted for each of the data points, along with the associated regression line and a 95% confidence interval for the data. Thus, with relatively little computational effort, one can obtain well-correlated canonical variables using a trial function with the CCAMC methodology. These new variables facilitate experimental planning, as described below.

Application 2: Electrical and optoelectronic properties of thin films

Thin-film solar cells, such as those based on, for example polycrystalline CuInSe₂ absorbers, are of considerable technological interest given their relatively high conversion efficiencies.³¹ Indeed, many workers have sought to describe the role of grain boundaries in determining device performance in these systems.³² In a recent paper, Abou-Ras et al. examined the grain-boundary character of these films in relation to their electrical and optoelectronic properties for a relatively large number of grain boundaries.³³ More specifically, EBSD was employed to acquire microstructural data that, in conjunction with an electron-beam-induced current (EBIC) analysis and a cathodoluminescence (CL) image study, permitted the correlation of boundary type with local electrical and optoelectronic properties. Figure3 shows the distribution of disorientation axes for grain boundaries in polycrystalline CuInSe₂, in multiples of a random distribution (MRD), presented as stereographic projections.

Several parameters were used to quantify the GBCD as extracted from EBSD maps, including, for a given boundary: the grain-boundary disorientation angle ψ, the three Rodrigues vector \(\overrightarrow R \) components of the disorientation and scalar measures of closeness k to the 〈100〉, 〈110〉, and 〈111〉 crystallographic directions. Our aim here is, as above, to to identify which (possibly non-linear) combination(s) of GBCD variables is (are) most correlated with electrical and optoelectronic properties. For this purpose, it is sensible to define sets of input and output variables, respectively, as

where the subscripts on k denote crystallographic directions and where ebic and cl are the EBIC and CL signals, respectively.

To identify those input and output variables that are most important in connecting microstructure to properties, one first performs a CCA using the variables given in Eq.7. For concreteness, we consider a data set comprising 104 samples for non-twin boundaries. From this analysis, one finds a single pair of canonical variates that are maximally, though relatively weakly correlated, having a correlation coefficient of 0.26 and a large p-value of 0.88 for the hypothesis test using Rao’s F-test.¹⁴ (Given this p-value, one cannot reject the possibility that these variables are uncorrelated.) The canonical weights for the input and output variables are given in Table1. In Fig.4a, the corresponding canonical variates, V and W, are plotted for each of the data points, along with the associated regression line.

From Table1 it is evident that several of the input variables, given their relatively small canonical weights, make little contribution to the input canonical variate. This situation is, in fact, ideally suited to CCAMC as the trial function need only depend on a few input variables. We identify three input variables, namely the components of \(\vec R\), that have significant canonical weights and ask whether the addition of some (possibly) non-linear function of these variables would lead to a stronger correlation with the output variables. Consider, for example, the Padé functional form given in Eq.5 with the variable set {R ₁, R ₂, R ₃} and expansion exponent m = 1. After approximately 150 CCAMC simulations, each comprising 15,000 iterations, expansion coefficients were found yielding a maximum observed correlation coefficient of 0.56 with an associated p-value < 10⁻³. In Fig.4b, the corresponding canonical variates are displayed for each of the data points, along with the regression line and a 95% confidence interval for the data. Clearly, this procedure has produced canonical variates with a substantially increased degree of correlation; moreover, the corresponding p-value gives one confidence that this correlation is real.