From Papers and Equations to FlagQuantum Code

From Papers and Equations to FlagQuantum Code

Turn a paper's mathematical model into a runnable FlagQuantum experiment, with a traceable correspondence between equations, tensors, circuit operations, and numerical checks.

Version: FlagQuantum 0.3.1, reference commit 5714b569cfc94e212b6d175eb7c87e1f36d8cf9d.

Verified: October 20, 2026, macOS Apple Silicon, Python 3.12.14, PyTorch 2.13.0, CPU, float64 parameters and complex128 statevectors. All three Python files below were executed. They use deterministic, noiseless statevector expectations; no hardware, shots, or accelerator access is required.

Scope: The worked model is a deliberately small teaching example, not a reproduction of a particular paper. Its equations are derived here. Reproducing a published experiment additionally requires its actual architecture, data, initialization, measurement protocol, and evaluation procedure.

Before starting, read Core Concepts and Qubit Ordering and Result Conventions. For environment setup, use Installation and Compatibility.

Prepare the Environment

Use a Python 3.12 virtual environment. For example, on macOS or Linux:

python3.12 -m venv .venv
source .venv/bin/activate
python -m pip install "flagquantum==0.3.1" "torch==2.13.0"
python -c "import flagquantum as fq, torch; print(fq.__version__, torch.__version__)"

The version check should report FlagQuantum 0.3.1 and PyTorch 2.13.0 (a platform-specific PyTorch build suffix may appear). Windows activation and alternative supported environments are covered in Installation and Compatibility. Save the three files below in one working directory and run them with this environment's Python.

1. Build an Equation-to-Code Map

Start with the paper's definitions and assumptions, then derive the computation before selecting API calls:

Paper definitions → mathematical derivation → tensor shapes and conventions → FlagQuantum implementation → independent numerical checks.

Read each derivation together with the equation/code table immediately below it. The tables show implementation excerpts; the complete runnable files follow the derivations. Matching equation numbers also appear in the script comments.

For each equation, record its equation number, the symbols it introduces, its inputs and outputs, and the code that implements it. Keep fixed constants, trainable parameters, classical inputs, and random variables distinct.

Mathematical object	Implementation in this guide	What to verify
Initial state $\lvert 00\rangle$	Default initial state of fq.Circuit(2)	Basis ordering and normalization
Variational parameters $\boldsymbol\theta$	torch.nn.Parameter, shape (2,)	Units, initialization, and gradient connection
Unitary $U(\boldsymbol\theta)$	.ry(...).ry(...).cx(...)	Gate order, wires, angle convention
Pauli products $P_k$	fq.Z(0), fq.X(0) @ fq.X(1)	Tensor-product convention
Expectations $\langle P_k\rangle$	Named fq.expectation(...) requests	Result names and batch dimension
Energy $E=\sum_k h_k\langle P_k\rangle$	Weighted sum of expectation tensors in PyTorch	Coefficients, signs, real scalar objective
Gradient $\nabla_\theta E$	torch.autograd.grad or loss.backward()	Closed-form, matrix, and parameter-shift checks
Update $\theta\leftarrow\mathrm{Optimizer}(\theta,\nabla E)$	torch.optim.Adam([theta], ...)	Fresh forward pass and post-update evaluation

This example assembles the Hamiltonian's weighted sum in PyTorch. It does not require an additional Hamiltonian-sum API.

Conventions to Set Before Coding

Convention	Used here	Translation risk
Basis	$\lvert q_0q_1\rangle$: 00, 01, 10, 11; qubit 0 is most significant	Another framework may reverse wire-to-bit mapping
Tensor product	$Z_0=Z\otimes I$, $Z_1=I\otimes Z$	Swapping factors changes observables
Angles	Radians; $R_Y(\phi)=e^{-i\phi Y/2}$	Some papers use $e^{-i\theta Y}$
Circuit order	Appended gates execute from left to right	Matrix products act on a ket from right to left
Parameter shape	(2,) means two circuit parameters in this example	Do not assume this is a batch of two circuits
Result shapes	State (1, 4), each expectation (1,), stacked terms (1, 4), energy (1,)	Preserve the leading circuit-batch dimension
Numerical types	Real parameters float64, complex states complex128	Lower precision needs different tolerances
Measurement	Exact simulator expectation	Finite-shot estimates have sampling error

