A Globally Optimal Alternative to MLP

1. Introduction

Current research suggests that in sufficiently large neural networks, most local minima in associated loss functions have relatively small values. Consequently, the requirement of locating a global minimum becomes less crucial, as it is more significant to identify a point within the parameter space with a low loss, though not necessarily minimal [1,2,3,4]. This perspective has been substantiated through various studies on large language models. These studies demonstrate that independent training of an identical model often results in final cost function values that display minimal variation. This observation suggests that the differences between most local minima, including the global minimum, are potentially negligible. Additionally, scaling laws [5,6] establish a relationship between the number of model parameters and the (minimum) value of cost function.

Nonetheless, these regular patterns have yet to be thoroughly investigated, while they have been predominantly inferred based on experimental results. For example, in the case of scaling laws, the constants in the empirical formula can be approximated solely through experimental data. Furthermore, practitioners face significant challenges in determining whether a model’s cost function has attained the global minimum during training, as there is a lack of research providing methodologies to calculate the lower bound of a given cost function or to visualize a model’s state upon reaching the global minimum.

In this study, we begin by establishing the equivalence between MLPs and Lagrange linear interpolation functions. Subsequently, we explore the performance of MLPs in regression tasks, with a visualization of the loss function’s states at both local minima (Figure 1a) and global minima (Figure 1d). Our analysis demonstrates that conventional training methods, which depend on continuous parameter updates, are fundamentally inadequate for guiding a loss function from a local minimum to a global minimum. To address this issue, we propose a Mesh-Refinement-Coarsening (MRC) loop (Figure 1 b & c) for updating model parameters, inspired by Finite Element Method (FEM) theory [7]. This algorithm reformulates the MLP as an equivalent Lagrange linear interpolation function, and we refer to the alternative model of MLP as Lagrange Regressor (LReg).

LReg serves as an alternative to MLPs, with the ability to ensure a rapid convergence of a loss function towards its global minimum. It is applicable to a range of recognition domains, including regression, image classification, image super-resolution, text classification, and language modeling. Furthermore, our experimental results on large datasets, such as ImageNet, demonstrate that LReg, when explained as a Parameter-Efficient Fine-Tuning (PEFT) method, it promotes further convergence of the validation loss function in large pre-trained models towards global minima (see Figure 5) and enhances the model performance. In contrast to other PEFT methods, LReg has a versatile structure that is not limited to transformer-based architectures, thereby increasing its applicability. Our contributions are as follows:

• Investigating and visualizing the challenge encountered by the loss function of an MLP in overcoming to escape local minima to attain the global minimum.
• Deriving the explicit scaling law for MLPs through error bounds and conducting experimental validation on regression tasks.
• Overcoming the limitations of neural networks in fitting multi-frequency data, thereby making LReg more versatile for complex applications.
• Applying LReg as a PEFT method for large models (not confined to transformer architectures) to investigate the potential of guiding their loss functions toward global minima, consequently further reducing test loss and improving the performance of pre-trained models.

2. Methods

2.1. Preliminary: Function Fitting via MLP

Mathematically, the composition of any piecewise linear functions remains piecewise linear. Furthermore, any piecewise linear function $f(x)$ can be represented as a linear combination of Lagrange basis functions:

$$f(x)=\sum _i{c}_i{\psi }_i(x),$$

(1)

where $c_i$ is a coefficient of the basis function ${\psi }_i(x)$.

It is noteworthy that both linear layers and ReLU activations are piecewise linear functions. Consequently, an MLP composed of these elements inherently constitutes a piecewise linear mapping. For such an MLP, it is always possible to determine a set of coefficients $\{c_1,c_2,\dots \}$ and basis functions $\{{\psi }_1,{\psi }_2,\dots \}$ such that the aforementioned equation can represent the MLP.

The geometry of piecewise linear functions is straightforward: in one-dimensional space, these functions appear as segmented lines, while in two-dimensional space, they form triangulated surfaces. Figure 1a illustrates the performance of an MLP on a function-fitting task. In this example, the MLP fits the function by the dashed segmented line, while the ground truth is represented by a solid curve. The area of the shaded region between them represents the MAE. Notably, any minor alteration to any node on the segmented line will result in an increase in the MAE, because the MAE is at a local minimum in this configuration.

2.2. Optimization: Discrete vs. Continuous Processes

Upon examining Figure 1 closely, the transitioning from the piecewise linear function depicted in Figure 1a to that in Figure 1d requires two key modifications: (1) relocation of nodes ${\mathit{p}}_1$ and ${\mathit{p}}_2$, and (2) the shift of nodes ${\mathit{p}}_7$ and ${\mathit{p}}_8$ into the interval between ${\mathit{p}}_4$ and ${\mathit{p}}_5$. This is a discrete updating process as it requires reordering $({\mathit{p}}_4,{\mathit{p}}_5,{\mathit{p}}_6,{\mathit{p}}_7,{\mathit{p}}_8,{\mathit{p}}_9)$ to $({\mathit{p}}_4,{\mathit{p}}_7,{\mathit{p}}_8,{\mathit{p}}_5,{\mathit{p}}_6,{\mathit{p}}_9)$, followed by adjusting the positions of ${\mathit{p}}_7$ and ${\mathit{p}}_8$. The former adjustment might be achievable through gradient-based methods, whereas the latter is almost impossible. This difficulty arises while updating parameters continuously with first-order adaptive learning rate methods. With the learning rate inevitably diminishing to very small values, the parameter updates become negligible, reducing the likelihood of the model undergoing the significant changes required for the transition from Figure 1a to Figure 1d.

2.3. Lagrange Basis Function

Our MRC algorithm is designed for discrete updates of model parameters. However, it is not directly applicable to the training of MLPs. In order to employ it effectively, we first need to reformulate the MLP as a linear combination of basis functions, as shown in Eqn. (1). The Lagrange basis functions are defined as tent-shaped linear functions that take a value of 1 at their specified nodes and 0 at all other nodes. In Appendix A.2, we derive the formula for the Lagrange basis function in detail. For the one-dimensional case, when fitting a continuous function $F:\mathbb{R}\to \mathbb{R}$ using an MLP, the Lagrange basis function can be expressed as:

$${\psi }_i(\mathit{x})=min\left({\displaystyle \frac{\mathrm{ReLU}(\mathit{x}-{\mathit{p}}_{i-1})}{{\mathit{p}}_i-{\mathit{p}}_{i-1}}},{\displaystyle \frac{\mathrm{ReLU}({\mathit{p}}_{i+1}-\mathit{x})}{{\mathit{p}}_{i+1}-{\mathit{p}}_i}}\right).$$

