Kernel Principal Component Analysis for Allen–Cahn Equation

1. Introduction

Differential equations, particularly partial differential equations (PDEs), are commonly used to mathematically represent most real-life problems in different scientific fields. Among them the Allen-Chan (AC) equation is a popular nonlinear PDE for modeling phase transition issues. The AC equation was firstly presented in [1] to simulate the movement of antiphase boundaries in crystalline solids. The AC equation is commonly applied to represent different natural phenomena and has become a fundamental equation for the diffuse interface method used to investigate phase transitions and interfacial dynamics in various fields such as materials science [2], image analysis [3], fluid dynamics [4] and mean curvature flow [5]. In materials science, it is used to simulate phase transformations in polymer, metal and ceramic materials [6]. AC equation is further used in multi-reconstruction [7] and predicting the evolution of microstructures in various materials and also in biophysics [8].

The AC equation define by

$$\begin{array}{ccc}\hfill u_t& ={\mu }^2\Delta u-f\left(u\right),\hfill & \hfill (x,t)\in \mathsf{\Omega }\times (0,T],\end{array}$$

(1)

explains the movement of the anti-phase borders of a dual metal mixture at a set temperature. In equation (1), $u(x,t)$ is used to denote the amount of one metal component in the alloy while the positive parameter $\mu$ represents the narrow inter-facial width that is smaller than the lab scale length. We examine the AC equation (1) on a domain $\mathsf{\Omega }\subset {\mathbb{R}}^d$, ($d=1,2,3$), enforced with homogeneous Neumann or periodic boundary conditions, and with the double-well potential $F\left(u\right)={(u^2-1)}^2/4$. The nonlinear functional in (1) is given as $f\left(u\right)=F^{\prime }\left(u\right)=u-u^3$, and the free energy functional of the AC equation is defined by

$$\mathcal{E}\left(u\right)={\int }_{\mathsf{\Omega }}\left({\displaystyle \frac{{\mu }^2}{2}}{|\nabla u|}^2+F\left(u\right)\right)dx.$$

(2)

In this form, the AC equation (1) is a dissipative gradient system

$$u_t=-{\displaystyle \frac{\delta \mathcal{E}}{\delta u}},$$

where $\delta \mathcal{E}/\delta u$ represents the variational derivative of the free energy. In other words, the AC equation meets a non-linear stability condition, which means that the free energy functional $\mathcal{E}\left(u\right)$ decreases over time, i.e., $\mathcal{E}\left(u(x,t_n)\right)\le \mathcal{E}\left(u(x,t_m)\right)$ for $t_n>t_m$.

Dimensionality reduction for dynamical systems is a crucial process aimed at simplifying complex models while retaining the essential features of their behavior. In many real-world applications, such as fluid dynamics, climate modeling, and biological systems, the state space can be incredibly high-dimensional, making simulations computationally expensive and difficult to analyze. Intrusive and data-driven reduced order modeling (ROM) are two distinct approaches used for dimensionality reduction. Intrusive ROM relies on integrating reduced models directly into the numerical algorithms of the original system, often requiring detailed knowledge of the governing equations. In contrast, data-driven ROM leverages discrete data to construct models, utilizing techniques like machine learning and statistical analysis to capture system behavior without explicit knowledge of the underlying physics. Among them, principal component analysis (PCA) is a robust statistical method used for reducing dimensionality and visualizing data across fields like image processing, genetics, finance, and social sciences. It transforms original variables into principal components that maximize variance, aiming to create a lower-dimensional representation while retaining original variability [9,10]. PCA’s effectiveness diminishes when data distributions exhibit low variance, making it less suitable for linear subspace decomposition.

Recently, kernel principal component analysis (KPCA) has emerged as an advanced alternative, which is an extension of traditional PCA that enables the analysis of nonlinear relationships within high-dimensional data. By employing a kernel function, KPCA implicitly maps the original data into a higher-dimensional feature space where linear separability is often achieved, allowing for more effective dimensionality reduction and feature extraction. This method excels in capturing complex patterns that conventional PCA might overlook, making it particularly useful in fields such as exploratory data analysis, pattern recognition, face recognition, manifold learning across fields like computer vision, bio-informatics, broad learning systems, and signal processing [11,12,13,14,15,16,17,18]. However, KPCA faces challenges in interpreting results in the original input space and the computational cost associated with kernel evaluations, with multidimensional scaling (MDS) being a notable approach for this purpose [14,19,20]. MDS allows for the iterative computation of the pre-image by employing various distance metrics between the data vectors. Despite these challenges, KPCA remains a powerful tool for exploratory data analysis and is widely utilized in applications

This paper explores the application of KPCA to data that derived from discrete solution vectors of the AC equation, a novel approach in ROM that preserves solution accuracy and conservation of energy dissipation property of the AC equation. In the literature, there are rare papers dealing with the ROM for AC equation. An energy stable ROM using finite differences and convex splitting, is derived for AC equation in [21], and a non-linear POD/Galerkin ROM is applied in [22]. In [23], the authors derive an energy stable ROM for AC equation using discontinuous Galerkin and average vector field methods, utilizing proper orthogonal decomposition-greedy adaptive sampling and discrete empirical interpolation method.

In the sequel, the paper outlines the AC equation’s mathematical formulation with space/time discretization in Section 2, details the KPCA process in Section 3, showcases results related to reduced solution accuracy and computational efficiency for the one-, two- and three-dimensional AC equations in Section 4, and concluding with key findings.

2. Discrete Data from AC Equation

Here, we present the spatio-temporal (full) discrete representation of the AC equation based on its free energy (2), where the solution vectors are utilized to construct the data set for further analysis.

