Laplacian-weighted Gradient Descent

• Laplacian-weighted Gradient Descent is a family of first-order optimization methods that replaces standard gradients with Laplacian-derived transformations for smoother updates.
• Its methodology employs the inverse of a circulant Laplacian matrix, which reduces gradient variance, improves conditioning, and permits larger stable learning rates.
• Empirical results across logistic regression, deep neural networks, and GANs show faster convergence, better generalization, and reduced sensitivity to noise.

Searching arXiv for recent and foundational papers on Laplacian-weighted gradient methods. Laplacian-weighted Gradient Descent can be understood as a family of first-order optimization methods in which the usual gradient, or the loss terms that determine it, is transformed by a Laplacian-derived operator or by Laplacian-based weights. The canonical instance is Laplacian Smoothing Gradient Descent (LSGD), which replaces the usual gradient in GD or SGD by the inverse of a positive definitive matrix coming from a one-dimensional discrete Laplacian or its high order generalizations; when applied to problems ranging from logistic regression to deep neural nets, the method is reported to dramatically reduce the variance, allow to take a larger step size, and improve the generalization accuracy (Osher et al., 2018). Related papers extend this construction to deterministic saddle-point avoidance, stochastic-gradient Langevin sampling, and differentiable reformulations of linearly-constrained $\ell_1$ minimization (Kreusser et al., 2019, Wang et al., 2019, Bonifaci, 2019).

1. Core formulation

For the empirical risk minimization problem

$$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$

standard GD or SGD uses the update

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$

LSGD instead uses

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$

where

$$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$

$\sigma>0$, and $\mathbf{L}$ is the discrete one-dimensional Laplacian operator on parameter indices. In the formulation emphasized in the original paper, $\mathbf{H}_\sigma$ is a circulant convolution matrix with periodic boundary conditions (Osher et al., 2018).

The matrix form of $\mathbf{H}_\sigma$ is circulant tridiagonal, with $1+2\sigma$ on the diagonal and $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$0 on the first off-diagonals and wrap-around corners. This structure is central: it makes the transformation local in coordinate space and diagonalizable by the discrete Fourier transform. The same paper also introduces higher-order generalizations,

$$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$1

including biharmonic smoothing for $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$2 (Osher et al., 2018).

A common conceptual point in this literature is that the one-dimensional Laplacian is used in parameter-space, not in data-space. The method therefore modifies the geometry of the update direction without changing the dataset or introducing graph structure over examples (Osher et al., 2018).

2. Smoothing mechanism and geometric effect

The discrete Laplacian acts as a local averaging or smoothing operator, and multiplication by $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$3 smooths high-frequency components of the gradient while preserving the mean. The paper states that Laplacian smoothing preserves the mean and has the effect of increasing the smaller components and decreasing the larger ones in the gradient, thereby compressing extremes caused by noisy gradient estimates (Osher et al., 2018).

This smoothing has three linked consequences. First, it reduces the variance of stochastic gradients. Second, it improves conditioning by damping oscillatory updates in stiff directions. Third, it permits larger stable learning rates in high dimensions because the smoothed operator effectively reduces the Lipschitz constant. In the step-size analysis given for an $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$4-Lipschitz objective, the maximal safe step size changes from $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$5 for GD to

$$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$6

where $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$7 quickly approaches $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$8 for large dimension $$\min_{\mathbf{w}} F(\mathbf{w}) := \frac{1}{n}\sum_{i=1}^n f_i(\mathbf{w}),$$9 (Osher et al., 2018).

The mechanism is therefore not merely heuristic noise filtering. It is a structured change of metric on the update direction. This suggests a useful interpretation: Laplacian weighting acts as an explicit smoothing prior on parameter increments, but does so through the optimizer rather than through an added regularization term.

3. Analytical properties

The variance reduction claim is formalized by modeling stochastic-gradient noise as Gaussian, $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$0. Under Laplacian smoothing,

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$1

and the paper gives a bound showing that the total variance is reduced relative to the unsmoothed case; the reported numerical evidence states that smoothing can reduce the variance of stochastic gradients by over an order of magnitude (Osher et al., 2018).