For example, ${\psi }_1$ and ${\psi }_9$ are shaded blue in Figure 1a for illustrative purposes. In FEM, a mesh refers to a network of interconnected points, called nodes (e.g., ten points $\{{\mathit{p}}_0,{\mathit{p}}_1,\dots \}$ in Figure 1a), which divide a complex geometry into smaller, simpler simplices (e.g., nine intervals $\{({\mathit{p}}_0,{\mathit{p}}_1),({\mathit{p}}_1,{\mathit{p}}_2),\dots \}$ in Figure 1a).

When fitting a continuous function $F:{\mathbb{R}}^d\to \mathbb{R}$ using an MLP, the formulation for the Lagrange basis becomes substantially more intricate. This necessitates an understanding of the terminologies associated with Delaunay Triangulation. In general scenarios, the mathematical expressions of the Lagrange basis are described in Appendix A.2, and the visualizations are provided in Appendix A.5.

2.4. Mesh-Refinement-Coarsening Loop

Inspired by the FEM, we designed the MRC loop to update model parameters. The process begins with identifying underfitting regions, where the fitted value significantly deviates from the target. As shown in Figure 1b, the model exhibits MAE in the elements $({\mathit{p}}_1,{\mathit{p}}_2)$, $({\mathit{p}}_2,{\mathit{p}}_3)$, and $({\mathit{p}}_4,{\mathit{p}}_5)$; Next, we perform mesh refinement by inserting fine nodes ${\mathit{p}}_{1.5}$, ${\mathit{p}}_{2.5}$, and ${\mathit{p}}_{4.5}$ at the center of these regions (For the two-dimensional case $F:{\mathbb{R}}^2\to \mathbb{R}$, we provide a visual example of fine node insertion in Appendix A.5) and denote others as coarse nodes. The corresponding coefficients $c_{1.5}$, $c_{2.5}$, and $c_{4.5}$ are computed via interpolation to construct the scaled-up model; Subsequently, we further update these new coefficients to continuously reduce MAE using model training or the least squares method (see Appendix A.5 for details); The scaled-up model will have a reduced MAE value but with a bigger model size. For example, there are three more basis functions involved. In order to keep model size unchanged, we perform coarse node removal to scale the model down after maintaining the MAE value as much as possible. In Figure 1c, three nodes are chosen to be removed, the adjacent simplices are merged to form a scaled Lagrange basis function (see Figure 1c); We then reindex the nodes, denoting them as coarse nodes to start the next iteration of the loop. The following transformations illustrate the MRC loop in Figure 1, detailing the key steps in mesh refinement, coefficient interpolation, mesh coarsening, and reindexing.

These transformations ensure the MRC loop efficiently adapts the model structure, reducing underfitted regions while maintaining computational efficiency.

Algorithm 1:Mesh-Refinement-Coarsening Loop

while global minimum not reached do     Identify top-k underfitted regions.     Refine the mesh by adding midpoints.     Compute coefficients of new basis via interpolation.     Adjust coefficients using least squares or training.     Remove redundant coarse nodes.     Reindex all nodes. end while

Figure 1 shows that as this loop progresses, the underfitting regions gradually shrink, and the MAE consistently decreases. This process significantly deviates from first-order adaptive learning rate methods. Our methodology utilizes discrete operations to update model parameters, check Algorithm 1 for details.

2.5. Scaling Law and Error-Bound Formula

Our approach demonstrates significant interpretability, supported by the universal approximation theorem and our quantitative experimental results. In the FEM theory, the number of parameters N is proportional to the number of simplices $n_t$ in the triangular mesh (refer to Appendix A.2). Since ${(n_t/d!)}^{1/d}=O(h^{-1})$, the error bound formula for fitting function $F:{\mathbb{R}}^d\to \mathbb{R}$ is derived as:

$$L(N)={|f(\mathit{x})-F(\mathit{x})|}_{\mathit{x}\in \Omega }=O(\underset{\xi \in \Omega }{max}\parallel \nabla F(\xi )\parallel ·n_t^{-1/d})=O(N^{-1/d}),$$

where $h={max}_i{max}_j{\mathbb{1}}_{\begin{array}{c}\mathrm{dist}({\mathit{p}}^{(i)},{\mathit{p}}^{(j)})=1\end{array}}·\parallel {\mathit{p}}^{(i)}-{\mathit{p}}^{(j)}\parallel$ represents the maximum length of mesh edges (see Appendix A), and $\Omega$ is the domain over which the approximation occurs.

Our experiment shown in Figure 5 demonstrates that we perfectly predicted the explicit scaling law.

2.6. LReg Applied to Large Models

Based on the interpolation formula Eqn. (1), we propose a Lagrange Regressor (LReg) which is equivalent to MLP. LReg adheres to the universal approximation theorem and exhibits strong interpretability for low-dimensional recognition tasks (see Sec. Section 3.1). However, as the data dimension d increases, the output dimension of LReg grows factorially ($O(d!)$) as described in Sec. Section 2.5. This exponential growth makes it computationally expensive for high-dimensional data, such as large images, thereby imposing significant challenges in environments with constrained computational resources.

To address this issue, LReg can be incorporated as a module within PEFT methodologies. As shown in Figure 2, consider the target $F(\mathit{x})$ and a pre-trained model $f(\mathit{x})$. When $\tilde{f}(\mathit{x})$ is well trained, the residual $F(\mathit{x})-\tilde{f}(\mathit{x})$ remains nearly flat around the value of 0 throughout the input space. This sparsity in the residual’s support set suggests that dimensionality reduction on features could serve as an effective strategy for applying LReg to high-dimensional inputs.

In Figure 1c, we use a SigmoidNorm layer to reduce the input dimensionality:

$$\mathit{x}\leftarrow \mathrm{Sigmoid}(\mathrm{PCA}(\mathit{x})).$$

The built-in PCA module can be obtained through standard unsupervised methods or by training. LReg is then applied to learn the residual:

$$F(\mathit{x})-\tilde{f}(\mathit{x})\approx \sum _i{c}_i{\psi }_i(\mathit{x}).$$

(2)

Subsequently, the composed $\mathit{x}$ is converted to ${\left\{{\psi }_i(\mathit{x})\right\}}_i$ by the encoder of LReg. Finally, the linear combination ${\sum }_i{c}_i{\psi }_i(\mathit{x})$ is completed via the Linear layer.

In this scenario, LReg enables the loss function of large pre-trained models to move from local minima. Experiments demonstrate that LReg connection further reduces cross-entropy loss and significantly improves the performance of pre-trained models within a single training epoch (see Table 2). This approach effectively fine-tunes the model while preserving computational efficiency.

3. Results

