A Symmetric Relative-Error Loss Function for Intermittent Multiscale Signal Modelling

1 Introduction

When performing regression tasks, regardless of context or modelling technique, the most common way to determine the fitness of a model is to calculate the Mean Squared Error (MSE). First described in [Legendre, 1806] as an algebraic procedure for fitting linear equations to data and later proven to be optimal for the lowest-sampling-variance estimator under unbiased assumptions [Henderson, 1975], it has persisted as the go-to loss function in the field of Machine Learning (ML) due to its context-free nature and, more importantly in the context of back-propagation, function smoothness. Unfortunately, this generality renders the MSE unsuitable for measuring model fitness when the sample space is not normally distributed or signal outliers are significant [Liano, 1996]. Such is the case of multiscale modelling problems, amongst which physical models such as the Navier-Stokes equation [Steinhauser, 2016] or intermittent signal forecasts like that of Photo-Voltaic (PV) energy generation can be mentioned.

To illustrate the problems associated with minimizing the MSE for multiscale spaces, a synthetic example is presented: let y be the target signal defined by

where each βi is a sample drawn from a single-trial Binomial distribution with a 0.15 probability of success. Two dynamical system estimators, namely ˆy(1) and ˆy(2) , are defined as

The first estimator introduces multiplicative noise to the reference signal. In contrast, the second one perfectly captures the peaks introduced by the intermittent signal β while holding the constant that minimizes the squared error in its absence. Simulating these conditions yields the curves shown in Figure 1a and MSEs of 0.215 and 0.107 for estimators 1 and 2 respectively, which empirically evidences how this metric mathematically neglects small-scale dynamics.

In ML regression, the traditional way to cope with these variations in scale is to fit the models to minimize the Mean Absolute Percentage Error (MAPE) [De Myttenaere et al., 2016], defined as

where ¯ δi := yi-ˆ yi yi is a term henceforth referred to as relative sample error.

The MAPE divides each sample's error by its target value, making errors at lower scales numerically comparable to those of larger reference magnitude. Figure 1b illustrates the equalization effect of this approach, with the MAPE of the

• (b) Squared and Relative Errors as perceived by the loss functions.

Figure 1: Synthetic multiscale example (100-sample window).

noisy signal now halving that of the small-scale-smoothing one (0.113 and 0.199 respectively). This effectively means a model trained to minimize the MAPE would better capture the behaviour of the signal. Unfortunately, the instability of this approach quickly becomes self-evident, since ∣∣¯ δ∣∣ → ∞ as |y| → 0 save for a perfect prediction, and even then stability cannot be guaranteed due to quantization.

Division by 0 is not uncommon in informatics and mathematical analysis, and some standard approaches have been proposed to handle this exception for regression applications [Hyndman and Koehler, 2006; Hodson et al., 2021]. However, [Kim and Kim, 2016] took an unconventional approach by interpreting ¯ δ as the hypotenuse slope of the "right triangle" defined by the perpendicular sides |y - ˆ y| and |y| to introduce Mean Arctangent Absolute Percentage Error (MAAPE), which measures the average angle of this slope as

Albeit a clever solution to the asymptotic behaviour at |yi| → 0+, choosing to preserve the absolute value of this ratio as a loss function introduces unwanted gradient behaviour at ˆ yi → yi. A smoother gradient can thus be attained by relaxing this one-to-one correspondence with the MAPE slope. Furthermore, by exploiting enforceable bounds, we can extend the definition to better capture the peak values of the signal without losing the small-scale dynamics.

This work modifies MAAPE and uses this extension as a functional unit to formulate the Symmetric Mean Average Squared Percentage Error (SMASPE). The specific contributions presented in this paper are the following:

Figure 2: Sample-wise behaviour of the MAAPE and its first two derivatives as a function of the error ratio.

1. The formulation of Mean Arctangent Squared Percentage Error (MASPE), a quadratic-relative-error variation of the MAAPE, and an analytical study of the mathematical properties that render it a more suitable loss function for gradient-based optimization.
2. SMASPE, a symmetrical extension of the MASPE that exploits heuristic signal bounds to allow models to capture multiscale dynamics, and a graphical study of its behaviour as a function of the predicted and reference signal values.
3. A numerical validation of SMASPE's adequacy for deterministic and stochastic optimization algorithms.