A distinctive theoretical feature of LSGD is its connection to Hamilton-Jacobi partial differential equations. The explicit LSGD iteration is described as an approximation of implicit gradient descent on the solution of

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$2

with Hopf-Lax representation

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$3

In this viewpoint, the implicit version of the new algorithm is almost the same as doing gradient descent on a new function that has the same global minima as the original function and is “more convex” (Osher et al., 2018).

The same paper also states that optimization algorithms with these surrogates converge uniformly in the discrete Sobolev $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$4 sense. The norm is written as

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$5

where $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$6 is the forward difference operator. Proposition-level results establish convergence in this norm and reduced optimality gap for convex optimization problems (Osher et al., 2018).

A recurrent empirical interpretation is that smoothing flattens sharp minima and convexifies the loss landscape while preserving global minimizers. In the terminology of the paper, this is associated with improved generalization because the method helps optimization avoid sharp minima (Osher et al., 2018).

4. Variants and algorithmic extensions

Several later papers adapt the Laplacian-weighted update principle to settings in which variance reduction is not the only concern.

Method	Update or objective	Reported distinguishing property
LSGD	$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$7	Variance reduction, larger step sizes, improved generalization
mLSGD	$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$8	Smaller attraction region for saddle points
LS-SGLD	$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \nabla f_{i_k}(\mathbf{w}^k).$$9	Strictly smaller discretization error in $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$0-Wasserstein distance
Laplacian approach to $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$1 minimization	$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$2	Differentiable convex reformulation leading to gradient-based IRLS-type methods

Modified Laplacian smoothing gradient descent (mLSGD) replaces the constant $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$3 by an iteration-dependent, monotonic, bounded function $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$4:

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$5

For the canonical quadratic

$$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$6

the dimension of the attraction region is stated to be $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$7 for GD and $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$8 for mLSGD; in $$\mathbf{w}^{k+1} = \mathbf{w}^k - \eta_k \mathbf{H}_\sigma^{-1}\nabla f_{i_k}(\mathbf{w}^k),$$9D, the attraction region for mLSGD reduces to the trivial set $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$0 (Kreusser et al., 2019).

Laplacian Smoothing Stochastic Gradient Markov Chain Monte Carlo transfers the same operator to stochastic-gradient Langevin dynamics. Its update is

$$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$1

The paper proves that for sampling from both log-concave and non-log-concave densities, LS-SGLD achieves strictly smaller discretization error in $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$2-Wasserstein distance, although its mixing rate can be slightly slower (Wang et al., 2019).

A distinct but related line of work uses Laplacian-weighted quantities to obtain gradient-based solvers for basis pursuit. The differentiable reformulation

$$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$3

induces gradient components

$$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$4

where $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$5 and $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$6. The resulting multiplicative update scheme

$$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$7

is analyzed as an IRLS-type method, with $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$8 iterations for the primal gradient scheme and $$\mathbf{H}_\sigma = \mathbf{I} - \sigma \mathbf{L},$$9 for the accelerated gradient scheme (Bonifaci, 2019).

5. Computational realization and empirical behavior

A major practical reason for the popularity of the original method is that $\sigma>0$0 is circulant. Applying $\sigma>0$1 reduces to convolution with a kernel and can be diagonalized by the discrete Fourier transform. For a gradient vector $\sigma>0$2, the smoothed direction can be computed by FFT and inverse FFT, and the reported complexity is $\sigma>0$3, described as negligible compared to one gradient step in modern ML models (Osher et al., 2018). In neural networks, the paper notes that gradients from each layer are smoothed layerwise, maintaining the physical structure (Osher et al., 2018).

The empirical results reported for LSGD cover quadratic test functions, MNIST logistic regression, deep neural networks, Wasserstein GANs, reinforcement learning, and non-convex toy problems. On quadratic objectives, LSGD and higher-order variants converge much faster and achieve substantially lower optimality gaps under gradient noise than GD or SGD, even with larger learning rates. On MNIST logistic regression with extreme mini-batches such as size $\sigma>0$4, standard SGD is described as struggling, while LSGD enables successful training and high test accuracy greater than $\sigma>0$5; histograms in the paper show a shift and narrowing of the accuracy distribution in favor of LSGD over SGD and SVRG baselines (Osher et al., 2018).