We begin by training LReg from scratch to validate the explicit scaling law in Sec. Section 2.5 and identify the limitations of conventional models in fitting multi-frequency data (Sec. Section 3.2). The exploration of model fine-tuning for achieving reduced loss function values and enhanced performance is discussed in Sec. Section 3.3. For extending the PEFT method to non-transformer-based models via LReg connection, along with comparison results, refer to Sec. Section 3.4. Implementation details are provided in Appendix A, while additional applications are discussed in Appendix C.

3.1. Scaling Law and Addressing Regression Challenges

Figure 3 illustrates the performance of the LReg method with varying numbers of model parameters. In these experiments, LReg is trained to approach convergence on fitting the mapping function $F:{\mathbb{R}}^d\to \mathbb{R}$. As shown, the $L(N)$, denoted as test MAE loss, is observed to be proportional to the power of the number of model parameters, $N^{-1/d}$, and the inverse power of the number of simplices, $n_t^{-1/d}$, on the triangulation mesh. This observation is in perfect agreement with the derived error-bound formula in Sec. Section 2.5.

Notably, the data learned by the model in Figure 3, represented as $\left\{(x,y)\right|{\displaystyle \frac{1}{X}}\sim U(0.02,1.0),$

$Y\sim \mathcal{N}(sin{\displaystyle \frac{1}{x}},0.01)\}$, exhibits multi-frequency characteristics that present significant challenges for conventional neural networks to learn effectively. It is worth noting that although distribution $\left\{(x,y)\right|{\displaystyle \frac{1}{X}}\sim U(0.02,1.0),Y\sim \mathcal{N}(sin{\displaystyle \frac{1}{x}},0.01)\}$ and distribution $\{(x,y)|X$

$\sim U(1,50),Y\sim \mathcal{N}(sinx,0.01)\}$ are equivalent in mathematics, the numerical results of function fitting are completely different. In our open-source GitHub repository, we provide such a demo: for the neural network without LReg connection, the model fits the latter well, but cannot fit the former anyway.

Traditional regressors often struggle with overfitting, while neural networks address this but perform poorly when fitting multi-frequency functions, a phenomenon known as the frequency principle [8]. As shown in Figure 4, Support Vector Regression [9] fits dataset $\mathbb{B}=\left\{(x,y)\right|{\displaystyle \frac{1}{X}}\sim u(0.02,0.5),y=sin{\displaystyle \frac{1}{x}}\}$ but overfits on dataset $\mathbb{A}$, whereas neural networks perform well on $\mathbb{A}=\left\{(x,y)\right|X\sim u(0,{\displaystyle \frac{\pi }{4}}),Y\sim N(sin8x,|cos8x|)\}$ but underfit on $\mathbb{B}$. Neither MLPs, CNNs, nor Transformers can effectively fit sharp transitions in dataset $\mathbb{B}$, often treating these regions as noise to avoid overfitting. Importantly, no existing model can fit the challenging dataset $\mathbb{C}=\left\{(x,y)\right|{\displaystyle \frac{1}{X}}\sim u(0.02,0.5),Y\sim N(sin{\displaystyle \frac{1}{x}},0.5x^2)\}$.

LReg addresses this gap by combining the strengths of traditional regressors and neural networks. It leverages a multiscale mesh to fit high-frequency regions using fine simplices adaptively. Figure 4 (right) demonstrates the LReg’s unique ability to handle dataset $\mathbb{C}$, a task no other model can achieve. In the next section, we further illustrate its gradual fitting process during training. Comparison experiments in Table 1 show that the LReg consistently achieves high $R^2$ scores across all test sets, demonstrating its robustness and effectiveness, even in high-noise and multi-frequency scenarios.

3.2. Comparison with Regressors: Training from Scratch

As mentioned earlier, neural networks often struggle to fit multi-frequency datasets effectively. Therefore, our primary focus is comparing the LReg with traditional regressors. To evaluate the effectiveness and generalization of the LReg, we have devised four diverse datasets, each generated from distinct probability distributions:

• ${\mathbb{A}}^1$: Generated from the distribution $\left\{(x,y)\right|X\sim U(-\pi ,\pi ),Y\sim \mathcal{N}(sinx,{\displaystyle \frac{1}{5}}{cos}^{2}x)\}$, with 1000 training examples and 200 test examples. The LReg was trained with a learning rate of 0.1.
• ${\mathbb{B}}^1$: Generated from the distribution $\left\{(x,y)\right|{\displaystyle \frac{1}{X}}\sim U(0.02,1.0),Y\sim \mathcal{N}(sin{\displaystyle \frac{1}{x}},0.01)\}$, with 1000 training examples and 200 test examples. We trained the LReg with a learning rate of 0.9.
• ${\mathbb{A}}^2$: Generated from the distribution $\left\{(\mathit{x},y)\right|{\mathit{x}}_i\sim U(-\pi ,\pi ),Y_i\sim \mathcal{N}(sin{\mathit{x}}_i,{\displaystyle \frac{1}{10}}{cos}^2{\mathit{x}}_i),$
  $Y={\displaystyle \frac{1}{2}}(Y_1+Y_2)\}$, with 7,500 training examples and 1,500 test examples. The LReg was trained with a learning rate of 0.1.
• ${\mathbb{B}}^2$: Generated from the distribution $\left\{(\mathit{x},y)\right|{\displaystyle \frac{1}{{\mathit{x}}_i}}\sim U(0.05,0.5),Y_i\sim \mathcal{N}(sin{\displaystyle \frac{1}{{\mathit{x}}_i}},0.01),$
  $Y={\displaystyle \frac{1}{2}}(Y_1+Y_2)\}$, with 50,000 training examples and 10,000 test examples. Training utilized a learning rate of 0.9.

Table 1 displays the coefficient of determination ($R^2$) scores for the LReg and traditional regressors [9,10,11,12,13,14,15,16,17,18] across fitting the four datasets. The LReg consistently achieves high $R^2$ scores across all test sets, demonstrating the effectiveness of the InterpolationNet on both high-noise and multi-frequency datasets. Furthermore, the minimal gap between training and test set evaluations underscores the robustness of the LReg, indicating its capability of generalization.

We also explore training LReg from scratch on non-regression tasks. Here, we delve into the practical application of LReg for text feature extraction using the AG News dataset [19] for classification tasks. Our experimental findings reveal that LReg can directly extract features from raw text. However, recognizing that tokens are unordered categorical variables, to enhance performance, we add a preprocessing layer that converts each token to a four-dimensional vector, the proportion of the token appearing in four categories.

Our experiment showcases the LReg achieving 90.01% test accuracy after the first epoch, with 90.4% test accuracy reached within just five epochs. It possesses a unique characteristic in text classification tasks: the number of its parameters remains independent of the token count. Compared to word2vec-based networks, our model performs equally well in classification but boasts only 256 model parameters, a significant reduction from the word2vec-based network with more than 6.13 million parameters.