The AC equation is first discretized in space. With this goal in mind, we employ finite difference discretization to the second order differentiation (Laplace) operator $\Delta$ in AC equation utilizing a tensor framework. Let the matrix $A\in {\mathbb{R}}^{m\times m}$ denotes the discrete Laplace operator, where the number m stands for the total degrees of freedom on the d-dimensional spatial mesh ($d=1,2,3)$. Let also the matrix $D_i$ stands, under periodic boundary conditions, for the matrix corresponding to the discretization of the Laplace operator by finite differences on a one-dimensional spatial interval consisting of $m_i$ grid points with the mesh size $\Delta x_i$, given by

$$D_i:={\displaystyle \frac{1}{\Delta x_i^2}}\left(\begin{array}{ccccc}-2& 1& & & 1\\ 1& -2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2& 1\\ 1& & & 1& -2\end{array}\right)\in {\mathbb{R}}^{m_i\times m_i}.$$

Firstly, consider a one-dimensional domain $\mathsf{\Omega }\subset \mathbb{R}$ that partitioned into $m_1$ uniform subintervals with the mesh size $\Delta x_1$. On this grid, simply, we have that $A=D_1$ and $m=m_1$. On a two-dimensional rectangular domain $\mathsf{\Omega }\subset {\mathbb{R}}^2$ that partitioned into $m_1$ and $m_2$ uniform subintervals with the mesh sizes $\Delta x_1$ and $\Delta x_2$ in either $x_1$ and $x_2$ coordinates, respectively, we can take the matrix A of discrete Laplace operator as

$$A=I_{m_2}\otimes D_1+D_2\otimes I_{m_1}.$$

On a three-dimensional cuboid $\mathsf{\Omega }\subset {\mathbb{R}}^3$, on the other hand, let the domain is partitioned into $m_1$, $m_2$ and $m_3$ uniform subintervals with the mesh sizes $\Delta x_1$, $\Delta x_2$ and $\Delta x_3$ in either $x_1$, $x_2$ and $x_3$ coordinates, respectively. Then, the matrix A of discrete Laplace operator is of the form

$$A=I_{m_3}\otimes I_{m_2}\otimes D_1+I_{m_3}\otimes D_2\otimes I_{m_1}+D_3\otimes I_{m_2}\otimes I_{m_1}.$$

In the above setting, ⊗ denotes the Kronecker product and $I_{m_i}$ is the identity matrix of size $m_i$, $i=1,2,3$. We note that the matrix $D_i$, and hence the matrix A, can be easily formed under homogeneous Neumann boundary condition; all we need is setting the most upper-right and the most bottom-left corner entries of the matrix $D_i$ to zero, and take $m_i:=m_i+1$.

After, we establish time-varying semi-discrete solution vector $\mathit{u}\left(t\right):[0,T]\mapsto {\mathbb{R}}^m$ as

$$\mathit{u}\left(t\right)={(u_1\left(t\right),\dots ,u_m\left(t\right))}^T,$$

with $u_i\left(t\right)$ representing the approximate solution at time t and i-th grid node of the partitioned mesh of the domain $\mathsf{\Omega }\subset {\mathbb{R}}^d$, in an appropriate order. After spatial discretization, we obtained the m-dimensional semi-discrete (dynamical) system

$$\dot{\mathit{u}}={\mu }^{2}A\mathit{u}-\mathit{f}\left(\mathit{u}\right),$$

(3)

where $\dot{\mathit{u}}$ is the ordinary differentiation with respect to time variable t, and the m-dimensional vector $\mathit{f}$ is such that ${\mathit{f}}_i\left(\mathit{u}\right)=u_i-u_i^3$, $i=1,\dots ,m$.

The full discrete system for the AC equation is derived by applying a temporal integration technique to the semi-discrete system (3). Since the effective non-intrusive dimension reduction methods rely on an accurate and reliable data, and the data set used in this paper is composed of the discrete solution vectors of the full discrete system of the AC equation, it is crucial to use a temporal integration that is both accurate and preserves the energy dissipation property of the AC equation. An energy-stable scheme is a numerical method that maintains the energy dissipation/conservation of a gradient/Hamiltonian flow at the discrete level. Extensive literature exists on gradient stable schemes for the AC equation. In this instance, we utilize Kahan’s method [24], which is an advanced numerical method designed to enhance the accuracy of time integration for dynamic systems. Kahan time integrator achieves higher precision compared to traditional integration methods, such as the standard Runge-Kutta or Euler methods. It is a second order, time-reversal, and linearly implicit technique for ODEs [24].

For the full discrete system, we partition the time interval $[0,T]$ into $n-1$ equal parts with time step-size $\Delta t=T/(n-1)$ to create discrete time points $t_j=j\Delta t$, $j=0,\dots ,n-1$. The solution vector at time $t_j$ is denoted by ${\mathit{u}}^j=\mathit{u}\left(t_j\right)\in {\mathbb{R}}^m$. Then, the application of Kahan’s method leads to the full discrete problem: given ${\mathit{u}}^0\in {\mathbb{R}}^m$ by the initial condition, compute ${\mathit{u}}^{j+1}$ from the linear system

$$\begin{array}{ccc}\hfill \left(I_m-{\displaystyle \frac{\Delta t}{2}}{J}_{\mathit{r}}\left({\mathit{u}}^j\right)\right)({\mathit{u}}^{j+1}-{\mathit{u}}^j)& =\Delta t\mathit{r}\left({\mathit{u}}^j\right),\hfill & \hfill j=0,1,\cdots ,n-1,\end{array}$$

(4)

where $\mathit{r}\left(\mathit{u}\right):={\mu }^{2}A\mathit{u}-\mathit{u}+{\mathit{u}}^3$ is the right-hand side vector of the dynamical system (3), and $J_{\mathit{r}}\left(\mathit{u}\right)$ is its diagonal Jacobian matrix. Under the given discretization, the discrete form of the energy $\mathcal{E}$ reads as

