The computation of Black-Scholes implied volatility is a fundamental task in quantitative finance, underpinning option valuation, model calibration and risk management. Although implied volatility is routinely used in practice, the inversion of the Black-Scholes pricing formula remains a challenging numerical problem, particularly in asymptotic regimes corresponding to extreme option prices, strikes or maturities, where the inverse map becomes highly sensitive to perturbations of the price. In this paper, we introduce a new family of asymptotically-informed neural-network architectures for implied-volatility computation. Exploiting the distinct behaviours of the Black-Scholes pricing function in different volatility regimes, we propose a family of architectures that learn a trainable partition of the price-log-moneyness domain through a system of gating functions and combines specialised local approximations of the implied-volatility function within each region. Extensive numerical experiments demonstrate that the proposed models consistently outperform standard feed-forward neural networks across a wide range of parameter domains, often by several orders of magnitude in relative accuracy while maintaining excellent generalisation properties. Furthermore, the neural-network outputs provide highly accurate initial guesses for a third-order Householder scheme, allowing near machine-precision implied-volatility computations after only two refinement iterations.
Figures & tables
Figure 1. (A): The partition of the space (B,C) into three regions defined by the tangent to the curve CBS(A,⋅) at its inflection point, for A=3 . (B): Inverse-volatility map Binv(A,)˙ for A=3 , and the corresponding partition of the (C,B) -domain into three regions defined by the tangent to the curve CBS(A,⋅) at its inflection point.
Figure 2. Partition of the (A,C) -domain induced by the inequalities ( 5 ). The central region corresponds to the approximately linear part of the pricing function, while the upper and lower regions correspond to the two asymptotic regimes.
Table 1. Families of gating functions f0 and f1 considered in this work.
Figure 3. Evolution of the training MSRE averaged over 25 independent runs for the selected architectures on the dataset (A,B)∈[0,16]×[10−5,7.07] . The shaded regions represent one standard deviation around the mean.
Figure 4. Pointwise relative error, on a logarithmic scale, of the chosen architectures over the domain (A,B)∈[0,16]×[10−5,7.07] . The evaluation grid contains 1000×1000 points.
Figure 5. Evolution of the logarithmic relative error ( 11 ) for a fixed value A=5 over the interval B∈[10−5,2.00] . The curves correspond to the neural-network output ( B0 ), after one Householder iteration, and after two Householder iterations. B is sampled uniformly over a grid of 500 points.
Model
Avg Log Error
Std Log Error
Max Log Error
SimpleExp
−33.56
7.38
−19.70
PolyACInvGenInter
−38.85
2.61
−35.35
PolyACInvGenFree
−38.85
2.60
−35.35
GaussACInvGenInter
−38.72
2.59
−35.35
PolyCInvExpInter
−38.70
2.62
−35.35
Table 3. Average, standard deviation and the maximum of the logarithm of the relative error ( 11 ) for the chosen architectures after 2 Householder iterations with A=5 and B∈[10−5,2.00] . B is sampled uniformly over a grid of 500 points.
Figure 6. Heat map of the logarithm of the relative error ( 11 ) of the chosen neural networks after 2 iterations of the Householder root finding method over the range of values (A,B)∈[0,16]×[10−5,7.07] , with 1000 grid values for each variable.
Figure 7. Comparison of the logarithm of the relative error ( 11 ) between Brent’s method (maximum of 50 iterations) and the proposed neural-network-plus-Householder approach on the domain (A,B)∈[0,0.5]×[10−5,2] . The evaluation grid contains 200×200 points.
Figure 8. Evolution of the training MSRE averaged over 25 independent runs for the selected architectures on the dataset (A,B)∈[0,3]×[10−7,1.22] . The shaded regions represent one standard deviation around the mean.
Figure 9. Pointwise relative error, on a logarithmic scale, of the chosen architectures over the domain (A,B)∈[0,3]×[10−7,1.22] . The evaluation grid contains 1000×1000 points.
Figure 10. Evolution of the logarithmic relative error ( 11 ) for a fixed value A=1.5 over the interval B∈[10−7,1.22] . The curves correspond to the neural-network output ( B0 ), after one Householder iteration, and after two Householder iterations. B is sampled uniformly over a grid of 500 points.
Model
Avg Log Error
Std Log Error
Max Log Error
SimpleExp
−34.39
6.00
−18.78
GaussACInvExpInter
−37.26
2.42
−34.65
GaussACInvGenInter
−37.25
2.43
−34.94
PolyACepsInvExpFree
−37.33
2.42
−34.79
PolyACSigExpFree
−37.13
2.34
−34.65
Table 5. Average, standard deviation and the maximum of the logarithm of the relative error ( 11 ) for the chosen architectures after 2 Householder iterations with A=1.5 and B∈[10−7,1.22] . B is sampled uniformly over a grid of 500 points.
Figure 11. Comparison of the logarithm of the relative error ( 11 ) between Brent’s method (maximum of 50 iterations) and the proposed neural-network-plus-Householder approach on the domain (A,B)∈[0,3]×[10−7,1.22] . The evaluation grid contains 200×200 points.
Figure 12. Evolution of the training MSRE averaged over 25 independent runs for the selected architectures on the dataset (A,B)∈[0,10−5]×[10−5,0.18] . The shaded regions represent one standard deviation around the mean.
Figure 13. Pointwise relative error, on a logarithmic scale, of the chosen architectures over the domain (A,B)∈[0,10−5]×[10−5,0.18] . The evaluation grid contains 500×500 points.
Figure 14. Evolution of the logarithmic relative error ( 11 ) for a fixed value A=5.10−6 over the interval B∈[10−5,0.18] . The curves correspond to the neural-network output ( B0 ), after one Householder iteration, and after two Householder iterations. B is sampled uniformly over a grid of 500 points.
Model
Avg Log Error
Std Log Error
Max Log Error
SimpleExp
−38.70
2.92
−10.948
PolyACepsInvExpInter
−39.06
2.88
−10.925
PolyCSigExpInter
−39.07
2.92
−10.926
GaussACInvExpInter
−39.40
2.85
−10.932
HardACInvExpInter
−38.62
2.91
−10.957
GaussCInvExpInter
−38.94
2.91
−10.953
Table 7. Average, standard deviation and the maximum of the logarithm of the relative error ( 11 ) for the chosen architectures after 2 Householder iterations with A=5.10−6 and B∈[10−5,0.18] . B is sampled uniformly over a grid of 500 points.
Figure 15. Comparison of the logarithm of the relative error ( 11 ) between Brent’s method (maximum of 50 iterations) and the proposed neural-network-plus-Householder approach on the domain (A,B)∈[0,10−5]×[10−5,0.18] . The evaluation grid contains 200×200 points.
Appendix figures & tables2 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 16. The logarithm of the ratio between the direct evaluation of the Black-Scholes formula ( 1 ) (denoted here Cdefault(h,t) ) and the Mills-ratio-based implementation ( 12 )(denoted Cexpanded(h,t) ), with varying h and fixed t . (A) t=0.001 , (B) t=0.01 , (C) t=0.1 and (D) t=1 . The magnitudes of the oscillations caused by subtractive cancellation decay as t increases.
Figure 17. The logarithm of the ratio between the direct evaluation of the Black-Scholes formula ( 1 ) (denoted here Cdefault(h,t) ) and the Mills-ratio-based implementation ( 12 )(denoted Cexpanded(h,t) ), with varying t and fixed h . (A) h=0.01 , (B) h=0.1 , (C) h=1 and (D) h=10 . The magnitudes of the oscillations decay as t increases. Larger values of h lead to an earlier disappearance of the oscillations.
Quantitative finance uses mathematical models, statistics, and machine learning to value assets, measure risk, and guide trading and portfolio decisions.
We present a novel application of Neural Networks with Local Converging Inputs (NNLCI) to improve the efficiency of existing numerical methods for pricing multi-asset options. The most concise input format for NNLCI has been introduced, offering substantial convenience and efficiency. NNLCI uses a neural network to locally correct solutions from a coarse mesh and a refined mesh (relative to the coarse one), requiring only a minimal amount of high-fidelity training data. We demonstrate this approach on cash-or-nothing options under the Black-Scholes equation in one, two, and three spatial dimensions, and on single-asset down-and-out barrier call options under the Heston stochastic-volatility model (whose pricing PDE is two-dimensional in the spot price S and the instantaneous variance v). In each case, NNLCI reduces the root-mean-square error (RMSE) of the refined-mesh numerical solution by a factor of approximately 4-12 on test sets, even when the neural network is trained on only a small subset of parameter combinations. These results demonstrate that NNLCI significantly reduces computational requirements for high-dimensional problems in real-time options trading and risk management, offering low training costs and strong generalization ability.
Harris Cobb, Wenbo Hao, Yingjie Liu
School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332
Modern option-learning systems operate in two coordinates: price space, where markets quote and no-arbitrage constraints are most naturally enforced, and implied volatility (IV) space, where volatility surfaces are smoothed, regularized, and evaluated. The bottleneck is interface, not approximation: Jäckel's seminal "Let's Be Rational" (LBR) solver already inverts the Black-Scholes price to machine precision efficiently. What is missing is a differentiable layer that preserves LBR in the forward pass and avoids backpropagating through its branch logic. Such a layer must also confront the unavoidable singularity of the inverse map in the low-vega regime, where the sensitivity 1/vega diverges as vega -> 0. We close this gap with PIVOT, the Price-Implied-Volatility Objective Translator. PIVOT keeps the LBR forward pass intact and supplies the backward pass by implicit differentiation through the smooth Black-Scholes/Black-76 price map, with an explicit gating contract: invalid domains return NaN, well-conditioned rows receive the exact 1/vega gradient, and low-vega rows are attenuated rather than silently regularized. On a single H100, a fused Triton kernel reaches 1.79e9 IV/s at machine precision (9.3e-14 max relative error vs. the reference C solver); end-to-end label generation sustains 48.9M/s on synthetic chains and 16.6M/s on SPX OptionMetrics. In a HyperIV-style one-day reproduction on SPX, PIVOT-augmented objectives Pareto-dominate the baselines, reducing held-out price MAE by up to 43.4% and the strongest three-seed gated objective improving price MAE by 38.8% and IV MAE by 21.3% jointly; cross-asset results on RUT, VIX, and NDX show directional price-MAE gains of 40.1%, 24.2%, and 16.7%, while an ungated IV-roundtrip control collapses to a degenerate near-zero surface, confirming the gate as a correctness contract rather than a tuning knob.
Neural network-based approaches have emerged as efficient alternatives to traditional optimization-based procedures for the calibration of stochastic volatility models. However, existing work has focused primarily on predictive accuracy, with comparatively little attention devoted to understanding the structure of the learned inverse calibration mappings. In this work, we analyze neural calibration mappings for the Heston and rough Heston models across multilayer perceptron, highway, and softmax-parametrized highway architectures, using complementary Shapley-based methods from explainable AI. Specifically, we consider SHAP and νSHAP explanations, which capture distinct, complementary notions of feature relevance, corresponding to sensitivity and sufficiency of feature subsets, respectively. Short maturities and smile wings consistently dominate parameter inference, and the dominant attribution structure remains qualitatively stable across architectures despite differences in predictive accuracy and parameter count. Parameter-specific differences between SHAP and νSHAP further reveal how distinct regions of the implied volatility surface contribute to parameter recovery and expose substantial redundancy in the calibration input. Building on this redundancy, we show that νSHAP explanations can guide a significant reduction in input dimensionality for the rough Heston model while matching calibration accuracy relative to the full implied volatility surface. These findings demonstrate that complementary Shapley-based methods provide structural insight into learned inverse calibration mappings beyond predictive error metrics, and offer a practical route to feature selection in neural calibration problems.
Shaïn Afzali, Serena Della Corte, Antonis Papapantoleon