The remainder of this paper is organized as follows. In Section 2, the benefits and drawbacks introduced using the above ratio's square instead of its absolute value are explored as a function of the auxiliary variable ¯ δ. Section 3 presents the symmetric bound-enforcement and studies its implications as a function of the reference and predicted values y and ˆ y. In Section 4, two real-world examples are used to showcase the capabilities of SMASPE. Finally, in Section 5, a summary of the obtained experimental results and possible extensions to the proposed methodology are laid out.

The code and complementary material is made available as a GitHub repository at [Vanegas Arias, 2025].

2 Unidimensional Analysis: Improvements over MAAPE

Figure 2 illustrates the MAAPE and its derivatives as functions of ¯ δ, with the first derivative being multiplied by -1 to highlight the direction of minimization. As mentioned in Section 1, optimizing the absolute value of this ratio makes this function C 0 (not continuously differentiable) over the real domain, which means it presents asymptotic behaviour at the origin and thus renders the optimization problem numerically unstable. This is the same phenomenon encountered when using the Mean Absolute Error (MAE) as a loss function, which makes the MSE better suited as a minimization target even if the MAE is being used as the main performance metric.

By loosening the trigonometrical correspondence with the MAPE and focusing on the utility of MAAPE as a minimization target, the equivalent MASPE can be derived from Equation (5), yielding

Since the arctangent is monotonically increasing, the predictor ˆy that minimizes MASPE for a given reference y also minimizes MAAPE. However, unlike its absolute-valued counterpart, Equation (6) belongs to C∞ (proof in complementary material) over the real domain. Thus, its gradient lacks MAAPE's discontinuity at the origin and, consequently, its second derivative is no longer asymptotic at 0 error, as shown in Figure 3.

It is worth noting that the price the MAAPE (and by extension the MASPE) pays for its boundedness is the vanishing gradient at infinity. Indeed, in both Figures 2 and 3 it is evident how the loss functions converge to π/2 as ¯ δ → ±∞, making its derivative converge to 0 by consequence. However, as shown in Section 3, this issue can be mitigated by enforcing loss-function symmetry w.r.t. known signal bounds.

3 Bidimensional Analysis: Symmetric Bounds

If by minimizing the expected squared error the MSE neglects small-scale system dynamics, then the relative-error approaches reviewed so far suffer an equivalent loss of accuracy at larger scales. Referring back to the synthetic example presented in Equation (1), the estimator defined by Equation (2) would be perceived to be twice as accurate as the one detailed in Equation (3) by the ¯ δ-based loss functions (MASPE included) despite consistently failing to capture the peak values of the intermittent signal. Therefore, an extension of the metrics explored so far that considers the behaviour of the loss as a function of both y and ˆ y is required.

Even if a system's dynamics are unknown, its domain can often be either characterized or empirically determined from data and thus be exploited to equalize measurable outputs reliably. Looking at the MASPE as a bidimensional function of y and ˆ y, let [y-, y+] ⊂ R be a real, finite interval bounding the target signal; i.e., -∞< y- ≤ yi ≤ y+ <∞ ∀i. Then, loss symmetry is enforced w.r.t. these bounds as

thus defining the SMASPE. Since this modification is merely the sum of two smooth functions evaluated on shifted (inverted) versions of the original signals y and ˆ y, the differentiability is preserved.

The second term of the above equation inverts the signal w.r.t. y+ so that large-scale dynamics are perceived as smallscale ones and vice versa. This effectively yields a numerically stable loss function within the user-defined interval, as illustrated by Figure 4, preventing the vanishing error gradients towards the upper bound of the modelled signal.

The [y-, y+] interval can be loosened as

where ∆y := y+ - y- and γ is a relaxation coefficient. This mitigates SMASPE's asymptotic behaviour near the symmetry bounds, further increasing the robustness of the approach. However, if overdone, the function may lose track of the difference in scale over the signal range, forfeiting the advantages of the relative-error approach. As evidenced by the vector field intensity in Figure 4, the trade-off SMASPE makes for its stability and domain coverage is its reliance on carefully defined bounds, which should follow a user-defined heuristic.

4 Applications

