Power-transform. Gaussion Smoothing, Homotopy for Optimizations, Stochastic Gradient Optimization

1 Introduction↩︎

In this work we aim to solve the following stochastic optimization problem \[\label{stochastic-objective} \max_{\boldsymbol{w}\in\mathbb{R}^d} G(\boldsymbol{w}):= \mathbb{E}_{\boldsymbol{x}\sim \mathcal{D}}[f_{\boldsymbol{w}}(\boldsymbol{x})],\tag{1}\] where \(\mathcal{D}\) is a known distribution, \(f\) is a known fitness function with accessible derivatives, and the objective \(G\) has a unique global maximum point \(\boldsymbol{w}^*\) and possibly multiple local maximums. If \(\boldsymbol{x}\) represents data and \(f_{\boldsymbol{w}}\) equals negative of the individual loss function associated with a model with \(\boldsymbol{w}\) being parameters, then (1 ) becomes the loss-minimization problem for model training in machine learning. The most popular methods for the loss-minimization problem are sample-based gradient methods, such as the full-batch gradient descent (GD), stochastic gradient ascent (SGD, [1]), and Adaptive Moment Estimation (Adam, [2]). While the full-batch GD is guaranteed to approach the stationary points of any sample estimate of \(G(\boldsymbol{w})\) under an appropriate learning schedule, it may get stuck near shallow local minima and saddle points. The randomness introduced by mini-batch selection helps SGD mitigate this issue and boosts performances. While the performance of SGD significantly depends on the learning rate schedule, Adam is less susceptible and converges much faster. However, these optimizers are still often trapped near local optimums. Hence, the optimization problem of locating the global optimum still remains open, which is crucial for training under-parameterized models.

A well-accepted explanation on SGD’s ability to escape from narrow local minimums is that the mini-batch selection randomness makes SGD actually optimize the smoothed objective \(\mathbb{E}_{\boldsymbol{w}}[G(\boldsymbol{w})]\) of \(G(\boldsymbol{w})\), which eliminates the narrow local minimums [3]. This relates to the smoothing technique ([4]) used by another type of popular derivative-free non-convex optimization methods—homotopy for optimization (e.g., [5], [6]). Its standard version solves the deterministic problem \(\min_{\boldsymbol{w}}h(\boldsymbol{w})\) by minimizing the Gaussian-smoothed surrogate \(H_{\sigma}(\boldsymbol{\mu}):=\mathbb{E}_{\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0}, I_d)}[h(\mu+\sigma \boldsymbol{\varepsilon})]\). Specifically, it incrementally decreases \(\sigma\) in an outer loop, and optimizes the surrogate in the inner loop under the current \(\sigma\) value, taking the solution found in the previous inner loop as the initial candidate for the current inner loop. Homotopy works because the smoothing technique removes sharp local minimums and the decreasing-\(\sigma\) mechanism drags the global optimum of the surrogate towards that of the original objective \(h(\boldsymbol{w})\) under certain conditions ([6]).

Traditional homotopy methods are often hindered by the inefficiency of their double-loop structure and a notorious sensitivity to the \(\sigma\)-decay schedule, which is notoriously difficult to tune. To address these challenges, recent research has pivoted toward more time-efficient single-loop mechanisms and strategies to minimize the displacement of the optimum caused by the smoothing process. For instance, the Single-Loop Gaussian Homotopy (SLGH) algorithm [7] introduces an automated \(\sigma\)-decay mechanism. GS-PowerOpt [8] utilizes a power transformation prior to Gaussian smoothing to effectively reduce the distance between the pre- and post-smoothed optima, while GS-PowerHP [9] integrates both the power transformation and a \(\sigma\)-decaying schedule. Despite these advancements, significant gaps remain for solving (1 ): SLGH lacks guarantees beyond convergence to a standard local optimum, while GS-PowerOpt and GS-PowerHP are restricted to zeroth-order, deterministic settings.

In this paper, we introduce Progressive Power Homotopy (Prog-PowerHP), a novel smoothing-based framework designed to solve the stochastic optimization problem defined in (1 ). It offers several critical advantages over existing optimization frameworks. Compared with the standard first-order methods of SGD and Adam and the homotopy method of SLGH, which are frequently susceptible to local optima, our method provides a rigorous theoretical pathway to approximating the global optimum under mild conditions (see Corollary 3). By utilizing a single-loop structure, it bypasses the excessive computational overhead and complex sub-problem convergence requirements of traditional double-loop homotopy. Compared with GS-PowerOpt and GS-PowerHP, our algorithm is designed for the stochastic optimization problem and adopts a gradually increasing power for better performance. A key innovation of our approach is the progressive power mechanism, which dynamically increases the power transformation as the smoothing parameter diminishes, effectively refining the search space and prioritizing deeper basins of attraction as the landscape’s fine-grained features emerge.
Contribution. The primary contributions of this work are three-fold. To the best of our knowledge, this is the first work to combine homotopy continuation with a progressive power-transformation of the objective for stochastic first-order optimization. Second, we establish the theoretical foundations of this approach by providing a comprehensive convergence analysis for Prog-PowerHP in the stochastic setting. Third, through extensive empirical evaluation on two-layer ReLU networks and phase retrieval, we demonstrate the algorithm’s superior landscape navigation capabilities. Our results show that in under-parameterized and limited-sample regimes—where standard optimizers consistently settle into higher-loss spurious minima—Prog-PowerHP achieves significantly higher approximation efficiency by locating deeper, more optimal local basins.
Related Work. There are existing studies that also exponentially transform the loss. For example, in tilted empirical risk minimization, this technique is used to serve multiple purposes, such as to decrease the influence of sample outliers and to improve model generalizations ([10], [11]). In the field of safe reinforcement learning, this technique is used to emphasize the worst-case scenarios ([12]–[14]). The main difference between these works and ours is, they treat the transformed objective as a risk-sensitive criterion, intentionally redefining optimality to emphasize tail performance, robustness, or safety. In contrast, our work uses a Gaussian-smoothed power-transformed objective purely as an optimization surrogate, designed to preserve the maximizer of the original fitness while improving concentration and optimization geometry.

