Probability Distribution Learning Framework and Its Application in Deep Learning

We begin by offering a concise summary of the classical mathematical and statistical theory underlying machine learning tasks and algorithms, which, in their most general form, can be formulated as follows [22].

To analyze the influencing factors of learning, the classical statistical machine learning framework introduces the optimization error ($e_{\text{opt}}$), generalization error ($e_{\text{gen}}$), and approximation error ($e_{\text{approx}}$).

The optimization error $e_{\rm opt}$ is primarily influenced by the algorithm $\mathcal{A}$ that is used to find the model $f_{s}$ in the hypothesis set $\mathcal{F}$ for given training data $s$. By far, the gradient-based optimization methods are the most prevalent methods to address this issue, leveraging their ability to precisely and efficiently compute the derivatives by automatic differentiation in a lot of widely used models. Owing to the hierarchical structure of neural networks (NNs), the back-propagation algorithm has been proposed to efficiently calculate the gradient with respect to the loss function [23, 24, 25, 26]. Stochastic gradient descent (SGD) [27] is a typical gradient-based optimization method used to identify the unique minimum of a convex function [28]. In scenarios where the objective function is non-convex, SGD generally converges to a stationary point, which may be either a local minimum, a local maximum, or a saddle point. Nonetheless, there are indeed conditions under which convergence to local minima is guaranteed [29, 30, 31, 32].

The approximation error, $e_{\text{approx}}$, describes the model’s fitting capability. The universal approximation theorem [33, 34, 35, 1] states that neural networks with only 2 layers can approximate any continuous function defined on a compact set up to arbitrary precision. This implies that the $e_{\rm approx}$, decreases as the model complexity increases.

The generalization error, $e_{\rm gen}(\mathcal{F})$, is defined as the difference between the true risk and the empirical risk. To bound the generalization error, Hoeffding’s inequality [36] alongside the union bound directly provides the PAC (Probably Approximately Correct) bound [37]. The uniform convergence bound represents a classic form of PAC generalization bound, which relies on metrics of complexity for the set $\mathcal{F}$, including the VC-dimension [4] and the Rademacher complexity [5]. Such bounds are applicable to all the hypotheses within $\mathcal{F}$. A potential drawback of the uniform convergence bounds is that they are independent of the learning algorithm, i.e., they do not take into account the way the hypothesis space is explored. To tackle this issue, algorithm-dependent bounds have been proposed to take advantage of some particularities of the learning algorithm, such as its uniform stability [6] or robustness [38].

However, the available body of evidence suggests that classical tools from statistical learning theory alone, such as uniform convergence, algorithmic stability, or bias-variance trade-off may be unable to fully capture the behavior of DNNs [12, 39]. The classical trade-off highlights the need to balance model complexity and generalization. While simpler models may suffer from high bias, overly complex ones can overfit noise (high variance). Reconciling the classical bias-variance trade-off with the empirical observation of good generalization, despite achieving zero empirical risk on noisy data using DNNs with a regression loss, presents a significant challenge [40]. Generally speaking, in deep learning, over-parameterization and perfectly fitting noisy data do not exclude good generalization performance. Moreover, empirical evidence across various modeling techniques, including DNNs, decision trees, random features, and linear models, indicates that test error can decrease even below the sweet-spot in the U-shaped bias-variance curve when further increasing the number of parameters [15, 41, 42]. Apart from inconsistencies with experiments, under the existing theoretical framework, the upper bound of the expected risk may be too loose. In other words, empirical risk and generalization ability may not completely determine model performance. A typical example is in highly imbalanced binary classification, such as when the ratio of positive to negative samples is 99:1. Suppose a model predicts all inputs as positive; in this case, the overall accuracy of the model is 99%, and the generalization error is 0, but clearly, the model is ineffective. Within the Definition 2.1, the result learned by the model $f_{s}$ is not completely deterministic and unique; it depends on the artificially chosen loss function, rather than being entirely determined by the data. This raises deeper questions: If the learning outcome depends on the artificial selection of the loss function, to what objects in the real physical world does the learned result correspond?

Deep learning poses several mysteries within the conventional learning theory framework, including the outstanding performance of overparameterized neural networks, the impact of depth in deep architectures, the apparent absence of the curse of dimensionality, and the surprisingly successful optimization performance despite the non-convexity of the problem. Theoretical insights into these problems from prior research have predominantly followed two main avenues: expressivity (fitting ability) and generalization ability.

Essentially, expressivity is characterized by the capability of DNNs to approximate any function. The universal approximation theorem [35, 43] states that a feedforward network with any ”squashing” activation function $\sigma(\cdot)$, such as the logistic sigmoid function, can approximate any Borel measurable function $f(x)$ with any desired non-zero amount of error, provided that the network is given enough hidden layer size $n$. Universal approximation theorems have also been proved for a wider class of activation functions, which includes the now commonly used ReLU [1]. Generally speaking, the greater the number of parameters a model possesses, the stronger its fitting capability. Moreover, research has demonstrated that extremely narrow neural networks can retain universal approximativity with an appropriate choice of depth [44, 45]. More precisely, enhancing the depth of a neural network architecture exponentially increases its expressiveness [46, 47].

