Solving 2D Black Scholes Equation via Hermitian Block Embedding and Generalised Quantum Signal Processing

Abstract: The Black Scholes equation provides a fundamental model for the no arbitrage pricing of financial derivatives. After finite difference discretisation, the pricing problem can be formulated as a finite dimensional linear algebra problem involving the inverse of a non Hermitian time step matrix. Recent advances in quantum linear algebra algorithms, particularly the generalised quantum signal processing (GQSP)algorithm, enable matrix functions to be implemented through polynomial transformations of a suitable unitary or Hermitian form. In this paper, we develop a Hermitian block embedding method that enables GQSP to be applied to the two dimensional Black Scholes equation. Numerical simulations for two asset European call options are performed to evaluate the proposed approach. GQSP based solutions are benchmarked against the classical polynomial approximation with backward Euler finite difference method, showing close agreement. This indicates that the Hermitian block embedding construction accurately captures the dynamics of the original non Hermitian operator. These results demonstrate the feasibility of combining Hermitian block embeddings with GQSP for multidimensional Black Scholes problems and provide a proof of principle for applying modern quantum linear algebra techniques to option pricing.

• The paper presents a Hermitian block embedding framework to enable efficient quantum simulation of the non-Hermitian 2D Black-Scholes PDE for multi-asset option pricing.
• It employs Generalised Quantum Signal Processing with odd polynomial approximations to accurately approximate the inverse time-step operator, validated against classical methods.
• Numerical simulations for basket call options show close agreement with classical backward Euler schemes, highlighting both the robustness and challenges of polynomial accuracy and discretisation stability.

Hermitian Block Embedding and GQSP for 2D Black-Scholes Option Pricing

Introduction

This work addresses the quantum solution of the two-dimensional Black-Scholes equation, a fundamental PDE in financial mathematics for multi-asset European option pricing. The paper introduces a Hermitian block embedding framework, permitting efficient application of Generalised Quantum Signal Processing (GQSP) to the inherently non-Hermitian discretised time-step operator arising in implicit finite-difference methods. This resolves the major obstacle in quantum linear algebra for multi-asset derivative pricing, extending prior one-dimensional approaches reliant on diagonal similarity transforms that break down in higher dimensions. Throughout, numerical validation is provided for basket call options, benchmarking quantum outputs against classical polynomial and backward Euler numerical schemes.

Black-Scholes Formulation and Discretisation

The two-dimensional Black-Scholes PDE describes the backward-in-time evolution of an option's value $V(S_1, S_2, t)$ dependent on two correlated assets with parameters $(\sigma_1, \sigma_2, r, \rho)$. After discretising asset domains using a lexicographic index and central stencils for first, second, and mixed derivatives, the backward Euler implicit scheme yields a linear system of the form

$$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$

with $\tilde{M}$ real, sparse, and significantly non-Hermitian due to drift and correlation terms. Unlike the tri-diagonal case in $1$D, the block structure cannot be Hermitianised via a diagonal similarity transform, as closed coupling loops challenge compatibility of the necessary ratio conditions.

Hermitian Block Embedding Construction

The paper introduces the Hermitian block embedding

$$H = \begin{pmatrix} 0 & \tilde{M} \ \tilde{M}^T & 0 \end{pmatrix}$$

whose spectrum is symmetric (eigenvalues $\pm s_j(\tilde{M})$). The inverse

$$H^{-1} = \begin{pmatrix} 0 & \tilde{M}^{-T} \ \tilde{M}^{-1} & 0 \end{pmatrix}$$

enables odd polynomial approximations to $1/x$ to target $\tilde{M}^{-1}$ in the relevant off-diagonal block, with application to embedded vectors of the form $(\sigma_1, \sigma_2, r, \rho)$0 yielding the evolved solution in its lower half.

Figure 1: Polynomial approximation of the inverse function.

This embedding allows recasting the inverse application as a matrix function synthesis problem suitable for Hermitian-GQSP, bypassing block-encoding and ancillary registers.

Odd Polynomial Approximation and Quantum Implementation

The Hermitian-GQSP protocol requires spectral scaling of $(\sigma_1, \sigma_2, r, \rho)$1:

$(\sigma_1, \sigma_2, r, \rho)$2