Our method is also fundamentally distinct from existing homotopy-based approaches. Unlike gradient-based standard homotopy methods [6] and SLGH [7], our approach employs a power-transformation mechanism. Moreover, in contrast to GS-PowerOpt [8], GS-PowerHP GS-PowerHP?, and the path-learning homotopy method [15]—all of which are developed for deterministic optimization—our framework is tailored to stochastic first-order optimization, where gradient noise and limited sample access fundamentally alter the optimization dynamics.

2 Prog-PowerHP: The Proposed Method↩︎

2.1 Motivations↩︎

The design of our algorithm is motivated by a simple geometric principle: emphasize dominant extrema while suppressing spurious ones. Specifically, the power transformation amplifies higher local maxima, making them more prominent relative to lower ones. When combined with Gaussian smoothing, this effect is preserved for dominant peaks, while smaller, spurious local maxima are attenuated or eliminated by the smoothing operation. Specifically, we construct the following surrogate objective \[\label{surrogate} \max_{\boldsymbol{\mu}\in\mathbb{R}^d} F_{N,\sigma}(\boldsymbol{\mu}):=\mathbb{E}_{\boldsymbol{x}\sim \mathcal{D},\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2I_d)}[e^{Nf_{\boldsymbol{w}}(\boldsymbol{x})}].\tag{2}\]

During optimization, we adopt a progressive homotopy strategy. The power parameter \(N\) is gradually increased so as to align the global optimum of the smoothed objective with that of the original, un-smoothed problem. Simultaneously, the smoothing scale \(\sigma\) is progressively decreased to maintain sufficient curvature for fast solution improvement. The effects of these two mechanism are illustrated with a simple example in Figure 1. Specifically, (a) shows that a larger \(N\) aligns \(\boldsymbol{w}^*=\arg\max_{\boldsymbol{w}}G(\boldsymbol{w})\) and \(\boldsymbol{\mu}^*=\arg\max_{\boldsymbol{\mu}}F_{N,\sigma}(\boldsymbol{\mu})\) better, and (b) shows that a larger \(\sigma\) leads to larger slopes when \(\boldsymbol{\mu}\) is away from \(\boldsymbol{w}^*\) and a smaller \(\sigma\) leads to larger slopes when \(\boldsymbol{\mu}\) is near \(\boldsymbol{w}^*\).

2.2 Gradients for Solving (2 )↩︎

Suppose the first-order gradient of \(f\) is accessible. It can be derived that: \[\label{F-gradient} \nabla_{\boldsymbol{\mu}} F_{N,\sigma}(\boldsymbol{\mu}) = \mathbb{E}_{\boldsymbol{w}\sim \mathcal{N}(\boldsymbol{\mu}, \sigma^2I_d),\boldsymbol{x}\sim \mathcal{D}}\left[\frac{\partial e^{Nf_{\boldsymbol{w}}(\boldsymbol{x})}}{\partial \boldsymbol{w}}\right].\tag{3}\] Therefore, an un-biased sample estimate for \(\nabla F_{N,\sigma}(\boldsymbol{\mu})\) is \[\label{ygmkhnuc} \hat{\nabla} F_{N,\sigma}(\boldsymbol{\mu}) = \frac{1}{KJ} \sum_{k=1}^K\sum_{j=1}^J \left.\frac{ \partial e^{Nf_{\boldsymbol{w}}(\boldsymbol{x}_j)}}{\partial \boldsymbol{w}}\right|_{\boldsymbol{w}=\boldsymbol{w}_k},\tag{4}\] where \(\{\boldsymbol{x}_j\}_{j=1}^J\) is a batch of independently independently distributed (i.i.d.) samples drawn from \(\mathcal{D}\), and \(\{\boldsymbol{w}_k\}_{k=1}^K\) is a batch of i.i.d. samples from \(\mathcal{N}(\boldsymbol{\mu},\sigma^2I_d)\).

2.3 Our Algorithm↩︎