Generalization ability represents a model’s predictive ability on unseen data. Increasing the depth of a neural network architecture naturally yields a highly overparameterized model, whose loss landscape is known for a proliferation of local minima and saddle points [48]. The fact that sometimes (for some specific models and datasets), the generalization ability of a model paradoxically increases with the number of parameters, even surpassing the interpolation threshold—where the model fits the training data perfectly. This observation is highly counterintuitive since this regimen usually corresponds to overfitting. DNNs show empirically that after the threshold, sometimes the generalization error tends to descend again, showing the so called Double-Descent generalization curve, and reaching a better minimum. In order to address this question, researchers have tried to focus on different methods and studied the phenomena under different names, e.g., over-parametrization [49, 50, 51, 52], double descent [53, 54, 55, 56], sharp/flat minima [57, 58, 59, 60, 61, 62], and benign overfit [63, 64, 65, 66]. In this paper, we review the related literature and organize it into the following six categories.

1. 1.

   Complexity-based methods. Conventional statistical learning theory has established a series of upper bounds on the generalization error (generalization bounds) based on the complexity of the hypothesis space, such as VC-dimension [67, 4] and Rademacher complexity [5]. Usually, these generalization bounds suggest that controlling the model size can help models generalize better. However, overparameterized DNNs make the generalization bounds vacuous.

2. 2.

   Algorithmic Stability. Algorithmic Stability (e.g., uniform stability, hypothesis stability) is straightforward: if an algorithm fits the dataset without being overly dependent on particular data points, then achieving small empirical error implies better generalization [6, 7, 8, 9, 10, 11]. Algorithmic Stability simply cares about the fact that the algorithm choice should not change too much and is stable enough if we slightly change the dataset.

3. 3.

   PAC-Bayesian. PAC-Bayesian (or simply “PAC-Bayes”) theory, which integrates Bayesian statistic theories and PAC theory for stochastic classifiers, has emerged as a method for obtaining some of the tightest generalization bounds [68, 69]. It emerged in the late 1990s, catalyzed by an early paper [70], and shepherded forward by McAllester [71, 72]. McAllester’s PAC-Bayesian bounds are empirical bounds, in the sense that the upper bound only depends on known computable quantities linked to the data. Over the years, many PAC-Bayes bounds were developed in many ever tighter variants, such as the ones for Bernoulli losses [73, 74, 75], exhibiting fast-rates given small loss variances [76, 77], data-dependent priors [78], and so on. We refer the reader to [79] for excellent surveys.

4. 4.

   The stochasticity of stochastic optimization algorithms as an implicit regularization. Due to the non-linear dependencies within the loss function relative to the parameters, the training of DNNs is a quintessential non-convex optimization problem, fraught with the complexities of numerous local minima. First-order optimization methods, such as SGD, have become the workhorse behind the training of DNNs. Despite its simplicity, SGD also enables high efficiency in complex and non-convex optimization problems [80]. It is widely believed that the inherent stochasticity of stochastic optimization algorithms performs implicit regularization, guiding parameters towards flat minima with lower generalization errors [81, 82, 48, 83]. To identify the flat minima, numerous mathematical definitions of flatness exist, with a prevalent one being the largest eigenvalue of the Hessian matrix of the loss function with respect to the model parameters on the training set, denoted as $\lambda_{\max}$ [84, 85, 58].

5. 5.

   The role of over-parameterization in non-convex optimization. As for the impact of parameter quantity on optimization algorithms, it is already known in the literature that DNNs in the infinite width limit are equivalent to a Gaussian process [86, 87, 88, 89, 90]. The work [91] elucidates that the evolution of the trainable parameters in continuous-width DNNs during training can be captured by the neural tangent kernel (NTK). Recent work has shown that with a specialized scaling and random initialization, the parameters of continuous width two-layer DNNs are restricted in an infinitesimal region around the initialization and can be regarded as a linear model with infinite dimensional features [49, 92, 93, 94]. Since the system becomes linear, the dynamics of gradient descent (GD) within this region can be tracked via properties of the associated NTK and the convergence to the global optima with a linear rate can be proved. Later, the NTK analysis of global convergence is extended to multi-layer neural nets [95, 96, 97]. An alternative research direction attempts to examine the infinite-width neural network from a mean-field perspective [98, 99, 100]. The key idea is to characterize the learning dynamics of noisy SGD as the gradient flows over the space of probability distributions of neural network parameters. When the time goes to infinity, the noisy SGD converges to the unique global optimum of the convex objective function for the two-layer neural network with continuous width. NTK and mean-field methods provide a fundamental understanding on training dynamics and convergence properties of non-convex optimization in infinite-width neural networks.

6. 6.

   Information-based methods. Information-theoretic generalization bounds have been empirically demonstrated to be effective in capturing the generalization behavior of Deep Neural Networks (DNNs) [101, 102]. The original information-theoretic bound introduced by [103] has been extended or improved in various ways, such as the chaining method [104, 105, 106], the random subset and individual technique [101, 107, 108, 109, 110]. The information bottleneck principle [111, 112, 113], from the perspective of the information change in neural networks, posits that for better prediction, the neural network learns useful information and filters out redundant information. Beyond supervised learning, information-theoretic bounds are also useful in a variety of learning scenarios, including meta-learning [114], semi-supervised learning [115], transfer learning [116], among others.

While there is yet no clarity in understanding why, for some algorithms and datasets, increasing the number of parameters actually improves generalization [117]. The accumulating evidence and existing theories indicate that the exceptional performance of DNNs is attributed to a confluence of factors, including their architectural design, optimization algorithms, and even training data. Elucidating the mechanisms underlying deep learning remains an unresolved issue at present.