3.3. Improving pre-trained Model Performance and Continuously Reducing Loss

As described in Sec. Section 2.6, LReg can serve as a connection to enable the loss function of large pre-trained models to escape local minima (see Figure 5). In this module, the SigmoidNorm layer and linear head contain trainable parameters (see Figure 2c). In this section, we illustrate the effectiveness of LReg connection through ablation studies. To ensure a fair comparison with baseline models, we adopt a stricter experimental setup: the original training recipe of the pre-trained model remains unchanged, the entire pre-trained model is frozen, and no additional training data is used.

This restriction is crucial to ablation studies because the training recipes for pre-trained models often have room for optimization [20]. Changes to batch size, optimizer, weight decay rate, or the order of applying data augmentations could potentially improve the performance of the pre-trained model [21,22]. To avoid such effects, we selected pre-trained models from TorchVision with publicly available training recipes as our baselines. The weights of these pre-trained models closely reproduce the results from the original papers on the ImageNet-1K dataset [23], with recipes available at [24].

Our method offers significant practical value. Unlike model scaling and adaptation approaches, it requires no modifications to the raw training recipe, converges within a single epoch, and is applicable to any model, not limited to the transformer-based architecture. For instance, on 4x A6000 GPUs, it adds fewer than one million parameters to a ResNet-50 model and achieves a 0.2% accuracy gain in just 40 minutes of training. Due to space constraints, we present some of the comparison results in Table 2.

In Appendix B, we highlight the strong stability of LReg, showing ±0.07% fluctuation in validation accuracy across independent model training (see Appendix B, Table A1). Trainable parameters in our method originate from the SigmoidNorm and linear layer (see Figure 2c). The total number of trainable parameters can be calculated as $N=(C_{\mathrm{in}}+1)\times d+d\times n\times C$, where N is the number of model parameters, $C_{\mathrm{in}}$ is the SigmoidNorm input dimension, d the PCA output dimension, n is the mesh node count (degrees of freedom), and C the total output dimension.

Increasing n and d improves model performance (see Appendix B, Table A2). We also study empirical scaling laws for model performance on cross-entropy loss to demonstrate how our method further the loss value of pre-trained models. [5] proposed an empirical formula where the loss scales as a power-law with model size:

$$L(N)={\left({\displaystyle \frac{c}{N}}\right)}^{\alpha }$$

where $L(N)$ is the cross-entropy loss on the validation dataset, c and $\alpha$ are constants related to the model type. Figure 5 shows that our method is more effective than simply scaling the model size. Our method increases the $\alpha$ for MobileNet [25] from 0.27 to 0.47, for ResNet [26] from 0.25 to 2.1, for ResNeXt [27] from 0.01 to 3.1, and for ViT [28] from -0.06 to 1.2.

3.4. Transfer Learning via LReg Connection

Through experiments, we identified three key limitations of existing PEFT methods: 1) Dependence on Transfer Learning: Methods like LoRA [29] struggle to outperform pre-trained models on the same or similar datasets, when no domain adaptation is needed. This means they cannot continuously reduce the loss function on the raw training set. Additional training data is often required for these methods to be effective. Therefore, their strength lies mainly in transfer learning scenarios (see Appendix B); 2) Sensitive Training Requirements: PEFT methods require specific training recipes, such as small learning rates. Without these, performance often deteriorates from the first epoch; 3) Task and Architecture Limits: PEFT methods, like LoRA, LoHA [30] and IA3 [31], are restricted to transformer-based architectures.

We compare fine-tuning (FT), BitFit [32], LoRA, and LReg connection on the GLUE dataset [33]; subsequently, we conduct additional experiments on WikiText benchmarks [34] using GPT2 [35] for causal language modeling and Roberta [36] for masked language modeling with a default random seed (42). Pre-trained models of the HuggingFace Transformers library [37] are used. A reduced learning rate ($\le {10}^{-5}$) is applied for LoRA to prevent performance degradation in these experiments. As shown in Table 3, our method overcomes the mentioned challenges, performing well even without additional training data. Therefore, it is also versatile for transfer learning.

Table 3. GPT2 and RoBERTa with Different Adaptation Methods on the GLUE dataset and WikiText Benchmark. We report the overall (both matched and mismatched) accuracy for MNLI and the accuracy for other tasks. For each model, we conducted multiple experiments and show the accuracy as “${\mathrm{mean}}_{\pm \mathrm{std}}$”. For a fair comparison, we adjusted the number of parameters of our method to be similar to those of LoRA. * indicates numbers published in prior works. We also report the number of trainable parameters, perplexity (PPL, with lower values being preferable), and training throughput (measured in sequences per second) for language modeling tasks.

MNLI

SST-2

MRPC

QNLI

QQP

RTE

Method

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

FT*

125M

87.6

125M

94.8

125M

90.2

125M

92.8

125M

91.9

125M

78.7

BitFit*

0.1M

84.7

0.1M

93.7

0.1M

92.7

0.1M

91.8

0.1M

84.0

0.1M

81.5

LoRA*

0.295M

$87.5_{\pm .3}$

0.295M

$95.1_{\pm .2}$

0.295M

$89.7_{\pm .7}$

0.295M

$93.3_{\pm .3}$

0.295M

$90.8_{\pm .1}$

0.295M

$86.6_{\pm .7}$

LReg

0.259M

$87.9_{\pm .5}$

0.259M

$95.3_{\pm .6}$

0.259M

$90.5_{\pm .4}$

0.259M

$93.5_{\pm .2}$

0.259M

$89.7_{\pm .4}$

0.259M

$83.5_{\pm .3}$

	GPT2			GPT2-Medium			GPT2-Large			Roberta-Base
Method	#Params	PPL	Seq/s	#Params	PPL	Seq/s	#Params	PPL	Seq/s	#Params	PPL	Seq/s
FT	124.4M	21.4234	6.98	354.8M	15.8900	3.50	774.0M	13.8468	1.58	124.7M	3.6415	10.05
LoRA	0.295M	21.4188	12.90	0.393M	15.8873	6.29	0.737M	13.8486	1.67	0.295M	3.6425	9.85
LReg	0.203M	21.3401	9.58	0.404M	15.8674	6.66	0.607M	13.8443	2.52	0.303M	3.6414	12.37

Table 3. GPT2 and RoBERTa with Different Adaptation Methods on the GLUE dataset and WikiText Benchmark. We report the overall (both matched and mismatched) accuracy for MNLI and the accuracy for other tasks. For each model, we conducted multiple experiments and show the accuracy as “${\mathrm{mean}}_{\pm \mathrm{std}}$”. For a fair comparison, we adjusted the number of parameters of our method to be similar to those of LoRA. * indicates numbers published in prior works. We also report the number of trainable parameters, perplexity (PPL, with lower values being preferable), and training throughput (measured in sequences per second) for language modeling tasks.