Our algorithm follows the step of stochastic gradient ascent, with an adaptive power \(N\) and an adaptive \(\sigma\). Specifically, the updating equations at step \(t\) are as follows. \[\begin{align} \label{update-rule} &N_t = N_{t-1} +\phi_t\Delta;\\ &\sigma_t = b + \sigma_{0} \beta^t;\\ &\boldsymbol{\mu}_{t} = \boldsymbol{\mu}_{t-1} + \alpha_t \hat{\nabla}F_{N_{t},\sigma_{t}}(\boldsymbol{\mu}_{t-1}),\\ \end{align}\tag{5}\] where \(\Delta\) denotes the pre-slected total increment of \(N\) in the updating process, \(\beta\in (0,1)\) is the decay factor for \(\sigma\), \(\phi_t\) is the increasing rate for \(N\), and \(\alpha_t\) is the learning rate for \(\boldsymbol{\mu}\). At the end of step \(T\), we output \(\boldsymbol{\mu}^* = \arg\max_{\boldsymbol{\mu} \in \{\boldsymbol{\mu}_1, \dots, \boldsymbol{\mu}_T\}} \hat{G}(\boldsymbol{\mu})\), where \(\hat{G}(\boldsymbol{\mu}):=\sum_{j=1}^{J}f_{\boldsymbol{\mu}_t}(\mathbf{x}_j)/J\) and \(\{\mathbf{x}_j\}\) is a set of validation samples drawn from \(\mathcal{D}\). The complete algorithm is in Algorithm 2.

3 Convergence Analysis↩︎

In this section, we prove that, under mild assumptions, Prog-PowerHP generates a sequence \(\{\boldsymbol{\mu}_t\}\) that converges in expectation to a point within an arbitrarily small neighborhood of the globally optimal solution \(\boldsymbol{w}^*\) of (1 ), as long as the total power increment \(\Delta\) is sufficiently large. The proof is divided to two parts. The first one (Section 3.1) shows that \(\{\boldsymbol{\mu}_t\}\) from (5 ) contains a candidate that converges in expectation to the stationary points of \(F_{N_0+\Delta,b}\) (see Remark 1), while the second part (Section 3.2) shows that, for any \(\delta>0\), there exists a threshold \(N_{\delta}\) such that for all \(N_0+\Delta>N_{\delta}\), every stationary point of \(F_{N_0+\Delta,b}\) lie within a \(\delta\)-neighborhood of \(\boldsymbol{w}^*\). These results are summarized in Corollary 3.

3.1 Convergence in Expectation to Stationary Points of \(F_{N_0+\Delta,b}\)↩︎

We derive the properties of \(F_{N,\sigma}\) in Lemma 1 – 3, which are used to derive a bound of \(\sum_{t=1}^{T}\alpha_t \mathbb{E}[\|\nabla F_{N_{t},\sigma_{t}}(\boldsymbol{\mu}_{t-1})\|^2]\) in Theorem 1. This bound is used by Corollary 1 to derive an iteration complexity for \(\lim_{T\rightarrow \infty}\min_{t\in\{1,...,T\}}\mathbb{E}[\|\nabla F_{N_{t},\sigma_{t}}(\boldsymbol{\mu}_{t-1})\|^2]=0\).

The assumption on the learning rate \(\alpha_t\) is the classic Robbins–Monro stochastic approximation conditions (e.g., [16]), and it appears in a wide range of convergence analyses for stochastic gradient–based algorithms. For example, [17] considers learning rates satisfying \(\lim_{t\rightarrow \infty} \alpha_t/t^{-\rho} \in (0,\infty)\) with \(\rho\in (0,1),\) which satisfies the Robbins–Monro conditions if \(\rho\in(1/2,1)\). Similarly, [18] established almost-sure convergence results in non-convex settings under the same step-size conditions.

The \(f\)-boundedness, \(\boldsymbol{x}\)-uniform Lipschitz continuity, and \(\boldsymbol{x}\)-uniform Lipschitz smoothness assumptions stated in Assumption 2 and Assumption 3 are also common in optimization for machine learning. Such regularity conditions are routinely adopted to control the sensitivity of the objective function with respect to the model parameters and to enable stability analyses. Representative examples include [7], [19], and [20]. In fact, Lipschitz continuity and Lipschitz smoothness are standard assumptions in theoretical works on non-convex optimization for machine learning models ([21], [22]).

3.2 Convergence to the Neighborhood of Global Optimum↩︎

The main task in this section is to show in Theorem 2 that, if \(N\) is greater than some threshold, all the stationary points of the surrogate \(F_{N,\sigma}(\boldsymbol{\mu})\) defined in (2 ) lie in a small neighborhood of \(\boldsymbol{w}^*:=\arg\max_{\boldsymbol{w}}G(\boldsymbol{w})\).

Intuitively, Part (1) assumes an alignment between average fitness and best-case fitness, namely that a parameter value \(\boldsymbol{w}\) which maximizes the expected fitness \(\mathbb{E}_{\boldsymbol{x}\sim\mathcal{D}}[f_{\boldsymbol{w}}(\boldsymbol{x})]\) also attains near-optimal best-case performance, measured by the essential supremum of \(f_{\boldsymbol{w}}(\boldsymbol{x})\) under \(\boldsymbol{x}\sim\mathcal{D}\). This assumption rules out pathological situations in which a model achieves a very high fitness on a vanishingly small subset of inputs while performing poorly on the bulk of the data.