This section showcases the performance and versatility of the SMASPE with two real-life datasets. In Section 4.1, a comparative study for the real-life application presented in [Kim and Kim, 2016], Stock-Keeping Unit (SKU)-demand forecasting, is performed to illustrate the characteristics of SMASPE w.r.t. the State-of-the-Art (SotA). Then, in Section 4.2, a Deep Neural Network (DNN) is trained to forecast PV-generation using the SMASPE to highlight the advantages of the proposed loss function in backpropagation. In both experiments, the SMASPE is tested using tight and loose bounds (γ = 0.1, sensitivity analysis in complementary material).

The experiments were implemented in Python using the Jax [Bradbury et al., 2018] and Keras 3 [Chollet and others, 2015] stacks for ML model construction and training.

4.1 BFGS Optimization: Stock-Level Next-Sample Forecasting

As the first experiment, the performance of the loss functions proposed in this paper was tested for the same use-case studied in [Kim and Kim, 2016]. The experiments in this section were carried out using a Kaggle dataset [Veera, 2020] which, to the authors' knowledge, is the only open-access dataset sampled with the same frequency (i.e., weekly) and containing approximately the same number of samples per sequence as the one used for the original MAAPE paper.

The dataset comprises 6000 time series, each containing 103 samples, representing the weekly sales volume per product across five retail stores, with one example from each of them being shown in Figure 5. Although not confirmed by the dataset description, it is reasonable to infer that the two sales peaks, approximately one year apart from each other, correspond to the sales near the holiday season, which explains the signal synchronization. Following the methodology of the MAAPE proposers [Kim and Kim, 2016], the first 95 samples were used to fit the models and the remaining 8 to evaluate their out-of-sample performance.

The considered models were the worst and best performers in [Kim and Kim, 2016]; namely, non-seasonal AutoRegressive Integrated Moving Average (ARIMA) and HoltWinters (HW) additive method (Triple Exponential Smoothing), respectively. Considering the length of the dataset, and since the authors did not provide the model parameters, a (4, 1, 4) order ARIMA and a HW with seasonality 4 were fitted separately to each time series, initializing the ARIMA model weights to 0 and the convex sum coefficients of the HW method to 0.5. These weights were adjusted by minimizing the MSE, MAAPE, MASPE, and SMASPE 1 . The performance of the fitted models was evaluated using all loss functions tested, the MAE, and their median convergence time (measured in iterations) in all sequences.

Table 1 contains the numerical results of the described experimental pipeline evaluated over the saturated model outputs, since negative demand values are not valid. From these results, the following can be stated:

• Despite the HWmodel converging generally faster than the ARIMA due to its lower parameter count (3 for HW vs. 10 for ARIMA), the previously worst-performing

1 Tight bounds set to each sequence's extrema.

(a) ARIMA model.

(b) HW model.

Figure 6: SKU demand model error distribution boxplots.

architecture ends up being the most accurate according to all error metrics when fitted using the SMASPE with moderately loose bounds. This, however, comes at an evident cost in convergence speed, requiring the most iterations out of all evaluated model-loss configurations.

• Figure 6 provides some insight on MAAPE's poor performance for training the HWmodels. Indeed, despite all loss functions producing a similar distribution box, the MAAPE yields significantly more outliers of larger magnitude than the MSE or the loose SMASPE, explaining the difference in order of magnitude for the metrics in Table 1. This phenomenon could remain hidden for a small sample size, which might explain why the results in this paper contradict those of [Kim and Kim, 2016].
• Finally, even when performing numerically similar (on average) to its loose counterpart while converging significantly faster, the worst-case scenarios in Figure 7 illustrate how tight bounds for the SMASPE often lead to unstable model weights. This pathological behaviour makes itself numerically evident through the MAE and MSE metrics for the HWmodel.

4.2 Backpropagation: PVGeneration 3-Day-Ahead Forecasting

The second dataset is the PV energy generation. It comprises 821 full days (after clean-up) of hourly-sampled floating-point data including timestamps, weather measurements, and average generated power in kilowatts. Furthermore, the output signal was normalized from 0 to ∼ 1 by dividing each daily curve over the median peak of the 15 preceding days to avoid numerically neglecting generation values during low-sunlight seasons such as winter and prolonged storms, thus setting a heuristic for the SMASPE bounds.

Missing samples in the dataset impose a natural division that was used to determine the training, validation, and test partitions, resulting in a 59.1/29.7/11.2% scheme. In addition to the historical PV generation signal, global radiation and air temperature measurements were used as exogenous features, the forecasts of which can be retrieved during deployment and used as additional inputs. Since only the historical measurements of weather variables were logged, forecast uncertainty was simulated by introducing additive noise with exponentially increasing variance; namely,