MNLI

SST-2

MRPC

QNLI

QQP

RTE

Method

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

#Params

Accuracy

FT*

125M

87.6

125M

94.8

125M

90.2

125M

92.8

125M

91.9

125M

78.7

BitFit*

0.1M

84.7

0.1M

93.7

0.1M

92.7

0.1M

91.8

0.1M

84.0

0.1M

81.5

LoRA*

0.295M

$87.5_{\pm .3}$

0.295M

$95.1_{\pm .2}$

0.295M

$89.7_{\pm .7}$

0.295M

$93.3_{\pm .3}$

0.295M

$90.8_{\pm .1}$

0.295M

$86.6_{\pm .7}$

LReg

0.259M

$87.9_{\pm .5}$

0.259M

$95.3_{\pm .6}$

0.259M

$90.5_{\pm .4}$

0.259M

$93.5_{\pm .2}$

0.259M

$89.7_{\pm .4}$

0.259M

$83.5_{\pm .3}$

	GPT2			GPT2-Medium			GPT2-Large			Roberta-Base
Method	#Params	PPL	Seq/s	#Params	PPL	Seq/s	#Params	PPL	Seq/s	#Params	PPL	Seq/s
FT	124.4M	21.4234	6.98	354.8M	15.8900	3.50	774.0M	13.8468	1.58	124.7M	3.6415	10.05
LoRA	0.295M	21.4188	12.90	0.393M	15.8873	6.29	0.737M	13.8486	1.67	0.295M	3.6425	9.85
LReg	0.203M	21.3401	9.58	0.404M	15.8674	6.66	0.607M	13.8443	2.52	0.303M	3.6414	12.37

4. Discussion

It is important to acknowledge the limitations associated with LReg. As described in Sec. Section 2.6, the substantial computational cost of extracting features from high-dimensional data restricts its direct application to large-scale datasets. Presently, we address this by incorporating LReg as an adaptation method, but we aim to extract features directly from large-scale data in the future. Additionally, the experiments in Sec. Section 3.3 indicate that, although our adaptation method significantly reduces cross-entropy on the ImageNet-1K validation set, the extent of accuracy improvement is relatively modest, warranting further investigation. We didn’t compare LReg with models such as Support Vector Machine (SVM) known for their robust interpretability, due to fundamentally different underlying principles, and the usual missing experimental results from interpretable models on large datasets.

5. Conclusions

In this paper, we reformulate MLPs and introduce a mathematically equivalent framework named LReg. Furthermore, we propose a discrete method, referred to as the MRC loop, which serves as an alternative to traditional continuous optimization techniques, thereby enabling the loss function to reach global minima. For general pre-trained models, we demonstrate that our approach can significantly reduce the loss function value.

From a theoretical perspective, we construct a visual example to elucidate the intrinsic challenge encountered by traditional optimization methods in reducing the loss function from local minima to global minima, and highlight the benefits of discrete parameter updates in addressing this challenge. Furthermore, we analyze the error-bound formula and, for the first time, derive an explicit scaling law, while we empirically validate through regression tasks. From an application perspective, LReg effectively mitigates the core limitation of traditional neural networks in learning multi-frequency data. Additionally, LReg serves as a plug-and-play module, and our analysis of the scaling law shows that it substantially reduces the loss function value of pre-trained models while providing moderate improvement in performance. In contrast to scaling up a model and training it from scratch, our approach is more efficient. Furthermore, LReg broadens the application scope of PEFT methods in transfer learning, without requiring the pre-trained model to be transformer-based.

This work opens new directions for research in model training and establishes a foundation for further advancements in discrete optimization.

Appendix A.1. Triangulated Mesh

In the context of FEM, elements often serve as the fundamental building blocks of the triangulation mesh, taking the form of simplices created by connecting nodes. For instance, in 1D FEM, simplices are intervals (see Figure 1a); in 2D FEM, triangles with three nodes are commonly used (see Figure A1, right), while 3D FEM often employs tetrahedra with four nodes. This concept is visually depicted in Figure A1 (left), where the mesh consists of eight nodes and seven triangles. This type of mesh is established by specifying the coordinates of discrete nodes and the vertex indices of simplices. Let d represent the dimensions, ${\left\{{\mathit{p}}^{(i)}\right\}}_{i=0}^{n-1}$ denote the grid nodes, and introduce a matrix $\mathit{p}$ to store the node coordinates:

$${\mathit{p}}_{i,j}={\mathit{p}}_j^{(i)}.$$

Additionally, utilize a matrix $\mathit{T}$ to store the indices of nodes constituting the simplices within the triangulation. Specifically, access the j-th sorted vertex of the i-th simplex in this mesh as ${\mathit{p}}_{{\mathit{T}}_{i,j},:}$. Figure A1 (left) illustrate a matrix $\mathit{T}$ takes the following form:

$$\mathit{T}={\left[\begin{array}{ccccccc}6& 0& 3& 2& 7& 3& 4\\ 7& 6& 7& 3& 4& 4& 3\\ 5& 5& 1& 1& 5& 7& 2\end{array}\right]}^T.$$

This matrix serves to describe all seven simplices within the mesh, such as the first simplex $△{\mathit{p}}^{(6)}{\mathit{p}}^{(7)}{\mathit{p}}^{(5)}$ and the last simplex $△{\mathit{p}}^{(4)}{\mathit{p}}^{(3)}{\mathit{p}}^{(2)}$.

Figure A1. Left - Mesh with eight nodes and seven triangles. Right - Contours of eight Lagrange basis functions, linear variation of ${\psi }_i$ associated with node ${\mathit{p}}^{(i)}$ across all triangles.

Figure A1. Left - Mesh with eight nodes and seven triangles. Right - Contours of eight Lagrange basis functions, linear variation of ${\psi }_i$ associated with node ${\mathit{p}}^{(i)}$ across all triangles.