Such an alignment is reasonable when the fitness \(f_{\boldsymbol{w}}(\boldsymbol{x})\) varies smoothly over the data domain and its distribution under \(\mathcal{D}\) is well behaved. For instance, if \(f_{\boldsymbol{w}}(\boldsymbol{x})\) is uniformly Lipschitz in \(\boldsymbol{x}\) and the data distribution assigns non-negligible probability mass to neighborhoods in the input space, then high fitness values cannot be achieved on isolated points without also improving the fitness on a set of inputs with positive probability, which in turn increases the expectation. Similarly, if the random variable \(f_{\boldsymbol{w}}(\boldsymbol{x})\) has light tails (e.g., is sub-Gaussian), then extreme values are tightly controlled by the mean, preventing sharp discrepancies between best-case and average performance.

This type of average-best alignment appears in intensively studied subfields of machine learning. In partial-label learning¹, for example, [26] show in their Theorem 1 that, under appropriate structural assumptions on label ambiguity, minimizing the expected true loss coincides with minimizing an infimum-based loss that captures the best achievable performance over admissible labels. While problem-specific, this result illustrates that, in structured learning settings, optimizing an average criterion can be equivalent to optimizing a best-case criterion.

More broadly, regularity assumptions such as Lipschitz continuity and light-tailed distributions are standard in optimization and generalization analyses in machine learning (e.g., [27], [28]). Under such conditions, it is natural to expect that models achieving superior average fitness also exhibit superior best-case behavior.

Part (3) in Assumption 4 assumes that, for each \(\boldsymbol{x}\), the near-optimal values of \(f_{\boldsymbol{w}}(\boldsymbol{x})\) are supported by a non-negligible fraction of the data distribution \(\mathcal{D}\), ruling out pathological cases where optimal fitness is achieved only on sets of vanishing probability. This belongs to the kind of assumptions that assume lower bounds on the probability of near-optimal sets, which are common among research works in machine learning (e.g., [29], [30]).

Below, we summarize our theoretical results.

4 Experiment↩︎

In this section, we evaluate our methods and focus on synthetic benchmarks—specifically phase retrieval and two-layer ReLU networks training—where the landscape geometry is well-characterized by existing theoretical literature ([31]–[34]). Standard benchmarks like MNIST or CIFAR-10 often obscure optimization failure modes due to the inherent redundancies of high-dimensional over-parameterization. By targeting limited-sample regimes in phase retrieval and under-parameterized settings in neural network training, we subject our algorithm to a rigorous stress test within landscapes known to harbor spurious local minima. This controlled setting more rigorously tests the effectiveness of power transformation and Gaussian smoothing by isolating their optimization behavior — something large-scale, generalization-oriented datasets cannot do.

We compare our performance with several widely used first-order optimizers, including the mini-batch SGD ([1]), SGD with Momentum ([35]), Adam ([2]), AdamW ([36]), and Sharpness-Aware Minimization (SAM, [37]). Additionally, we include the homotopy-based method of SLGHr ([7]) as a benchmark for ablation study.

4.1 Phase Retrieval↩︎

Given the unseen phase \(\boldsymbol{x}_0\in\mathbb{R}^d\) and measurements of the form \(y=\langle \boldsymbol{a}, \boldsymbol{x}_0 \rangle\), the phase retrieval (PR) problem seeks to recover \(\boldsymbol{x}_0\in\mathbb{R}^d\) from these observations. Beyond this classical formulation, modern phase retrieval research has increasingly incorporated sophisticated artificial intelligence paradigms, including deep generative priors ([38]), unsupervised teacher-student distillation protocols ([39]), and conditional generative models ([40]).

The classical phase retrieval problem is typically formulated as the following optimization task (e.g., [34]), which is notorious for exhibiting multiple spurious local minima: \[\label{PR-objective} \min_{\boldsymbol{x}\in\mathbb{R}^d} \mathcal{L}(\boldsymbol{x}):= \mathbb{E}_{\boldsymbol{a}\sim\mathcal{N}(\boldsymbol{0},I_d)} \left[( \langle \boldsymbol{a}, \boldsymbol{x} \rangle^2-\langle\boldsymbol{a},\boldsymbol{x}_0\rangle^2 )^2\right].\tag{6}\] In practice, a batch of observed signals, \(\{\boldsymbol{a}_n\}_{n=1}^N\) and \(\{y_n\}_{n=1}^N\) with \(y_n = \langle \boldsymbol{a}_n, \boldsymbol{x} \rangle\), are usually available. The sample-based objective is constructed as \[\label{PR-objective-est} \min_{\boldsymbol{x}\in\mathbb{R}^d} \hat{\mathcal{L}}(\boldsymbol{x}):=\frac{1}{N}\sum_{n=1}^N ( \langle \boldsymbol{a}_n, \boldsymbol{x} \rangle^2-y_n^2 )^2.\tag{7}\] In our simulated experiments, we use the following standard metric ([38], [41]) for measuring the quality of a given solution \(\boldsymbol{x}\): \[\label{PR-metric} \mathcal{M}(\boldsymbol{x}):= \min(\|\boldsymbol{x}-\boldsymbol{x}_0\|_2, \|\boldsymbol{x}+\boldsymbol{x}_0\|_2)/\|\boldsymbol{x}_0\|_2,\tag{8}\] where \(\|\cdot\|_2\) denotes the \(L_2\)-norm.