The one-element batch is reduced to a scalar with .squeeze(0) for checking or .mean() for training. Extending to multiple circuits requires a deliberate batching and loss-reduction policy. The examples below accept a single two-parameter vector.

2. Derive and Check a Single Rotation

For the computational basis and the usual Pauli $Y$,

$$R_Y(\theta)=\begin{pmatrix} \cos(\theta/2)&-\sin(\theta/2)\\\ \sin(\theta/2)&\cos(\theta/2) \end{pmatrix},\qquad R_Y(\theta)\lvert0\rangle= \cos(\theta/2)\lvert0\rangle+\sin(\theta/2)\lvert1\rangle.$$

Therefore,

$$p(0)=\cos^2(\theta/2),\quad p(1)=\sin^2(\theta/2),\quad E(\theta)=\langle Z\rangle=p(0)-p(1)=\cos\theta, \qquad E'(\theta)=-\sin\theta.$$

Equation	Corresponding Python	Meaning
$R_Y(\theta)\lvert0\rangle$	circuit = fq.Circuit(1, dtype=torch.complex128).ry(0, theta)	Apply the paper's rotation to the default zero state.
$E(\theta)=\langle Z\rangle$	request = fq.expectation(fq.Z(0), name="z") energy = fq.run(circuit, outputs=request, options=options).expectation("z").squeeze(0)	Request the observable, then read a scalar expectation.
$E(\theta)=\cos\theta$	torch.testing.assert_close(energy, theta.cos())	Compare the simulator value with the derived expression.
$dE/d\theta=-\sin\theta$	gradient, = torch.autograd.grad(energy, theta) torch.testing.assert_close(gradient, -theta.sin())	Differentiate the implementation and check the derivative.

If a paper instead defines $U(\theta)=e^{-i\theta Y}$, implement ry(0, 2 * theta). The chain rule gives

$$E_{\mathrm{paper}}(\theta)=\cos(2\theta),\qquad \frac{d E_{\mathrm{paper}}}{d\theta}=-2\sin(2\theta).$$

Paper convention	Code that preserves the chain rule
$U(\theta)=e^{-i\theta Y}=R_Y(2\theta)$	paper_energy = z_expectation(2 * theta)
$dE_{\mathrm{paper}}/d\theta=-2\sin(2\theta)$	paper_gradient, = torch.autograd.grad(paper_energy, theta) torch.testing.assert_close(paper_gradient, -2 * (2 * theta).sin())

Passing theta directly to ry would implement a different model. Converting 2 * theta to a Python number would disconnect its gradient.

Save the following as rotation_translation.py:

import torch
import flagquantum as fq

options = fq.ExecutionOptions(
    mode="statevector", backend="pytorch", device="cpu",
    precision="complex128", allow_backend_fallback=False,
    allow_approximate=False,
)

def z_expectation(angle):
    circuit = fq.Circuit(1, dtype=torch.complex128).ry(0, angle)
    return fq.run(
        circuit, outputs=fq.expectation(fq.Z(0), name="z"), options=options,
    ).expectation("z").squeeze(0)

theta = torch.tensor(0.6, dtype=torch.float64, requires_grad=True)
energy = z_expectation(theta)
gradient, = torch.autograd.grad(energy, theta)
torch.testing.assert_close(energy, theta.cos())
torch.testing.assert_close(gradient, -theta.sin())
print(f"RY(theta): expectation={energy.item():.6f}, gradient={gradient.item():.6f}")

# The paper defines U(theta) = exp(-i * theta * Y), without the factor 1/2.
# Preserve 2 * theta inside the computation graph.
paper_energy = z_expectation(2 * theta)
paper_gradient, = torch.autograd.grad(paper_energy, theta)
torch.testing.assert_close(paper_energy, (2 * theta).cos())
torch.testing.assert_close(paper_gradient, -2 * (2 * theta).sin())
print(f"exp(-i*theta*Y): expectation={paper_energy.item():.6f}, "
      f"gradient={paper_gradient.item():.6f}")