$$E\left(\mathit{u}\right)=\sum _{i=1}^m\left[{\displaystyle \frac{{\mu }^2}{2}}{\left(D_f\mathit{u}\right)}_i^2+{\displaystyle \frac{1}{2}}{(u_i^2-1)}^2\right],$$

(5)

where the matrix $D_f\in {\mathbb{R}}^{m\times m}$ mimics the first order forward finite difference differentiation.

3. Reduced Order Model

In this part, we will cover the linear/nonlinear ROM formulation of the AC equation. First, we briefly outline the conventional PCA method, followed by an explanation of KPCA, the nonlinear counterpart of the PCA. The process involves two steps: starting with vectors from the column space of a data matrix, known as the input space, moving to the reduced space via projection, and then mapping back from the reduced space to the input space to get the reduced approximation. Because the ROM approach we are examining relies on a given data, we will focus in the following sections on the data matrix

$$U=[{\mathit{u}}^1\cdots {\mathit{u}}^n]\in {\mathbb{R}}^{m\times n},$$

(6)

where the column vectors represent the solution vectors of the AC equation derived from the full discrete system (4). We make a point to simplify notation by beginning the initial superscript of the snapshot vectors in equation (6) at 1 instead of 0. This means that the column vector ${\mathit{u}}^i\in {\mathbb{R}}^m$ in equation (6) corresponds to the solution vector ${\mathit{u}}^{i-1}\in {\mathbb{R}}^m$, $i=1,\dots ,n$, of the AC equation derived from the system (4).

3.1. Linear Dimension Reduction (PCA)

In the case of a data matrix (6) containing the solution vectors for AC equation, PCA aims to find a linear model of a reduced dimension $k\ll m$, which can accurately represent the variance in the vectors ${\mathit{u}}^i$. We assume that the matrix U has zero column sum, otherwise, we can achieve this by subtracting the average column value $\overline{\mathit{u}}=(1/n)\sum {\mathit{u}}^i$ from each column. The covariance matrix C of the matrix U is defined by

$$C=\sum _{i=1}^n{\mathit{u}}^i{\left({\mathit{u}}^i\right)}^T=U{U}^T\in {\mathbb{R}}^{m\times m},$$

whose diagonalization is given by

$$C=P\mathsf{\Lambda }{P}^T,$$

where the column vectors of the orthogonal matrix $P=[{\mathit{p}}^1\cdots {\mathit{p}}^m]\in {\mathbb{R}}^{m\times m}$ represent the eigenvectors related to the (sorted) eigenvalues ${\lambda }_1>\cdots >{\lambda }_m$ lined up on the diagonal entries of diagonal matrix $\mathsf{\Lambda }\in {\mathbb{R}}^{m\times m}$. Next, the eigenvectors $\{{\mathit{p}}^1,\dots ,{\mathit{p}}^k\}\subset {\mathbb{R}}^m$ of the covariance matrix C with the k largest eigenvalues ${\lambda }_1>\cdots >{\lambda }_k$ are considered as the basis of the k-dimensional linear subspace in a straightforward manner [17]. Ultimately, any random vector ${\mathit{u}}^{*}\in {\mathbb{R}}^m$ in the input space can be roughly expressed by the pre-image ${\widehat{\mathit{u}}}^{*}\in {\mathbb{R}}^m$ as a linear combination of the eigenvectors

$${\mathit{u}}^{*}\approx {\widehat{\mathit{u}}}^{*}=\sum _{i=1}^k{\mathit{z}}_i^{*}{\mathit{u}}^i=P_k{\mathit{z}}^{*},$$

(7)

with the matrix $P_k=[{\mathit{p}}^1\cdots {\mathit{p}}^k]\in {\mathbb{R}}^{m\times k}$ being composed of the initial k columns of the orthogonal matrix P. The reel coefficients ${\mathit{z}}_i^{*}$ represent the components of the projection vector ${\mathit{z}}^{*}=[{\mathit{z}}_1^{*},\dots ,{\mathit{z}}_k^{*}]\in {\mathbb{R}}^k$, obtained by projecting ${\mathit{u}}^{*}$ onto the smaller linear subspace. The $k\ll m$ dimensional vector ${\mathit{z}}^{*}=P_k^T{\mathit{u}}^{*}$ lies in the reduced space, which is the projection of the m-dimensional vector ${\mathit{u}}^{*}$ from the input space [25].

3.2. Nonlinear Dimension Reduction (KPCA)

PCA is restricted to reducing dimensionality in a linear manner. However, standard PCA may be inefficient when confronted with data featuring intricate structures that cannot be accurately represented within a linear subspace. Yet, KPCA provides a resolution by allowing us to expand the linear PCA for nonlinear dimensionality reduction [12,14,16].

By utilizing KPCA, one maps the vectors from the m-dimensional input space to a larger $M\gg m$ dimensional (potentially infinite-dimensional) space known as the feature space, through a nonlinear map $\mathsf{\Phi }(·):{\mathbb{R}}^m\mapsto {\mathbb{R}}^M$. Then, the traditional PCA is subsequently implemented on the vectors within this feature space. In order to achieve this, we label as $\tilde{U}$ the converted data matrix created as

$$\tilde{U}=[\mathsf{\Phi }\left({\mathit{u}}^1\right)\cdots \mathsf{\Phi }\left({\mathit{u}}^n\right)]\in {\mathbb{R}}^{M\times n},$$