Under each set of selected values of \(d\) and \(N\), we perform 100 experiments. In each experiment, we randomly generate a true solution \(\boldsymbol{x}_0\) from \(\mathcal{N}(\boldsymbol{0},0.1I_d)\), an initial solution candidate \(\boldsymbol{x}^{(0)}\) from \(\mathcal{N}(\boldsymbol{0},I_d/\sqrt{d})\), and \(N\) independent samples \(\{\boldsymbol{a}_n\}_{n=1}^N\) from \(\mathcal{N}(\boldsymbol{0},I_d)\). Then, with these values, we apply each compared optimization method to solve (7 ), where each solution update uses a mini-batch of \(\{\boldsymbol{a}_n\}_{n=1}^N\). Let \(\{\boldsymbol{x}^{(t)}\}_{t=1}^T\) denote the solution candidates produced by the performed. If \(\min_{t} \mathcal{M}(\boldsymbol{x}^{(t)})\leq 0.001\), we consider that the true solution \(\boldsymbol{x}_0\) is recovered and call the method to be successful for this experiment. Define \(\boldsymbol{x}_i^*=\arg\min_{t}\mathcal{M}(\boldsymbol{x}^{(t)})\), where \(i\) denotes the experiment index. Average results obtained over the 100 experiments are reported in Table 1, where we also report the two-sided \(t\)-test result to check if our algorithm produces a significantly lower \(\mathcal{M}(\boldsymbol{x}^*)\) than other algorithms.

It is well-known that the PR problem exhibits a transition into difficulty when the measurement-to-dimension ratio \(N/d\) decreases. In particular, prior work shows that the recovery is impossible when \(N/d < 1\), while above this threshold spectral and first-order methods begin to succeed [33]. However, in the regime where \(N/d\in [2,4]\), recovery is theoretically possible but empirically challenging, as the non-convex loss landscape is poorly conditioned [32], [34]. Motivated by these observations, we focus our experiments on this information-scarce regime by fixing \(N=400\) and varying \(d\) so that \(N/d \in [2,4]\), where algorithmic performance is most sensitive to optimization strategy.

The experiment results in Table 1 show that our algorithm outperforms all the compared algorithms. The \(t\) statistics indicate that the results are statistically significant.

Our results align with the established theory of phase retrieval: when \(N/d\) is small (e.g, \(N/d\in [2,4)\)), the optimization geometry is fragile and recovery is sensitive to the choice of algorithm. Beyond this regime, the non-convex loss landscape becomes benign—characterized by the absence of spurious local minima and the presence of negative curvature at all saddle points ([42]), and hence, the performance differences between various optimization methods largely vanish.

4.2 Two-layer Neural Network Training↩︎

We evaluate our algorithm using the experimental setup for training two-layer neural networks (without bias terms) as described in [31]. Specifically, given a set of orthonormal ground-truth weight vectors \(\{\boldsymbol{v}_i\}_{i=1}^k\), the objective is to solve \(\min_{\boldsymbol{W}\in\mathbb{R}^{d\times n}} \mathcal{L}(\boldsymbol{W})\), where \[\label{twolayer-loss} \mathcal{L}(\boldsymbol{W})=\mathbb{E}_{\boldsymbol{x}\sim \mathcal{N}(\boldsymbol{0},I_k)}\left[\frac{1}{2}\left(\sum_{j=1}^n[\boldsymbol{w}_j'\boldsymbol{x}]_+ - \sum_{i=1}^k[\boldsymbol{v}_i'\boldsymbol{x}]_+\right)^2\right],\tag{9}\] \(\boldsymbol{W}:=[\boldsymbol{w}_1,...,\boldsymbol{w}_n]\) represents the learnable weights and \([z]_+=\max(0,z)\) is the ReLU activation. [31] demonstrate that locating the global optimum is particularly challenging in the exactly-parameterized case \(k=n\in[16,20]\), especially for \(n\in[16,20]\), due to the prevalence of spurious local minima. While over-parameterization (\(n>k\)) is known to significantly reduce this difficulty, we choose to address the more challenging regimes. We evaluate our algorithms across four settings: the exactly-parameterized cases where \(k=n=20\) and \(k=n=25\), as well as the under-parameterized cases where \((k,n)=(20,19)\) and \((k,n)=(25,24)\).

With each pair of selected values for \((k,n)\), we perform 100 experiments on each of the compared optimization algorithms. In each experiment, we randomly generate a set \(\{\boldsymbol{v}_i\}_{i=1}^k\) of orthonormal vectors in \(\mathbb{R}^k\), and sample a batch of 30 \(\boldsymbol{x}\)’s from \(\mathcal{N}(\boldsymbol{0},I_k)\) for each update of the algorithm’s solution. All the algorithms share the same initial solution candidate. Following [31], we use the test error as the evaluation metric, which is computed with a set of 5,000 test samples of \(\boldsymbol{x}\): \[\label{test-error} e(\boldsymbol{W}) = \frac{1}{5000} \sum_{l=1}^{5000}\frac{1}{2} \left(\sum_{j=1}^n[\boldsymbol{w}_j'\boldsymbol{x}_l]_+ - \sum_{i=1}^k[\boldsymbol{v}_i'\boldsymbol{x}_l]_+\right)^2.\tag{10}\] For each algorithm, we report statistics about the optimal solution \(\boldsymbol{W}^*:=\arg\min_{\tau=1}^T e(\boldsymbol{W}_\tau)\) found in the solution-update process. The pre-selected total number \(T\) of iterations is set as \(35,000\) each method (we choose a population size of 1 for our method). Again, we consider an algorithm successful (i.e., the global minimum is found) if \(e(\boldsymbol{W}^*)\leq 0.001\).