In statistical data analysis, determining the probability density function (PDF) or probability mass function (PMF) from a limited number of random samples is a common task, especially when dealing with complex systems where accurate modeling is challenging. Generally, there are three main methods for this problem: parametric, semi-parametric, and non-parametric models.

Parametric models assume that the data follows a specific distribution with unknown parameters. Examples include the Gaussian model, Gaussian mixture models (GMMs), and distributions from the generalized exponential family. A major limitation of these methods is that they make quite simple parametric assumptions, which may not be sufficient or verifiable in practice [118]. Non-parametric models estimate the PDF/PMF based on data points without making assumptions about the underlying distribution [18]. Examples include the histogram [119, 120, 121] and kernel density estimation (KDE) [122, 123, 124]. While these methods are more flexible, they may suffer from overfitting or poor performance in high-dimensional spaces due to the curse of dimensionality. Semi-parametric models [19] strike a balance between parametric and non-parametric methods. They do not assume the a priori-shape of the PDF/PMF to estimate as non-parametric techniques. However, unlike the non-parametric methods, the complexity of the model is fixed in advance, in order to avoid a prohibitive increase of the number of parameters with the size of the dataset, and to limit the risk of overfitting. Finite mixture models are commonly used to serve this purpose.

Statistical methods for determining the parameters of a model conditional on the data, including the standard likelihood method [17, 125, 20], M-estimators [126, 127], and Bayesian methods [128, 129, 130], have been widely employed. However, traditional estimators struggle in high-dimensional space: because the required sample size or computational complexity often explodes geometrically with increasing dimensionality, they have almost no practical value for high-dimensional data.

Recently, neural network-based methods have been proposed for PDF/PMF estimation and have shown promising results in problems with high-dimensional data points such as images. There are mainly two families of such neural density estimators: auto-regressive models [131, 132] and normalizing flows [133, 134]. Autoregression-based neural density estimators decompose the density into the product of conditional densities based on probability chain rule $p(x)=\prod_{i}p(x_{i}|x_{1:i-1})$. Each conditional probability $p(x_{i}|x_{1:i-1})$ is modeled by a parametric density (e.g., Gaussian or mixture of Gaussians), and the parameters are learned by neural networks. Density estimators based on normalizing flows represent $x$ as an invertible transformation of a latent variable $z$ with a known density, where the invertible transformation is a composition of a series of simple functions whose Jacobian is easy to compute. The parameters of these component functions are then learned by neural networks.

• 1.

  The set of real numbers is denoted by $\mathbb{R}$, and the cardinality of set $\mathcal{Z}$ is denoted by $|\mathcal{Z}|$.

• 2.

  Random variables are denoted using upper case letters such as $X,Y$, which take values in sets $\mathcal{X}$ and $\mathcal{Y}$, respectively.

• 3.

  The expectation of a random variable $X$ is denoted by $\mathbb{E}[X]$ or $\mathbb{E}_{X}[X]$, with $\mathbb{E}_{X\sim p}[X]$ when the distribution of $X$ is specified. We define $\mathbb{E}_{X\sim p}[f(X)]=\sum_{x\in\mathcal{X}}p(X=x)f(x)$, where $\mathcal{X}$ is the value space of $X$.

• 4.

  The space of models, which is a set of functions endowed with some structure, is represented by $\mathcal{F}_{\Theta}=\{f_{\theta}:\theta\in\Theta\}$, where $\Theta$ denotes the parameter space. Similarly, we define the hypothesis space with feature $x$ as $\mathcal{F}_{\Theta,x}:=\{f_{\theta,x}:\theta\in\Theta\}$.

• 5.

  The sub-differential is defined as $\partial f(x)=\{s\in\mathbb{R}^{n}\mid f(y)-f(x)-\langle s,y-x\rangle\geq 0,\forall y\in\mathrm{dom}(f)\}$ (see e.g., [135], Page 27). If $f(x)$ is smooth and convex, the sub-differential $\partial f(x)$ reduces to the usual gradient, i.e., $\partial f(x)=\{\nabla f(x)\}$.

• 6.

  The Legendre-Fenchel conjugate of a function $\Omega$ is denoted by $\Omega^{*}(f):=\sup_{p\in\mathrm{dom}(\Omega)}\langle p,f\rangle-\Omega(p)$. By default, $\Omega$ is a continuous strictly convex function, and its gradient with respect to $p$ is denoted by $p_{\Omega}^{*}$. When $\Omega(\cdot)=\frac{1}{2}\|\cdot\|_{2}^{2}$, we have $p=p_{\Omega}^{*}$.

• 7.

  The Fenchel-Young loss $d_{\Omega}\colon\mathrm{dom}(\Omega)\times\mathrm{dom}(\Omega^{*})\to\mathbb{R}^{+}$ generated by $\Omega$ is defined as [136]:

  $$d_{\Omega}(q;f):=\Omega(q)+\Omega^{*}(f)-\langle q,f\rangle,$$

  (1)

  where $\Omega^{*}$ denotes the conjugate of $\Omega$.

• 8.

  The maximum and minimum eigenvalues of matrix $A$ are denoted by $\lambda_{\max}(A)$ and $\lambda_{\min}(A)$, respectively.