where the columns consist of the converted vectors $\mathsf{\Phi }\left({\mathit{u}}^l\right)\in {\mathbb{R}}^M$ corresponding to the input space vectors ${\mathit{u}}^l$. In most cases, the matrix $\tilde{U}$, associated with the arbitrary map $\mathsf{\Phi }(·)$, does not necessarily possess zero column sum for PCA. Subtracting the mean $\overline{\mathsf{\Phi }}=(1/n)\sum \mathsf{\Phi }\left({\mathit{u}}^i\right)$ from each column results in a data matrix

$$\tilde{\overline{U}}=[\tilde{\mathsf{\Phi }}\left({\mathit{u}}^1\right)\cdots \tilde{\mathsf{\Phi }}\left({\mathit{u}}^n\right)]=\tilde{U}H\in {\mathbb{R}}^{M\times n},$$

with zero column sum, where the columns are $\tilde{\mathsf{\Phi }}\left({\mathit{u}}^l\right)=\mathsf{\Phi }\left({\mathit{u}}^l\right)-\overline{\mathsf{\Phi }}$, and H is the centering matrix defined by

$$H=I_n-{\displaystyle \frac{1}{n}}\mathbf{1}{\mathbf{1}}^T\in {\mathbb{R}}^{n\times n},$$

(8)

where $I_n$ is the n-dimensional identity matrix and $\mathbf{1}={[1,\dots ,1]}^T\in {\mathbb{R}}^n$ is the vector of ones. Next, we apply the standard PCA steps described earlier, to the data matrix $\tilde{\overline{U}}$ which has columns that represent the feature space. This involves identifying the eigenvectors of the covariance matrix

$$\tilde{C}=\sum _{l=1}^n\tilde{\mathsf{\Phi }}\left({\mathit{u}}^l\right)\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^l\right)}^T=\tilde{\overline{U}}{\tilde{\overline{U}}}^T\in {\mathbb{R}}^{M\times M}.$$

(9)

Currently, the KPCA has two significant disadvantages. In the first place, the size M could be so large that it becomes nearly impossible to calculate the eigenvectors of the M-dimensional covariance matrix $\tilde{C}$. Second, the random nonlinear function $\mathsf{\Phi }(·)$ is often inaccessible. To tackle these problems, a kernel trick is utilized [16,17]. To comprehend this technique, let’s examine the eigenvalue problem

$$\begin{array}{ccc}\hfill \tilde{C}{\mathit{v}}^i& ={\tilde{\lambda }}_i{\mathit{v}}^i,\hfill & \hfill i=1,\dots ,M,\end{array}$$

(10)

where $\{{\tilde{\lambda }}_i,{\mathit{v}}^i\}$ represent the eigenpair of the covariance matrix $\tilde{C}$. It is important to understand that, by definition, the eigenvectors cover the feature space, which is equivalent to the column space of the altered data matrix $\tilde{\overline{U}}=[\tilde{\mathsf{\Phi }}\left({\mathit{u}}^1\right)\cdots \tilde{\mathsf{\Phi }}\left({\mathit{u}}^n\right)]$. Consequently, there are real coefficients $a_{ij}$ for each eigenvector ${\mathit{v}}^i$ that satisfy the linear combination

$$\begin{array}{ccc}\hfill {\mathit{v}}^i& =\sum _{j=1}^n{a}_{ij}\tilde{\mathsf{\Phi }}\left({\mathit{u}}^j\right),\hfill & \hfill i=1,\dots ,M.\end{array}$$

(11)

By replacing the connection (11) and the identity (9) into the equation (10), we get that

$$\sum _{l=1}^n\tilde{\mathsf{\Phi }}\left({\mathit{u}}^l\right)\sum _{j=1}^n{a}_{ij}\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^l\right)}^T\tilde{\mathsf{\Phi }}\left({\mathit{u}}^j\right)={\tilde{\lambda }}_i\sum _{j=1}^n{a}_{ij}\tilde{\mathsf{\Phi }}\left({\mathit{u}}^j\right).$$

All eigenvectors are within the range of the altered vectors ${\left\{\tilde{\mathsf{\Phi }}\left({\mathit{u}}^s\right)\right\}}_{s=1}^n$, so for $s=1,\dots ,n$, we can examine the equivalent equations by projecting onto the vectors $\tilde{\mathsf{\Phi }}\left({\mathit{u}}^s\right)$, resulting in

$$\begin{array}{c}\hfill \sum _{l=1}^n\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^s\right)}^T\tilde{\mathsf{\Phi }}\left({\mathit{u}}^l\right)\sum _{j=1}^n{a}_{ij}\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^l\right)}^T\tilde{\mathsf{\Phi }}\left({\mathit{u}}^j\right)={\tilde{\lambda }}_i\sum _{j=1}^n{a}_{ij}\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^s\right)}^T\tilde{\mathsf{\Phi }}\left({\mathit{u}}^j\right).\end{array}$$

(12)

At this stage, a kernel function $\kappa (·,·):{\mathbb{R}}^m\times {\mathbb{R}}^m\mapsto \mathbb{R}$ is defined so that

$$\begin{array}{ccc}\hfill \kappa ({\mathit{u}}^s,{\mathit{u}}^l)& =\langle \mathsf{\Phi }\left({\mathit{u}}^s\right),\mathsf{\Phi }\left({\mathit{u}}^l\right)\rangle =\mathsf{\Phi }{\left({\mathit{u}}^s\right)}^T\mathsf{\Phi }\left({\mathit{u}}^l\right),\hfill & \hfill s,l=1,\dots ,n,\end{array}$$

with the goal of depicting the Euclidean inner products $\mathsf{\Phi }{\left({\mathit{u}}^s\right)}^T\mathsf{\Phi }\left({\mathit{u}}^l\right)$ of non-centered transformed vectors in the feature space using input space vectors. Different types of kernel functions are utilized in various studies, including linear, polynomial, and Gaussian kernels [16,25]. Here, we employ the Gaussian kernel defined for any vector $\mathit{x},\mathit{y}$ by