The first-order Lagrange basis studied in this article, denoted as $\{{\psi }_0(\mathit{x}),\dots ,{\psi }_{n-1}(\mathit{x})\}\subset P_1({\mathbb{R}}^d)$, are piecewise linear polynomials associated with nodes $\{{\mathit{p}}^{(0)},\dots ,{\mathit{p}}^{(n-1)}\}$. These functions are defined such that ${\psi }_i({\mathit{p}}^{(j)})={\mathbb{1}}_{\mathrm{i}}=\mathrm{j}$. Figure A1 (right) illustrates this: ${\psi }_i(\mathit{x})$ corresponds to node ${\mathit{p}}^{(i)}$, exhibiting linear variation across all elements. Its support encompasses the union of all neighboring elements of node ${\mathit{p}}^{(i)}$ (refer to Appendix B for a 3-dimensional visualization). For example, $\mathrm{supp}({\psi }_3)=△{\mathit{p}}^{(3)}{\mathit{p}}^{(4)}{\mathit{p}}^{(7)}\cup △{\mathit{p}}^{(3)}{\mathit{p}}^{(7)}{\mathit{p}}^{(1)}\cup △{\mathit{p}}^{(3)}{\mathit{p}}^{(1)}{\mathit{p}}^{(2)}\cup △{\mathit{p}}^{(3)}{\mathit{p}}^{(2)}{\mathit{p}}^{(4)}$.

Now, we formulate the Lagrange basis from its original definition to establish the foundational architecture of LReg. It is important to highlight that the traditional Lagrange basis involves unbalanced computing of barycentric coordinates, which may not be well-suited for parallel deep learning platforms (see Appendix C for details on the traditional definition of the Lagrange basis). Consequently, in this subsection, we re-derive the Lagrange basis to enhance parallel computing.

Appendix A.2. Lagrange Basis Function

Let $n_t$ represent the number of simplices in the multiscale mesh. We introduce the Parameters Tensor $\mathbf{S}$ defined as:

$${\mathbf{S}}_{j,:,:}={\left[\begin{array}{cccc}{\mathit{p}}_0^{({\mathit{T}}_{j,0})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,0})}& 1\\ ⋮& \ddots & ⋮& ⋮\\ {\mathit{p}}_0^{({\mathit{T}}_{j,d-1})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,d-1})}& 1\\ {\mathit{p}}_0^{({\mathit{T}}_{j,d})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,d})}& 1\end{array}\right]}^{-1},j=0,\dots ,n_t-1.$$

Additionally, we introduce the Node Membership tensor $\mathbf{M}$ defined as:

$${\mathbf{M}}_{i,j,k}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\mathrm{the}i-\mathrm{th}\mathrm{node}\mathrm{matches}\mathrm{the}k-\mathrm{th}\mathrm{vertex}\mathrm{of}\mathrm{the}j-\mathrm{th}\mathrm{simplex},\hfill \\ 0,\hfill & \mathrm{other}\mathrm{cases}.\hfill \end{array}\right.$$

By defining:

$${\mathit{U}}_{j,k}(\mathit{x})=\sum _{\tau =0}^{d-1}{\mathbf{S}}_{j,\tau ,k}·{\mathit{x}}_{\tau }+{\mathbf{S}}_{j,d,k},j=0,\dots ,n_t-1,k=0,\dots ,d.$$

We will demonstrate in Appendix A that the following function qualifies the definition of Lagrange basis:

$${\psi }_i(\mathit{x})={\displaystyle \frac{{\sum }_{j=0}^{n_t-1}{\sum }_{k=0}^d{\mathbb{1}}_{{min}_{\tau }{\mathit{U}}_{j,\tau }(\mathit{x})\ge 0}·{\mathbf{M}}_{i,j,k}·{\mathit{U}}_{j,k}(\mathit{x})}{max({\sum }_{j=0}^{n_t-1}{\sum }_{k=0}^d{\mathbb{1}}_{{min}_{\tau }{\mathit{U}}_{j,\tau }(\mathit{x})\ge 0}·{\mathbf{M}}_{i,j,k},1)}},i=0,\dots ,n-1.$$

(A1)

So far, we have successfully constructed :

$$\begin{array}{cc}\hfill \mathrm{Encoder}:& {\mathbb{R}}^d\to {[0,1]}^n,\hfill \\ \hfill & \mathit{x}\mapsto ({\psi }_0(\mathit{x}),\dots ,{\psi }_{n-1}(\mathit{x})).\hfill \end{array}$$

The above basis is in the format of $P_1({\mathbb{R}}^d)$, which is very useful for low-dimensional regression tasks. In the one-dimensional case, Eq. (A1) simplifies to the Lagrange basis $P_1({\mathbb{R}}^1)$:

$${\psi }_i(\mathit{x})=min\left({\displaystyle \frac{\mathrm{ReLU}(\mathit{x}-{\mathit{p}}_{i-1})}{{\mathit{p}}_i-{\mathit{p}}_{i-1}}},{\displaystyle \frac{\mathrm{ReLU}({\mathit{p}}_{i+1}-\mathit{x})}{{\mathit{p}}_{i+1}-{\mathit{p}}_i}}\right),i=0,1,\dots ,n-1.$$

The $P_1({\mathbb{R}}^1)$ basis has exceptionally low computational complexity, making its direct sum well-suited for recognition tasks on large-scale datasets. In Sec. Section 2.6, we introduce the specific approach in detail. We use the $P_1({\mathbb{R}}^1)$ basis to implement an adaptive method for learning a decoupled residual system. This method can enhance the performance of pre-trained models on large datasets.

Appendix A.3. Proof of Lagrange Basis Expression

We will now demonstrate, in three concise steps, that Eqn. (A1) qualifies as a Lagrange basis function.

Piecewise Linear: Since $\{{\mathit{U}}_{j,0},\dots ,{\mathit{U}}_{j,d}\}$ are piecewise linear functions, their linear combination $L_i$ is also piecewise linear.

Kronecker Delta: From the definition of $\mathit{U}$ and $\mathbf{S}$, we have the following equation:

$$\left[\begin{array}{cccc}{\mathit{x}}_0& \dots & {\mathit{x}}_{d-1}& 1\end{array}\right]=\left[\begin{array}{ccc}{\mathit{U}}_{j,0}(\mathit{x})& \dots & {\mathit{U}}_{j,d}(\mathit{x})\end{array}\right]\left[\begin{array}{cccc}{\mathit{p}}_0^{({\mathit{T}}_{j,0})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,0})}& 1\\ ⋮& \ddots & ⋮& ⋮\\ {\mathit{p}}_0^{({\mathit{T}}_{j,d-1})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,d-1})}& 1\\ {\mathit{p}}_0^{({\mathit{T}}_{j,d})}& \dots & {\mathit{p}}_{d-1}^{({\mathit{T}}_{j,d})}& 1\end{array}\right].$$

Decomposing this equation, we obtain $\mathit{x}={\sum }_{k=0}^d{\mathit{U}}_{j,k}(\mathit{x}){\mathit{p}}^{({\mathit{T}}_{j,k})}$ and ${\sum }_{k=0}^d{\mathit{U}}_{j,k}(\mathit{x})=1$. This implies two important conclusions: ${\mathit{U}}_{j,k^{\prime }}({\mathit{p}}^{({\mathit{T}}_{j,k^{\prime \prime }})})={\mathbb{1}}_{k^{\prime }=k^{\prime \prime }}$ and ${min}_{\tau }{\mathit{U}}_{j,\tau }(\mathit{x})\ge 0$ is true if and only if $\mathit{x}$ belongs to the j-th simplex. Therefore, we have