The results are summarized in Table 2. While our algorithm exhibits comparable or slightly improved performance relative to most baseline methods in the exactly-parameterized regime (\(n=k\)), these gains are generally not statistically significant. In contrast, in under-parameterized settings (\(n<k\)), our method demonstrates a clear and consistent advantage, achieving substantially lower loss values with statistically significant margins. These results highlight the superior landscape navigation capabilities of our approach; specifically, our algorithm consistently identifies deeper local minima in scenarios where capacity constraints preclude a zero-loss solution. This empirical behavior aligns with our theoretical framework, which suggests that the power transformation effectively reweights the optimization landscape to favor more optimal local basins.

4.3 Explanation of Results↩︎

We hypothesize that the consistent outperformance of our approach stems from a dual-mechanism synergy between the power transform and landscape smoothing. Specifically, the power transformation reshapes the remaining topography to reweight the ‘attractiveness’ of various basins, allowing the iterates to prioritize deeper, more optimal valleys. Complementing this, the smoothing component serves to regularize the optimization surface by effectively marginalizing shallow, spurious local minima that typically trap first-order methods in data-sparse or under-parameterized regimes. Together, these mechanisms enable our algorithm to bypass high-loss noise and converge toward superior local basins that are often inaccessible to standard optimizers. Interestingly, the convergence to identical success rates in the exactly-parameterized case (\(n=k\)) suggests that while this dual strategy is highly proficient at navigating ‘cluttered’ environments, it still encounters the fundamental symmetry-breaking hurdles inherent in perfect recovery tasks. Ultimately, these findings position our method as a robust solution for non-convex problems where model capacity is constrained or sampling is limited.

5 Conclusion↩︎

In this work, we propose a novel global optimization framework that integrates a progressive power-transformation mechanism within a homotopy scheme, and we establish its convergence to the global optimum with provable iteration complexity. We evaluate the proposed method on constrained nonconvex landscapes arising from limited data availability—such as phase retrieval in the regime (\(n/d \in [2,4]\))—and limited model capacity, as in under-parameterized neural networks with (\(n < k\)).

In the under-parameterized neural network setting, where exact recovery is mathematically impossible, our method consistently converges to solutions with significantly lower loss than other first-order optimizers. A similar advantage is observed in phase retrieval, particularly near the information-theoretic limit (n/d = 2) [43]. Together, these results suggest that the proposed combination of power transformation and Gaussian smoothing provides an effective mechanism for navigating degenerate non-convex landscapes, with potential implications for reducing model size and computational demands in large-scale deep learning systems.

6 Appendix↩︎

6.1 Single Loop Gaussian Homotopy (SLGH)↩︎

[7] proposes two single-loop Gaussian homotopy algorithms that applies \(\nabla f(\boldsymbol{x})\), SLGHd and SLGHr. These two algorithms solves the surrogate problem of \[\label{slgh-goal} \min_{\boldsymbol{\mu}\in\mathbb{R}^d,\sigma\in\mathbb{R}^+} F(\boldsymbol{\mu},\sigma),\tag{11}\] where \[\begin{align} F(\boldsymbol{\mu},\sigma) &= \mathbb{E}_{\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2 I_d)}\left[\mathbb{E}_{\boldsymbol{x}\sim \mathcal{D}}[f_{\boldsymbol{w}}(\boldsymbol{x})]\right] = \mathbb{E}_{\boldsymbol{x}\sim \mathcal{D}}[ \mathbb{E}_{\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2 I_d)}[f_{\boldsymbol{w}}(\boldsymbol{x})]]. \end{align}\]

We let \[\begin{align} \bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})&= \mathbb{E}_{\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2 I_d)}[f_{\boldsymbol{w}}(\boldsymbol{x})]\\ &= \mathbb{E}_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)} [f_{\boldsymbol{\mu}+\sigma \boldsymbol{\epsilon}}(\boldsymbol{x})]\\ &= \frac{1}{\sqrt{2d}}\int_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)} f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x}) e^{-\| \boldsymbol{\epsilon} \|^2} d\boldsymbol{\epsilon}. \end{align}\] Under the assumption of interchangeable differentiation and expectation, we have \[\label{slgh-derivative} \begin{align} \nabla_{\boldsymbol{\mu}}\bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})&=\frac{1}{\sqrt{2d}}\int_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)} \nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x}) e^{-\| \boldsymbol{\epsilon} \|^2} d\boldsymbol{\epsilon} \\ & = \mathbb{E}_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)}[\nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x})] \in\mathbb{R}^d. \end{align}\tag{12}\] \[\begin{align} \nabla_{\sigma}\bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})&= \frac{1}{\sqrt{2d}}\int_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)} \nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x})\cdot \boldsymbol{\epsilon} e^{-\| \boldsymbol{\epsilon} \|^2} d\boldsymbol{\epsilon}\\ &= \mathbb{E}_{\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)}[\nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x})\cdot\boldsymbol{\epsilon}]\in\mathbb{R}, \end{align}\] where \(\nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x}):= \left.\frac{\partial f_{\boldsymbol{w}}(\boldsymbol{x})}{\partial\boldsymbol{w}}\right|_{\boldsymbol{w}=\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}\in\mathbb{R}^d\), and \(\cdot\) denotes the dot product between vectors. These two results imply that \[\begin{align} \nabla_{\boldsymbol{\mu}} F(\boldsymbol{\mu},\sigma) = \mathbb{E}_{\boldsymbol{x}\sim \mathcal{D}}[\nabla_{\boldsymbol{\mu}}\bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})] = \mathbb{E}_{\boldsymbol{x}\sim \mathcal{D},\boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)}[\nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x})] \in\mathbb{R}^d.\\ \nabla_{\sigma} F(\boldsymbol{\mu},\sigma) = \mathbb{E}_{\boldsymbol{x}\sim \mathcal{D}}[\nabla_{\sigma}\bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})] = \mathbb{E}_{\boldsymbol{x}\sim\mathcal{D}, \boldsymbol{\epsilon}\sim \mathcal{N}(\boldsymbol{0},I_d)}[\nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}}(\boldsymbol{x})\cdot\boldsymbol{\epsilon}]\in\mathbb{R}. \end{align}\]