print("Rotation convention checks passed.")

Run:

python rotation_translation.py

Expected output, rounded to six decimals:

RY(theta): expectation=0.825336, gradient=-0.564642
exp(-i*theta*Y): expectation=0.362358, gradient=-1.864078
Rotation convention checks passed.

3. Translate a Two-Qubit Variational Model

Consider the dimensionless Hamiltonian

Equation (1).

$$H=0.7Z_0+0.2Z_1+0.4X_0X_1-0.3Z_0Z_1.$$

Mathematical term	Observable request in evaluate	Coefficient
$0.7Z_0$	fq.expectation(fq.Z(0), name="z0")	0.7
$0.2Z_1$	fq.expectation(fq.Z(1), name="z1")	0.2
$0.4X_0X_1$	fq.expectation(fq.X(0) @ fq.X(1), name="xx")	0.4
$-0.3Z_0Z_1$	fq.expectation(fq.Z(0) @ fq.Z(1), name="zz")	-0.3

The requests and COEFFICIENTS = torch.tensor([0.7, 0.2, 0.4, -0.3], dtype=REAL) use the same term order. Swapping result order without swapping coefficients changes the Hamiltonian.

We minimize its energy over the trial state

Equation (2).

$$\lvert\psi(\boldsymbol\theta)\rangle= \mathrm{CX}_{0\rightarrow1} \bigl(R_Y(\theta_0)\otimes R_Y(\theta_1)\bigr)\lvert00\rangle.$$

Formula operation, in execution order	Circuit code
Prepare $\lvert00\rangle$	circuit = fq.Circuit(2, dtype=COMPLEX)
Apply $R_Y(\theta_0)\otimes R_Y(\theta_1)$	circuit.ry(0, theta[0]).ry(1, theta[1])
Apply $\mathrm{CX}_{0\rightarrow1}$	circuit.cx(0, 1)

The last two code rows are chained together in the complete script. Parameter $\theta_0$ maps to theta[0], and $\theta_1$ maps to theta[1].

Equation (2) becomes .ry(0, theta[0]).ry(1, theta[1]).cx(0, 1): rotations execute first, then controlled-X. The first two gates commute because they act on different qubits; moving controlled-X ahead of them would generally change the model.

Expand the State and Observables

Let $c_j=\cos(\theta_j/2)$ and $s_j=\sin(\theta_j/2)$. Before controlled-X, the amplitudes are $(c_0c_1,c_0s_1,s_0c_1,s_0s_1)$. Controlled-X exchanges basis states 10 and 11, giving

Equation (3).

$$\lvert\psi\rangle= \begin{pmatrix}c_0c_1\\c_0s_1\\s_0s_1\\s_0c_1\end{pmatrix}.$$

State entry	Statevector position	Analytical value
Amplitude of 00	state[0]	c0 * c1
Amplitude of 01	state[1]	c0 * s1
Amplitude of 10	state[2]	s0 * s1
Amplitude of 11	state[3]	s0 * c1

Here state = result.statevector().squeeze(0), c0, c1 = (theta / 2).cos().unbind(), and s0, s1 = (theta / 2).sin().unbind(). The validation below constructs this analytical state and compares all four amplitudes, then checks an independent dense matrix implementation.

For real amplitudes $a_{00},a_{01},a_{10},a_{11}$, diagonal Pauli products are signed probability sums, while $X_0X_1$ exchanges both bits:

Equation (4).

$$\begin{aligned} \langle Z_0\rangle &= a_{00}^2+a_{01}^2-a_{10}^2-a_{11}^2=\cos\theta_0,\\\ \langle Z_1\rangle &= a_{00}^2-a_{01}^2+a_{10}^2-a_{11}^2=\cos\theta_0\cos\theta_1,\\\ \langle X_0X_1\rangle &=2(a_{00}a_{11}+a_{01}a_{10})=\sin\theta_0,\\\ \langle Z_0Z_1\rangle &=a_{00}^2-a_{01}^2-a_{10}^2+a_{11}^2=\cos\theta_1. \end{aligned}$$