ensuring $(\sigma_1, \sigma_2, r, \rho)$3 (with $(\sigma_1, \sigma_2, r, \rho)$4 controlled by the spectral gap). The odd polynomial $(\sigma_1, \sigma_2, r, \rho)$5 approximates $(\sigma_1, \sigma_2, r, \rho)$6 on this domain, with higher polynomial degree dictated by the gap $(\sigma_1, \sigma_2, r, \rho)$7.

Figure 2: Comparison of the two-dimensional terminal payoff, classical backward Euler solution, polynomial approximation and the scaled quantum/GQSP output. The terminal payoff is given by the basket call payoff at maturity, while the classical and quantum surfaces show the result after one backward time step. The close agreement between the scaled quantum output, polynomial approximation and the classical backward Euler solution indicates that the GQSP block-embedding construction reproduces the intended inverse time-step action.

Implementation proceeds via construction of a unitary embedding $(\sigma_1, \sigma_2, r, \rho)$8, enabling GQSP circuits to enact polynomial transformations, ultimately extracting the evolved solution via post-selecting on the embedded register.

Numerical Simulations and Empirical Validation

Simulations with PennyLane and classical polynomial fitting demonstrate that the quantum GQSP approach achieves close agreement to classical backward Euler results, capturing the option surface after a backward time-step. The polynomial approximation errors are minimal except for spectral gaps near zero, where degree and numerical instability escalate.

Figure 3: Absolute error surfaces for the two-dimensional simulation. The left panel shows $(\sigma_1, \sigma_2, r, \rho)$9, comparing the polynomial approximation with the same classical benchmark. The right panel shows $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$0, comparing the scaled quantum/GQSP output with the classical backward Euler solution.

Slices fixing one asset and varying the other confirm the validity for each asset condition.

Figure 4: Fixed $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$1 solution comparison.

Figure 5: Classical fixed $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$2 solution comparison.

The main sources of error are: quantum noise, polynomial fit limitations as $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$3 (ill-conditioned cases), and discretisation artifacts from large time steps. Comparative analyses using small versus large time steps quantify the underlying discretisation error, which dominates in some parameter regimes.

Figure 6: Polynomial approximation of the inverse function.

Figure 7: Comparison of the two-dimensional terminal payoff, classical backward Euler solution, classical polynomial approximation and the scaled quantum/GQSP output.

Figure 8: Absolute error surfaces for the two-dimensional simulation.

Figure 9: Fixed $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$4 solution comparison.

Figure 10: Classical fixed $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$5 solution comparison.

Alternate parameter sets further corroborate robustness and parameter universality, but highlight the trade-off between polynomial accuracy and discretisation stability.

Implications and Future Directions

The Hermitian block embedding with GQSP enables a block-encoding-free extension from $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$6D to $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$7D option pricing in quantum linear algebra. The method demonstrates strong numerical fidelity for two-asset basket options and lays foundational methodology for tackling higher-dimensional derivatives, where classical computation is intractable. The primary limitations stem from spectral gap dependencies, degree of polynomial required, and stability of discretisation with large time steps.

Theoretical implications include the general applicability of Hermitian embedding for non-Hermitian sparse matrices and parity-preserving polynomial construction. Practically, the quantum resource cost for multi-step time evolution remains to be optimised, as odd-polynomial powers cannot replicate repeated time steps. There is considerable scope for improved polynomial fitting over disjoint spectral domains, adaptive discretisation, non-uniform grids, and integration into portfolio valuation.

Future developments may include algorithmic refinements for higher-order time evolution, resource scaling in $$\tilde{M} V^{n+1} = V^n, \qquad V^{n+1} = \tilde{M}^{-1} V^n$$8-asset settings, and advances in quantum circuit stability to diminish noise floors, which will be crucial for eventual financial applications.

Conclusion

This work rigorously develops and empirically validates a Hermitian block embedding approach for the quantum solution of the two-dimensional Black-Scholes equation via GQSP. The method circumvents structural limitations of prior approaches, ensuring direct access to the necessary inverse operation through parity-preserving polynomial approximations. Numerical evidence confirms quantum and classical solution consistency within discretisation and approximation limits. The framework is poised to generalise for multi-asset derivative pricing, with implications for quantum supremacy in high-dimensional financial computation (2606.00458).