• 9.

  To improve clarity and conciseness, we provide the following definitions of symbols:

  		$\displaystyle q_{Y|x}(y):=q_{Y|X}(y|x),$		(2)
  		$\displaystyle q_{Y|x}:=(q_{Y|x}(y^{(1)}),\cdots,q_{Y|x}(y^{(|\mathcal{Y}|)}))^{\top},$
  		$\displaystyle\bar{q}:=\mathbb{E}_{q\sim\mathcal{Q}}[q],$
  		$\displaystyle\bar{q}_{X}:=\mathbb{E}_{q_{X}\sim\mathcal{Q}_{X}}[q_{X}],$
  		$\displaystyle\bar{q}_{Y|x}:=\mathbb{E}_{q_{Y|x}\sim\mathcal{Q}_{Y|x}}[q_{Y|x}],$
  		$\displaystyle\mathrm{Var}(q_{Y|x}):=\mathbb{E}_{q_{Y|x}\sim\mathcal{Q}_{Y|x}}[\|q_{Y|x}-\bar{q}_{Y|x}\|_{2}^{2}],$

  where $q_{X}(x)$ and $q_{Y|X}(y|x)$ represent the marginal and conditional PMFs/PDFs, respectively. The conditional distribution $q_{Y|x}$ is also modeled as a random variable.

The theorems in the subsequent sections rely on the following lemmas.

Throughout this paper, in line with practical applications and without affecting the conclusions, we assume that $\Omega$ is a proper, lower semi-continuous (l.s.c.), and strictly convex function or functional. Consequently, the Fenchel-Young loss attains zero if and only if $q=f_{\Omega^{*}}^{*}$.

The random pair $Z=(X,Y)$ follows the distribution $q$ and takes values in the product space $\mathcal{Z}=\mathcal{X}\times\mathcal{Y}$, where $\mathcal{X}$ represents the input feature space and $\mathcal{Y}$ denotes a finite set of labels. The training data, denoted by an $n$-tuple $s^{n}:=\{z^{(i)}\}_{i=1}^{n}=\{(x^{(i)},y^{(i)})\}_{i=1}^{n}$, is composed of i.i.d. samples drawn from the unknown true distribution $q$. Given the available information, the true underlying distribution is undetermined. Consequently, we model $q$ as a random variable in the space $\Delta$. Let $\mathcal{Q}^{(s^{n},\mathcal{P})}$ (abbreviated as $\mathcal{Q}$) be the posterior distribution of $q$, conditional on the i.i.d. sample set $s^{n}$ and the prior distribution $\mathcal{P}$. We define $\bar{q}=\mathbb{E}_{q\sim\mathcal{Q}}q$. Furthermore, let $\tilde{q}$ denote an unbiased and consistent estimator of $\bar{q}$. In the absence of prior knowledge about $q$, the Monte Carlo method (MC) is commonly utilized. This method estimates $q$ based on the event frequencies, denoted by $\hat{q}(s^{n}):=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\{z\}}(z^{(i)})$. In this article, we refer to $\hat{q}$ as the frequency distribution estimator to distinguish it from the general estimator $\tilde{q}$. We also assume that $\tilde{q}$ is a better estimator relative to $\hat{q}$. Real-world distributions are usually smooth, which inspires the common practice of smoothing $\hat{q}$ to construct $\tilde{q}$ to approximate the true distribution $q$. These smoothing techniques effectively integrate prior knowledge and are often advantageous, as exemplified by classical Laplacian smoothing. For simplicity, when the sample set is given, we abbreviate the results estimated by the estimator $\tilde{q}(s^{n})$ and $\hat{q}(s^{n})$ as $\tilde{q}$ and $\hat{q}$, respectively. It is important to note that, given the sample set $s^{n}$ and the prior distribution $\mathcal{P}$, calculating $\bar{q}$ and $\mathcal{Q}$ is feasible. The specific calculation methods are detailed in Theorem 4.2. In this article, we utilize the notation $f_{\theta}$ (hereinafter abbreviated as $f$) to denote the model characterized by the parameter vector $\theta$. Additionally, $f_{\theta,x}$ (abbreviated as $f_{x}$) represents the model with a parameter vector $\theta$ and a specific input $x$. For ease of reference, we list the symbol variables defined above in Table 1.

This subsection introduces the optimization objectives, related definitions, and the necessary conditions that loss functions, models, and optimization algorithms must satisfy within the PD learning framework.

We refer to the tasks with learning objectives of $q$ and $q_{Y|X}$ as native PD learning and conditional PD learning, respectively. Since the learning target is modeled as a random variable, the learning errors in native PD learning and conditional PD learning are defined as follows:

where $q_{X}$, $q_{Y|X}$, and $q$ are modeled as random variables. We thus formulate PD learning in two situations:

• 1.

  Native PD learning

  $$\min_{\theta}L_{n}(q,p)\,\ \mathrm{s.t.}\,\ s=\{z_{i}\}_{i=1}^{n},$$

  (8)

  where $p=(f_{\theta})_{\Omega^{*}}^{*}$, $\epsilon\in\mathbb{R}^{+}$.

• 2.

  Conditional PD learning

  $$\min_{\theta}L_{c}(q,p)\,\ \mathrm{s.t.}\,\ s=\{z_{i}\}_{i=1}^{n},$$

  (9)

  where $p_{Y|x}=(f_{\theta,x})_{\Omega^{*}}^{*}$, $\epsilon\in\mathbb{R}^{+}$.