$$\kappa (\mathit{x},\mathit{y})=exp\left(-{\displaystyle \frac{{\parallel \mathit{x}-\mathit{y}\parallel }^2}{2{\sigma }^2}}\right),$$

where $\parallel ·\parallel$ represents the Euclidean norm and $\sigma$ being a parameter.

In addition, to depict the Euclidean inner products $\tilde{\mathsf{\Phi }}{(·)}^T\tilde{\mathsf{\Phi }}(·)$ of transformed vectors that are centered in the feature space, we utilize the framework [14]

$$\tilde{\kappa }({\mathit{u}}^s,{\mathit{u}}^l)=\kappa ({\mathit{u}}^s,{\mathit{u}}^l)-{\displaystyle \frac{1}{n}}{\mathbf{1}}^T{\mathit{k}}_{{\mathit{u}}^s}-{\displaystyle \frac{1}{n}}{\mathbf{1}}^T{\mathit{k}}_{{\mathit{u}}^l}+{\displaystyle \frac{1}{n^2}}{\mathbf{1}}^{T}K\mathbf{1},$$

where the kernel matrix $K\in {\mathbb{R}}^{n\times n}$, and the vector ${\mathit{k}}_{\mathit{u}}\in {\mathbb{R}}^n$ are defined as

$$K_{ij}=\kappa ({\mathit{u}}^i,{\mathit{u}}^j),{\mathit{k}}_{\mathit{u}}={(\kappa (\mathit{u},{\mathit{u}}^1),\dots ,\kappa (\mathit{u},{\mathit{u}}^n))}^T.$$

We note that it is not necessary to compute all the vectors ${\mathit{k}}_{{\mathit{u}}^i}$, as they are simply the ith columns of the symmetric kernel matrix K. The use of the kernel function results in an equivalent equation to (12) being

$$\begin{array}{c}\hfill \sum _{l=1}^n\kappa ({\mathit{u}}^s,{\mathit{u}}^l)\sum _{j=1}^n{a}_{ij}\kappa ({\mathit{u}}^l,{\mathit{u}}^j)={\tilde{\lambda }}_i\sum _{j=1}^n{a}_{ij}\kappa ({\mathit{u}}^s,{\mathit{u}}^j).\end{array}$$

(13)

The equation (13) can be represented in matrix-vector form as

$$\begin{array}{ccc}\hfill {\tilde{K}}^2{\mathit{a}}^i& ={\tilde{\lambda }}_i\tilde{K}{\mathit{a}}^i,\mathrm{or}\tilde{K}{\mathit{a}}^i={\tilde{\lambda }}_i{\mathit{a}}^i,\hfill & \hfill i=1,\dots ,M,\end{array}$$

(14)

for the coefficient vector ${\mathit{a}}^i={(a_{i1},\dots ,a_{in})}^T\in {\mathbb{R}}^n$, where $\tilde{K}=HKH$ and H is the centering matrix given in (8). Ultimately, for any vector ${\mathit{u}}^{*}\in {\mathbb{R}}^m$ in the input space, its transformed vector $\tilde{\mathsf{\Phi }}\left({\mathit{u}}^{*}\right)\in {\mathbb{R}}^M$ in the feature space is found. The projection vector ${\mathit{z}}^{*}={(z_1^{*},\dots ,z_k^{*})}^T\in {\mathbb{R}}^k$ represents the projection of $\mathsf{\Phi }\left({\mathit{u}}^{*}\right)$ onto the reduced k-dimensional space ($k\ll m\ll M$) spanned by the eigenvectors $\{{\mathit{v}}^1,\dots ,{\mathit{v}}^k\}$ associated with the largest k eigenvalues $\{{\tilde{\lambda }}_1,\dots ,{\tilde{\lambda }}_k\}$. After calculating the coefficients ${\mathit{a}}^i$ from the eigenvalue problem (14) and utilizing the identity (11) of the eigenvectors ${\mathit{v}}^i$, the projection vectors’ entries $z_i^{*}$ can be determined using the kernel function as the following

$$\begin{array}{ccc}\hfill z_i^{*}& =\tilde{\mathsf{\Phi }}{\left({\mathit{u}}^{*}\right)}^T{\mathit{v}}^i=\sum _{j=1}^n{a}_{ij}\tilde{\kappa }({\mathit{u}}^{*},{\mathit{u}}^j),\hfill & \hfill i=1,\dots ,k.\end{array}$$

Reconstructing the original image of any random vector ${\mathit{u}}^{*}\in {\mathbb{R}}^m$ in the input space can be achieved by approximating its pre-image ${\widehat{\mathit{u}}}^{*}\in {\mathbb{R}}^m$ using standard PCA as shown in equation (7). Nonetheless, KPCA does not follow this pattern. Denote $P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)\in {\mathbb{R}}^k$ as the projection of $\mathsf{\Phi }\left({\mathit{u}}^{*}\right)$ onto the reduced space spanned by the eigenvectors $\{{\mathit{v}}^1,\dots ,{\mathit{v}}^k\}$, namely

$$P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)=\sum _{i=1}^k{z}_i^{*}{\mathit{v}}^i+\overline{\mathsf{\Phi }}.$$

Next, we can find an estimated pre-image ${\widehat{\mathit{u}}}^{*}$ such that the vector $\mathsf{\Phi }\left({\widehat{\mathit{u}}}^{*}\right)$ after transformation is the vector closest to the projected vector $P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)$. This involves determining the minimum of the objective functional $\rho \left({\widehat{\mathit{u}}}^{*}\right)={\parallel \mathsf{\Phi }\left({\widehat{\mathit{u}}}^{*}\right)-P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)\parallel }^2$ using a least-squares approach. By imposing ${\nabla }_{{\widehat{\mathit{u}}}^{*}}\rho =0$, we obtain the following equation to be solved