Derived expectation	Simulator result	Analytical reference in expected_terms
$\langle Z_0\rangle=\cos\theta_0$	result.expectation("z0")	a.cos()
$\langle Z_1\rangle=\cos\theta_0\cos\theta_1$	result.expectation("z1")	a.cos() * b.cos()
$\langle X_0X_1\rangle=\sin\theta_0$	result.expectation("xx")	a.sin()
$\langle Z_0Z_1\rangle=\cos\theta_1$	result.expectation("zz")	b.cos()

In the script, a, b = theta.unbind() gives the scalar tensors for $\theta_0,\theta_1$. torch.stack(..., dim=-1) preserves the result order and builds terms with shape (1, 4).

For general complex amplitudes, use absolute squares for probabilities and the real part of conjugate overlap terms for $X_0X_1$. The simplified real-amplitude expressions above apply to this ansatz.

Substituting (4) into (1) gives an independently checkable energy:

Equation (5).

$$E(\boldsymbol\theta)=0.7\cos\theta_0 +0.2\cos\theta_0\cos\theta_1 +0.4\sin\theta_0-0.3\cos\theta_1.$$

Formula step	Corresponding code	Shape
Collect $(\langle Z_0\rangle,\langle Z_1\rangle,\langle X_0X_1\rangle,\langle Z_0Z_1\rangle)$	terms = torch.stack([result.expectation(name) for name in NAMES], dim=-1)	(1, 4)
Multiply each term by $h_k$	terms * COEFFICIENTS	(1, 4)
Sum over terms: $E=\sum_k h_k\langle P_k\rangle$	energy = (terms * COEFFICIENTS).sum(dim=-1)	(1,)
Select the single-circuit scalar	scalar = energy.squeeze(0)	()
Check equation (5)	torch.testing.assert_close(scalar, closed_form(theta), rtol=1e-9, atol=1e-9)	Scalar comparison

closed_form translates the trigonometric expression literally; evaluate obtains the same quantity by executing the circuit. Keeping these two paths separate makes their agreement meaningful.

Differentiate (5) term by term:

Equation (6).

$$\begin{aligned} \partial_{\theta_0}E &=-0.7\sin\theta_0-0.2\sin\theta_0\cos\theta_1+0.4\cos\theta_0,\\\ \partial_{\theta_1}E &=-0.2\cos\theta_0\sin\theta_1+0.3\sin\theta_1. \end{aligned}$$

Derivative	Literal formula implementation in closed_gradient
$\partial_{\theta_0}E$	-0.7 * a.sin() - 0.2 * a.sin() * b.cos() + 0.4 * a.cos()
$\partial_{\theta_1}E$	-0.2 * a.cos() * b.sin() + 0.3 * b.sin()
Simulator gradient $\nabla_\theta E$	gradient, = torch.autograd.grad(scalar, theta)
Compare both derivatives	torch.testing.assert_close(gradient, closed_gradient(theta), rtol=1e-9, atol=1e-9)

The first two rows are stacked into a (2,) tensor. Autograd differentiates the executed circuit path, while closed_gradient implements the paper-style derivation explicitly.

Connect Each Equation to the Implementation

Equation	Code location below
(1): coefficients and Pauli products	COEFFICIENTS, named requests in evaluate, and hamiltonian in dense_reference
(2): ansatz	Circuit construction in evaluate; matrix product in dense_reference
(3): state amplitudes	result.statevector() compared with the independent matrix state
(4): individual expectations	expected_terms compared with the named results
(5): energy	Weighted sum in evaluate, closed_form, and dense torch.vdot
(6): derivatives	closed_gradient compared with autograd and parameter shift

Save the following as paper_model.py, in the same directory as the later training script. Running it performs the validation; importing it provides the reusable model functions.

import torch
import flagquantum as fq

REAL = torch.float64
COMPLEX = torch.complex128
OPTIONS = fq.ExecutionOptions(
    mode="statevector", backend="pytorch", device="cpu",
    precision="complex128", allow_backend_fallback=False,
    allow_approximate=False,
)
# Eq. (1): coefficients in z0, z1, xx, zz order.
COEFFICIENTS = torch.tensor([0.7, 0.2, 0.4, -0.3], dtype=REAL)
NAMES = ("z0", "z1", "xx", "zz")

