Extrinsic Finite-Difference Scheme

• Extrinsic finite-difference schemes are numerical methods that approximate derivatives using function evaluations in the ambient space, applicable to manifolds and fractal subsets.
• They enable derivative-free optimization and efficient PDE analysis by bypassing intrinsic computations and complex spectral decompositions.
• Adaptive algorithms leveraging extrinsic differences offer practical benefits, including dimension-linear cost and on-the-fly parameter tuning for improved accuracy.

An extrinsic finite-difference scheme is a class of numerical methods that approximate derivatives, gradients, or differential operators by evaluating function values at points in the ambient Euclidean space, potentially outside the geometric domain of interest, such as a manifold or fractal subset. This approach extends standard finite-difference techniques to contexts where the solution space is embedded in a higher-dimensional space or possesses a nontrivial geometric or topological structure. Extrinsic finite differences are particularly valuable in derivative-free optimization on Riemannian submanifolds and in the analysis of PDEs on self-similar fractal geometries, enabling efficient and accurate computations without recourse to intrinsic (on-domain-only) operations or spectral decompositions.

1. Defining the Extrinsic Finite-Difference Gradient on Manifolds

For a smooth, $d$-dimensional Riemannian submanifold $M \subset \mathbb{R}^n$ and a differentiable function $f:\mathbb{R}^n \to \mathbb{R}$ evaluable anywhere in the ambient space, the extrinsic finite-difference gradient at $x \in M$ is constructed as follows. An orthonormal basis $\{e_1(x),\dots,e_d(x)\}$ of the tangent space $T_xM$ is selected. With a small parameter $h > 0$, the finite-difference gradient is

$$g_h(x) := \sum_{i=1}^d \frac{f(x + h\,e_i(x)) - f(x)}{h} \, e_i(x) \in T_xM.$$

Equivalently, letting $$U(x) = [e_1(x), \dots, e_d(x)] \in \mathbb{R}^{n \times d}$$,

$$g_h(x) = \mathrm{Proj}_x\left( \frac{f(x + h\,U(x)) - f(x)}{h} \, U(x) \right),$$

where $\mathrm{Proj}_x$ denotes orthogonal projection onto $T_xM$ in the ambient Euclidean norm. Under standard smoothness assumptions, this construction approximates the Riemannian gradient $\mathrm{grad} f(x)$ with $O(h)$ consistency error. Unlike intrinsic schemes, extrinsic finite differences require $f$ to be defined in a neighborhood outside $M$ and may offer computational advantages by avoiding repeated retractions onto the manifold (Taminiau et al., 13 Jan 2026).

2. Algorithmic Structure and Adaptive Accuracy

A fundamental instance is the extrinsic derivative-free Riemannian optimization (DFRO) algorithm, designed to minimize $f(x)$ over $x \in M$ using only function evaluations. The method maintains estimates for the ambient smoothness constant (denoted $\tau_k$) and the step-size controller ($\sigma_k$). The finite-difference parameter is set adaptively at each iteration via $h_k = \frac{2\epsilon}{5\sqrt{d}\tau_k}$ for target accuracy $\epsilon > 0$. Armijo-type sufficient decrease conditions guide adaptive selection of $\sigma_k$ and $\tau_k$ until the Riemannian gradient norm is below the desired threshold. These updates allow the method to efficiently learn near-optimal step and accuracy parameters, with no prior knowledge of Lipschitz constants required (Taminiau et al., 13 Jan 2026).

Algorithm Outline (Ext-RFD):

Step	Description	Formula/Condition
1	Compute basis and $g_k$	$h_k = \frac{2\epsilon}{5\sqrt{d}\tau_k}$, $g_k$ as above
2	Check norm, refine if necessary	If $\|g_k\| < \frac{4\epsilon}{5}$: $\tau_k \leftarrow 2\tau_k$
3	Armijo-type line search for $\sigma_k$	$$f(x_k) - f(x_k^+(\sigma_k)) \geq \frac{1}{4\sigma_k}\|g_k\|^2$$
4	Make step, update parameters	$x_{k+1} = x_k^+(\sigma_k)$, $\sigma_{k+1} = \sigma_k/2$
5	Terminate if criticality reached	If $\|\mathrm{grad} f(x_{k+1})\| \leq \epsilon$

3. Theoretical Properties: Consistency, Complexity, and Assumptions