$${\widehat{\mathit{u}}}^{*}={\displaystyle \frac{{\sum }_{i=1}^k{\tilde{\gamma }}_{i}exp\left(-\frac{\parallel {\widehat{\mathit{u}}}^{*}-{\mathit{u}}^i{\parallel }^2}{2{\sigma }^2}\right){\mathit{u}}^i}{{\sum }_{i=1}^k{\tilde{\gamma }}_{i}exp\left(-\frac{\parallel {\widehat{\mathit{u}}}^{*}-{\mathit{u}}^i{\parallel }^2}{2{\sigma }^2}\right)}},$$

(15)

where we set

$${\tilde{\gamma }}_i={\gamma }_i+{\displaystyle \frac{1}{n}}\left(1-\sum _{j=1}^n{\gamma }_j\right),{\gamma }_i=\sum _{l=1}^k{z}_l^{*}{a}_{li}.$$

The implicit equation (15) can be solved by fixed-point iteration or other Newton type nonlinear iterative solvers. But, the nonlinear iteration technique is unstable and depends highly on the initial estimation [14,25]. In the KPCA method presented in this study, the pre-image is calculated using a non-iterative approach [14], which relies on the correlation between the distance of vectors in the input space and the distance of the vectors in the transformed space, employing the Gaussian kernel function as a distance metric function $\parallel ·\parallel$. By utilizing MDS concept [14,19], we can roughly discover a pre-image ${\widehat{\mathit{u}}}^{*}$ where the differences in distances between ${\mathit{u}}^{*}$ and each input vector ${\mathit{u}}^i$ and between projected vector $P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)$ and each transformed feature vector $\mathsf{\Phi }\left({\mathit{u}}^i\right)$ remain constant, see Figure 1. To achieve this goal, let

$$d_{ij}^2:=d^2({\mathit{u}}^i,{\mathit{u}}^j)=\parallel {\mathit{u}}^i-{\mathit{u}}^j{\parallel }^2,{\tilde{d}}_{ij}^2:={\tilde{d}}^2(\mathsf{\Phi }\left({\mathit{u}}^i\right),\mathsf{\Phi }\left({\mathit{u}}^j\right))={\parallel \mathsf{\Phi }\left({\mathit{u}}^i\right)-\mathsf{\Phi }\left({\mathit{u}}^j\right)\parallel }^2,$$

represent the distance metrics between input space vectors and their transformed feature space vectors. After manipulating the Gaussian kernel function as $f\left(d_{ij}^2\right)=exp(-d^2({\mathit{u}}^i,{\mathit{u}}^j)/2{\sigma }^2)$, we can establish a relation between input space distance and feature space distance, and also its inverse map, as [14]

$$f\left(d_{ij}^2\right)={\displaystyle \frac{1}{2}}(K_{ii}+K_{jj}-{\tilde{d}}_{ij}^2),d_{ij}^2=-2{\sigma }^{2}ln\left({\displaystyle \frac{K_{ii}+K_{jj}-{\tilde{d}}_{ij}^2}{2}}\right).$$

(16)

In addition, it can be demonstrated using the kernel matrix that the distance in feature space between the projected vector $P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)$ and a transformed vector $\mathsf{\Phi }\left({\mathit{u}}^i\right)$ is

$$\begin{array}{cc}\hfill & {\tilde{d}}^2(\mathsf{\Phi }\left({\mathit{u}}^i\right),P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)={\parallel \mathsf{\Phi }\left({\mathit{u}}^i\right)-P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)\parallel }^2\hfill \\ \hfill & =\parallel \mathsf{\Phi }\left({\mathit{u}}^i\right){\parallel }^2+{\parallel P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)\parallel }^2-2P_k\mathsf{\Phi }{\left({\mathit{u}}^{*}\right)}^T\mathsf{\Phi }\left({\mathit{u}}^i\right)\hfill \\ \hfill & =\left({\mathit{k}}_{{\mathit{u}}^{*}}+{\displaystyle \frac{1}{n}}K\mathbf{1}-2{\mathit{k}}_{{\mathit{u}}^i}\right)H^T{C}_{a}H\left({\mathit{k}}_{{\mathit{u}}^{*}}-{\displaystyle \frac{1}{n}}K\mathbf{1}\right)+{\displaystyle \frac{1}{n^2}}{\mathbf{1}}^{T}K\mathbf{1}+K_{ii}-{\displaystyle \frac{2}{n}}{\mathbf{1}}^T{\mathit{k}}_{{\mathit{u}}^{*}},\hfill \end{array}$$

(17)

for the matrix

$$C_a=\sum _{j=1}^k{\displaystyle \frac{1}{{\tilde{\lambda }}_j}}{\mathit{a}}^j{\left({\mathit{a}}^j\right)}^T.$$

Ultimately, by taking $P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)\approx \mathsf{\Phi }\left({\widehat{\mathit{u}}}^{*}\right)$ and utilizing the metric relation (16) within equation (15), we can derive the non-iterative solution formula [14]