In PD learning, we utilize the Legendre-Fenchel conjugates of $f_{\theta}$ and $f_{\theta,x}$, denoted as $p$ and $p_{Y|x}$, to approximate the $q$ and $q_{Y|x}$. The purpose of this treatment has two aspects: first, it provides a general method to map the model’s output from Euclidean space to the probability simplex; second, this approach aligns with the form of the standard loss function proposed in this paper. Native PD learning can be considered as a specific case of conditional PD learning when the feature variable $X$ is fixed, i.e., $X\sim\delta_{x}$.

We list the definitions used in the PD learning framework as follows:

1. 1.

   Fitting error. This error measures the performance of the model in fitting the target $\tilde{q}$, denoted by $\mathrm{Fit}(\tilde{q},p)=\|p-\tilde{q}\|_{2}^{2}$.

2. 2.

   Uncertainty. We refer to the variances $\mathrm{Var}(q)$ and $\mathrm{Var}(q_{Y|x})$ as the uncertainties of $q$ and $q_{Y|x}$, respectively. They are independent of the models $f$ and $f_{x}$ and depend solely on the sample set and prior knowledge.

3. 3.

   Estimation error. This error quantifies the distance between $\bar{q}$ and the corresponding provided estimate, given by $\mathrm{Est}(\tilde{q})=\|\tilde{q}-\bar{q}\|_{2}^{2}$. The estimation error arises from the finite sample size and serves as a model-independent metric that characterizes the estimator. Since conditional PD learning involves estimation errors in both the feature and label dimensions, we formally define the label estimation error as $\mathrm{Est}(\tilde{q}_{Y|x})$ and the feature estimation error as $\mathrm{Est}(\tilde{q}_{X})$.

4. 4.

   Standard loss function: A loss function $\ell(\tilde{q},p):\Delta\times\Delta\to\mathbb{R}$ is regarded as a standard loss if it meets the following conditions :

   • (a)

     Unique optimum. The loss function is strictly convex, ensuring that it has a unique minimum when predicting the correct probability $\tilde{q}$, i.e., $p=\tilde{q}$. The unique optimum condition guarantees that the standard loss function is effective in directing the model towards the desired fitting target.

   • (b)

     Sampling stability. The expected loss, $\mathbb{E}_{q}\ell(q,p)$, achieves its minimum value when $p=\mathbb{E}_{q}[q]$, expressed as $\mathbb{E}_{q}[q]=\arg\min_{p}\mathbb{E}_{q}[\ell(q,p)]$. The sampling stability condition ensures the consistency of optimization results across different optimization methods, such as SGD, batch-GD, and GD. In other words, minimizing the loss function on the entire dataset yields results consistent with those obtained by minimizing the expected value of the loss function under random sampling.

5. 5.

   Structural matrix. We define $A:=\nabla_{\theta}f^{\top}\nabla_{\theta}f$ and $A_{x}:=\nabla_{\theta}f_{x}^{\top}\nabla_{\theta}f_{x}$ as the structural matrices corresponding to the model $f$ and $f_{x}$.

6. 6.

   Structural error. We define the structural error as follows:

   $$S(\alpha,\beta,\gamma,A)=\alpha G(A)+\beta U(A)+\gamma L(A),$$

   (10)

   where

   	$\displaystyle U(A)$	$\displaystyle=-\log\lambda_{\min}(A),$
   	$\displaystyle L(A)$	$\displaystyle=-\log\lambda_{\max}(A),$
   	$\displaystyle G(A)$	$\displaystyle=U(A)-L(A).$

   Here, $\alpha$, $\beta$, and $\gamma$ are positive real numbers representing the weights associated with $G(A)$, $U(A)$, and $L(A)$, respectively.

7. 7.

   GSC algorithm. We refer to the optimization algorithm capable of simultaneously reducing both the gradient norm and the structural error as the GSC algorithm . The GSC algorithm for native PD learning is mathematically formulated as follows:

   $$\min_{\theta}\{\log\|\nabla_{\theta}\ell(\tilde{q},p)\|_{2}^{2}+S(\alpha,\beta,\gamma,A)\}.$$

   (11)

   Similarly, the GSC algorithm under conditional PD learning is formulated as

   $$\min_{\theta}\mathbb{E}_{X}[\log\|\nabla_{\theta}\ell(\tilde{q}_{Y|X},p_{Y|X})\|_{2}^{2}+S(\alpha,\beta,\gamma,A_{X})].$$

   (12)

The loss function, model, and optimization algorithm within the PD learning framework must satisfy the following conditions :

• 1.

  Loss function. The loss function $\ell(\tilde{q},p):\Delta\times\Delta\to\mathbb{R}$ is a standard loss function.

• 2.

  Model. For any $\theta\in\Theta$, $f_{\theta}$ and $f_{\theta,x}$ possess gradients or subgradients. This condition ensures that we can search for the model’s stationary points using gradient or subgradient optimization algorithms.

• 3.

  Optimization algorithm. The GSC algorithm is applied.

In Subsections 4.3 and 4.4, we will demonstrate that the conditions for loss functions and optimization algorithms are also met in real-world deep learning scenarios.

Based on the aforementioned definitions, employing Bias-Variance Decomposition, we obtain the following theorem.

The calculation method for $\bar{q}$ is outlined in the following theorem, and this method is equally applicable to $\bar{q}_{X}$ and $\bar{q}_{Y|x}$.

In the absence of further information regarding the i.i.d. sample $s^{n}$ for $q$, we invoke the equal-probability hypothesis, according to which $\mathcal{P}$ is assumed to be the uniform distribution across the space $\Delta$. The equal-probability hypothesis is one of the fundamental assumptions of statistical physics, positing that all probabilities are uniformly distributed over the space $\Delta$. Assuming the equal-probability hypothesis holds, it follows that

