# Customized Gaussian process for representing polycrystalline texture

Bingqian Li, Piotr Breitkopf, Ludovic Cauvin

► To cite this version:

Bingqian Li, Piotr Breitkopf, Ludovic Cauvin. Customized Gaussian process for representing polycrystalline texture. Computer Methods in Applied Mechanics and Engineering, 2025, 440, pp.117934. 10.1016/j.cma.2025.117934. hal-05021289

HAL Id: hal-05021289

https://utc.hal.science/hal-05021289v1

Submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Customized Gaussian process for representing polycrystalline texture

Bingqian Li<sup>ID</sup>, Piotr Breitkopf<sup>ID</sup>, Ludovic Cauvin<sup>ID *</sup>

Université de Technologie de Compiègne, Alliance Sorbonne Université, Centre de recherche Royallieu, Roberval, CS 60319, 60203 Compiègne
Cedex, France

## GRAPHICAL ABSTRACT

### Customized Gaussian Process for Representing Polycrystalline Texture

#### Highlights

Pole Density Function
reconstruction with customized
Gaussian Process Regression

Function's periodicity, non-negativity ensured

Study of the Confidence Interval
(CI) of plausible reconstructions

#### GPR customization

- Periodic distance measures: intrinsic/geodesic/Euclidean

- Stationary kernel functions: Squared exponential / De La Vallée Poussin / Matérn

- Hyperparameters optimization: Leave-One-Out / Marginal likelihood maximization

- Data trasnformation: Logarithmic / log-linear

#### Synthetic examples

Three synthetic pole figure
examples constructed with
quadratic splines (data provided)

Both top and under views of the
pole figures presented

<figure><img src="image_5.png" /></figure>

One of the three synthethic
reference pole figures

<figure><img src="image_6.png" /></figure>

1500 random
measurement points

#### Result (best configurations)

<figure><img src="image_7.png" /></figure>

Upper bound of 95% CI
e = 14.52%

<figure><img src="image_8.png" /></figure>

A posteriori mean of GPR,
as reconstruction result
e = 2.03%

<figure><img src="image_9.png" /></figure>

Lower bound of 95% CI
e = 6.37%

#### Comparisons

<figure><img src="image_10.png" /></figure>

Spherical harmonics expansion
approach (MTEX), order of 28
e = 8.04%

<figure><img src="image_11.png" /></figure>

Linear interpolation (MTEX)
e = 11.21%

#### Conclusion

Log-linear transformation ensures
non-negativity of mean and CI.

Choices of geodesic distance, Matérn kernel (ν=5/2), marginal likelihood maximization give smallest reconstruction error.

Bingqian Li, Piotr Breitkopf, Ludovic Cauvin (2025)

## ARTICLE INFO

### Keywords:

Microstructure

Texture

Pole figures

Digital material representation

Gaussian Process Regression

Kriging

## ABSTRACT

A customized Gaussian Process Regression (GPR) model is developed to reconstruct Pole Density Functions in texture analysis. The model integrates spherical-periodic distance measures with conventional stationary kernels, adapting the GPR framework to capture localized texture features. A key contribution is the introduction of a log-linear data transformation, which enforces the non-negativity of both interpolated function values and stochastic intervals, ensuring physically meaningful reconstructions. To assess the model's effectiveness, a systematic investigation examines the impact of distance measures, kernel selection, and hyperparameter

<footnote>* Corresponding author.
E-mail address: ludovic.cauvin@utc.fr (L. Cauvin).</footnote>
optimization using synthetic texture datasets, provided as part of this work, with evaluations focusing on reconstruction accuracy, feature preservation, and uncertainty quantification in comparison to the conventional spherical harmonics approach. GPR with a log-linear transformation, geodesic distance, and a Matérn $v=5/2$ kernel, shows promise in achieving higher accuracy than traditional spherical harmonics for reconstructing non-negative pole density functions, while additionally providing confidence intervals for uncertainty quantification.

# 1. Introduction

The paper introduces targeted, domain-specific adaptations of Gaussian Process Regression (GPR) [1–3], including a log-linear transformation to enforce non-negativity, geodesic distance for spherical periodicity, stationary kernel selection, and volume preservation—that collectively redefine its application to texture analysis, enabling accurate, uncertainty-aware reconstructions of Pole Density Functions beyond the capabilities of conventional GPR [4] or harmonic methods

The physical properties of polycrystalline materials, particularly their anisotropy, are closely linked to the crystallographic orientation distribution of individual grains, commonly referred to as *texture* [5–7]. Electron Backscatter Diffraction (EBSD) directly provides the 3D orientations of crystallites in the form of Euler angles. In the case of plastic deformation, which is highly dependent on crystal slip systems, the orientation distribution of slip planes is of particular interest and is often represented by *pole figures*. These figures can be derived by projecting the texture or obtained directly through X-ray or neutron diffraction methods.

However, these experimental techniques only capture the orientations of thousands of independent points on the sample surface (depending on resolution settings). For quantitative analysis, it is common to construct a continuous function model based on the experimental data. The Orientation Distribution Function (ODF) represents texture in Euler angles space, showing the volume fraction of crystallites with a specific 3D orientation. The function representing the pole figure, defined on a unit sphere, provides the volume fraction of crystallites with their slip planes oriented in a particular direction and is referred to as the Pole Density Function. The Pole Density Function can also be obtained by projecting the ODF, and various methods have been proposed for reconstructing the ODF using Pole Density Functions [8–10]. Directly expressing these functions is typically infeasible. Instead, for decades, the approach of series expansion of harmonic functions is used to approximate these functions: generalized spherical harmonics for the ODF and spherical harmonics for the Pole Density Function.

This harmonic expansion approach has several advantages. It effectively handles spherical symmetry and crystallographic symmetry, as spherical harmonics are naturally periodic on a sphere. With this method, the ODF or Pole Density Function can be represented by a series of coefficients and a pre-established spherical harmonics library, which requires minimal storage. Moreover, the precision of the function representation can be controlled by limiting the number of terms in the series, resulting in a compact mathematical form [11]. Additionally, integration calculations are relatively straightforward with this method.

However, harmonic expansion, while effective for symmetric representations, is less suitable for asymmetrical data or sharply localized distributions. Moreover, the need for a pre-defined function library restricts its adaptability, and this method does not allow the incorporation of prior knowledge about the data.

Present work considers GPR to construct the Pole Density Function from experimental results for three key reasons:

- mathematical convenience: GPR provides closed-form expressions for predictions and uncertainty, which offer practical advantages such as ease of interpretation and analysis;

- physical adaptability: The kernel function in GPR can be tailored to fit raw data characteristics, adapting to symmetric, asymmetric, or localized datasets. Experimental errors are also accounted for through a noise hyperparameter. [12,13];

- statistical potential: as a stochastic process, GPR enables statistical evaluations [14], providing insights into the variance and likelihood of the model as well as stochastic sampling of Pole Density Functions.

Additionally, compared to other machine learning methods, GPR provides both accurate predictions and inherent uncertainty estimates, which are essential for assessing the reliability of texture representation. While methods such as random forests and neural networks are powerful, they do not naturally offer the same level of uncertainty quantification. This makes GPR particularly suitable for our application, where uncertainty assessment is a key consideration.

The paper is organized as follows: Section 2 introduces the fundamental concepts of texture representation, focusing on the ODF and its relation to pole figure data [15]. This section provides the mathematical framework and normalization conditions essential for texture analysis. Section 3 presents GPR as an alternative to traditional reconstruction methods, explaining its mathematical formulation, kernel functions, and the optimization of hyperparameters, with a particular emphasis on enforcing physical constraints such as non-negativity and volume preservation. Section 4 describes the customization of GPR for pole density function reconstruction, detailing the selection of kernel functions, incorporation of different distance measures, and the role of data transformations to improve accuracy. Section 5 applies the proposed framework to synthetic textures, presenting the numerical protocol, error measures, and a thorough comparison of results under different configurations. This section also includes hyperparameter optimization results and an analysis of how various distance measures and transformations affect reconstruction performance, with a comparison to other techniques like spherical harmonics and linear reconstruction. Section 6 concludes the paper by summarizing the findings and highlighting the benefits of GPR for texture representation, along with suggestions for future work.
<figure><img src="image_12.png" /><figcaption>Fig. 1. Pole Figure representation of a crystalline orientation $g = \{\phi_1, \Phi, \phi_2\}$ for a pyramidal plane in Hexagonal Close-Packed (HCP) lattice.</figcaption></figure>

## 2. Texture and orientation distribution function representation

Considering a polycrystalline material, the ODF, $f(g)$, describes the volume fraction $dV$ of all crystallites possessing a specific crystalline orientation $g$ within an angular element $dg$, relative to the total sample volume $V$, also referred to the density

$$ \frac{dV_g}{V} = f(g)dg. \qquad (1) $$

In EBSD analysis, an orientation $g$ is defined as a triplet of Euler angles $\phi_1, \Phi$ and $\phi_2$, following the $zx'z''$ convention [8], and referring to three consecutive rotations of the crystal with respect to the laboratory reference frame in $\mathbb{E} = [0, 2\pi) \times [0, \pi] \times [0, 2\pi)$. Consequently, the ODF satisfies the normalization condition

$$ \oint_{\mathbb{E}} f(g)dg = 1. \qquad (2) $$

When studying the orientation distribution of a slip plane family $h$, the ODF defined in $\mathbb{E}$ can be projected to the Pole Density Function, defined on the unit spherical surface $S^2 = \{x \in \mathbb{R}^3 : \|x\| = 1\}$.

Consider for example a crystalline orientation $g = \{\phi_1, \Phi, \phi_2\}$ applied to a Hexagonal Close-Packed (HCP) lattice, and focus on one of its pyramidal slip planes. The orientation of the slip plane, represented by the normal vector $n$, is transformed by the rotations $g$, resulting in a new orientation $n'$ (Fig. 1). This transformed orientation can also be represented by a point $y = \{az, el\}$ on $S^2$, commonly referred to as a pole, expressed in terms of azimuth ($az$) and elevation ($el$) angles.

In fact, there is a region in the orientation space $\mathbb{E}$ that all rotations within this region correspond to the same pole $y$. Physically, this refers to as the crystal rotating about $y$. Denote this rotational angle as $\gamma$, then over a full rotation from $\gamma = 0$ to $2\pi$, all the crystalline orientations forming this region contribute to the pole. By integrating the ODF density $f(g)$ over this region, the Pole Density Function $P_h(y)$, defined on $S^2$, can be obtained. For a specific slip plane family $h$, $P_h(y)$ is expressed as