$$\begin{array}{cc}\hfill {\widehat{\mathit{u}}}^{*}& ={\displaystyle \frac{{\sum }_{i=1}^k{\tilde{\gamma }}_{i}exp\left(-\frac{\parallel {\widehat{\mathit{u}}}^{*}-{\mathit{u}}^i{\parallel }^2}{2{\sigma }^2}\right){\mathit{u}}^i}{{\sum }_{i=1}^k{\tilde{\gamma }}_{i}exp\left(-\frac{\parallel {\widehat{\mathit{u}}}^{*}-{\mathit{u}}^i{\parallel }^2}{2{\sigma }^2}\right)}}\hfill \\ \hfill & ={\displaystyle \frac{{\sum }_{i=1}^k{\tilde{\gamma }}_i\left(\frac{1}{2}(2-{\tilde{d}}^2(\mathsf{\Phi }\left({\widehat{\mathit{u}}}^{*}\right),\mathsf{\Phi }\left({\mathit{u}}^i\right)))\right){\mathit{u}}^i}{{\sum }_{i=1}^k{\tilde{\gamma }}_i\left(\frac{1}{2}(2-{\tilde{d}}^2(\mathsf{\Phi }\left({\widehat{\mathit{u}}}^{*}\right),\mathsf{\Phi }\left({\mathit{u}}^i\right)))\right)}}\hfill \\ \hfill & \approx {\displaystyle \frac{{\sum }_{i=1}^k{\tilde{\gamma }}_i\left(\frac{1}{2}(2-{\tilde{d}}^2(P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right),\mathsf{\Phi }\left({\mathit{u}}^i\right)))\right){\mathit{u}}^i}{{\sum }_{i=1}^k{\tilde{\gamma }}_i\left(\frac{1}{2}(2-{\tilde{d}}^2(P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right),\mathsf{\Phi }\left({\mathit{u}}^i\right)))\right)}},\hfill \end{array}$$

(18)

where the distance $\parallel P_k\mathsf{\Phi }\left({\mathit{u}}^{*}\right)-\mathsf{\Phi }\left({\mathit{u}}^i\right){\parallel }^2$ in the formula is determined through relation (17). Furthermore, the k vectors ${\mathit{u}}^i$ mentioned in (18) are chosen from the n input space vectors that are closest to the given vector ${\mathit{u}}^{*}$ [14,25].

4. Numerical Results

In this section, we investigate how well the KPCA nonlinear dimensionality reduction method performs on one-, two- and three-dimensional AC equations. In either example, we construct the data matrix U whose columns are formed by the discrete solution vectors of the given AC equation, obtained by the full discrete system discussed in Section 2. In order to measure the accuracy of a reduced approximation ${\widehat{\mathit{u}}}^{*}$ to a fixed full solution ${\mathit{u}}^{*}$, we use the relative absolute error defined by

$$\parallel {\mathit{u}}^{*}-{\widehat{\mathit{u}}}^{*}{\parallel }_R={\displaystyle \frac{|{\mathit{u}}^{*}-{\widehat{\mathit{u}}}^{*}|}{|{\mathit{u}}^{*}|}}.$$

4.1. One-Dimensional Problem

The AC equation (1) will be numerically modeled over the one-dimensional domain $\mathsf{\Omega }=[0,2\pi ]$ with periodic boundary conditions, and with the target time $T=600$. The inter-facial width is set as $\mu =0.16$, and the initial phase is

$$u(x,0)=0.8+sin\left(x\right).$$

Data matrix U is created using discrete solutions at each time $t_k$ taken from the instant information of the system with space-time discretization, with the mesh sizes $\Delta x=\pi /100$ and $\Delta t=0.5$ [6]. The discretized data matrix for space and time coordinates is written as $U=[{\mathit{u}}^1\cdots {\mathit{u}}^{1200}]\in {\mathbb{R}}^{200\times 1200}$.

Figure 2 displays the plots of the precise solutions at $t=0,250,400,500$, and also the solution trajectory for the time span $[0,600]$. The AC equation with the potential function $F\left(u\right)={(u^2-1)}^2/4$ has one stable ($u=0$) and two unstable ($u=\pm 1$) equilibrium points. This phenomenon, known as phase separation, involves solutions transitioning between the equilibrium points. The boundaries separating the two unstable equilibrium zones shift across areas over extended periods, known as meta-stability. Figure 2, left, shows this state of meta-stability.

Table 1 compares PCA and KPCA methods in terms of accuracy and computational efficiency. For this comparison, the input space solution vector ${\mathit{u}}^{*}={\mathit{u}}^{100}$ corresponding to the solution at time $t=50$ is used. Results are presented for reduced dimensions $k=1,\dots ,6$. The relative absolute errors between the full solution ${\mathit{u}}^{*}$ and the reduced approximate solution ${\widehat{\mathit{u}}}^{*}$ are shown when PCA and KPCA are used. When using the same number of reduced bases, more accurate results are obtained with KPCA.

On the other hand, in terms of computational efficiency, Table 2 shows the solution times required by the KPCA method using the formula (15), which is solved by fixed point iteration, and the method using the algebraic equation (18). The first two columns show the errors between the full and reduced solutions obtained with both formulae, while the last two columns show the processing time required to generate the preliminary images. It can be seen that the errors obtained with the solutions of the non-linear equation solved with iteration (2-3 iterations) with the formula (15) and without iteration with the formula (18) are almost the same, which shows that the solutions obtained without iteration are acceptable in terms of precision. It can be seen in the last two columns that the solution time required by the KPCA method without iteration is much less than the solution time required with iteration, i.e., the KPCA method is quite fast.

Figure 3 shows the reduced solutions together with the full solutions at times $t=250$ and $t=500$. For the reduced dimension $k=2$, it can be seen that the exact and reduced solutions coincide at both times. The full and reduced solution profiles are given in Figure 4. From the figures, it is observed that the full and reduced solutions coincide, and hence show the same phase transition behavior.

The energy graphs obtained with full solutions and preliminary images are given in Figure 5. It can be seen that the same energy reduction behavior is observed in both graphs.

4.2. Two-Dimensional Problem