def evaluate(theta):
    # Rebuild the circuit for every forward pass; theta stays in the graph.
    circuit = fq.Circuit(2, dtype=COMPLEX)
    # Eq. (2): rotations execute before controlled-X.
    circuit.ry(0, theta[0]).ry(1, theta[1]).cx(0, 1)
    result = fq.run(circuit, outputs=[
        fq.expectation(fq.Z(0), name="z0"),
        fq.expectation(fq.Z(1), name="z1"),
        fq.expectation(fq.X(0) @ fq.X(1), name="xx"),
        fq.expectation(fq.Z(0) @ fq.Z(1), name="zz"),
    ], options=OPTIONS)
    terms = torch.stack([result.expectation(name) for name in NAMES], dim=-1)
    # Eq. (5): sum the weighted expectation values.
    energy = (terms * COEFFICIENTS).sum(dim=-1)
    return energy, terms, result

def closed_form(theta):
    # Eq. (5): independently derived scalar expression.
    a, b = theta.unbind()
    return (0.7 * a.cos() + 0.2 * a.cos() * b.cos()
            + 0.4 * a.sin() - 0.3 * b.cos())

def closed_gradient(theta):
    # Eq. (6): explicit derivatives with respect to theta[0] and theta[1].
    a, b = theta.unbind()
    return torch.stack([
        -0.7 * a.sin() - 0.2 * a.sin() * b.cos() + 0.4 * a.cos(),
        -0.2 * a.cos() * b.sin() + 0.3 * b.sin(),
    ])

def dense_reference(theta):
    # Independent dense matrix implementation, using |q0 q1> ordering.
    identity = torch.eye(2, dtype=COMPLEX)
    x = torch.tensor([[0, 1], [1, 0]], dtype=COMPLEX)
    z = torch.diag(torch.tensor([1, -1], dtype=COMPLEX))
    hamiltonian = (0.7 * torch.kron(z, identity)
                   + 0.2 * torch.kron(identity, z)
                   + 0.4 * torch.kron(x, x) - 0.3 * torch.kron(z, z))
    def ry(angle):
        c, s = (angle / 2).cos(), (angle / 2).sin()
        return torch.stack([torch.stack([c, -s]),
                            torch.stack([s, c])]).to(COMPLEX)
    cx = torch.tensor([[1, 0, 0, 0], [0, 1, 0, 0],
                       [0, 0, 0, 1], [0, 0, 1, 0]], dtype=COMPLEX)
    zero = torch.tensor([1, 0, 0, 0], dtype=COMPLEX)
    state = cx @ torch.kron(ry(theta[0]), ry(theta[1])) @ zero
    energy = torch.vdot(state, hamiltonian @ state).real
    return energy, state, hamiltonian