Based on Theorem 4.1, optimizing the learning error is equivalent to optimizing the expected fitting error of the model for $\bar{q}$ and $\bar{q}_{Y|x}$ under $X\sim\bar{q}_{X}$, as shown in the following formula:

Clearly, $p=\bar{q}$, $p_{Y|x}=\bar{q}_{Y|x}$, and $p_{X}=\bar{q}_{X}$ are the optimal solutions for the optimization objectives. However, when the prior distribution $\mathcal{P}$ is unknown or the equal-probability hypothesis clearly does not hold, Theorem 4.2 is not applicable. Therefore, in practice, we use simpler unbiased and consistent estimators based on $s^{n}$ as substitutes, as shown below:

Thus, the model’s learning error comprises three parts: 1. Uncertainty that depends solely on the data, 2. The fitting error of the model for $\tilde{q}$ and $\tilde{q}_{Y|x}$, and 3. The estimation error generated by $\tilde{q}$, $\tilde{q}_{X}$, and $\tilde{q}_{Y|x}$.

In this subsection, we provide a comprehensive description of the PD learning framework and summarize the algorithm logic and components in PD learning. In addition, we compare PD learning with the traditional deep learning framework to clarify their similarities and differences, highlighting the advantages inherent in the PD learning framework.

We summarize the algorithm logics of native PD learning and conditional PD learning in Figure 1, where the white box denotes the intermediate quantity of reasoning, whereas the gray box signifies the final processable components. We illustrate the algorithmic logic of the PD learning framework using native PD learning as an example.

1. 1.

   The learning error can be decomposed into two components: model-independent uncertainty and the fitting error of the model with respect to $\bar{q}$.

2. 2.

   The uncertainty depends solely on the sampling data and represents the lower bound of the learning error, as shown in Theorem 4.1. As the amount of sampling data increases, the uncertainty gradually approaches zero, as described in Theorem 4.6.

3. 3.

   In high-dimensional scenarios, the computational cost of calculating $\bar{q}$ is substantial as shown in Theorem 4.2. Therefore, we use an unbiased and consistent estimator $\tilde{q}$ of $\bar{q}$. This estimator serves as a bridge between the optimal estimator $\bar{q}$ and the model output $p$. Thus, minimizing the fitting error $\mathrm{Fit}(\bar{q},p)$ is equivalent to minimizing both the fitting error $\mathrm{Fit}(\tilde{q},p)$ and the estimation error $\mathrm{Est}(\tilde{q})$.

4. 4.

   The estimation error $\mathrm{Est}(\tilde{q})$ does not depend on the model and can only be reduced by leveraging prior knowledge or increasing the sample size. For instance, $\hat{q}$ is theoretically feasible, in practical scenarios, especially in the case of insufficient samples, $\hat{q}$ may not be the best choice. Based on the prior knowledge, we can generate more accurate estimators. If the underlying true distribution is known to be smooth, we can construct $\tilde{q}$ by applying Laplacian smoothing to $\hat{q}$, which often yields smaller estimation error and superior performance. In the field of machine learning, techniques such as label smoothing [141] and synthetic minority oversampling technique (SMOTE) [142], which introduce a degree of smoothness to the existing data, also yield favorable outcomes.

5. 5.

   Optimization of the fitting error, given by $\min_{\theta}\mathrm{Fit}(\tilde{q},p)$, is a typical non-convex problem. Simple gradient-based optimization algorithms do not guarantee a global optimal solution. By using the standard loss function, the GSC algorithm facilitates the approximation of global optimal solutions in fitting error minimization by minimizing both $\|\nabla_{\theta}d_{\Omega}(\tilde{q},f)\|_{2}^{2}$ and the structural error as described in Theorem 4.4 and Corollary 4.1. Gradient-based optimization algorithms can reduce $\|\nabla_{\theta}d_{\Omega}(\tilde{q},f)\|_{2}^{2}$, while the effectiveness of optimizing the structural error is guaranteed by Theorem 4.5.

6. 6.

   Based on the above conclusions, both model-dependent and model-independent bounds on learning error can be derived, as shown in Theorem 4.7 and Theorem 4.8.

Corresponding to the algorithm logic above, the PD learning framework encompasses five principal components: standard loss, model, prior knowledge ($\mathcal{P}$), estimator ($\tilde{q}$), and the GSC algorithm (represented by $\mathcal{A}$), as shown in Figure 2. The estimator yields an approximation of $\bar{q}$, denoted as $\tilde{q}$, leveraging both the sample set $s^{n}$ and prior knowledge $\mathcal{P}$. The standard loss function proposed is a type of loss function that possesses a unique optimal solution and sampling stability. We have demonstrated that under these conditions, it has a unique mathematical expression and can be represented as a Young–Fenchel loss function, which encompasses various common losses as shown in Theorem 4.3. The model in PD learning is utilized to fit the target distribution $\tilde{q}$ via its Legendre-Fenchel conjugate function. The GSC algorithm is designed to minimize both the gradient norm and the structural error.

To elucidate PD learning more effectively, we compare it with the conventional deep learning theory framework.

• 1.

  Learning objective. The underlying true distribution $q$ or $q_{Y|x}$ represents the learning objectives in PD learning. Within the conventional deep learning theoretical framework, the primary objective is to ascertain the unknown mapping between features and labels. However, the existence and identifiability of these mappings are uncertain, and the mappings that are learned depend on the loss function chosen.