The extrinsic FD scheme achieves $O(h)$ accuracy in approximating the Riemannian gradient, contingent on the ambient Lipschitz-gradient property (smoothness constant $L_E$):

$$\|g_h(x) - \mathrm{grad} f(x)\| \leq C_f h, \quad C_f = \frac{L_E \sqrt{d}}{2}$$

With further safeguards on FD parameter choice and descent conditions, one obtains matching relative-error lower and upper bounds and guarantees sufficient function decrease per iteration. The complexity to reach an $\epsilon$-criticality is

• $O(d\epsilon^{-2})$ function evaluations
• $O(\epsilon^{-2})$ retractions

where constants depend on initial value and smoothness parameters but not the ambient dimension $n$ (Taminiau et al., 13 Jan 2026).

Critical assumptions include:

• A3 (Ambient smoothness): Ensures accuracy for extrinsic FD schemes.
• A1 (Manifold smoothness): Needed for descent estimates after (retracted) steps.
• Global retraction: Needed for every iterated step to remain on $M$.
• Evaluability outside $M$: Essential for extrinsic differences; absent this, only intrinsic approaches are feasible.

4. Extrinsic Finite-Differences on Fractals: The Sierpiński Simplex Case

In the context of self-similar sets such as the Sierpiński gasket and tetrahedron, extrinsic finite-difference schemes approximate the Laplacian by operating on recursive graph-approximations embedded in the Euclidean space (Riane et al., 2018). For $d \in \{3, 4\}$, a sequence of finite vertex sets $V_m$ is constructed via contraction maps from the corners of a regular simplex. The unweighted graph Laplacian on this structure is

$$\Delta_m u(X) = \sum_{Y \sim_m X} (u(Y) - u(X)), \quad X \in V_m \setminus V_0$$

With appropriate renormalization constants $C_m$ derived from Kigami–Strichartz theory, the continuous Laplacian is approximated by

$$\Delta u(X) \approx C_m \Delta_m u(X), \quad C_m = \begin{cases} \frac{3}{2} 5^m, & d=3 \ 2 \cdot 6^m, & d=4 \end{cases}$$

This approach enables explicit construction of the Laplacian and subsequent time-stepping (Euler, implicit, Crank–Nicolson) for the heat equation on these fractal domains, with all computations performed extrinsically via the embedding (Riane et al., 2018).

5. Error, Stability, and Convergence Analyses

For PDEs on Sierpiński simplices, assuming the exact solution is Hölder-continuous with exponent $\alpha > 0$ in space and $C^1$ in time, local truncation errors decompose as

$\varepsilon^m_{k,i} = O(h) + O(2^{-m\alpha})$

where $h$ is the time step and $m$ defines the graph approximation level. Explicit schemes are stable under a CFL-type constraint:

$h \leq \frac{d^2}{2(d+2)^m}$

Implicit and Crank–Nicolson schemes are unconditionally stable. By synchronizing $h$ and $m$ so that $h^{-1}2^{-m\alpha}=O(1)$, the overall error converges at $O(2^{-m\alpha})$ (Riane et al., 2018).

6. Computational Advantages and Broader Implications

Extrinsic finite-difference schemes on manifolds and fractals bypass the need for spectral decompositions, diagonalizations, or explicit eigenvector/eigenvalue computations. On manifolds, they yield dimension-linear evaluation cost independent of ambient embedding, and on fractals, they leverage sparse Laplacian construction from Euclidean embedding and recursive subdivision. The enabling assumption in both contexts is evaluability of $f$ or the discrete Laplacian at off-domain or embedded points; where this is unavailable, only less computationally advantageous intrinsic formulations are viable. For optimization, adaptivity to unknown smoothness via on-the-fly parameter updates further reduces manual tuning (Taminiau et al., 13 Jan 2026, Riane et al., 2018).

7. Avoidance of Eigenvalue Approximations and Comparative Position

The extrinsic finite-difference approach distinguishes itself—particularly in the context of fractal domains—by avoiding reliance on spectral decimation and approximate spectral calculation. Through explicit recursive graph Laplacian construction, precise scaling, and elementary time-discretization, it provides a fully constructive algorithmic pathway for both solution and theoretical analysis. This is in contrast to methods necessitating prior computation or estimation of Laplacian spectra, thereby obviating major computational bottlenecks for large or complex structures (Riane et al., 2018).