if __name__ == "__main__":
    torch.set_printoptions(precision=6, sci_mode=False)
    theta = torch.tensor([0.6, -0.4], dtype=REAL, requires_grad=True)
    energy, terms, result = evaluate(theta)
    scalar = energy.squeeze(0)  # One circuit: batch size is one.
    gradient, = torch.autograd.grad(scalar, theta)
    reference_energy, reference_state, _ = dense_reference(theta)
    dense_gradient, = torch.autograd.grad(reference_energy, theta)
    a, b = theta.unbind()
    # Eq. (4): check each observable before checking their weighted sum.
    expected_terms = torch.stack([a.cos(), a.cos() * b.cos(),
                                  a.sin(), b.cos()]).unsqueeze(0)
    torch.testing.assert_close(terms, expected_terms, rtol=1e-9, atol=1e-9)
    # Eq. (3): compare all four amplitudes in |q0 q1> order.
    c0, c1 = (theta / 2).cos().unbind()
    s0, s1 = (theta / 2).sin().unbind()
    analytical_state = torch.stack([c0*c1, c0*s1, s0*s1, s0*c1]).to(COMPLEX)
    torch.testing.assert_close(result.statevector().squeeze(0), analytical_state,
                               rtol=1e-9, atol=1e-9)
    torch.testing.assert_close(result.statevector().squeeze(0), reference_state,
                               rtol=1e-9, atol=1e-9)
    torch.testing.assert_close(scalar, closed_form(theta), rtol=1e-9, atol=1e-9)
    torch.testing.assert_close(scalar, reference_energy, rtol=1e-9, atol=1e-9)
    torch.testing.assert_close(gradient, closed_gradient(theta), rtol=1e-9, atol=1e-9)
    torch.testing.assert_close(gradient, dense_gradient, rtol=1e-9, atol=1e-9)
    with torch.no_grad():
        # Eq. (7): independent plus/minus gate-angle evaluations.
        shifted = []
        for index in range(theta.numel()):
            shift = torch.zeros_like(theta)
            shift[index] = torch.pi / 2
            plus = evaluate(theta + shift)[0].squeeze(0)
            minus = evaluate(theta - shift)[0].squeeze(0)
            shifted.append((plus - minus) / 2)
        shift_gradient = torch.stack(shifted)
    torch.testing.assert_close(gradient, shift_gradient, rtol=1e-9, atol=1e-9)
    # Additional points avoid validating only one accidental agreement.
    for angles in ([0.0, 0.0], [1.1, 0.8], [-0.7, 2.0]):
        point = torch.tensor(angles, dtype=REAL, requires_grad=True)
        value = evaluate(point)[0].squeeze(0)
        derivative, = torch.autograd.grad(value, point)
        torch.testing.assert_close(value, dense_reference(point)[0], rtol=1e-9, atol=1e-9)
        torch.testing.assert_close(derivative, closed_gradient(point), rtol=1e-9, atol=1e-9)
    print("terms:", terms.detach())
    print(f"energy: {scalar.item():.6f}")
    print("autograd:", gradient)
    print("parameter shift:", shift_gradient)
    print("State, observable, energy, and gradient checks passed.")

Run:

python paper_model.py

Expected output:

terms: tensor([[0.825336, 0.760184, 0.564642, 0.921061]], dtype=torch.float64)
energy: 0.679311
autograd: tensor([-0.169130, -0.052545], dtype=torch.float64)
parameter shift: tensor([-0.169130, -0.052545], dtype=torch.float64)
State, observable, energy, and gradient checks passed.

The assertions use absolute and relative tolerances of 1e-9, rather than testing rounded printed values. State, individual observables, energy, and derivatives are checked at the initial point; energy and derivatives are also checked at three additional points. These are focused consistency checks, not a proof that an arbitrary paper implementation is correct.

Why the Parameter-Shift Check Applies Here

Each parameter occurs in exactly one gate $e^{-i\theta_jY/2}$, and the objective is a linear combination of exact expectations. Thus

Equation (7).

$$\partial_{\theta_j}E(\boldsymbol\theta)=\frac{ E(\boldsymbol\theta+\tfrac{\pi}{2}\mathbf e_j) -E(\boldsymbol\theta-\tfrac{\pi}{2}\mathbf e_j)}{2}.$$

Formula component	Code inside the parameter-shift loop
$(\pi/2)\mathbf e_j$	shift = torch.zeros_like(theta) shift[index] = torch.pi / 2
$E(\theta+(\pi/2)\mathbf e_j)$	plus = evaluate(theta + shift)[0].squeeze(0)
$E(\theta-(\pi/2)\mathbf e_j)$	minus = evaluate(theta - shift)[0].squeeze(0)
$(E_+-E_-)/2$	shifted.append((plus - minus) / 2)

The script evaluates these shifted circuits explicitly. This check does not imply that the selected simulator's autograd internally uses parameter shift.

If a paper reuses one parameter in several gates, differentiate each occurrence and sum the chain-rule contributions. Shifting the shared parameter globally by $\pi/2$ generally does not give (7). Rescaled angles need their scale factors too. For a nonlinear loss $L=f(E)$, apply the classical chain rule $dL/d\theta=f'(E)\,dE/d\theta$; substituting the full nonlinear loss into (7) is generally invalid. PyTorch autograd tracks these dependencies when all angle and loss operations remain differentiable tensors.

4. Optimize and Validate the Final Result

For a normalized pure state and Hermitian $H$, the variational principle gives

Equation (8).

$$E(\boldsymbol\theta)\geq\lambda_{\min}(H).$$