• 2.

  Modeling method. In PD learning, the learning objective is conceptualized as a random variable within the probability simplex. Conversely, within the conventional deep learning theoretical framework, the unknown mapping is represented as a function that processes input features to generate corresponding labels within a defined function space $\mathcal{F}_{\Theta}$. When the unknown mapping is not a functional mapping, it cannot be described by $\mathcal{F}_{\Theta}$.

• 3.

  Optimization objective and its reformulation. The optimization objective in PD learning is the minimization of learning error, as defined in (9). This objective is achieved by controlling the estimation error and the fitting error. In contrast, within the conventional deep learning theoretical framework, the optimization objective is the minimization of expected risk. Expected risk minimization is reformulated as Empirical Risk Minimization (ERM) and generalization error minimization, collectively referred to as Structural Risk Minimization (SRM).

• 4.

  Fitting method of model. In PD learning, the target is fitted by employing the Legendre-Fenchel conjugate function of the model’s output. Conversely, in the conventional deep learning theoretical framework, neural networks are typically used to directly fit the target without the intermediary step of the conjugate function.

• 5.

  Loss function and its role. In PD learning, the standard loss function is used to facilitate the minimization of the gradient norm. In contrast, within the conventional deep learning theoretical framework, a variety of loss functions may be employed. The minimization of the loss function, in conjunction with the regularization term corresponding to SRM, serves as the optimization objective.

• 6.

  Optimization algorithm. The GSC algorithm is used to optimize the fitting error, which can ensure obtaining an approximate global optimal solution for fitting error minimization in the PD learning framework. Within the conventional deep learning theoretical framework, gradient-based iterative algorithms, including GD, SGD, and Adaptive Moment Estimation(Adam), are utilized to achieve SRM. While these algorithms strive to identify the global optimal solution, they do not ensure its attainment.

• 7.

  Bounds on Learning Performance. In the PD learning framework, both model-dependent and model-independent lower and upper bounds on learning error are provided. These bounds characterize the impact of model fitting capacity, data uncertainty, and mutual information between features and labels on the learning error. Within the conventional deep learning theoretical framework, existing evidence suggests that classical bounds from statistical learning theory alone may be insufficient to fully capture the behavior of DNNs.

In this subsection, the equivalence between standard losses and the Fenchel-Young loss is established. Common loss functions are all standard loss functions, which also demonstrate the generality of the concept of standard loss functions.

In our study, we represent the standard loss as $d_{\Omega}(q,f_{\theta})$, where $f_{\theta}$ denotes the model parameterized by $\theta$. To maintain generality without compromising clarity, we assume that $d_{\Omega}(q,f_{\theta})$ is continuous and differentiable with respect to $\theta$. The standard loss function is a broad concept that encompasses a wide range of commonly used loss functions in statistics and machine learning, as illustrated in Table 2.

This subsection demonstrates that under the conditions specified in 4.1, the model can converge to the global optimal solution of the fitting error minimization problem as follows:

Specifically, we provide upper and lower bounds for the fitting error based on the gradient norm of the loss function and the structural error through the following theorem and corollary. These upper and lower bounds can converge to zero under the conditions specified in 4.1, thereby ensuring that the optimization objective (LABEL:eq:fitting_obj) can achieve the global optimal solution.

First, we provide the upper and lower bounds of the fitting error under the condition of using a standard loss function.

This theorem leads to the following significant conclusions:

1. 1.

   Since $\mathrm{Fit}(\tilde{q},p)$ is equivalent to $\|\nabla_{\theta}d_{\Omega}(\tilde{q},f)\|_{2}$ when $\lambda_{\min}(A)\neq 0$, in the context of the aforementioned non-convex optimization problem LABEL:eq:fitting_obj, the global optima can be approximated by concurrently minimizing $S(\alpha,\beta,\gamma,A)$ and $\|\nabla_{\theta}d_{\Omega}(\tilde{q},f)\|$.

2. 2.

   The gradient-based/sub-gradient-based iterative algorithms, such as GD and SGD, allow the model parameters to converge to stationary points, where $\|\nabla_{\theta}d_{\Omega}(\tilde{q},f)\|_{2}=0$.

3. 3.

   The above theorem also provides the condition under which the gradient-based iterative algorithms can converge to a globally optimal solution, i.e., $S(\alpha,\beta,\gamma,A)=0$.

Based on the established Theorem 4.4, we can derive the upper bound on the fitting error in conditional PD learning.

Similarly, in conditional PD learning, the gradient-based iterative algorithms are able to converge to a globally optimal solution, provided that the following condition is satisfied:

Theorem 4.4 and Corollary 4.1 provide a theoretical foundation for optimizing fitting error using the GSC algorithm. Next, we will analyze and explain how to implement the GSC algorithm to minimizing fitting error. According to Definition 7, the GSC algorithm aims to simultaneously reduce both the gradient norm and the structural error. Gradient-based optimization methods, including SGD, Batch-GD, and Newton’s method, primarily focus on minimizing the squared $L_{2}$ norm of the gradient, $\|\nabla_{\theta}d(\tilde{q},f)\|_{2}^{2}$. Alternative methods, such as simulated annealing and genetic algorithms, can also be employed as complementary strategies to achieve this objective. Regarding the optimization of $S(\alpha,\beta,\gamma,A)$, we next demonstrate that the architecture of the model, the number of parameters, and the independence among parameter gradients can be leveraged to diminish structural error.