$$\left\{\begin{array}{cc}{\mathbb{1}}_{{min}_{\tau }{\mathit{U}}_{j,\tau }({\mathit{p}}^{(i)})\ge 0}={\mathbf{M}}_{i,j,k}={\mathit{U}}_{j,k}({\mathit{p}}^{(i)})=1,\hfill & \mathrm{if}i={\mathit{T}}_{j,k}\hfill \\ {\mathbf{M}}_{i,j,k}={\mathit{U}}_{j,k}({\mathit{p}}^{(i)})=0,\hfill & \mathrm{if}i\ne {\mathit{T}}_{j,k}\hfill \end{array}\right.$$

This proves that ${\psi }_i({\mathit{p}}^{(j)})={\mathbb{1}}_{i=j}$.

Globally Continuity: Lastly, since ${\psi }_i$ is inherently linear within all simplices and exhibits continuity across all grid nodes, we can conclude that ${\psi }_i$ is globally continuous.

Appendix A.4. Validation of Lagrange Basis Functions’ Mathematical Expression

This section focuses on the mathematical validation of Lagrange basis function expressions to ensure they satisfy our definition in Sec. Section 2.3, the Lagrange basis functions are defined as tent-shaped linear functions that take a value of 1 at their specified nodes and 0 at all other nodes. We start with a two-dimensional case and a triangulated mesh. Let triangle △ be defined by nodes $\{{\mathit{p}}^{(i)},{\mathit{p}}^{(j)},{\mathit{p}}^{(k)}\}$. The following barycentric coordinates $\{{\lambda }_{△,i},{\lambda }_{△,j},{\lambda }_{△,k}\}$ are three first-degree polynomials of $\mathit{x}$

$$\left[\begin{array}{c}{\lambda }_{△,i}(\mathit{x})\\ {\lambda }_{△,j}(\mathit{x})\\ {\lambda }_{△,k}(\mathit{x})\end{array}\right]={\left[\begin{array}{ccc}{\mathit{p}}_0^{(i)}& {\mathit{p}}_0^{(j)}& {\mathit{p}}_0^{(k)}\\ {\mathit{p}}_1^{(i)}& {\mathit{p}}_1^{(j)}& {\mathit{p}}_1^{(k)}\\ 1& 1& 1\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathit{x}}_0\\ {\mathit{x}}_1\\ 1\end{array}\right].$$

Referring to the instance depicted in Figure A1 (left), the mesh consists of eight nodes and seven triangles. Specifically, let $\{{△}^{(j)}=△{\mathit{p}}^{({\mathit{T}}_{j,0})}{\mathit{p}}^{({\mathit{T}}_{j,1})}{\mathit{p}}^{({\mathit{T}}_{j,2})}|j=0,\dots ,6\}$. We will now verify that the following ${\psi }_3$ corresponds to the third basis function in this mesh

$${\psi }_3(\mathit{x})={\displaystyle \frac{1}{max({\sum }_{i\in \{2,3,6,5\}}{\mathbb{1}}_{\mathit{x}\in {△}^{(i)}},1)}}\sum _{i\in \{2,3,6,5\}}{\mathbb{1}}_{\mathit{x}\in {△}^{(i)}}{\lambda }_{{△}^{(i)},3}(\mathit{x}).$$

First, ${\psi }_3$ possesses values of Kronecker Delta:

$${\psi }_3(\mathit{x})=\left\{\begin{array}{c}\begin{array}{ccccccccc}\hfill 35{\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(0)},\\ \hfill {\displaystyle \frac{1}{max(2,1)}}(& 1·0\hfill & \hfill +& 1·0\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(1)},\\ \hfill {\displaystyle \frac{1}{max(2,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 1·0\hfill & \hfill +& 1·0\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(2)},\\ \hfill {\displaystyle \frac{1}{max(4,0)}}(& 1·1\hfill & \hfill +& 1·1\hfill & \hfill +& 1·1\hfill & \hfill +& 1·1)\hfill & \hfill =1,\mathrm{if}\mathit{x}={\mathit{p}}^{(3)},\\ \hfill {\displaystyle \frac{1}{max(2,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 1·0\hfill & \hfill +& 1·0)\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(4)},\\ \hfill {\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(5)},\\ \hfill {\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(6)},\\ \hfill {\displaystyle \frac{1}{max(2,1)}}(& 1·0\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 1·0)\hfill & \hfill =0,\mathrm{if}\mathit{x}={\mathit{p}}^{(7)},\end{array}\hfill \end{array}\right.$$

Then, ${\psi }_3$ is a first-degree polynomial in every triangle and $\mathrm{supp}{\psi }_3=\mathrm{inn}{\cup }_{i\in \{2,3,6,5\}}{△}^{(i)}$:

$${\psi }_3(\mathit{x})=\left\{\begin{array}{c}\begin{array}{ccccccccccc}\hfill 36{\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& 0,\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_0,\\ \hfill {\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& 0,\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_1,\\ \hfill {\displaystyle \frac{1}{max(1,1)}}(& 1·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& {\lambda }_{{△}^{(2)},3}(\mathit{x}),\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_2,\\ \hfill {\displaystyle \frac{1}{max(1,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 1·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& {\lambda }_{{△}^{(3)},3}(\mathit{x}),\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_3,\\ \hfill {\displaystyle \frac{1}{max(0,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& 0,\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_4,\\ \hfill {\displaystyle \frac{1}{max(1,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 1·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& {\lambda }_{{△}^{(5)},3}(\mathit{x}),\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_5,\\ \hfill {\displaystyle \frac{1}{max(1,1)}}(& 0·{\lambda }_{{△}^{(2)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(3)},3}(\mathit{x})\hfill & \hfill +& 1·{\lambda }_{{△}^{(6)},3}(\mathit{x})\hfill & \hfill +& 0·{\lambda }_{{△}^{(5)},3}(\mathit{x}))\hfill & \hfill =& {\lambda }_{{△}^{(6)},3}(\mathit{x}),\hfill & \mathrm{if}\mathit{x}\in \mathrm{inn}{\mathit{T}}_6,\end{array}\hfill \end{array}\right.$$

Finally, since ${\psi }_3$ is continuous in all nodes and first-degree in all triangles, it is globally continuous.

Appendix A.5. Visualization of Lagrangian basis

In Sec. Section 2.3 and Appendix A.2, we introduced first-order Lagrange basis functions, a set of piecewise linear functions defined on a mesh. Each basis function corresponds to a node.

Consider the grid depicted in Figure A2 (left). Taking the inserted fine node ${\mathit{p}}^{(20)}$ as an example, it has a total of four neighboring nodes: ${\mathit{p}}^{(5)}$, ${\mathit{p}}^{(0)}$, ${\mathit{p}}^{(7)}$, and ${\mathit{p}}^{(4)}$. By connecting these nodes, we can determine the support of the basis function ${\psi }_{20}$. In Figure A2 (middle), we present the function graphs of ${\psi }_{20}$ and ${\psi }_7$. It can be observed that these functions exhibit linear variations on each mesh triangle. Taking ${\psi }_{20}$ as an example, its function value at ${\mathit{p}}_{20}$ is 1, and 0 at all other nodes. Similarly, ${\psi }_7$ has a function value of 1 at ${\mathit{p}}_7$ and 0 at other nodes. In Figure A2 (right), the orange triangles represent the function graph of $f(\mathit{x})$ on the domain $△{\mathit{p}}^{(0)}{\mathit{p}}^{(7)}{\mathit{p}}^{(9)}$, the green dots represent a subset of the training set ${\left\{({\mathit{x}}^{(i)},y^{(i)})\right\}}_i$, where the projections (raw data) fall on $△{\mathit{p}}^{(0)}{\mathit{p}}^{(7)}{\mathit{p}}^{(9)}$. Base on the definition of Lagrange basis function, we have:

$${f(\mathit{x})|}_{\mathit{x}\in △{\mathit{p}}^{(0)}{\mathit{p}}^{(7)}{\mathit{p}}^{(9)}}=\sum _i{c}_i{\psi }_i{(\mathit{x})|}_{\mathit{x}\in △{\mathit{p}}^{(0)}{\mathit{p}}^{(7)}{\mathit{p}}^{(9)}}=c_0{\psi }_0(\mathit{x})+c_7{\psi }_7(\mathit{x})+c_9{\psi }_9(\mathit{x}).$$

The training process mentioned in Sec. Section 2.4 uses the least squares method to adjust coefficients, minimizing the mean gap (MAE) between the green points and the triangular surface.

$$\left\{\begin{array}{c}\underset{(c_0,c_7,c_9)}{argmin}\sum _i|c_0{\psi }_0({\mathit{x}}^{(i)})+c_7{\psi }_0({\mathit{x}}^{(i)})+c_9{\psi }_0({\mathit{x}}^{(i)})-y^{(i)}|,\mathrm{Using}\mathrm{model}\mathrm{training}\mathrm{method};\hfill \\ \underset{(c_0,c_7,c_9)}{argmin}\sum _i{|c_0{\psi }_0({\mathit{x}}^{(i)})+c_7{\psi }_0({\mathit{x}}^{(i)})+c_9{\psi }_0({\mathit{x}}^{(i)})-y^{(i)}|}^2,\mathrm{Using}\mathrm{least}\mathrm{squares}\mathrm{method}.\hfill \end{array}\right.$$

Figure A2. Data visualization. Left - an example of refining a two-dimensional mesh. Middle - the graphs of basis functions ${\psi }_{20}$ and ${\psi }_7$. Right - the graphs of the function $f(\mathit{x})$ and a subset of the training set, adjusting coefficients for minimizing MAE Loss.

Figure A2. Data visualization. Left - an example of refining a two-dimensional mesh. Middle - the graphs of basis functions ${\psi }_{20}$ and ${\psi }_7$. Right - the graphs of the function $f(\mathit{x})$ and a subset of the training set, adjusting coefficients for minimizing MAE Loss.

Appendix C.1. Solve PDEs

In this section, we utilize the LReg to address the following partial differential equations (PDEs):

$$\left\{\begin{array}{cc}\Delta u+{(u-\beta )}^2={(\alpha cosxsiny-1)}^2+1,\hfill & (x,y)\in \Omega ;\hfill \\ u=\beta ,\hfill & (x,y)\in \partial \Omega .\hfill \end{array}\right.$$

Here, $\Omega =[0,1]\times [0,1]$. We construct a dataset that takes $(\alpha ,\beta )$ as input data and assigns the corresponding numerical solution of the PDEs as the target output. This dataset comprises 12,000 examples, with $\alpha$ randomly selected from the distribution $U(-\pi /2,\pi /2)$ and $\beta$ randomly chosen from the distribution $U(0,2\pi )$. We then split the dataset into two parts: 10,000 for training and 2,000 for testing. Figure A3 illustrates how well the network predicts the exact solution.

Figure A3. Residual - The gap between the exact solution and the model output.

Figure A3. Residual - The gap between the exact solution and the model output.

Appendix C.2. fitting high-noise data

In this section, we conduct the LReg on learning the dataset $\mathbb{A}=\left\{(x,y)\right|X\sim U(-4,4),Y\sim \mathcal{N}(sinx,0.2{cos}^{2}x)\}$. A comprises 6000 examples, with 5000 for training and 1000 for testing. Figure A4 shows the training progress.

Figure A4. Blue dots represent the training set, while the orange curve represents the network.

Figure A4. Blue dots represent the training set, while the orange curve represents the network.

Appendix C.3. fitting multi-frequency data

In this section, we conduct the LReg on learning the dataset $\mathbb{A}=\left\{(x,y)\right|{\displaystyle \frac{1}{x}}\sim U(0.02,0.5),y=sin{\displaystyle \frac{1}{x}}\}$. A comprises 6000 examples, with 5000 for training and 1000 for testing. Figure A5 illustrates the training progress. Remarkably, after just 4 epochs of training, the neural network outputs closely approximate the target values. By the 32nd epoch’s conclusion, the neural network outputs and target values are nearly indistinguishable.

Figure A5. Blue dots represent the training set, while the orange curve represents the network.

Figure A5. Blue dots represent the training set, while the orange curve represents the network.

Appendix C.4. Fit a Vector-Valued Function

In this instance, we utilize the LReg to fit spherical harmonics. Our dataset denoted as $\mathbb{A}=\left\{(\mathit{x},\mathit{y})\right|\mathit{x}=(\theta ,\varphi ),\mathit{y}=(\mathrm{Real}(Y_4^2(\theta ,\varphi )),\mathrm{Imag}(Y_4^2(\theta ,\varphi ))),\Theta \sim U(0,2\pi ),\Phi \sim U(0,\pi )\}$ , comprises 48,000 examples, with 40,000 allocated for training and an additional 8,000 for testing. Figure A6 shows the training progress.

Figure A6. In each block, the left panel represents the real part of our model output, while the right panel represents the imaginary part of the model output.

Figure A6. In each block, the left panel represents the real part of our model output, while the right panel represents the imaginary part of the model output.