Here $H$ preserves even and odd computational-basis parity. In the ordered even basis (00, 11) and odd basis (01, 10), its two blocks are

$$H_{\mathrm{even}}=\begin{pmatrix}0.6&0.4\\0.4&-1.2\end{pmatrix},\qquad H_{\mathrm{odd}}=\begin{pmatrix}0.8&0.4\\0.4&-0.2\end{pmatrix}.$$

The smaller eigenvalues are $-0.3-\sqrt{0.97}$ and $0.3-\sqrt{0.41}$, respectively. Therefore,

Equation (9).

$$\lambda_{\min}(H)=-0.3-\sqrt{0.97}\approx-1.284885780.$$

Mathematical check	Corresponding code in the training script
Find $\lambda_{\min}(H)$	exact_ground = torch.linalg.eigvalsh(hamiltonian).min()
Check equation (8), allowing roundoff	assert final_energy >= exact_ground - 1e-9
Check that training reaches equation (9)	torch.testing.assert_close(final_energy, exact_ground, rtol=0, atol=1e-6)

The ansatz can reach this ground state: setting $\theta_1=0$ restricts it to the even subspace, and $\theta_0=\pi+\arctan(0.4/0.9)$ minimizes its energy. This is a special property of this small example; a general ansatz may not contain the true ground state.

Match the Optimization Equations to the Training Loop

Mathematical operation	PyTorch code	Interpretation
Initialize $\theta^{(0)}=(0.6,-0.4)$	theta = torch.nn.Parameter(torch.tensor([0.6, -0.4], dtype=REAL))	Register the two trainable angles.
Define $L(\theta)=E(\theta)$	loss = evaluate(theta)[0].mean()	For the one-circuit batch, the mean is its energy.
Clear the previous gradient	optimizer.zero_grad(set_to_none=True)	Each iteration uses its own gradient.
Compute $g_t=\nabla_\theta L(\theta^{(t)})$	loss.backward()	The gradient is accumulated into theta.grad.
Apply the Adam update	optimizer.step()	Use Adam's moment estimates and learning rate to update theta.
Evaluate $E(\theta^{(t+1)})$	final_energy = evaluate(theta)[0].mean()	Recompute after the last parameter update.

A formula such as $\theta^{(t+1)}=\theta^{(t)}-\eta g_t$ describes plain gradient descent. To implement that update directly, select torch.optim.SGD([theta], lr=eta) with its default zero momentum. The complete example uses Adam, whose update also depends on gradient history; matching a paper requires matching its actual optimizer.

Save as train_paper_model.py, alongside paper_model.py:

import torch
from paper_model import REAL, evaluate, dense_reference

theta = torch.nn.Parameter(torch.tensor([0.6, -0.4], dtype=REAL))
optimizer = torch.optim.Adam([theta], lr=0.05)
initial_energy = evaluate(theta)[0].mean().item()
for _ in range(300):
    optimizer.zero_grad(set_to_none=True)
    loss = evaluate(theta)[0].mean()  # Mean over the one-element batch.
    loss.backward()
    optimizer.step()

# Recompute after the last update, rather than printing a pre-update loss.
with torch.no_grad():
    final_energy = evaluate(theta)[0].mean()
    _, _, hamiltonian = dense_reference(theta)
    exact_ground = torch.linalg.eigvalsh(hamiltonian).min()
    assert final_energy >= exact_ground - 1e-9
    torch.testing.assert_close(final_energy, exact_ground, rtol=0, atol=1e-6)
    print(f"initial energy: {initial_energy:.6f}")
    print(f"final energy: {final_energy.item():.6f}")
    print(f"exact ground energy: {exact_ground.item():.6f}")
    print(f"energy gap below 1e-6: {abs((final_energy - exact_ground).item()) < 1e-6}")
    print("Training and variational-bound checks passed.")

Run from the directory containing both files:

python train_paper_model.py

Expected output:

initial energy: 0.679311
final energy: -1.284886
exact ground energy: -1.284886
energy gap below 1e-6: True
Training and variational-bound checks passed.

The final energy is checked against dense exact diagonalization with an absolute tolerance of 1e-6. Small numerical differences across platforms are possible. The fixed starting point and 300 Adam steps were verified in the environment above; they are not a convergence guarantee for other models or settings. Equivalent angles can represent the same physical state, so compare energies and observables before requiring parameters to match.