[7] only states that \(\nabla_{\boldsymbol{\mu}} F(\boldsymbol{\mu},\sigma)\) and \(\nabla_{\sigma} F(\boldsymbol{\mu},\sigma)\) can be estimated by \(\nabla_{\boldsymbol{\mu}} \bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})\) and \(\nabla_{\sigma} \bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x})\), using a fixed sample of \(\boldsymbol{x}\). But they did not specify how the two estimators are computed or estimated given \(\boldsymbol{x}\). Hence, it seems that they assume the two estimators can be computed exactly. Below, we provide a natural formula for approximating the two estimators.

\[\begin{align} & \hat{\nabla}_{\boldsymbol{\mu}}\bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x}) = \frac{1}{K}\sum_{k=1}^K \nabla f_{\boldsymbol{\mu}+\sigma \boldsymbol{\epsilon}_k}(\boldsymbol{x})\\ & \hat{\nabla}_{\sigma} \bar{F}(\boldsymbol{\mu},\sigma;\boldsymbol{x}) = \frac{1}{K}\sum_{k=1}^K \nabla f_{\boldsymbol{\mu}+\sigma\boldsymbol{\epsilon}_k}(\boldsymbol{x})\cdot \boldsymbol{\epsilon}_k\\ &\hat{\nabla}_{\boldsymbol{\mu}} F(\boldsymbol{\mu},\sigma) = \frac{1}{KJ}\sum_{j=1}^J \sum_{k=1}^K \nabla f_{\boldsymbol{\mu}+\sigma \boldsymbol{\epsilon}_k}(\boldsymbol{x}_j)\\ &\hat{\nabla}_{\boldsymbol{\mu}} F(\boldsymbol{\mu},\sigma) = \frac{1}{KJ}\sum_{j=1}^J \sum_{k=1}^K \nabla f_{\boldsymbol{\mu}+\sigma \boldsymbol{\epsilon}_k}(\boldsymbol{x}_j)\cdot\boldsymbol{\epsilon}_k \end{align}\]

6.2 SLGHd↩︎