We establish the following condition as the foundation for our analysis.

This condition essentially requires that the derivative values of the model’s output with respect to the parameters are finite and remain approximately independent of each other. In practical scenarios, gradients are typically finite, allowing for an approximate equivalence between the independence of parameter gradients and gradient independence condition. It should be noted that the aforementioned condition is not inherently and unconditionally valid. In section 5, we will explain how numerous existing deep learning techniques contribute to the substantiation of this condition. Based on gradient independence condition, we deduce the following theorem to provide a theoretical basis for reducing structural error.

The derivation for this theorem is presented in Appendix A.3. The aforementioned theorem indicates that, under the gradient independence condition, the structural error can be controlled through the parameter count $|\theta|$ and the ball radius $\epsilon$. Consequently, based on the theorem, the following are the specific methodologies for reducing structural error:

• 1.

  Since $Z(|\theta|,|\tilde{q}|)$ decreases as $|\theta|$ increases, increasing the number of parameters can help reduce the $G(A)$ under the gradient independence condition. The proposed strategies include:

  1. (a)

     Mitigating the dependencies among parameters to facilitate the validity of the gradient independence condition;

  2. (b)

     Increasing the number of model parameters when the model aligns with the gradient independence condition.

  3. (c)

     Utilizing a specialized model architecture to ensure that $G(A)\approx 0$ ($A\approx aI$), and this approximation is maintained throughout the training process.

  In section 5, we will elucidate that numerous existing deep learning techniques essentially correspond to the above strategies.

• 2.

  Given that both $U(A)$ and $L(A)$ demonstrate a negative correlation with $\epsilon^{2}$, while $G(A)$ remains independent of $\epsilon$, it follows that augmenting $\|\nabla_{\theta}f_{i}\|_{2}^{2}$ can effectively contribute to the reduction of structural error.

  1. (a)

     Gradient-based iteration methods with stochasticity. Following random parameter initialization, the gradient norm $\|\nabla_{\theta}f_{i}\|_{2}^{2}$ is typically small. The application of gradient-based iteration methods subsequently increases $\|\nabla_{\theta}f_{i}\|_{2}^{2}$. This increase affects the independence of gradients, leading to an increasement in $G(A)$. However, as $\|\nabla_{\theta}f_{i}\|_{2}^{2}$ stabilizes over iterations, the stochasticity introduced during the process gradually re-establish the gradient independence condition, causing $G(A)$ to converge. This analytical conclusion is supported by the experimental results presented in Figures 6 and 7.

  2. (b)

     Application of residual blocks. Since the gradient norm $\|\nabla_{\theta}f_{i}\|_{2}^{2}$ is initially small by random parameter initialization, incorporating residual blocks can effectively enhance $\|\nabla_{\theta}f_{i}\|_{2}^{2}$, thereby reducing $U(A)$ and $L(A)$. However, this strategy leads to the gradient independence condition no longer holding, thereby causing $G(A)$ to increase. This finding has been experimentally validated, as illustrated in Figure 7.

Using the GSC algorithm to optimize fitting error means that the gradient norm can control the fitting error only when the structural error remains constant. Therefore, when employing random parameter initialization and SGD, model training can be divided into the following two stages:

1. 1.

   Structural error optimization stage. Random parameter initialization methods ensure that the gradient independence condition is approximately satisfied. Random parameter initialization methods, along with feature normalization techniques, typically center input features and parameter values around 0, resulting in a relatively small $\|\nabla_{\theta}f_{i}\|_{2}^{2}$ and a correspondingly smaller ball radius $\epsilon$. At the onset of training, SGD leads to an increase in $\|\nabla_{\theta}f\|_{2}^{2}$, thereby increasing $\epsilon$. Since both $U(A_{x})$ and $L(A_{x})$ exhibit a negative correlation with $\epsilon$, they consequently decrease as $\epsilon$ increases. As training progresses, the backpropagated error decreases to a level insufficient to significantly alter $\|\nabla_{\theta}f\|_{2}^{2}$, leading to the stabilization of $U(A_{x})$, $L(A_{x})$, and $G(A_{x})$.

2. 2.

   Gradient norm control stage. After the structural error converges to a stable value, the fitting error begins to be controlled by the gradient norm and subsequently decreases only after the structural error has approached stability.

In other words, from the perspective of the GSC algorithm in analyzing the training process of deep models, the convergence of structural error is a prerequisite for the convergence of fitting error. This theoretical prediction is validated through experimental results 6.

This subsection examines the properties of the estimation error and uncertainty. Since both $\tilde{q}$ and $\hat{q}$ serve as unbiased and consistent estimators of $\bar{q}$, we can derive the following theorem.

This theorem suggests that, given a sufficiently large sample size, the uncertainty and estimation error associated with an arbitrary unbiased and consistent estimator are negligible. It is evident that this theorem also holds true in the context of conditional PD learning. Beyond the sample set, the prior distribution $\mathcal{P}$ also exerts a significant influence on the estimation error. For a typical deterministic single-label classification with negligible noise, the label probability distribution conditioned on feature $x$ is given by $q_{Y|x}=\delta_{y}$, where $y$ is the label associated with $x$. Under these conditions, it follows that $\forall x\in\mathcal{X}$, $\mathrm{Est}(\tilde{q}_{Y|x})=0$ and $\mathrm{Var}(q_{Y|x})=0$.