where σx denotes the variance of the original variable x, ς represents the standard deviation of the white noise for the first sample over the prediction horizon, and β sets the exponential growth factor of the noise variance as a function of the sample index over the prediction horizon. These parameters were calibrated to ensure that the noise had a standard deviation of 0.1% of each feature's standard deviation at the beginning of the forecast horizon and increased tenfold by the end of the forecast window (i.e. β = ln 10 71 ).

The chosen architecture was an encoder/decoder Recursive Neural Network (RNN) with Gated-Recurrent Unit (GRU) layers similar to the scheme proposed by [Cho, 2014] with its hyperparameters fixed (1737 weights in total) for 72-hour context and forecast windows. The weights were randomly initialized five times and used as common starting points for all loss functions considered in Section 4.1, keeping the weights that achieved the median validation loss. These were minimized using the Adam optimizer (Learning Rate of 5 × 10-4 , 64 sequences per batch) for a maximum of 2000 epochs (500 patience epochs). Finally, rather than the actual extrema of the training partition, the expected [0, 1] interval was used to define the tight bounds.

Table 2 contains the numerical results of the described experimental pipeline evaluated on the saturated model outputs, since negative generation values are not possible in this scenario. Figure 8 showcases the best and worst-case PVgeneration forecasts based on the IAE for the displayed pre-

• (a) Loss history for the training partition (unsaturated output signal to avoid gradient clipping).
• (b) Loss history for the validation partition (saturated output signal).

Figure 10: PV-generation model error distribution boxplots.

diction horizon. From these, the following can be stated:

• Despite the MSE and MAAPE converging the fastest out of all considered loss functions, Figure 9b reveals that this is due to the divergence of the validation loss. Since the MASPE history does not exhibit the same behaviour, it can be inferred that this issue stems frommodel overfitting and gradients around zero-relative-error respectively. This makes the models trained using the fast-converging loss functions perform significantly worse than the proposed ones for all error metrics evaluated over the test partition, as illustrated by Figure 8.
• On the opposite end of the spectrum, and as shown in Figure 9, MASPE models take the longest to converge due to the vanishing gradients brought up in Section 2.

Furthermore, Figure 10 reveals all proposed loss functions converge towards similar error distributions, yet the loose SMASPE remains the clear winner in terms of convergence speed and stability.

• The numerical results in Table 2 validate the increase in accuracy achieved by SMASPE, yielding relative improvements of 42% and 13% in terms of MAE w.r.t. MAAPE and MASPE-trained models respectively. Moreover, Figure 11 provides further insight into the signederror behaviour: while the MAAPE-trained model fails at capturing PV-generation peaks, the relaxed SMASPEtrained one achieves a near-symmetrical distribution at these critical points over the prediction horizon.

5 Conclusion

We have proposed the Mean Arctangent Squared Percentage Error (MASPE) as an improvement over MAAPE to prevent the asymptotic behaviour of the gradient at zero relative error. Moreover, we have proposed a novel loss function, the Symmetric Mean Average Squared Percentage Error (SMASPE), using the MASPE as a building block to exploit a priori knowledge on the bounds of the target signal and capture its dynamics across all numerical scales. An analysis of the mathematical properties that render these metrics more suitable as loss functions for gradient-based optimization was laid out, followed by a numerical validation of the theorized advantages for both deterministic and stochastic optimization applications. The results not only show that a SMASPE with properly calibrated bounds surpasses both the ubiquitous MSE and the SotA MAAPE for deterministic optimization in terms of accuracy, but it also significantly outperforms the latter when targetted by stochastic optimization algorithms (e.g., ADAM). As possible extensions, we propose multiple symmetry intervals for scalar output models and the generalization of the SMASPE to vectorial outputs, either using hypercubes for independent dimensions or hyperspheres for vectors with known magnitude bounds.

Acknowledgements

Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Executive Agency (REA). Neither the European Union nor the granting authority can be held responsible for them.

This programme has received funding from the European Union through Marie Skłodowska-Curie actions under project number 101081466: SEED - Systems and Engineering Science Doctorate.

Special thanks to Dr. Jose Rodrigo Rojo Garc´ıa, who selflessly lent his expertise during the formalization of the proof for the smoothness of the SMASPE.