At every time step \(t\), the updating equation is \[\begin{align} \label{SLGH-d} \boldsymbol{\mu}_{t+1} &= \boldsymbol{\mu}_{t} + \beta \hat{\nabla}_{\boldsymbol{\mu}} F(\boldsymbol{\mu}_t,\sigma_t);\\ \sigma_{t+1} &= \max\{\min\{\sigma_t-\eta \hat{\nabla}_{\sigma} F(\boldsymbol{\mu}_t,\sigma_t),\gamma \sigma_t\}, \epsilon'\}. \end{align}\tag{13}\]

6.3 SLGHr↩︎

At every time step \(t\), the updating equation is \[\begin{align} \label{qkcyeght} \boldsymbol{\mu}_{t+1} &= \boldsymbol{\mu}_{t} + \beta \hat{\nabla}_{\boldsymbol{\mu}} F(\boldsymbol{\mu}_t,\sigma_t);\\ \sigma_{t+1} &= \gamma \sigma_t. \end{align}\tag{14}\]

6.4 Proof to Theoretical Results↩︎

6.4.1 Proof to Lemma 1↩︎

Lemma 1 Under Assumption 2, given any \(N>0\) and \(\sigma>0\), both \(F_{N,\sigma}(\boldsymbol{\mu})\) and \(\nabla_{\boldsymbol{\mu}}F_{N,\sigma}(\boldsymbol{\mu})\) are well-defined for all \(\boldsymbol{\mu}\in\mathbb{R}^d\).

6.4.2 Proof to Lemma 2↩︎

Lemma 2 Under Assumption 2 and 3, for any \(\boldsymbol{\mu}_1\), \(\boldsymbol{\mu}_2\in\mathbb{R}^d\), \(\| \nabla F_{N,\sigma}(\boldsymbol{\mu}_1) - \nabla F_{N,\sigma}(\boldsymbol{\mu}_2) \|\leq L_F\|\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\|\), where \(L_F=e^{NM_f}N(NL_{0f}^2+L_{1f})\).

6.4.3 Proof to Lemma 3↩︎

Lemma 3 Under Assumption 2, \(\| \hat{\nabla} F_{N,\sigma}(\boldsymbol{\mu}) \|^2 \leq N^2e^{2NM_f} L_{0f}^2\).

6.4.4 Proof to Theorem 1↩︎

Theorem 1. Let \(\{\boldsymbol{\mu}_t\}_{t=1}^T\) be the series produced by (5 ) with a pre-selected deterministic \(\boldsymbol{\mu}_0\). Under Assumption 2 and 3, we have \[\sum_{t=1}^{T}\alpha_t \mathbb{E}[\|\nabla F_{N_{t},\sigma_{t}}(\boldsymbol{\mu}_{t-1})\|^2] \leq e^{(N_0+\Delta)M_f} -F_{N_{0},b+\sigma_{0}}(\boldsymbol{\mu}_{0}) + H_1 \sum_{t=0}^{T-1}\alpha_t^2 +H_2d \sum_{t=0}^{T-1}\beta^t+H_3 \sum_{t=1}^{T} \phi_{t},\] where \(H_1=(N_0+\Delta)^3 e^{3(N_0+\Delta)M_f}((N_0+\Delta)L^2_{0f}+L_{1f})L_{0f}^2\), \(H_2 = e^{(N_0+\Delta)M_f}(N_0+\Delta)L_{0f}\sigma_0(1-\beta)\sqrt{2}\), \(H_3:=\max\{ (N_0e)^{-1}, e^{(N_0+\Delta)M_f} M_f \}\Delta\).

6.5 Proof to Corollary 1↩︎

Corollary 1. Let \(\{\boldsymbol{\mu}_t\}_{t=1}^T\) be the series produced by (5 ) with a pre-selected deterministic \(\boldsymbol{\mu}_0\) and a learning rate \(\alpha_t=t^{-(1/2+\gamma)}\) where \(\gamma\in (0,1/2)\) and \(T\geq 2\). Under Assumption 1, 2, and 3, for any \(\varepsilon>0\), whenever \(T> (\frac{C_{0,\Delta}(1-2\gamma)}{2-2^{1/2+\gamma}}d\varepsilon^{-1})^{2/(1-2\gamma)}=O_{\Delta}((d \varepsilon^{-1})^{2/(1-2\gamma)})\), we have that \[\min_{t\in\{\lfloor T/2\rfloor,...,T\}} \mathbb{E}[\|\nabla F_{N_{t},\sigma_{t}}(\boldsymbol{\mu}_{t-1})\|^2] \leq \varepsilon.\] Here, \(C_{0,\Delta}:=e^{(N_0+\Delta)M_f} - F_{N_{0},b+\sigma_{0}}(\boldsymbol{\mu}_{0}) + H_1 \sum_{t=1}^{\infty}t^{-(1/2+\gamma)} +H_2(1-\beta)^{-1}+H_3\), \(C_{1,\Delta}= \max\{1, 2/C_{0,\Delta}\}\), \(\lfloor T/2 \rfloor\) denotes the largest integer no greater than \(T/2\), and \(\{H_1,H_2,H_3\}\) are as those defined in Theorem 1.

6.6 Proof for Corollary 2↩︎

Corollary 2 The dependence of the \(O\)-term on \(\Delta\) can be removed if we further assume that \(M_f\leq 0\).

6.7 Proof to Lemma 4↩︎

Lemma 4 Under Assumption 2 and 4, for any \(\delta>0\), we have (1) \(\lim_{N\rightarrow +\infty}s_N/v_N=0\), where \(s_N:=\sup_{\|\boldsymbol{w}-\boldsymbol{w}^*\|\geq\delta}G_N(\boldsymbol{w})\), \(v_N:=G_N(\boldsymbol{w}^*)e^{-N\eta/4}\), and \(\eta:= \Psi(\boldsymbol{w}^*)-\sup_{\|\boldsymbol{w}-\boldsymbol{w}^*\|\geq \delta} \Psi(\boldsymbol{w})\), and (2) there exists \(\delta'\in(0,\delta)\) and \(N_{\delta'}\) such that \(G_N(\boldsymbol{w})>v_N\) whenever \(\|\boldsymbol{w}-\boldsymbol{w}^*\|<\delta'\) and \(N>\frac{8}{\eta}\log c_{\epsilon}^{-1}\). Here, \(\epsilon=\eta/16\) and \(c_{\epsilon}\) is specified in Assumption 4.

6.7.1 Proof to Theorem 2↩︎

Theorem 2 Suppose Assumption 2 and 4 hold. Denote \(\boldsymbol{\mu}=[\mu_1,...,\mu_d]\) and \(\boldsymbol{w}^*=[w_1^*,...,w_d^*]\). For any \(\delta>0\), \(\sigma>0\), and \(M>0\), there exists a threshold \(N_{\delta,\sigma,M}>0\) such that, whenever \(N>N_{\delta,\sigma,M}\) and \(\|\boldsymbol{\mu}\|\leq M\), for all \(i\in\{1,2,...,d\}\) we have \(\frac{\partial F_{N,\sigma}(\boldsymbol{\mu})}{\partial \mu_i}<0\) if \(\mu_i>w_i^*+\delta\), and \(\frac{\partial F_{N,\sigma}(\boldsymbol{\mu})}{\partial \mu_i}>0\) if \(\mu_i<w_i^*-\delta\).

References↩︎