Dense reference matrices require storage exponential in qubit count. Use them for small correctness cases, not as a scalable implementation strategy. For larger experiments, validate small reductions first and then follow Choosing a Simulation Method.

For a complete source-based case study, follow Reproducing Molecular Hydrogen VQE, including corrected paper data, a 54-geometry energy curve, and independently checked gradients.

5. Reuse This Workflow for a Real Paper

First identify exactly what you intend to reproduce: a mathematical identity, an equivalent circuit, a reported experiment, or a performance result. Agreement at one level does not establish agreement at the others.

Use this template when documenting each translation:

# Paper-to-Code Experiment: <Short Title>

## Source and Scope
- Paper: <title, authors, persistent URL, version/date>
- Equations/figures implemented: <numbers and page references>
- Target claim: <specific observable, accuracy, energy, or runtime result>
- Assumptions: <pure/mixed state, noise, boundary conditions, approximations>
- Not reproduced: <explicit omissions or deviations>

## Equation-to-Code Map
| Source equation | Meaning and assumptions | Shapes/conventions | Code function | Independent check |
| --- | --- | --- | --- | --- |
| <Eq. N> | <definition> | <basis, units, dimensions> | <function> | <identity/reference> |

## Experiment Specification
- FlagQuantum: <release and/or full commit>
- Environment: <Python, PyTorch, OS, hardware, backend, precision>
- State preparation and circuit: <wires, gate order, depth>
- Parameters: <initialization, sharing, random seeds>
- Inputs/data: <source, preprocessing, train/test split, access requirements>
- Observables/objective: <coefficients, normalization, reduction>
- Gradients/optimizer: <method, learning rate, steps, stopping rule>
- Measurement: <exact expectations or shots, noise, grouping, readout>

## Runnable Implementation
<Complete files, dependency setup, and exact execution commands.>

## Validation and Results
- Small analytical or dense reference: <expected and observed>
- State/observable/energy/gradient checks: <tolerances and observed errors>
- Paper comparison: <reported value, reproduced value, uncertainty, deviations>
- Performance claims: <warmup, synchronization, repeats, timing boundary>

## Limitations and Evidence
<Unsupported features, failed checks, approximation effects, and scope of evidence.>

For sampled measurements, compare statistical estimates with suitable uncertainty and repeat counts; deterministic 1e-9 assertions are unsuitable. For noisy or open systems, use the density-matrix or trajectory formulation appropriate to the model rather than assuming the pure-state equations above still apply. A matching energy alone can hide incorrect states or observables, which is why this guide checks several layers.

6. Common Translation Errors

Symptom	Check
Values look plausible but disagree with the paper	Angle half-factor, coefficient signs, identity offsets, units, and wire order
Correct expectations but wrong probabilities	Expectation versus bit probability; for one Pauli-Z measurement, $p(0)=(1+\langle Z\rangle)/2$
Correct energy but wrong state	Degeneracy, global phase, hidden basis permutations, and incomplete observable checks
Gradient is missing or unexpectedly zero	Tensor detachment, .item() in the forward path, parameter sharing, and stationary initialization
Gradient differs by a constant factor	Rescaled angles and sum-versus-mean loss reduction
Parameter-shift and autograd disagree	Generator spectrum, repeated parameters, nonlinear loss, shifts, or shot noise
Printed final loss is stale	Recompute the forward pass after the last optimizer update
Tiny example agrees but a larger run does not	Backend, precision, approximation tolerance, measurement protocol, and model scaling

Physical pure states are equivalent up to global phase. The direct statevector equality in this guide is appropriate because both implementations use the same explicit gate matrices and initial state; when comparing other implementations, align global phase or compare fidelity and observables.

Sources and Further Reading

Release-pinned implementation references:

• Circuit construction
• Observable construction
• Execution options
• Execution results

Continue with Training Parameterized Circuits, Hybrid Classical-Quantum Models, Noise and Open-System Simulation, and Reproducibility and Benchmarks. Use Documentation Guide to maintain source references, runnable examples, and version labels.