Numerical approximation of the first $p$-Laplace eigenpair

1. Problem Statement

The paper addresses the numerical approximation of the first Dirichlet eigenpair $(\lambda_p, u_p)$ for the $p$-Laplace operator on both Euclidean domains and surface manifolds embedded in $\mathbb{R}^3$. The operator is defined as:
$$\Delta_p u = \text{div}(|\nabla u|^{p-2}\nabla u)$$
The study focuses specifically on the regime where $p$ is large ($2 \le p < \infty$), approaching the limit $p \to \infty$.

Key Challenges:

• Nonlinearity and Degeneracy: As $p$ increases, the operator becomes highly nonlinear and degenerate. The eigenfunctions lose smoothness, belonging to the regularity class $C^{1,\sigma}$ where $\sigma \to 0$ as $p \to \infty$. This makes standard numerical methods unstable or inaccurate.
• Numerical Conditioning: The eigenvalues $\lambda_p$ grow or decay exponentially with $p$ depending on the domain size, leading to numerical overflow or underflow.
• Surface Complexity: Existing literature largely focuses on Euclidean domains or low $p$ values. Extending these methods to curved surfaces (manifolds) for large $p$ is an open challenge.
• Geometric Connection: The limit $p \to \infty$ connects the eigenpair to the geometry of the domain, specifically the distance-to-boundary function. However, bridging the gap between finite $p$ computations and the asymptotic limit has been difficult.

2. Methodology

The authors propose a robust numerical framework combining Surface Finite Element Methods (SFEM) with a Newton inverse-power iteration scheme and a novel domain rescaling strategy.

A. Discretization (SFEM)

• Space: The continuous domain $\Omega$ is discretized using a surface mesh $\Omega_h$. The authors use bilinear ($Q_1$) Lagrange elements on quadrilateral meshes.
• Rationale: Since eigenfunctions for $p > 2$ are not $C^2$, higher-order elements are not beneficial. $Q_1$ elements are sufficient and computationally efficient.
• Weak Formulation: The problem is solved in the weak sense, seeking $(\lambda_p, u_p) \in \mathbb{R}^+ \times W^{1,p}_0(\Omega)$ such that:
  $$\int_\Omega |\nabla_S u_p|^{p-2} \nabla_S u_p \cdot \nabla_S \phi \, dA = \lambda_p \int_\Omega |u_p|^{p-2} u_p \phi \, dA$$

B. Newton Inverse-Power Iteration

To solve the nonlinear system, the authors employ an alternating update scheme:

1. Eigenvalue Update: Computed via the discrete nonlinear Rayleigh quotient.
2. Eigenfunction Update: Solved using Newton's method. The nonlinear term is regularized with a small parameter $\eta$ ($\gamma(u;\eta) = (\eta^2 + |\nabla u|^2)^{(p-2)/2}$) to ensure the stiffness matrix remains positive definite.
3. Linear Solver: The resulting linear systems are solved using the Conjugate Gradient (CG) method preconditioned with Algebraic Multigrid (AMG). This is crucial for handling the ill-conditioning that arises at large $p$.

C. Continuation in $p$ and Domain Rescaling

• Continuation: Instead of solving for a large $p$ directly, the algorithm starts at $p=2$ (linear Laplacian) and incrementally increases $p$ (e.g., $\delta p = 1$), using the solution from $p_{k}$ as the initial guess for $p_{k+1}$.
• Domain Rescaling Strategy: This is a key innovation. The eigenvalues scale as $\lambda_p \sim (\max \text{dist})^{-p}$. To prevent numerical blowup (if the domain is large) or underflow (if small), the domain is dynamically rescaled.
  • A scaling factor $\alpha$ is chosen such that the maximum distance to the boundary is close to 1.
  • An adaptive thresholding mechanism monitors the computed eigenvalue. If $\lambda_p$ exceeds a threshold (e.g., $>20$ or $<2$), the domain is rescaled, and the eigenvalue/eigenfunction are adjusted accordingly. This allows stable computation for $p$ up to 100.

3. Key Contributions

1. Extension to Surface Domains: The paper presents the first systematic treatment of the $p$-Laplacian eigenvalue problem on surface manifolds (codimension 1 in $\mathbb{R}^3$) for a wide range of $p$.
2. Large $p$ Capability: The method successfully computes eigenpairs for $p$ as large as 100, significantly exceeding previous limits (typically $p \le 10$ or restricted to unit balls).
3. Stabilization via Rescaling: The introduction of the adaptive domain rescaling strategy effectively mitigates the extreme growth/decay of eigenvalues, enabling stable convergence where standard methods fail.
4. Verification of Asymptotics: The numerical results confirm theoretical predictions regarding the $p \to \infty$ limit, specifically the convergence of the eigenfunction to the normalized distance-to-boundary function (for symmetric domains) and the convergence of $\lambda_p^{1/p}$ to the inverse of the inradius.

4. Numerical Results

The authors validate their method through extensive experiments:

• 1D Case (Benchmark): Compared against known analytical solutions involving generalized trigonometric functions ($\sin_p$).
  • Convergence: Eigenvalues converge at $\approx$ second order for small $p$, dropping slightly as $p$ increases due to reduced regularity. Eigenfunctions converge at $\approx$ first order.
  • Asymptotics: Confirmed that $\lambda_p^{1/p} \to 1/R$ and the rate of convergence is $\mathcal{O}((\log p)/p)$.
• Euclidean Planar Domains:
  • Disk: The limiting eigenfunction converges to the distance-to-boundary function. Results show good agreement with existing literature for $p \le 10$ and extend to $p=100$.
  • Square: The limiting eigenfunction is not the distance function (due to the ridge set $R$ not equaling the maximal set $M$). The method successfully captures the nonsmooth behavior along the diagonals.
• Surface Domains:
  • Simple Surfaces (Hemisphere, Half-Torus): For the hemisphere (where $R=M$), the eigenfunction converges to the geodesic distance function. The half-torus shows behavior similar to 1D due to its geometry.
  • Complex Surfaces (Duck, Hand, Mug, Pig): The method handles complex, non-symmetric geometries. For domains where the boundary is a point or small region, the eigenfunctions show high sensitivity at low $p$ but stabilize into a distance-like structure as $p$ increases.
• Monotonicity: The results verify the theoretical inequality $p(\lambda_p)^{1/p} < s(\lambda_s)^{1/s}$ for $p < s$, confirming the stability of the numerical scheme.

5. Significance and Applications

• Geometric Analysis: The work provides a computational tool to explore the deep connection between nonlinear spectral theory and the intrinsic geometry of domains (specifically the distance function).
• Shape Optimization: The first eigenpair is a critical component in shape optimization algorithms. Accurate computation for large $p$ is a prerequisite for optimizing shapes based on $p$-Laplacian criteria.
• Data Science: The $p$-Laplacian is used in graph-based clustering and dimensionality reduction. While current data applications often use $p \le 2$, this work lays the groundwork for exploring the non-convex, large-$p$ regime in geometric data analysis.
• Algorithmic Robustness: The combination of Newton iteration, AMG preconditioning, and adaptive rescaling offers a template for solving other degenerate nonlinear PDEs on complex manifolds.

In summary, the paper successfully bridges the gap between finite $p$ computations and the $p \to \infty$ limit, providing a robust, scalable numerical framework for analyzing the first eigenpair of the $p$-Laplacian on both flat and curved domains.