$$ P_h(y) = \frac{1}{2\pi} \cdot \int_{y=0}^{2\pi} f(g)d\gamma, \qquad (3) $$

with the normalization condition

$$ \frac{1}{4\pi} \cdot \oint_{S^2} P_h(y)dy = 1. \qquad (4) $$

The visualization of $P_h(y)$ is commonly referred to as a pole figure.

The value of both the ODF and the Pole Density Function indicating a density, their value has to be non-negative

$$ \left\{ \begin{array}{l} f(g) \ge 0, \quad \forall g \in \mathbb{E}, \\ P_h(y) \ge 0, \quad \forall y \in S^2. \end{array} \right. \qquad (5) $$

To allow for quantitative analysis and comparison of texture, the ODF and the Pole Density Function are constructed from the discrete experimental datasets $D = \{(y_i, d_i) | i = 1, \dots, N\}$ of $N$ orientations $y_i$, corresponding to densities $d_i$. $N$ is determined by the precision settings of the measurement equipment, usually at the order of thousands. The Pole Density Function is conventionally constructed by spherical harmonic series expansion [8].

Fig. 2(a) gives an example of an authentic EBSD observation of a zinc alloy, in which different colors represent different crystalline orientations. Fig. 2(b) shows the projected pole figures of basal and prismatic slip plane families.
<figure><img src="image_13.png" /><figcaption>Fig. 2. EBSD observation and pole figures of zinc alloy [16]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Such EBSD observations give directly the ODF. Based on plasticity theories [17,18], this texture information is essential for mean-field homogenization, as typically applied in the Visco-Plastic Self-Consistent (VPSC) model [19–21] or in other analytical models used for homogenization and damage studies [22,23]. Alternatively, the ODF can be determined through the inversion of Pole Density Functions, as proposed in [8,10].

On the other hand, as based on statistical information, ODF can serve as a numerical basis for digital material representation [24,25]. By attributing texture characteristics to artificial microstructures [26–28], it is possible to generate physically admissible Representative Volume Elements (RVEs) [29]. The size of the RVE must be determined [30,31] before performing full-field homogenization of polycrystalline materials [32–34].

The aim of present work is to interpolate a continuous Pole Density Function respecting Eqs. (4) and (5), associated with Confidence Intervals (CI).

### 3. Gaussian process regression

This work is grounded in the hypothesis that the pole density function follows a Gaussian process [1]

$$P_h(y) \sim GP(m(y), k_{\theta}(y, y')) \quad (6)$$

where $m(y)$ and $k_{\theta}(y, y')$ are the mean and kernel functions, respectively. The subscript $\theta$ denotes the set of hyperparameters.

Application of the GPR considers the a posteriori mean density, denoted as $\tilde{d}(y)$, depending on the observed data $D$, as the regression result,

$$\tilde{d}(y) = c^{\top} C^{-1} d \quad (7)$$

with the a posteriori variance

$$\mathbb{V}[d(y)] = k_{\theta}(y, y) - c^{\top} C^{-1} c \quad (8)$$

Here, $d$ is the vector of observed density values, $C$ is the covariance matrix of data points, and $c$ records the covariance values between the predictor point $y$ and all observations $y = (y_1, y_2, \dots, y_N)^{\top}$,

$$C_{ij} = k_{\theta}(y_i, y_j) + \delta_{ij} \sigma_n^2 \quad i, j = 1, \dots, N; \quad (9)$$

$$c_i = k_{\theta}(y_i, y_i) \quad i = 1, \dots, N, \quad (10)$$

with $\delta_{ij}$ the Kronecker delta and $\sigma_n$ the measurement error.

For a set of observations $D$, the log marginal likelihood of the model is given by

$$\mathcal{L}(D, \theta) = -\frac{1}{2} d^{\top} C^{-1} d - \frac{1}{2} \ln |C| - \frac{N}{2} \ln 2\pi. \quad (11)$$
## 3.1. Hyperparameters' optimization methods

The optimization of the kernel function's hyperparameters is a critical step in GPR. This study considers two hyperparameter optimization methods: Leave-One-Out (LOO) and the maximization of marginal likelihood.

The LOO method excludes one data point and uses the remaining points to predict the value at the missing location. This prediction is recorded as

$$ \mathbf{d}_{i}^{\text{LOO}} = \tilde{d}_{-i}(y_{i}), \qquad (12) $$

where $\mathbf{d}_{i}^{\text{LOO}}$ is the ith component of the LOO prediction vector $\mathbf{d}^{\text{LOO}}$, and $\tilde{d}_{-i}(y)$ represents the a posteriori mean (see Eq. (7)) at the location y, computed without the data point $(y_i, d_i)$ in D.

This procedure is repeated for all N points in D, while keeping the hyperparameters $\theta$ fixed. The LOO reconstruction error is then calculated by comparing the predicted values with the actual values

$$ e_{\text{LOO}}(D, \theta) = \frac{\|\mathbf{d}^{\text{LOO}} - \mathbf{d}\|}{\|\mathbf{d}\|}. \qquad (13) $$

The hyperparameters $\theta$ are optimized by minimizing this error, expressed as

$$ \theta_{\text{opt}} = \arg \min_{\theta} (e_{\text{LOO}}(D, \theta)) \qquad (14) $$

The marginal likelihood, depending on the experimental dataset D and the hyperparameters $\theta$, represents the probability of the model given the experimental data.

To optimize $\theta$, the negative log marginal likelihood Eq. (11) is minimized, a process equivalent to maximizing the marginal likelihood.

$$ \theta_{\text{opt}} = \arg \min_{\theta} (-\mathcal{L}(D, \theta)) \qquad (15) $$

While marginal likelihood is consistent with mathematical grounds of GPR, it relies heavily on prior assumptions and may perform suboptimally if the model is misspecified. In contrast, LOO directly estimates the model's predictive performance by approximating how it generalizes to unseen data. However, it can be computationally intensive and sensitive to data sparsity.

Given the complexity of our model and the potential sensitivity to prior choices, it is crucial to assess both methods to ensure robust hyperparameter optimization. By comparing them, we aim to balance computational efficiency and predictive accuracy while mitigating the risk of model excessive smoothing or overfitting.

## 3.2. Volume and non-negativity constraints

Denoting $\alpha = C^{-1}\mathbf{d}$, Eq. (7) becomes a linear combination of kernels

$$ \tilde{d}(y) = c^{\top}(y)\alpha. $$

Satisfying the conditions Eqs. (4) and (5), the goal is to find a coefficient vector $\tilde{\alpha}$ that is as close to $\alpha$ as possible under the constraints of volume conservation Eq. (4) and of non-negativity Eq. (5)

$$ \begin{array}{ll} \text{minimize} & \|\tilde{\alpha} - \alpha\| \\ \text{subject to} & \int_{S^2} c(y)^{\top} \tilde{\alpha} \, dy = 4\pi \\ & c(y)^{\top} \tilde{\alpha} \ge 0, \forall y \in S^2 \end{array} \qquad (16) $$

To adapt GPR to the problem of Pole Density Function representation, following section proposes successively kernel selections, choice of appropriate distance measure, variable transformation to ensure non-negativity, and scaling for texture volume preservation.

## 4. GPR customization

The key points of the GPR customization are presented in this section.

Gaussian Process Eq. (6) requires the choice of the a priori mean function $m(y)$ and of the kernel function $k_{\theta}(y, y')$. In the case of pole density functions, the choice of $m(y) \equiv 0$ is motivated by the spherical periodicity, and by the hypothesis that densities should extrapolate to zero far from the data.

While a periodic kernel (ex. a sinusoidal) could be considered, Section 4.1 presents a range of isotropic and stationary kernels with periodicity to be ensured by the choice of an appropriate distance measure among the functions introduced in Section 4.2.

Section 4.3 ensures the non-negativity of density values as well as their CI. The non-negative nature of a density function is guaranteed by a specific data transformation of the density values.

Finally, Section 4.4 provides the scaling method to preserve the texture volume.
<table><caption>**Table 1**
Four isotropic, stationary kernel functions, $r = dist(y, y')$.</caption><thead><tr><th>Kernel</th><th>Function</th></tr></thead><tbody><tr><td>Squared exponential (SE)</td><td>$k_{SE}(y, y') = \sigma_f^2 \exp(-\frac{r_f^2}{\ell^2})$</td></tr><tr><td>De La Vallée Poussin (DLVP)</td><td>$k_{DLVP}(y, y') = \sigma_f^2 (\frac{1+\cos(r)}{2})^\kappa$, where $\kappa = \frac{\ln(\frac{1}{\ell})}{2\ln(\cos(r))}$</td></tr><tr><td>Matérn, $\nu = 3/2$</td><td>$k_{\nu=3/2}(y, y') = \sigma_f^2 (1 + \frac{\sqrt{3\ell}}{\ell}) \exp(-\frac{\sqrt{3\ell}}{\ell})$</td></tr><tr><td>Matérn, $\nu = 5/2$</td><td>$k_{\nu=5/2}(y, y') = \sigma_f^2 (1 + \frac{\sqrt{5\ell}}{\ell} + \frac{5\ell^2}{3\ell^2}) \exp(-\frac{\sqrt{5\ell}}{\ell})$</td></tr></tbody></table>

### 4.1. Kernel choice

Four isotropic, stationary kernels are tested and compared in Table 1, with the kernel depending on two orientations $y$ and $y'$ expressed in terms of a distance function $r = dist(y, y')$.

The hyperparameters $\theta = \{\sigma_f, \ell\}$, representing the standard deviation and length scale, must be optimized to reflect the specific nature of the data.

The choice of kernel function depends on the characteristics of the data. For example, the Squared Exponential (SE) kernel is frequently used for its simplicity and the smoothness it provides in reconstructed functions. However, Stein [35] suggests that this level of smoothness may not be realistic for many physical processes.

The Matérn kernel family allows adjustable smoothness through the parameter $\nu$, where smaller values produce rougher reconstructions. When $\nu \to +\infty$, the Matérn kernel becomes equivalent to the SE kernel. Two commonly used configurations are $\nu = 3/2$ and $\nu = 5/2$.

The De La Vallée Poussin (DLVP) kernel is distinct from the other kernels in that it can be expressed exactly using a finite number of Fourier coefficients. This links it to the spherical harmonics approach, which is commonly applied in texture analysis. Additionally, the DLVP kernel is periodic due to its inclusion of cosine, which makes it suitable for use with angular coordinates in one-dimensional problems. In this study, it can also be used in the $S^2$ domain, as a spherical distance $r$ is defined and limited to $\pi$. In this case, the kernel functions within a range where its periodicity is not apparent.

The stationarity of the kernel functions relies on the definition of the distance function, which is introduced below.

### 4.2. Distance measures on $S^2$

Three different measures of the distance between two orientations [36], $y = \{az, el\}$ and $y' = \{az', el'\}$ are studied,

$$ \operatorname{dist}(y, y') = \begin{cases} \sqrt{\Delta az^2 + \Delta el^2} & , \text{ intrinsic distance;} \\ \sqrt{|x_1 - x'_1|^2 + |x_2 - x'_2|^2 + |x_3 - x'_3|^2} & , \text{ Euclidean distance;} \\ \arccos((x_1 \ x_2 \ x_3) \cdot (x'_1 \ x'_2 \ x'_3)^\mathrm{T}) & , \text{ geodesic distance.} \end{cases} \quad (17) $$

Intrinsic distance uses the angular distance $\Delta\varphi$, between angles $\varphi, \varphi' \in [0, 2\pi]$, given by

$$ \Delta\varphi = \min\{|\varphi - \varphi'|, 2\pi - |\varphi - \varphi'|\} \quad (18) $$

Euclidean and geodesic distance, $\{x_1, x_2, x_3\}$ are the Cartesian coordinates of the orientation $y$ on the unit sphere.

$$ \langle x_1 \ x_2 \ x_3 \rangle^\mathrm{T} = \langle \cos(el) \ \cos(az) \ \cos(el) \ \sin(az) \ \sin(el) \rangle^\mathrm{T} \quad (19) $$

The Euclidean distance is measured along the straight line (chord) connecting the two points, representing the shortest path in the 3D Euclidean space. In contrast, the geodesic distance is defined on the spherical surface, representing the shortest arc connecting the two points along the surface of the sphere.

The intrinsic distance, on the other hand, is the Cartesian distance between two points on a cylindrical projection of the sphere. A comparative analysis of the impact of these three distance definitions on the GPR results is discussed in the application section.

### 4.3. Density non-negativity

To ensure Eq. (5), non-negativity constraints [37–40] or simply a logarithmic data transformation can be applied to the density values

$$ d_{\log} = \log(d + \epsilon). \quad (20) $$

This transformation allows GPR to be performed on $d_{\log}$, the logarithmic values of the density, avoiding negative values. This approach is also named as the log-Gaussian process model [41]. For experimental data where the density $d$ is non-negative, a
small positive constant $\epsilon = 10^{-6}$ is added to handle cases where $d$ might be zero. The inverse transformation of the GPR prediction $\tilde{d}_{\text{log}}$ is

$$ \tilde{d} = \exp(\tilde{d}_{\text{log}}). \qquad (21) $$

However, Eq. (21) tends to compress larger values, and enlarges the difference for the values close to zero. This study proposes a generalization, limiting the logarithmic transformation to values below a chosen threshold $d_0 > 0$ while applying a linear transformation for values above

$$ d_{\text{lli}} = \begin{cases} d_{\text{log}} & \text{if } d \le d_0, \\ \frac{1}{d_0} \cdot d + \log(d_0) - 1 & \text{if } d > d_0. \end{cases} \qquad (22) $$

The slope remains continuous at the threshold point. When $d_0$ increases, the log-linear tends to the logarithmic transformation. Section 5.4 discusses the choice of $d_0$ value. The inverse transformation is given as

$$ \tilde{d} = \begin{cases} \exp(\tilde{d}_{\text{lli}}) & \text{if } \tilde{d}_{\text{lli}} \le \log(d_0), \\ d_0 \cdot (\tilde{d}_{\text{lli}} - \log(d_0) + 1) & \text{if } \tilde{d}_{\text{lli}} > \log(d_0). \end{cases} \qquad (23) $$

Particular attention must be given to the variance when reconstructing the pole density function. Not only must the reconstruction itself remain non-negative, but its CI must also respect this constraint

$$ CI = [\tilde{d}^L, \tilde{d}^U] \subset [0, +\infty], \qquad (24) $$

with $\tilde{d}^L$ and $\tilde{d}^U$ the lower and upper bounds of the CI.

Under the logarithmic transformation, the CI of $n$ standard deviations ($\pm n\sigma_{\text{log}}$), is given as

$$ CI_{\text{log}} = [\tilde{d}_{\text{log}}^L, \tilde{d}_{\text{log}}^U] \qquad (25) $$

$$ = [\exp(\tilde{d}_{\text{log}} - n\sigma_{\text{log}}), \exp(\tilde{d}_{\text{log}}^U + n\sigma_{\text{log}})], \quad n \in \mathbb{R}^+, \qquad (26) $$

and $CI_{\text{lli}}$ is defined in the same manner but using Eq. (23). Thanks to these transformations, Eq. (24) is ensured.

However, with this non-linear or semi-linear data transformation, the symmetry of CI around the mean value is no more guaranteed.

Another important consideration is the effect of the data transformation on the kernel's variance hyperparameter, $\sigma_f^2$. This hyperparameter depends on the scale of the transformed data. To facilitate comparisons of $\sigma_f^2$ across different data transformations, the Coefficient of Variation (CV) provides a normalized measure of the kernel's variance relative to the magnitude of the (transformed) data,

$$ CV = \frac{\sigma_f}{d_{\text{max}} - d_{\text{min}}}. \qquad (27) $$

offering a qualitative interpretation of the optimized hyperparameters.

## 4.4. Texture volume preservation

Numerical integration of the volume conservation condition Eq. (4) yields

$$ \tilde{d}^{\text{T}} w = 4\pi, \qquad (28) $$

where $w$ is a weight vector depending on the size of the area element associated to the predictor grid nodes.

Suppose $M$ nodes uniformly distributed on the unit sphere surface, the weight for each node equals to $4\pi/M$, and Eq. (28) becomes

$$ \sum_{i}^{M} \tilde{d}_{i} = M \qquad (29) $$

The scaled density satisfying Eq. (4) becomes

$$ \tilde{d} \leftarrow \frac{M}{\sum_{i} \tilde{d}_{i}} \tilde{d}. \qquad (30) $$

# 5. Application to synthetic textures

This section applies GPR to three synthetic experimental pole figure measurements. A Hexagonal Close-Packed (HCP) lattice is assumed, with the basal plane (0001) activated as the slip plane. The normal vector of the basal plane corresponds to the central axis of the hexagonal lattice. In this case, the projection of a single basal plane leads to only two points on the sphere, symmetric to the sphere center, the reason why the upper view and the lower view of the pole figure are mirrored. Both views are given in following examples as the samplings of the upper and lower spheres are independent, resulting in different grids.
<table><caption>**Table 2**
Synthetic pole figures and random sampling sets.</caption><thead><tr><th></th><th>Analytical</th><th>Sampling set, 1500 points</th></tr></thead><tbody><tr><th>A</th><td></td><td></td></tr><tr><th>B</th><td></td><td></td></tr><tr><th>C</th><td></td><td></td></tr></tbody></table>

## 5.1. Numerical experimental protocol

The reference pole figures (Table 2 left column) are constructed analytically using quadratic splines (c.f. Appendix), providing a known “exact” function for evaluating reconstruction accuracy. The sampling datasets, representing experimental measurements, are independent on each sphere. One given dataset consists of 1500 randomly selected points, and the measured density is assumed to be error-free. The sampling sets are shown in the right columns of Table 2.

The following sections begin by introducing the regularly distributed Fibonacci predictor grid for evaluating the reconstructed Pole Density Function. Next, two error measures are presented: the $L_\infty$ reconstruction error on experimental data points, and $L_2$ error across the entire sphere. Before delving into the customization of GPR, the influence of data transformation on the reconstruction is analyzed.

The study then proceeds to validate the components of GPR customization. First, three distance measures are compared to determine the most suitable one. Then, two hyperparameter optimization methods are computed to identify the most efficient approach for further use. After selecting the appropriate distance measure and optimization method, reconstructions of the Pole Density Function are performed using four isotropic, stationary kernels to identify the best-performing kernel.

Finally, leveraging the statistical nature of GPR, the variance of the reconstructions is analyzed to provide additional insights into the reliability and uncertainty of the reconstructed Pole Density Function.

## 5.2. Fibonacci predictor grid

On the unit sphere, the coordinates of the ith node of the Fibonacci grid are given as

$$ 
\begin{aligned}
z_i &= (2i - 1)/M - 1, \\
x_i &= \sqrt{1 - z_i^2} \cos(2\pi i \phi), \\
y_i &= \sqrt{1 - z_i^2} \sin(2\pi i \phi),
\end{aligned}
$$

with $\phi = (\sqrt{5} + 1)/2$, the golden ratio.

The reconstruction of the pole density function is achieved by interpolating the density given at the random sampling points to the 10 000 points predictor grid (Fig. 3) by the GPR.
<figure><img src="image_14.png" /><figcaption>Fig. 3. Reference grid: Fibonacci grid with 10 000 points.</figcaption></figure>

### 5.3. Error measures

To assess the quality of the reconstruction, two error measures are employed: the interpolation error, $e_e$, between the reconstructed Pole Density Function and the experimental observations at the sampled data points. The interpolation error

$$e_e = \frac{\max(|\tilde{d} - d|)}{\max(d)} \qquad (31)$$

is expected to vanish due to the interpolation property of GPR in absence of measurement error.

The global reconstruction error $e_g$, given by

$$e_g = \frac{\sum |\tilde{d} - d_a|}{4\pi}, \qquad (32)$$

measures the difference between the reconstructed Pole Density Function and the analytical reference function across the densely sampled predictor grid.

### 5.4. Data transformation

Pole figures often contain regions with zero density. When transitioning from positive density values to zero, GPR naturally strives to maintain continuity. This behavior can lead to the prediction of negative density values, which are not physically meaningful. To mitigate this issue, a logarithmic data transformation Eq. (21), is a natural choice. A more general log-linear transformation, Eqs. (22) and (23), is introduced also to address this limitation.

In Table 3, reconstructions are carried out under three conditions: without data transformation, with logarithmic transformation, and with log-linear transformation. For these reconstructions, GPR employs geodesic distance and a squared exponential kernel. The hyperparameters of the GPR are optimized by maximizing the marginal likelihood. Subsequent discussion details the rationale for these choices.

From the figures in the left column, negative densities are obtained using GPR without data transformation, whereas both logarithmic (central column) and log-linear (right column) transformations ensure non-negative values. However, the logarithmic transformation results present significant reconstruction error, of the order from $10^3$ to $10^5$. Indeed, logarithmic function expands small values while compressing larger ones leading to overfitting.

In contrast, the log-linear transformation in Eqs. (22) and (23) considers a logarithmic function for values less than $d_0$ and a linear function for larger values. This approach preserves slope continuity at the threshold point, $d_0 = 3.656 \times 10^{-4}$, which, in this case, is optimized by minimizing the global reconstruction error $e_g$.

As the measurement noise $\sigma_n$ is assumed to be zero, the three considered GPR (without/logarithmic/log-linear transformation) do not present interpolation error. However, for the 3 considered samples (A, B and C), the GPR using the log-linear transformation gives the smallest reconstruction error, $e_g$, while the logarithmic transformation gives the largest. Moreover, comparing to GPR without data transformation, using the log-linear transformation guarantees density non-negativity.

Finally, considering the optimized hyperparameters' values presented in Table 4, it is noticeable that the variance hyperparameter $\sigma_n^2$ depends significantly on the data transformation applied. In fact, this variation is directly related to the magnitude of the data, which can be illustrated by the value of the normalized standard deviation CV Eq. (27), since the latter is around 0.17 in both cases of raw and log-linear GPR. In the case of logarithmic transformation, CV is between 0.45 to 0.49, in the three samples, A, B and C, meaning that the uncertainty region of the reconstruction is almost three times larger in the other two cases.

Regarding the characteristic length hyperparameter $\ell$, its optimization aligns with the intrinsic characteristics of the data. A smaller length scale $\ell$ corresponds to a more oscillatory regression function. The optimized $\ell$ values for the Raw GPR and the log-linear transformation are consistent. However, for the logarithmic transformation, a decrease in the optimized length scale parameter $\ell$, even though effectively preventing negative values in the reconstruction, leads to overfitting. Consequently, all subsequent reconstructions in this study are processed using the log-linear data transformation.
<caption>**Table 3**</caption>

Comparison of reconstructions with different data transformations — without transformation/logarithmic/log-linear for Samples A, B, and C. (geodesic distance, SE kernel, hyperparameters optimized by maximizing $\mathcal{L}(D, \theta)$).

<figure><img src="image_15.png" /><figcaption>**Table 4**</figcaption></figure>

Comparison of optimized hyperparameters with different data transformations — without data transformation/logarithmic transformation/log-linear transformation for Samples A, B and C. (geodesic distance, SE kernel, hyperparameters optimized by maximizing $\mathcal{L}(D, \theta)$).

<table><thead><tr><th>Sample</th><th>Raw GPR</th><th>Logarithmic</th><th>Log-linear</th></tr></thead><tbody><tr><td>A</td><td>$\sigma_f^2 = 1.2584$, $CV = 0.17585$, $\ell = 0.14659$</td><td>$\sigma_f^2 = 56.721$, $CV = 0.48067$, $\ell = 0.09595$</td><td>$\sigma_f^2 = 9.4223e6$, $CV = 0.17587$, $\ell = 0.14671$</td></tr><tr><td>B</td><td>$\sigma_f^2 = 1.0252$, $CV = 0.17186$, $\ell = 0.14829$</td><td>$\sigma_f^2 = 56.5390$, $CV = 0.48234$, $\ell = 0.09329$</td><td>$\sigma_f^2 = 7.6778e6$, $CV = 0.17190$, $\ell = 0.14818$</td></tr><tr><td>C</td><td>$\sigma_f^2 = 0.91073$, $CV = 0.17448$, $\ell = 0.14954$</td><td>$\sigma_f^2 = 50.689$, $CV = 0.45889$, $\ell = 0.09644$</td><td>$\sigma_f^2 = 6.8113e6$, $CV = 0.17439$, $\ell = 0.14965$</td></tr></tbody></table>

### 5.5. Choice of distance measure

Considering the characteristics of the Pole Density Function, its definition on the unit spherical surface $S^2$ imposes a periodicity constraint on the GPR kernel. To account for this, this work proposes employing a distance measure inherently respecting spherical periodicity and applying it within isotropic, stationary kernels.

Table 5 presents the reconstruction results obtained using three different distance measures: intrinsic, geodesic, and Euclidean Eq. (17). To isolate the influence of the distance measure, all other GPR configurations remain consistent: the Squared Exponential kernel is used, and the hyperparameters are arbitrarily set to {$\sigma_f = 1000$, $\ell = 0.4$, $\sigma_n = 1 \times 10^{-5}$}. A small noise term $\sigma_n$ is introduced to ensure numerical stability.

The reconstructed figures reveal a noticeable polar singularity when using the intrinsic distance (Table 5 “intrinsic” column, rows “A”, “B” and “C”), with regions in dark red exceeding the limits of the color bar. This occurs because the intrinsic distance
<caption>**Table 5**</caption>

Comparison of reconstructions with different distance measures for Samples A, B and C — intrinsic/geodesic/Euclidean distance. (log-linear transformation, SE kernel, arbitrarily chosen hyperparameters: $\sigma_f = 1000$, $\ell = 0.4$, $\sigma_n = 1 \times 10^{-5}$).

<table><thead><tr><th></th><th>intrinsic</th><th>geodesic</th><th>Euclidean</th></tr></thead><tbody><tr><td>A</td><td><img src="intrinsic_A.png" alt="A intrinsic distance measure for Sample A"></td><td><img src="geodesic_A.png" alt="A geodesic distance measure for Sample A"></td><td><img src="euclidean_A.png" alt="A euclidean distance measure for Sample A"></td></tr><tr><td>B</td><td><img src="intrinsic_B.png" alt="B intrinsic distance measure for Sample B"></td><td><img src="geodesic_B.png" alt="B geodesic distance measure for Sample B"></td><td><img src="euclidean_B.png" alt="B euclidean distance measure for Sample B"></td></tr><tr><td>C</td><td><img src="intrinsic_C.png" alt="C intrinsic distance measure for Sample C"></td><td><img src="geodesic_C.png" alt="C geodesic distance measure for Sample C"></td><td><img src="euclidean_C.png" alt="C euclidean distance measure for Sample C"></td></tr></tbody></table>

does not account for spherical periodicity. For instance, at $el = \pi/2$, all azimuth angles correspond to the same polar point, yet the intrinsic distance formula produces a non-zero value. Generally, as el approaches $\pi/2$, the distance estimation error increases, resulting in the observed polar singularities. In contrast, reconstructions using geodesic (Table 5 second column) and Euclidean distances (Table 5 third column) do not exhibit such singularities.

Consequently, regarding the reconstruction error $e_g$, the intrinsic distance yields a significant error. In contrast, geodesic and Euclidean distances achieve comparable errors of less than 5%, with Euclidean distance performing slightly better.

This behavior is not expected as Euclidean distance introduces distortions with increasing separation between orientations [42]. While intrinsic distance requires a periodic kernel to eliminate the influence of the polar singularity, the geodesic and Euclidean distance measures are employed in the following section with stationary kernels.

## 5.6. Hyperparameters optimization

For squared exponential kernel, two hyperparameters, $\theta = \{\sigma_f, \ell\}$, are optimized, where $\sigma_f$ represents the a priori variance and $\ell$ denotes the characteristic length scale.

Table 6 presents the reconstruction results for the three experimental datasets using the Squared Exponential kernel with geodesic and Euclidean distance measures, allowing comparison between LOO optimization and marginal likelihood maximization.

Firstly, for all reconstructions, with the log-linear data transformation applied, none of them show negative values, regardless of the chosen optimization method. From a visual perspective, the reconstructions using hyperparameters optimized by maximizing the marginal likelihood (Table 6, the right two columns) are smoother. These reconstructions successfully reproduce the features of the reference pole figures. On the other hand, the reconstructions using $\theta$ optimized by the LOO approach (the left two columns) appear more fragmented and irregular, with most of the pole figures showing low-density regions, failing to capture the distinctive characteristics of the reference figures.

The error measures support these observations. The interpolation errors $e_e$ remain numerically negligible in all cases, meaning the reconstruction values at the experimental points are accurate. However, the reconstruction errors $e_g$ are much larger for the LOO approach, reaching 89% (Table 6, first row, left two columns), while the errors for the marginal likelihood method are around 3% (all reconstructions in the right two columns), regardless of whether geodesic or Euclidean distance is used.

Additionally, the marginal likelihood values of the models optimized by the two methods differ significantly. For hyperparameters optimized by marginal likelihood (Table 6, the right two columns), the maximized values are on the order of $-1 \times 10^{-4}$, while for
<caption>**Table 6**</caption>

Comparison of reconstructions with hyperparameter optimized by LOO and Maximization of $\mathcal{L}(D, \theta)$ for Samples A, B and C. (log-linear transformation, geodesic/Euclidean distance, SE kernel).

<table><caption>**Table 7**</caption><thead><tr><th></th><th colspan="3">LOO + geo.</th><th colspan="3">LOO + Euc.</th><th colspan="3">L + geo.</th><th colspan="3">L + Euc.</th></tr></thead><tbody><tr><td>A</td><td>e<sub>e</sub> =<br>4e-16</td><td>e<sub>g</sub> =<br>0.891</td><td>L =<br>-1.17e11</td><td>e<sub>e</sub> =<br>4e-16</td><td>e<sub>g</sub> =<br>0.891</td><td>L =<br>-1.31e11</td><td>e<sub>e</sub> =<br>2e-12</td><td>e<sub>g</sub> =<br>0.036</td><td>L =<br>-1.17e4</td><td>e<sub>e</sub> =<br>1e-12</td><td>e<sub>g</sub> =<br>0.035</td><td>L =<br>-1.17e4</td></tr><tr><td>B</td><td>e<sub>e</sub> =<br>5e-16</td><td>e<sub>g</sub> =<br>0.777</td><td>L =<br>-1.19e11</td><td>e<sub>e</sub> =<br>6e-16</td><td>e<sub>g</sub> =<br>0.777</td><td>L =<br>-1.30e11</td><td>e<sub>e</sub> =<br>7e-12</td><td>e<sub>g</sub> =<br>0.032</td><td>L =<br>-1.15e4</td><td>e<sub>e</sub> =<br>3e-12</td><td>e<sub>g</sub> =<br>0.031</td><td>L =<br>-1.15e4</td></tr><tr><td>C</td><td>e<sub>e</sub> =<br>3e-16</td><td>e<sub>g</sub> =<br>0.875</td><td>L =<br>-1.31e11</td><td>e<sub>e</sub> =<br>5e-16</td><td>e<sub>g</sub> =<br>0.875</td><td>L =<br>-1.56e11</td><td>e<sub>e</sub> =<br>2e-11</td><td>e<sub>g</sub> =<br>0.030</td><td>L =<br>-1.14e4</td><td>e<sub>e</sub> =<br>1e-11</td><td>e<sub>g</sub> =<br>0.030</td><td>L =<br>-1.13e4</td></tr></tbody></table>

Comparison of optimized hyperparameters and CPU time *t* by LOO and maximization of $\mathcal{L}(D, \theta)$ for Samples A, B and C. (log-linear transformation, geodesic/Euclidean distance, SE kernel).

<table><thead><tr><th>Sample</th><th>LOO</th><th>max L(D, θ)</th><th>Sample</th><th>LOO</th><th>max L(D, θ)</th></tr></thead><tbody><tr><td>A, geo.</td><td>σ<sub>f</sub> = 0.168,<br>CV = 2.3e-5,<br>ℓ = 0.018,<br>t = 16672 s</td><td>σ<sub>f</sub> = 9.422e6,<br>CV = 0.1759,<br>ℓ = 0.147,<br>t = 50.620 s</td><td>A, Euc.</td><td>σ<sub>f</sub> = 0.150,<br>CV = 2.2e-5,<br>ℓ = 0.018,<br>t = 23931 s</td><td>σ<sub>f</sub> = 9.468e6,<br>CV = 0.1763,<br>ℓ = 0.148,<br>t = 49.470 s</td></tr><tr><td>B, geo.</td><td>σ<sub>f</sub> = 0.134,<br>CV = 2.3e-5,<br>ℓ = 0.027,<br>t = 25906 s</td><td>σ<sub>f</sub> = 7.678e6,<br>CV = 0.1719,<br>ℓ = 0.148,<br>t = 47.290 s</td><td>B, Euc.</td><td>σ<sub>f</sub> = 0.122,<br>CV = 2.2e-5,<br>ℓ = 0.027,<br>t = 18843 s</td><td>σ<sub>f</sub> = 7.684e6,<br>CV = 0.1720,<br>ℓ = 0.149,<br>t = 46.080 s</td></tr><tr><td>C, geo.</td><td>σ<sub>f</sub> = 0.142,<br>CV = 2.5e-5,<br>ℓ = 0.019,<br>t = 25552 s</td><td>σ<sub>f</sub> = 6.811e6,<br>CV = 0.1744,<br>ℓ = 0.150,<br>t = 46.790 s</td><td>C, Euc.</td><td>σ<sub>f</sub> = 0.119,<br>CV = 2.3e-5,<br>ℓ = 0.019,<br>t = 15148 s</td><td>σ<sub>f</sub> = 6.826e6,<br>CV = 0.1746,<br>ℓ = 0.151,<br>t = 48.220 s</td></tr></tbody></table>

those optimized by LOO (the left two columns), the values are significantly smaller, about $-1 \times 10^{-11}$. These results clearly show that reconstructions using hyperparameters optimized through marginal likelihood maximization provide much significantly better results compared to those optimized by the LOO approach.

Table 7 shows the values of the optimized hyperparameters and the CPU time, denoted as *t*. The hyperparameters optimized by the two methods show significant differences.

In the first place, the computational cost of LOO is significant (comparison of CPU time *t* in LOO columns and in max $\mathcal{L}(D, \theta)$ columns). In fact, each iteration requires *N* GPR predictions, involving inverting *N* matrices of size ($N-1 \times (N-1)$). This makes the method time-consuming and unsuitable for the current dataset.

Considering the hyperparameters' values, the CV obtained by the LOO optimization are on the order of $10^{-5}$ (Table 7, values of CV in both LOO columns), which is too small. Such small deviation values eliminate the variability of the data, which is also unacceptable.

Additionally, the characteristic length scales $ℓ$ from the LOO approach are approximately 0.02 (Table 7, values of $ℓ$ in both LOO columns), corresponding to a central angle of about 1.15°. This small value leads to overfitting, as the model focuses excessively on local features. This explains the appearance of the reconstructed pole figures (Table 6, the left two columns): in regions without
<caption>**Table 8**</caption>

Comparison of reconstructions with different kernel functions for Samples A, B and C — Squared Exponential/De La Vallée Poussin/Matérn $v = 3/2$ (log-linear transformation, geodesic distance, hyperparameters optimized by maximizing $\mathcal{L}(D, \theta)$).

<table><thead><tr><th></th><th>Squared exponential</th><th>De La Vallée Poussin</th><th>Matérn, v = 3/2</th><th>Matérn, v = 5/2</th></tr></thead><tbody><tr><td>A</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td>B</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

$$e_e = 2e - 12 \qquad e_g = 0.0358 \qquad \mathcal{L} = -1.1707e4$$

$$e_e = 2e - 12 \qquad e_g = 0.0360 \qquad \mathcal{L} = -1.1714e4$$

$$e_e = 5e - 13 \qquad e_g = 0.0226 \qquad \mathcal{L} = -1.1480e4$$

$$e_e = 5e - 13 \qquad e_g = 0.0203 \qquad \mathcal{L} = -1.1237e4$$

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

$$e_e = 7e - 12 \qquad e_g = 0.0320 \qquad \mathcal{L} = -1.1476e4$$

$$e_e = 4e - 12 \qquad e_g = 0.0324 \qquad \mathcal{L} = -1.1465e4$$

$$e_e = 3e - 13 \qquad e_g = 0.0183 \qquad \mathcal{L} = -1.0881e4$$

$$e_e = 6e - 12 \qquad e_g = 0.0143 \qquad \mathcal{L} = -1.0325e4$$

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

$$e_e = 2e - 11 \qquad e_g = 0.0298 \qquad \mathcal{L} = -1.1347e4$$

$$e_e = 8e - 12 \qquad e_g = 0.0300 \qquad \mathcal{L} = -1.1359e4$$

$$e_e = 4e - 13 \qquad e_g = 0.0147 \qquad \mathcal{L} = -1.0676e4$$

$$e_e = 1e - 11 \qquad e_g = 0.0124 \qquad \mathcal{L} = -1.0151e4$$

nearby experimental points, the reconstruction defaults to the *a priori* mean value of zero. At the sample points, it forms sharp, localized peaks that exactly reproduces experimental densities, leaving most areas of the pole figures with low density.

In contrast, hyperparameters optimized by maximizing marginal likelihood yield CV of about 0.17 (Table 7). These values allow the reconstruction to incorporate more variability. The characteristic lengths $\ell$ are around 0.15, corresponding to a central angle of approximately 8.59°. This larger length scale incorporates more information from the surrounding data points, resulting in smoother reconstructed pole figures (Table 6, the right two columns).

A possible explanation for the limitations of the LOO approach in this context is that texture data is inherently independent for each orientation. Each data point represents unique information that cannot be inferred from the presence of other data points. This lack of redundancy reduces the effectiveness of the LOO approach which relies on the assumption that leaving out one point does not significantly alter the predictive capacity of the model trained on the remaining data.

Consequently, hyperparameter optimization in this study proceeds by maximizing the marginal likelihood, which is both more computationally efficient and better suited to the nature of the data.

Regarding the choice of distance measures, the difference in reconstruction results between the geodesic and Euclidean distances is minor, with the Euclidean distance showing a slight advantage in some cases. The optimized hyperparameters are also comparable between the two measures. This similarity can be attributed to the small differences between the two distance measures when the central angle is small.

Given this minor difference and the fact that the geodesic distance scales linearly with the central angle between two points, the geodesic distance is selected for further reconstructions due to its consistency and alignment with the geometry of the spherical domain.

## 5.7. Kernel comparison

With the choice of the log-linear data transformation, geodesic distance measure, and hyperparameter optimization via marginal likelihood maximization, it becomes possible to employ any isotropic, stationary kernel function for the reconstruction of the Pole Density Function.

In this section, four kernel functions are compared: the Squared Exponential (SE), De La Vallée Poussin (DLVP), and Matérn kernels with $v = 3/2$ and $v = 5/2$, as summarized in Table 1. Table 8 presents the reconstruction results for samples A, B, and C using these four kernels.

A key point is that all four kernels produce plausible reconstructions of the pole figures. As expected, due to the interpolation property of GPR, these reconstructions accurately match experimental data at the measurement points.
<caption>Table 9</caption>

Comparison of optimized hyperparameters of different kernels for Samples A, B and C — Squared Exponential/De La Vallée Poussin/Matérn $v = 3/2$/Matérn $v = 5/2$ (log-linear transformation, geodesic distance, hyperparameters optimized by maximizing $L(D, \theta)$).

<table><thead><tr><th></th><th>SE</th><th>DLVP</th><th>Matérn, v = 3/2</th><th>Matérn, v = 5/2</th></tr></thead><tbody><tr><td>A</td><td>σ<sub>f</sub> = 9.422e6,<br>CV = 0.17587,<br>ℓ = 0.1467</td><td>σ<sub>f</sub> = 9.398e6,<br>CV = 0.17565,<br>ℓ = 0.1462</td><td>σ<sub>f</sub> = 4.346e7,<br>CV = 0.37770,<br>ℓ = 0.5280</td><td>σ<sub>f</sub> = 3.112e7,<br>CV = 0.31962,<br>ℓ = 0.2934</td></tr><tr><td>B</td><td>σ<sub>f</sub> = 7.677e6,<br>CV = 0.17190,<br>ℓ = 0.1482</td><td>σ<sub>f</sub> = 7.6734e6,<br>CV = 0.17185,<br>ℓ = 0.1476</td><td>σ<sub>f</sub> = 3.1755e7,<br>CV = 0.34960,<br>ℓ = 0.6233</td><td>σ<sub>f</sub> = 6.8935e7,<br>CV = 0.51509,<br>ℓ = 0.4627</td></tr><tr><td>C</td><td>σ<sub>f</sub> = 6.8113e6,<br>CV = 0.17439,<br>ℓ = 0.1497</td><td>σ<sub>f</sub> = 6.8027e6,<br>CV = 0.17428,<br>ℓ = 0.1490</td><td>σ<sub>f</sub> = 2.4523e7,<br>CV = 0.33090,<br>ℓ = 0.6276</td><td>σ<sub>f</sub> = 5.3171e7,<br>CV = 0.48725,<br>ℓ = 0.4605</td></tr></tbody></table>

<figure><img src="image_28.png" /><figcaption>Fig. 4. Visualization of four stationary kernel functions with optimized characteristic length $ℓ$ and an additional SE kernel function with optimized $ℓ$ from Matérn $v = 5/2$ kernel (Sample A, $r = dist(y, y')$). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In terms of visual results, the SE and DLVP kernels show similar reconstructions (Table 8, comparing between the left two columns), and the two Matérn kernels give comparable results (comparing between the right two columns). A notable difference appears in regions where the reference density is zero. The SE and DLVP kernels tend to predict non-zero positive densities (represented by light blue areas) more frequently than the Matérn kernels.

Quantitatively, the reconstruction error $e_g$ and the log marginal likelihood $L(D, \theta)$ for SE and DLVP kernels are almost identical (Table 8, comparing the values between the left two columns). Both Matérn kernels yield smaller reconstruction errors, with $v = 5/2$ achieving the smallest. The values of $L(D, \theta)$ for all four kernels are in the same order of magnitude, around $-1 \times 10^{-4}$. Higher the values of $L(D, \theta)$ lower reconstruction errors are. This consistency reinforces the validity of using marginal likelihood maximization as the hyperparameter optimization method.

The optimized hyperparameters' values of each kernel for samples A, B, and C are presented in Table 9. When comparing the CV values, the SE and DLVP kernels exhibit similar values, around 0.17 (the left two columns). In contrast, Matérn kernels, despite having smaller reconstruction errors $e_g$ (Table 8, $e_g$ in the fourth column), show a greater uncertainty in their reconstructions. The CV values for Matérn $v = 3/2$ are consistently around 0.35 across all samples (Table 9 third column), while for Matérn $v = 5/2$, the CV varies more, ranging from 0.32 to 0.52 (fourth column).

The characteristic length's values $ℓ$ cannot be directly compared due to the unique nature of each kernel. Fig. 4 illustrates the four kernel functions using the optimized $ℓ$ for Sample A. Additionally, an SE kernel is plotted with the $ℓ$ optimized for the Matérn $v = 5/2$ kernel. Notably, with the same $ℓ$, the Matérn $v = 5/2$ kernel covers a broader data range.

The SE and DLVP kernels overlap significantly (blue continuous and orange dotted lines), explaining the similarity in their optimized hyperparameters and reconstruction results. Given that the geodesic distance already accounts for spherical periodicity and that the DLVP kernel yields nearly identical results to the SE kernel while being more complex, the SE kernel is preferred.

In contrast, the Matérn kernels exhibit distinct behaviors, generally considering a broader range of data points during interpolation (green and red lines). This reduces the occurrence of unexpected positive densities in regions with large null densities (comparison of white regions between the left and the right half of Table 8) and leads to smaller reconstruction errors. However, the inclusion of more data also increases uncertainty, as indicated by higher CV values (Table 9).

In summary, the SE and DLVP kernels focus on data points closer to the measurement locations, reducing a priori uncertainty. However, their current optimized hyperparameters make them prone to oscillatory behavior during extrapolation, resulting in presence of positive densities in null density regions and approximately 50% higher $e_g$ compared to the Matérn kernels. Between the tested Matérn variants, the $v = 5/2$ achieves lower $e_g$ and higher marginal likelihoods than the $v = 3/2$ across all samples A, B and C. Thus, the Matérn $v = 5/2$ kernel is chosen for further analysis.
<caption>**Table 10**</caption>

Visualization of GPR reconstructions and 95% CI bounds for Samples A, B and C — Lower bound $\tilde{d}^L$/Reconstruction (a posteriori mean) $\tilde{d}$/Upper bound $\tilde{d}^U$ (log-linear transformation, geodesic distance, hyperparameters optimized by maximizing $L(D, \theta)$, Matérn $\nu = 5/2$ kernel).

<figure><img src="image_29.png" /></figure>

## 5.8. Confidence intervals

GPR, as a statistical method, enables the estimation of reconstruction confidence. The 95% CI in the space of log-linearly transformed data

$$95\%CI_{I_{\text{III}}} = [\tilde{d}_{\text{III}} - 1.96\sigma_{\text{III}}, \tilde{d}_{\text{III}} + 1.96\sigma_{\text{III}}], \qquad (33)$$

is symmetric around the a posteriori mean function $\tilde{d}_{\text{III}}$. Here, $\sigma_{\text{III}}$ is the standard deviation function Eq. (8), which depends on the data set D, the a priori variance $\sigma_f^2$, and the predictor point y.

When transformed back into the original data space Eq. (23), the 95% CI, denoted as $[\tilde{d}^L, \tilde{d}^U]$, ensures non-negativity. In this space, the interval bounds are no longer necessarily symmetric to the mean, reflecting the non-linear nature of the transformation.

Table 10 presents visualizations of the 95% CI lower bounds (left column), reconstructions (middle column), and 95% CI upper bounds (right column) for samples A, B, and C.

Visually, the most noticeable differences occur in the zero-density regions. The lower bounds $\tilde{d}^L$ rarely exhibit positive values, the mean $\tilde{d}$ shows only a few, and the upper bounds $\tilde{d}^U$ generally display small positive values. In these regions, the difference between $\tilde{d}$ and $\tilde{d}^U$ is larger than that between $\tilde{d}$ and $\tilde{d}^L$ because $\tilde{d}^L$ does not contain negative values. This confirms that the data transformation ensures the non-negativity of the CI. Apart from this, the visible differences are minor, indicating that within the 95% CI, the reconstructions exhibit minimal variation.

Regarding interpolation error $e_e$, although the values for $\tilde{d}^L$ and $\tilde{d}^U$ are higher than for $\tilde{d}$, they remain on the order of $10^{-7}$, which is negligible. This confirms that the variance at sample points is zero and the experimental data is fully respected.

For reconstruction error $e_g$, across all three samples, $\tilde{d}^L$ exhibits errors two to three times larger than $\tilde{d}$, while $\tilde{d}^U$ shows errors twice as large as $\tilde{d}^L$. However, with the exception of the largest error ($e_g = 0.1452$ for $\tilde{d}^U$ in sample A, possibly due to its highly localized texture), the maximum errors for other bounds are around 0.06 to 0.07, suggesting numerically that the reconstructions are robust.

## 5.9. Comparison with other reconstruction approaches

In this last section, the reconstructions obtained with the current GPR settings are compared against those produced by the conventionally used spherical harmonics expansion approach and additionally, by a linear reconstruction (both algorithms provided by [43]). The figures and corresponding error measures are presented in Table 11.
<caption>**Table 11**</caption>

Comparison of different reconstruction approaches for Sample A, B and C — GPR (log-linear transformation, geodesic distance, Matérn $v = 5/2$ kernel, hyperparameters optimized by maximizing $L(D, \theta)$)/Spherical harmonics approximation (Raw data)/Linear interpolation (Raw data).

<table><thead><tr><th></th><th>GPR</th><th>Harmonic</th><th>Linear</th></tr><tr><th></th><th>1500 functions</th><th>$m = 28, 841$ functions</th><th></th></tr></thead><tbody><tr><th>A</th><td><img src="image_A_gpr.png" alt="GPR reconstruction of a point with true density 5e-13"></td><td><img src="image_A_harmonic.png" alt="Harmonic reconstruction of a point with true density 0.0203"></td><td><img src="image_Ae</td><td>5e-13</td></tr><tr><th>B</th><td><img src="image_B_gpr.png" alt="GPR reconstruction of a point with true density 6e-12"></td><td><img src="image_B_harmonic.png" alt="Harmonic reconstruction of a point with true density 0.0143"></td><td><img src="image_B_harmonic.png" alt="Harmonic reconstruction of a point with true density 0.0158"></td><td>6e-13</td></tr><tr><th>C</th><td><img src="image_C_gpr.png" alt="GPR reconstruction of a point with true density 1e-11"></td><td><img src="image_C_harmonic.png" alt="Harmonic reconstruction of a point with true density 0.0124"></td><td><img src="image_C_harmonic.png" alt="Harmonic reconstruction of a point with true density 0.0253"></td><td>1e-11</td></tr></tbody></table>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

GPR and linear reconstructions avoid negative densities by construction, unlike the harmonic reconstructions. For GPR, this is achieved through the log-linear transformation, ensuring all reconstructed values are non-negative. For linear reconstruction, the lack of smoothness prevents interpolated values from exceeding the observed data range, thereby avoiding negative densities.

In contrast, harmonic reconstructions are inherently oscillatory, potentially resulting in both negative values and unexpected positive values in regions where the true density is zero. For example, in Sample A, the center zone exhibits a highly localized peak requiring high-order harmonics with large coefficients for accurate reproduction. However, with the order limited to $m = 28$ (as set by the default MTEX configuration [43]), the method struggles to capture the adjacent pit, leading to significant negative values.

Regarding interpolation error $e_g$, both GPR and linear reconstructions produce exact values at sample points, while the harmonic approach does not. This is because harmonic reconstruction is an approximation rather than an interpolation method. Increasing the order of the basis functions may improve accuracy in matching sample values but significantly increases computational and storage demands, as the total number of basis functions grows quadratically with $(m+1)^2$, $m$ being the order.

For the reconstruction error $e_g$ shown in Table 11, the linear interpolation exhibits the largest error among the three methods for each sample. This is because linear reconstruction cannot extrapolate higher density peaks not present in the range of the measurements. While the harmonic approach achieves lower $e_g$ values than the linear interpolation, its error is still at least three times larger than that of GPR, due to negative values, imprecise sample capture, and oscillatory behavior.

From these comparisons, the customized GPR approach achieves the smallest reconstruction error across all samples (A, B, and C) while preventing negative densities in the results. Also, GPR provides a CI which is also non-negative thanks to the data transformation (Section 5.8).

## 6. Conclusion and future work

The GPR approach presents an alternative to the traditional spherical harmonics series expansion method for reconstructing pole figures. The adaptability of GPR, enabled by the ability to design kernel functions based on prior knowledge, allows for greater adaptability to experimental data.

For the tested synthetic samples, the best performance was obtained using a log-linear data transformation to ensure the non-negativity of both reconstructions and CI, a geodesic distance measure to account for spherical periodicity, optimizing hyperparameters through marginal likelihood maximization for stability and efficiency, and employing a Matérn kernel with $v = 5/2$ for minimal reconstruction errors. This setup consistently outperformed the spherical harmonics approach with functions of order 28, in both interpolation and reconstruction accuracy. These results contribute to a deeper understanding of the trade-offs between accuracy, computational efficiency, and the representation of complex texture features.
<caption>**Table A.12**</caption>

Sample A construction: radial quadratic spline functions' center points coordinates, and corresponding radial parameter *r* and amplitude *a*.

<table><thead><tr><th>az</th><th>el</th><th>r</th><th>a</th><th>az</th><th>el</th><th>r</th><th>a</th></tr></thead><tbody><tr><td>-π/2</td><td>0</td><td>0.8</td><td>0.97</td><td>π/2</td><td>0</td><td>0.8</td><td>0.97</td></tr><tr><td>π + 0.6</td><td>0.6</td><td>0.96</td><td>0.9</td><td>0.6</td><td>-0.6</td><td>0.96</td><td>0.9</td></tr><tr><td>π + 0.6</td><td>-0.6</td><td>0.96</td><td>0.9</td><td>0.6</td><td>0.6</td><td>0.96</td><td>0.9</td></tr><tr><td>π - 0.5</td><td>-1.0</td><td>0.93</td><td>0.5</td><td>-0.5</td><td>-1.0</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.5</td><td>-0.5</td><td>0.93</td><td>0.5</td><td>-0.5</td><td>-0.5</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.5</td><td>0</td><td>0.93</td><td>0.5</td><td>-0.5</td><td>0</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.5</td><td>0.5</td><td>0.93</td><td>0.5</td><td>-0.5</td><td>0.5</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.5</td><td>1.0</td><td>0.93</td><td>0.5</td><td>-0.5</td><td>1.0</td><td>0.93</td><td>0.5</td></tr><tr><td>π + 0.05</td><td>0</td><td>0.99</td><td>0.97</td><td>0.05</td><td>0</td><td>0.99</td><td>0.97</td></tr></tbody></table>

The statistical foundation of GPR not only provides reconstructions with quantifiable CI but also accommodates measurement uncertainties and supports the sampling of alternative plausible reconstructions. These features enhance the utility and robustness of GPR for pole figure analysis and related applications.

In future work, GPR might be compared to other machine learning methods, such as random forests or neural networks, to further evaluate their relevance in the context of texture representation. Moreover, the EBSD measurement error could also be included in the reconstruction of Pole Density Function using proposed method. Finally, the ultimate goal is to integrate the proposed texture representation method in data-driven numerical simulations of polycrystalline materials.

### CRediT authorship contribution statement

**Bingqian Li:** Writing – review & editing, Writing – original draft, Investigation. **Piotr Breitkopf:** Writing – review & editing, Supervision, Funding acquisition. **Ludovic Cauvin:** Writing – review & editing, Validation, Supervision, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors gratefully acknowledge the financial support of the Région Hauts-de-France, France and the SAFE AI Chair of the Heudiasyc Laboratory (UMR-CNRS 7253), France, under the co-funding of the Ph.D. project.

## Appendix. Construction of reference pole figure samples on the unit sphere

Each reference pole figures sample in A, B and C is constructed by a summation of radial quadratic spline functions.

### Mathematical representation

An analytical Pole Density Function $P_h(y)$ defined on unit spherical surface $S^2 = \{x \in \mathbb{R}^3 : \|x\| = 1\}$ is given by:

$$P_h(y) = \sum_i a_i \cdot s(y \cdot y_i^c, r_i),$$

where $y \in S^2$ is an orientation on the sphere, $y_i^c \in S^2$ is the center of the $i$th radial function, $r_i$ is the radial parameter controlling the extent of the spline, $a_i$ is the amplitude of the $i$th radial function and finally, $y \cdot y_i^c$ is the dot product of the vectors in Cartesian coordinates, providing the cosine of the angle between them.

The radial quadratic spline function $s(z, r)$ is defined as

$$s(z, r) = \begin{cases} \frac{(z-r)^2}{(1-r)^2}, & \text{if } z > r, \\ 0, & \text{otherwise.} \end{cases}$$

### Data points and parameters

The centers $y_i^c$ are specified using their polar coordinates $(az_i, el_i)$, which are converted into Cartesian unit vectors Eq. (19). The corresponding radial parameters $r_i$ and amplitudes $a_i$ are as follows (see Tables A.12–A.14):
<caption>Table A.13</caption>

Sample B construction: radial quadratic spline functions' center points coordinates, and corresponding radial parameter *r* and amplitude *a*.

<table><caption>Table A.14</caption><thead><tr><th>az</th><th>el</th><th>r</th><th>a</th><th>az</th><th>el</th><th>r</th><th>a</th></tr></thead><tbody><tr><td>π/2</td><td>0</td><td>0.8</td><td>0.97</td><td>-π/2</td><td>0</td><td>0.8</td><td>0.97</td></tr><tr><td>π + 0.6</td><td>0.6</td><td>0.96</td><td>0.9</td><td>0.6</td><td>-0.6</td><td>0.96</td><td>0.9</td></tr><tr><td>π + 0.6</td><td>-0.6</td><td>0.96</td><td>0.9</td><td>0.6</td><td>0.6</td><td>0.96</td><td>0.9</td></tr><tr><td>π - 0.7</td><td>-0.7</td><td>0.95</td><td>0.5</td><td>-0.7</td><td>-0.7</td><td>0.95</td><td>0.5</td></tr><tr><td>π - 0.6</td><td>-0.5</td><td>0.93</td><td>0.4</td><td>-0.6</td><td>-0.5</td><td>0.93</td><td>0.4</td></tr><tr><td>π - 0.5</td><td>0</td><td>0.93</td><td>0.7</td><td>-0.5</td><td>0</td><td>0.93</td><td>0.7</td></tr><tr><td>π - 0.6</td><td>0.5</td><td>0.93</td><td>0.4</td><td>-0.6</td><td>0.5</td><td>0.93</td><td>0.4</td></tr><tr><td>π - 0.7</td><td>0.7</td><td>0.95</td><td>0.5</td><td>-0.7</td><td>0.7</td><td>0.95</td><td>0.5</td></tr><tr><td>0.05</td><td>0</td><td>0.93</td><td>0.97</td><td>0.05 - π</td><td>0</td><td>0.93</td><td>0.97</td></tr></tbody></table>

Sample C construction: radial quadratic spline functions' center points coordinates, and corresponding radial parameter *r* and amplitude *a*.

<table><thead><tr><th>az</th><th>el</th><th>r</th><th>a</th><th>az</th><th>el</th><th>r</th><th>a</th></tr></thead><tbody><tr><td>π/2</td><td>0</td><td>0.8</td><td>0.97</td><td>-π/2</td><td>0</td><td>0.8</td><td>0.97</td></tr><tr><td>π + 0.6</td><td>0.4</td><td>0.93</td><td>0.97</td><td>0.6</td><td>0.4</td><td>0.93</td><td>0.97</td></tr><tr><td>π + 0.6</td><td>-0.4</td><td>0.93</td><td>0.97</td><td>0.6</td><td>-0.4</td><td>0.93</td><td>0.97</td></tr><tr><td>π - 0.8</td><td>-0.6</td><td>0.95</td><td>0.6</td><td>-0.8</td><td>-0.6</td><td>0.95</td><td>0.6</td></tr><tr><td>π - 0.6</td><td>-0.5</td><td>0.93</td><td>0.5</td><td>-0.6</td><td>-0.5</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.5</td><td>0</td><td>0.93</td><td>0.7</td><td>-0.5</td><td>0</td><td>0.93</td><td>0.7</td></tr><tr><td>π - 0.6</td><td>0.5</td><td>0.93</td><td>0.5</td><td>-0.6</td><td>0.5</td><td>0.93</td><td>0.5</td></tr><tr><td>π - 0.8</td><td>0.6</td><td>0.95</td><td>0.6</td><td>-0.8</td><td>0.6</td><td>0.95</td><td>0.6</td></tr><tr><td>π + 0.05</td><td>0</td><td>0.97</td><td>0.7</td><td>0.05</td><td>0</td><td>0.97</td><td>0.7</td></tr></tbody></table>

## Data availability

Data will be made available on request.

## References

[1] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, in: Adaptive Computation and Machine Learning, MIT Press, Cambridge, Mass, 2006, OCLC: ocm61285753.

[2] L. Laurent, R. Le Riche, B. Soulier, P.-A. Boucard, An overview of gradient-enhanced metamodels with applications, Arch. Comput. Methods Eng. 26 (1) (2019) 61–106, http://dx.doi.org/10.1007/s11831-017-9226-3.

[3] D. Gaudrie, R. Le Riche, V. Picheny, B.t. Enaux, V. Herbert, Modeling and optimization with Gaussian processes in reduced eigenbases, Struct. Multidiscip. Optim. 61 (6) (2020) 2343–2361, http://dx.doi.org/10.1007/s00158-019-02458-6.

[4] H.B. Nielsen, S.N. Lophaven, J. Søndergaard, DACE - a matlab kriging toolbox, 2002, URL http://www2.compute.dtu.dk/pubdb/pubs/1460-full.html.

[5] F. Kocks, C. Tomé, H.-R. Wenk, Texture and anisotropy. preferred orientations in polycrystals and their effect on material properties, 2000.

[6] D. Mainprice, R. Hielscher, H. Schaeben, Calculating anisotropic physical properties from texture data using the MTEX open-source package, Geol. Soc. Lond. Spec. Publ. 360 (1) (2011) 175–192, http://dx.doi.org/10.1144/SP360.10, URL https://www.lyellcollection.org/doi/10.1144/SP360.10.

[7] J. Jin, L. Cauvin, B. Raghavan, P. Breitkopf, S. Dutta, M. Xiao, Towards a data-driven paradigm for characterizing plastic anisotropy using principal components analysis and manifold learning, Comput. Mater. Sci. 235 (2024) 112834, http://dx.doi.org/10.1016/j.commatsci.2024.112834, URL https://www.sciencedirect.com/science/article/pii/S0927025624000557.

[8] H.J. Bunge, Texture Analysis in Materials Science: Mathematical Methods, english ed., Butterworths, London Boston, 1982.

[9] R.-J. Roe, Description of crystallite orientation in polycrystalline materials. III. General solution to pole figure inversion, J. Appl. Phys. 36 (6) (1965) 2024–2031, http://dx.doi.org/10.1063/1.1714396, URL https://pubs.aip.org/jap/article/36/6/2024/166209/Description-of-Crystallite-Orientation-in.

[10] R. Hielscher, H. Schaeben, A novel pole figure inversion method: specification of the MTEX algorithm, J. Appl. Crystallogr. 41 (6) (2008) 1024–1037, http://dx.doi.org/10.1107/S0021889808030112, URL https://scripts.iucr.org/cgi-bin/paper?S0021889808030112.

[11] A. Eghtesad, T. Barrett, M. Knezevic, Compact reconstruction of orientation distributions using generalized spherical harmonics to advance large-scale crystal plasticity modeling: Verification using cubic, hexagonal, and orthorhombic polycrystals, Acta Mater. 155 (2018) 418–432, http://dx.doi.org/10.1016/j.actamat.2018.06.017.

[12] S. Ben Elhaj Salah, R. Le Riche, J. Bruchon, N. Moulin, Génération de champs aléatoires de perméabilité pour les composites par méthodes à noyaux, in: 16ème Colloque National en Calcul de Structures, CNRS, CSMA, ENS Paris-Saclay, CentraleSupélec, Ciens, France, 2024, URL https://hal.science/hal-04610944.

[13] O. Ozbayram, A. Olivier, L. Graham-Brady, Heteroscedastic Gaussian process regression for material structure–property relationship modeling, Comput. Methods Appl. Mech. Engrg. 431 (2024) 117326, http://dx.doi.org/10.1016/j.cma.2024.117326, URL https://linkinghub.elsevier.com/retrieve/pii/S0045782524005814.

[14] B. Zheng, Z. Zheng, G.X. Gu, Uncertainty quantification and prediction for mechanical properties of graphene aerogels via Gaussian process metamodels, Nano Futur. 5 (4) (2021) 045004, http://dx.doi.org/10.1088/2399-1984/ac3c8f, Publisher: IOP Publishing.

[15] O. Engler, V. Randle, Introduction to Texture Analysis: Macrotexture, Microtexture, and Orientation Mapping, Second Edition, second ed., CRC Press, Boca Raton, 2009, http://dx.doi.org/10.1201/9781420063660.

[16] X. Wang, Multi-Scale Characterization and Modeling of Highly Anisotropic Zinc Based Alloys for Forming Process Simulation (Ph.D. thesis), Université de Technologie de Compiègne, 2017.

[17] G.I. Taylor, Plastic strain in metals, J. Inst. Met. 62 (1938) 307–324, URL https://cir.nii.ac.jp/crid/1571135650487126144.
[18] R. Masson, M. Bornert, P. Suquet, A. Zaoui, An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals, J. Mech. Phys. Solids 48 (6) (2000) 1203-1227, http://dx.doi.org/10.1016/S0022-5096(99)00071-X, URL https://www.sciencedirect.com/science/article/pii/S002250969900071X.

[19] R.A. Lebensohn, C.N. Tomé, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys, Acta Met. Mater. 41 (9) (1993) 2611-2624, http://dx.doi.org/10.1016/0956-7151(93)90130-K, URL https://www.sciencedirect.com/science/article/pii/S095671519390130K.

[20] A. Molinari, S. Ahzi, R. Kouddane, On the self-consistent modeling of elastic-plastic behavior of polycrystals, Mech. Mater. 26 (1) (1997) 43-62, http://dx.doi.org/10.1016/S0167-6636(97)00017-3, URL https://www.sciencedirect.com/science/article/pii/S0167663697000173.

[21] L. Cauvin, B. Raghavan, S. Bouvier, X. Wang, F. Meraghni, Multi-scale investigation of highly anisotropic zinc alloys using crystal plasticity and inverse analysis, Mater. Sci. Eng.: A 729 (2018) 106-118, http://dx.doi.org/10.1016/j.msea.2018.05.038, URL https://linkinghub.elsevier.com/retrieve/pii/S0921509318306865.

[22] E. Herve, A. Zaoui, n-Layered inclusion-based micromechanical modelling, Internat. J. Engrg. Sci. 31 (1) (1993) 1-10, http://dx.doi.org/10.1016/0020-7225(93)90059-4, URL https://www.sciencedirect.com/science/article/pii/S0020722593900594.

[23] B. Panicaud, K. Saanouni, A. Baczmański, M. François, L. Cauvin, L. Le Joncour, Theoretical modelling of ductile damage in duplex stainless steels – Comparison between two micro-mechanical elasto-plastic approaches, Comput. Mater. Sci. 50 (6) (2011) 1908-1916, http://dx.doi.org/10.1016/j.commatsci.2011.01.038, URL https://www.sciencedirect.com/science/article/pii/S0927025611000620.

[24] L. Xia, B. Raghavan, P. Breitkopf, W. Zhang, Numerical material representation using proper orthogonal decomposition and diffuse approximation, Appl. Math. Comput. 224 (2013) 450-462, http://dx.doi.org/10.1016/j.amc.2013.08.052, URL https://linkinghub.elsevier.com/retrieve/pii/S009630031300920X.

[25] L. Madej, L. Rauch, K. Perzynski, P. Cybulka, Digital material representation as an efficient tool for strain inhomogeneities analysis at the micro scale level, Arch. Civ. Mech. Eng. 11 (3) (2011) 661-679, http://dx.doi.org/10.1016/S1644-9665(12)60108-3, URL https://www.sciencedirect.com/science/article/pii/S1644966512601083.

[26] R. Quey, P. Dawson, F. Barbe, Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing, Comput. Methods Appl. Mech. Engrg. 200 (17-20) (2011) 1729-1745, http://dx.doi.org/10.1016/j.cma.2011.01.002, URL https://hal.science/hal-00858028. Publisher: Elsevier.

[27] R. Quey, L. Renversade, Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron X-ray diffraction data, Comput. Methods Appl. Mech. Engrg. 330 (2018) 308-333, http://dx.doi.org/10.1016/j.cma.2017.10.029, URL https://www.sciencedirect.com/science/article/pii/S0027025617307028.

[28] P. Seibert, A. Raßloff, K.A. Kalina, J. Gussone, K. Bugelnig, M. Diehl, M. Kästner, Two-stage 2D-to-3D reconstruction of realistic microstructures: Implementation and numerical validation by effective properties, Comput. Methods Appl. Mech. Engrg. 412 (2023) 116098, http://dx.doi.org/10.1016/j.cma.2023.116098, URL https://linkinghub.elsevier.com/retrieve/pii/S0045782523002220.

[29] S. Bargmann, B. Klusemann, J. Markmann, J.E. Schnabel, K. Schneider, C. Soyarslan, J. Wilmers, Generation of 3D representative volume elements for heterogeneous materials: A review, Prog. Mater. Sci. 96 (2018) 322-384, http://dx.doi.org/10.1016/j.pmatsci.2018.02.003, URL https://linkinghub.elsevier.com/retrieve/pii/S0079642518300161.

[30] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, Int. J. Solids Struct. 40 (13) (2003) 3647-3679, http://dx.doi.org/10.1016/S0020-7683(03)00143-4, URL https://www.sciencedirect.com/science/article/pii/S0020768303001434.

[31] T. Kanit, F. N'Guyen, S. Forest, D. Jeulin, M. Reed, S. Singleton, Apparent and effective physical properties of heterogeneous materials: Representativity of samples of two materials from food industry, Comput. Methods Appl. Mech. Engrg. 195 (33) (2006) 3960-3982, http://dx.doi.org/10.1016/j.cma.2005.07.022, URL https://www.sciencedirect.com/science/article/pii/S0045782505003889.

[32] F. Roters, P. Eisenlohr, T.R. Bieler, D. Raabe, Crystal Plasticity Finite Element Methods: In Materials Science and Engineering, first ed., Wiley, 2010, http://dx.doi.org/10.1002/9783527631483, URL https://onlinelibrary.wiley.com/doi/book/10.1002/9783527631483.

[33] S.S. Vel, A.C. Cook, S.E. Johnson, C. Gerbi, Computational homogenization and micromechanical analysis of textured polycrystalline materials, Comput. Methods Appl. Mech. Engrg. 310 (2016) 749-779, http://dx.doi.org/10.1016/j.cma.2016.07.037, URL https://linkinghub.elsevier.com/retrieve/pii/S0045782516308192.

[34] F. Fritzen, S. Forest, T. Böhlke, D. Kondo, T. Kanit, Computational homogenization of elasto-plastic porous metals, Int. J. Plast. 29 (2012) 102-119, http://dx.doi.org/10.1016/j.ijplas.2011.08.005, URL https://www.sciencedirect.com/science/article/pii/S0749641911001458.

[35] M.L. Stein, Interpolation of Spatial Data, in: Springer Series in Statistics, Springer, New York, NY, 1999, http://dx.doi.org/10.1007/978-1-4612-1494-6, URL http://link.springer.com/10.1007/978-1-4612-1494-6.

[36] F. Sittel, T. Filk, G. Stock, Principal component analysis on a torus: Theory and application to protein dynamics, J. Chem. Phys. 147 (24) (2017) 244101, http://dx.doi.org/10.1063/1.4998259, URL https://pubs.aip.org/jcp/article/147/24/244101/195589/Principal-component-analysis-on-a-torus-Theory-and.

[37] M. Randon, B. Quost, N. Boudaoud, D.V. Wissel, Vehicle consumption estimation via calibrated Gaussian process regression, in: 2022 IEEE Intelligent Vehicles Symposium, IV, IEEE, Aachen, Germany, 2022, pp. 71-76, http://dx.doi.org/10.1109/IV51971.2022.9827326, URL https://ieeexplore.ieee.org/document/9827326/.

[38] L.P. Swiler, M. Gulian, A.L. Frankel, C. Safta, J.D. Jakeman, A survey of constrained gaussian process regression: Approaches and implementation challenges, J. Mach. Learn. Model. Comput. 1 (2) (2020) http://dx.doi.org/10.1615/JMachLearnModelComput.2020035155, URL https://www.dl.begellhouse.com/journals/558048804a15188a, 2cbcbe11139f18e5, 0776649265326db4.html. Publisher: Begel House Inc...

[39] S. Da Veiga, A. Marrel, Gaussian process modeling with inequality constraints, Ann. Faculté Sci. Toulouse : Math. 21 (3) (2012) 529-555, http://dx.doi.org/10.5802/afst.1344, URL https://afst.centre-mersenne.org/articles/10.5802/afst.1344/.

[40] H. Maatouk, X. Bay, Gaussian process emulators for computer experiments with inequality constraints, Math. Geosci. 49 (5) (2017) 557-582, http://dx.doi.org/10.1007/s11004-017-9673-2.

[41] J. Møller, A. Syversveen, R. Waagepetersen, Log Gaussian Cox processes, Scand. J. Stat. 25 (3) (1998) 451-482, http://dx.doi.org/10.1111/1467-9469.00115.

[42] K. Hinsen, Comment on: "Energy landscape of a small peptide revealed by dihedral angle principal component analysis", Proteins: Struct. Funct. Bioinform. 64 (3) (2006) 795-797, http://dx.doi.org/10.1002/prot.20900, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/prot.20900. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/prot.20900.

[43] F. Bachmann, R. Hielscher, H. Schaeben, Texture analysis with MTEX-Free and open source software toolbox, Solid State Phenom. 160 (2010) http://dx.doi.org/10.4028/www.scientific.net/SSP.160.63.