The AC equation (1) will be numerically modeled over the two-dimensional domain $\mathsf{\Omega }={[0,2\pi ]}^2$ with periodic boundary conditions, and with the target time $T=5$. The inter-facial width is set as $\mu =0.05$, and the initial condition is

$$u(x_1,x_2,0)=2e^{sin\left(x_1\right)+sin\left(x_2\right)-2}+2.2e^{-sin\left(x_1\right)-sin\left(x_2\right)-2}-1.$$

Data matrix U is created using discrete solutions with the mesh sizes $\Delta x_1=\Delta x_2=\pi /32$ and $\Delta t=0.01$ [6]. The discretized data matrix for this setting is formed as $U=[{\mathit{u}}^1\cdots {\mathit{u}}^{500}]\in {\mathbb{R}}^{4096\times 500}$.

Table 3 compares the PCA and KPCA methods in terms of accuracy and computational efficiency. For this comparison, the solution vector ${\mathit{u}}^{*}={\mathit{u}}^{100}$ input space corresponding to the solution at time $t=1$ is used. Results are presented for reduced dimensions $k=1,\dots ,6$. The relative absolute errors between the exact solution ${\mathit{u}}^{*}$ and the reduced approximate solution ${\widehat{\mathit{u}}}^{*}$ are shown when PCA and KPCA are used. When using the same number of bases, more accurate results are obtained with KPCA.

In terms of computational efficiency, Table 4 shows the solution times required by the KPCA method using the formula (15) solved by fixed point iteration and the method using the algebraic equation (18). The first two columns show the errors between the full and reduced solutions obtained with both formulas, while the last two columns show the processing time required to generate the preliminary images. It can be seen that the errors obtained with the solutions of the non-linear equation solved with iteration (3 iterations) and without iteration are almost the same, which shows that the solutions obtained without iteration are acceptable in terms of precision. It can be seen in the last two columns that the solution time required by the KPCA method without iteration is considerably less than the solution time required with iteration.

Figure 6 shows the initial profile and profiles of the reduced solution together with the full solution at the final time $t=5$. For the reduced dimension $k=2$, the full and reduced solutions are similar. This means that the reduced solutions have the same behavior as the full one.

On the other hand, the energy plots obtained with full solutions and preliminary images are given in Figure 7. Here, it is clearly seen that the same energy reduction behavior occurs in both plots, similar to the one-dimensional AC equation.

4.3. Three-Dimensional Problem

In this last example, the scaled AC equation [6,26,27]

$$u_t=\Delta u-{\displaystyle \frac{1}{{\mu }^2}}f\left(u\right),$$

will be numerically modeled over the three-dimensional domain $\mathsf{\Omega }={[0,1]}^3$ with homogenous Neumann boundary conditions, and with the target time $T=0.035$. The inter-facial width is set as $\mu =0.1$, and the initial condition is

$$u(x_1,x_2,x_3,0)=tanh\left({\displaystyle \frac{R_0-\sqrt{{(x_1-0.5)}^2+{(x_2-0.5)}^2+{(x_3-0.5)}^2}}{\sqrt{2}\mu }}\right).$$

Here, the given parameter $R_0$ represents the initial radius of the spherical zero isosurface, and at time t, the exact radius R of the sphere follows from the formula $R\left(t\right)=\sqrt{R_0^2-4t}$ [26]. This indicates a decrease in the radius of the spherical zero isosurface over time. In the numerical simulation, we take the initial radius as $R_0=0.4$.

For the discrete solution matrix U, we apply the mesh sizes $\Delta x_1=\Delta x_2=\Delta x_3=0.05$ and $\Delta t=0.0001$, i.e., the discretized data matrix is $U=[{\mathit{u}}^1\cdots {\mathit{u}}^{350}]\in {\mathbb{R}}^{9261\times 350}$. Table 5 compares the PCA and KPCA methods in terms of accuracy and computational efficiency. For this comparison, the solution vector ${\mathit{u}}^{*}={\mathit{u}}^{100}$ input space corresponding to the solution at time $t=0.01$ is used. Results are presented for reduced dimensions $k=1,\dots ,6$. The relative absolute errors between the full solution ${\mathit{u}}^{*}$ and the reduced approximate solution ${\widehat{\mathit{u}}}^{*}$ are shown when PCA and KPCA are used. When using the same number of bases, more accurate results are obtained with KPCA.

In terms of computational efficiency, Table 6 shows the solution times required by the KPCA method with and without iteration. It can be seen that the errors obtained with the solutions of the non-linear equation solved with iteration (2 iterations) and without iteration are almost the same, similar to the one- and two-dimensional cases, and that the KPCA method without iteration is fastest.

In Figure 8, left, the energy plots obtained with full and reduced solutions are given. Again, it is clearly seen that the same energy reduction behavior occurs by both full and reduced solutions. Figure 9 shows the profiles of the zero isosurface by full and reduced solutions at times $t=0.02,0.03,0.035$. For the reduced dimension $k=2$, the full and reduced profiles are highly similar, the radius of the spherical zero isosurface decrease as the time progresses. The decrease of the exact and numerical radii is also demonstrated in Figure 8, right.

5. Discussion

In this article, we suggest a nonlinear dimensionality reduction method that upholds precision and conserves the energy decay property of the one-, two-, and three-dimensional AC equations. The method for reducing dimensionality is non-intrusive and relies on the KPCA. The KPCA utilizes the MDS method and considers the input space and feature space distance metrics to develop a non-iterative algorithm for reconstructing pre-images. The numerical test examples show the accuracy of solutions and the preservation of energy dissipation structure. A comparison is made with traditional PCA and KPCA, as well as an iterative scheme with a non-iterative one. In our future research, we aim to explore various types of kernel functions that may offer more effective performance compared to the Gaussian kernel.