For deep neural networks, the same experiments report improved generalization accuracy and reduced experiment-to-experiment variability on ResNet and LeNet benchmarks. The paper further states that smoothing during early epochs is sufficient for benefit, with the example that only the first $\sigma>0$6 epochs are smoothed in ResNet56, and that the generalization gap is reduced (Osher et al., 2018). In Wasserstein GANs, training instability with large step sizes is reported to be resolved with LS-RMSProp; in reinforcement learning on Cartpole, LS-RMSProp reaches optimal policy more reliably and in fewer episodes. On Rosenbrock and related non-convex toy problems, LSGD is described as bypassing narrow or poor minima where GD would get stuck (Osher et al., 2018).

The sampling extension shows an analogous computational pattern. LS-SGLD also exploits the circulant structure of $\sigma>0$7 for FFT-based application of both $\sigma>0$8 and $\sigma>0$9, and its experiments are reported on posterior sampling, Bayesian logistic regression, and Bayesian convolutional neural networks (Wang et al., 2019).

6. Scope, distinctions, and related Laplacian weighting

A recurring misconception is to identify all Laplacian-based optimization with the same mechanism. The papers distinguish several meanings of “Laplacian weighting.” In LSGD and mLSGD, the gradient itself is premultiplied by the inverse of a one-dimensional discrete Laplacian-based matrix in parameter-space. In the basis-pursuit reformulation, Laplacian-like matrices appear in the objective and in the weighted least-squares subproblems rather than as a direct smoothing operator on neural-network gradients (Bonifaci, 2019). In LS-SGLD, both the drift and the injected Gaussian noise are transformed so that the stationary distribution is preserved (Wang et al., 2019).

A second distinction concerns where Laplacian effects enter the learning problem. In wide two-layer networks for mean squared error regression, gradient descent is shown to have an implicit bias toward minimizing a curvature regularization functional; in the univariate case, the learned solution is within $\mathbf{L}$0 of the interpolant whose difference from the initialization has the smallest $\mathbf{L}$1-norm of the second derivative weighted by a curvature penalty, and for constant penalty the solution is the natural cubic spline interpolation of the training data. In the multivariate case, the analogue involves the Radon transform of a fractional Laplacian, and constant penalty yields polyharmonic splines (Jin et al., 2020). This suggests that Laplacian weighting may arise either explicitly in the optimizer or implicitly in the function-space bias induced by gradient descent.

A third distinction concerns loss-level weighting. In Laplacian-Enhanced Gaussian Splatting, normalized second-order Laplacian responses are mapped into pixel-wise weights through nonlinear response-to-weight functions, and the weighted loss is used to improve structure-aware Gaussian optimization while keeping the original 3DGS rendering pipeline unchanged. The paper reports PSNR improvements of up to $\mathbf{L}$2 dB over 3DGS and up to $\mathbf{L}$3 dB over EGGS, with further gains when incorporated into FastGS and FasterGS (Guo et al., 6 Jun 2026). Although this is not LSGD, it exemplifies a broader pattern in which Laplacian-derived weights alter the descent direction indirectly by reweighting the objective.

In graph-based semi-supervised learning, properly weighting the graph Laplacian near labeled points is used to prevent degeneracy when the ratio of labeled to unlabeled data decreases. The proposed weight

$\mathbf{L}$4

yields a continuum variational problem whose smooth solution attains the labeled values continuously (Calder et al., 2018). This is not gradient smoothing in parameter-space, but it shows that Laplacian weighting can also serve as a stability mechanism in the regularization geometry of learning problems.

Taken together, these works support a broad but technically specific usage of the term: Laplacian-weighted Gradient Descent is most naturally centered on Laplacian smoothing of the update direction, yet it also has closely related forms in which Laplacian-derived operators or weights reshape the optimization landscape, the stochastic dynamics, or the effective regularization of the learned solution.