Quantum Monte Carlo (QMC) Simulation
Installation Required: This functionality requires MATLAB Support Package for Quantum Computing.
This example shows how to use Quantum Monte Carlo (QMC) simulation in MATLAB® to compute the mean of a function of a random variable. Running this example requires Statistics and Machine Learning Toolbox™ as well as MATLAB Support Package for Quantum Computing.
There are a broad range of tasks in finance and economics that depend on Monte Carlo simulation, from option pricing to macroeconomic stress testing. While this example does not explore computational efficiency, research shows that QMC offers a quadratic speed-up compared to classic Monte Carlo methods.
This example follows  in considering an application of QMC to compute the mean of a function of a random variable. Specifically, the mean of a trigonometric function of an underlying normal distribution is computed. This calculation generalizes to many practical applications. For example, if the probability distribution represents the price of an underlying asset, the function could be the price of an option on that asset.
Assume there is a random variable x that is normally distributed, and the function to be evaluated is
The expected value of f(x) is given analytically by
Calculate the value of μ.
AnalyticMean = sinh(1)/exp(1)
AnalyticMean = 0.4323
Classic Monte Carlo
For comparison, compute the mean μ with classic Monte Carlo. The steps to do this are:
Generate a sample of points from the underlying distribution of x
Compute the function value at each point in the sample.
Compute the sample mean of the function values.
func = @(x) sin(x).^2; numSamples = 1000; sampleData = randn(numSamples,1); funcValues = func(sampleData); MCMean = mean(funcValues)
MCMean = 0.4210
Plot the sample, function values, and calculated mean.
tiledlayout(1,2) nexttile h1 = histogram(sampleData); title("Sample Points") nexttile h2 = histogram(funcValues); hold on xline(MCMean,'--','MCMean'); hold off title("Classic Monte Carlo","Mean of Function Values")
First, define some parameters for the quantum simulation. Define the number of qubits for the probability distribution register as 5, and the number of qubits for the estimation register as 6. Then, discretize the probability distribution and function over a grid and compute the discrete mean of the function.
m = 5; % Number of probability qubits n = 6; % Number of estimation qubits M = 2^m; N = 2^n; probQubits = 1:m; x_max = pi; x = linspace(-x_max,x_max,M)'; p = normpdf(x); p = p./sum(p); DiscreteMean = sum(func(x).*p)
DiscreteMean = 0.4326
The QMC methodology covered in  consists of the steps
Load the probability distribution onto a quantum register of
Encode the value of the random variable onto a value qubit.
Use amplitude estimation, with a register of
nqubits, to estimate the probability of the value qubit being in the state .
While there are other approaches for amplitude estimation (such as iterative amplitude estimation), this example uses quantum phase estimation.
The circuit diagram below summarizes the approach. The F block loads the probability distribution and encodes the random variable onto a value qubit. The Q blocks and the inverse Quantum Fourier Transform (QFT) are used for quantum phase estimation.
1. Load Probability Distribution
Load the probability distribution onto
m qubits. While  uses a quantum
variational circuit to approximate the probability distribution, this example loads the
probabilities directly onto the circuit. This requires the use of some custom gate
functions, which you can save as functions in the current folder:
A = initStateToProbabilities(p); A.Name = 'Ainit';
Simulate the register to verify that the probability distribution loaded correctly.
circ = quantumCircuit(A); state = simulate(circ); histogram(state,probQubits) title("Simulated Probabilities")
2. Encode Random Variable
Encode the discretized random variable function onto a value qubit.
valueQubit = m + 1; % Value qubit index estQubits = (valueQubit + 1):(n + m + 1); % Estimation qubits index anglesR = 2*asin(sqrt(func(x))); R = inv(ucrGate(probQubits,valueQubit,anglesR,"y")); R.Name = "R";
Compute the probability of the value qubit being in the state . The purpose of the F circuit is to evaluate the expected value of the random variable function.
F = compositeGate([A; R],[probQubits valueQubit],Name="F"); sF = simulate(quantumCircuit(F)); probability(sF,valueQubit,"1")
ans = 0.4326
3. Quantum Phase Estimation
Estimate the value of the value qubit using quantum phase estimation.
The controlled-Q block in the circuit diagram is implemented with several composite gates as outlined in . Create a composite gate to represent this circuit element and plot the circuit to see the gates it contains.
controlQubit = valueQubit+1; % Re-mapped to estQubits in the loop below cV = czGate(controlQubit,valueQubit); % For simplicity, do not wrap in a named CompositeGate cZgates = [xGate([probQubits valueQubit]) hGate(probQubits(1)) mcxGate([probQubits(2:end) valueQubit controlQubit],probQubits(1),) hGate(probQubits(1)) xGate([probQubits valueQubit])]; cZ = compositeGate(cZgates,1:controlQubit,Name="cZ"); % Put the gates together to form controlled Q: cQ = compositeGate([cV; inv(F); cZ; F; cV; inv(F); cZ; F],[probQubits valueQubit controlQubit],Name='cQ'); QPE = ; for ii=1:n % Repeat the cQ gate as needed and map to the current control qubit: cQrep = compositeGate(repmat(cQ,2^(n-ii),1),[probQubits valueQubit estQubits(ii)],Name="cQ"+2^(n-ii)); QPE = [QPE; cQrep]; %#ok<AGROW> end plot(cQ)
Simulate Final Circuit
Assemble the final circuit using the various component circuits, and plot the result.
QubitBlocks name-value argument to display the two registers
and the single value qubit.
circ = quantumCircuit([F; hGate(estQubits); QPE; inv(qftGate(estQubits))]); plot(circ,"QubitBlocks",[m 1 n])
Simulate the circuit.
state = simulate(circ);
In terms of the phase θ, the analytical mean is given by
Solve for θ using the analytic value for μ.
AnalyticPhase = acos(-2*AnalyticMean + 1)/pi
AnalyticPhase = 0.4568
Compare the analytic phase with the quantum phase estimation.
[states,probabilities] = querystates(state,estQubits); plot(bin2dec(states)/N,probabilities) xlim([.35 .65]) xline(AnalyticPhase,'--',"AnalyticPhase") xlabel("\theta") ylabel("Probability") title("Phase Estimation with QMC")
The estimated phase is the first maximum in amplitude. Use the estimated phase to compute the mean value of the function, and compare the quantum value against the analytic value, classic Monte Carlo value, and discrete value.
y_ind = find(abs(probabilities - max(probabilities)) < 1e3*eps,1); phase_estimated = bin2dec(states(y_ind))/N; QMCMean = (1 - cos(pi * phase_estimated)) / 2; ResultTable = table(AnalyticMean,MCMean,DiscreteMean,QMCMean)
ResultTable = 1×4 table AnalyticMean MCMean DiscreteMean QMCMean ____________ _______ ____________ _______ 0.43233 0.42299 0.43264 0.42663
 Priazhkina, S., V. Skavysh, D. Guala, & T. Bromley. "Quantum Monte Carlo for Economics: Stress Testing and Macroeconomic Deep Learning." Bank of Canada working paper (2022). https://www.bankofcanada.ca/wp-content/uploads/2022/06/swp2022-29.pdf
 Brassard, G., P. Hoyer, M. Mosca, & A. Tapp. "Quantum Amplitude Amplification and Estimation." Contemporary Mathematics 305 (2002): 53-74. https://arxiv.org/pdf/quant-ph/0005055.pdf
 Rebentrost, P., B. Gupt, & T. R. Bromley. "Quantum Computational Finance: Monte Carlo Pricing of Financial Derivatives." Physical Review A, 98(2), 022321 (2018). https://arxiv.org/pdf/1805.00109.pdf
 Stamatopoulos, N., D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, & S. Woerner. "Option Pricing Using Quantum Computers." Quantum 4, 291 (2020). https://arxiv.org/pdf/1905.02666.pdf
 Woerner, S., & D. J. Egger. "Quantum Risk Analysis." npj Quantum Information 5(1), 15 (2019). https://arxiv.org/pdf/1806.06893.pdf