---
title: Multi-Product Formula (MPF) in Hamiltonian Simulation
url: https://www.emergentmind.com/topics/multi-product-formula-mpf
type: topic
---

# Multi-Product Formula (MPF) in Hamiltonian Simulation

Multi-product formula (MPF) usually denotes a construction in which a target object is approximated or expanded as a linear combination of lower-order product expressions. In quantum Hamiltonian simulation, the term refers to approximations of the time-evolution operator $U(t)=e^{-iHt}$ obtained by combining product-formula circuits evaluated at different effective step sizes so that successive terms in the Trotter error expansion cancel. In this setting, MPFs interpolate between product-formula methods and linear-combination-of-unitaries techniques: they retain the low-order circuit structure of Trotterization, but aim for near-linear dependence on evolution time and poly-logarithmic dependence on precision [2403.08922].

## 1. Formal construction

For a time-independent Hamiltonian decomposed as $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$, an MPF starts from a low-order product formula $U_p(s)$ or $S_\chi(s)$ and forms a weighted sum of repeated short-time evolutions. Two equivalent notational forms appear in the literature:
\[
U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j},
\qquad
M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).
\]
Here the $k_j$ are positive integers and the coefficients are chosen so that low-order terms in the expansion of the product formula cancel [2403.08922].

The cancellation conditions are Vandermonde-type linear constraints. In the hybrid hardware-oriented formulation, one imposes
\[
\sum_{j=1}^l a_j=1,
\qquad
\sum_{j=1}^l \frac{a_j}{k_j^\eta}=0
\]
for a prescribed range of exponents $\eta$, which removes the first $l-1$ error terms. For a second-order base formula, solving the linear system with rows $1,k_j^{-2},\ldots,k_j^{-2m+2}$ yields a short-interval global error of order $O(\Delta^{2m+1})$ [2207.11268].

This algebraic structure is closely related to Richardson extrapolation. In the notation of one later analysis, an MPF of “order” $m=2J$ is written as
\[
M_{p,m,J}(\tau)=\sum_{j=1}^J c_j\,[T_p(\tau/k_j)]^{k_j},
\qquad
M_{p,m,J}(\tau)=e^{-iH\tau}+O(\tau^{m+1}),
\]
with integer dilation factors $k_j$ and real weights $c_j$ fixed so as to cancel all terms through $\tau^m$ [2507.06557].

## 2. Error cancellation and commutator structure

A central question is whether MPF error bounds can reflect the same locality- and commutator-sensitive structure that makes product formulas attractive. In the rigorous analysis of the well-conditioned MPF, the relevant quantities are the nested-commutator sums
\[
\alpha_{\comm,j}
=
\sum_{\gamma_1,\dots,\gamma_j}
\bigl\|\,[H_{\gamma_1},[H_{\gamma_2},\dots,[H_{\gamma_{j-1}},H_{\gamma_j}]\dots]]\,\bigr\|.
\]
For one short step $\Delta$, the MPF error is bounded by a series whose terms are weighted by products of these nested-commutator sums. In the commuting case, all $\alpha_{\comm,j}$ vanish for $j\ge 2$, and the MPF is exact [2403.08922].

That result established explicit commutator scaling, but it also prompted a technical dispute about how system-size dependence should be interpreted. One later paper argues that the nested-commutator bound expressed through arbitrarily large $q$ does not by itself resolve size-efficient complexity, because locality benefits are absent when arbitrarily high nested commutators are required. The same work therefore introduces an alternative analysis based on a truncated BCH or Floquet-Magnus expansion, with truncation order
\[
p_0(N,\epsilon)=\bigl\lceil \log(3N/\epsilon)\bigr\rceil,
\]
and derives an MPF error bound involving only commutators up to order $p_0$ [2507.06557].

A common simplification is to describe MPF as merely “higher-order Trotter.” That characterization is incomplete. The defining feature is not a single higher-order product formula, but a linear combination of low-order product formulas whose error terms are arranged to cancel. This distinction matters because the resulting bounds and implementation strategies differ from those of ordinary Suzuki recursion.

## 3. Complexity and asymptotic scaling

The main asymptotic promise of MPF is a combination of near-optimal dependence on evolution time $T$ and target precision $\epsilon$ with commutator-sensitive dependence on the Hamiltonian decomposition. In the explicit-commutator analysis, one defines an MPF commutator parameter
\[
\mu_m
=
2\max_{\substack{j\in 2\mathbb Z^+,\,j\ge 2m\\1\le l\le m}}
\biggl(
\sum_{\substack{j_1+\cdots+j_l=j\\ j_k\in 2\mathbb Z^+}}
\prod_{\kappa=1}^l \alpha_{\comm,j_\kappa+1}
\biggr)^{1/(j+l)}.
\]
If these stabilize to $\mu=\sup_{m\gg 1}\mu_m$, then the MPF query cost is
\[
O\!\Bigl(\mu\,T\;\mathrm{poly}\log\!\frac{\mu T}{\epsilon}\Bigr),
\]
which is almost linear in $T$ and poly-logarithmic in $1/\epsilon$ [2403.08922].

This contrasts with conventional Trotter–Suzuki simulation. A $p$th-order product formula has cost
\[
O\bigl(\alpha_{\comm,p+1}^{1/p}\,T^{1+1/p}\,\epsilon^{-1/p}\bigr),
\]
which is super-linear in $T$ and polynomial in $1/\epsilon$. In the comparison drawn in the same analysis, post-Trotter methods such as qubitization and quantum signal processing have cost
\[
O(\|H\|\,T+\log(1/\epsilon)),
\]
which is optimal in $(T,\epsilon)$ but scales with the spectral norm $\|H\|$ and cannot exploit commutator cancellations [2403.08922].

The later truncated-commutator analysis seeks to recover size-efficient locality dependence rigorously. For a $k$-local, $g$-extensive Hamiltonian on $N$ qubits, it derives a total query cost
\[
O\!\Bigl(N^{\tfrac1{p+1}+(\log(Ngt/\varepsilon))^2}\,g\,t\;\times\;\mathrm{polylog}(Ngt/\varepsilon)\Bigr),
\]
and interprets this as polynomial in $N$ with exponent $1/(p+1)$ together with polylogarithmic dependence on $1/\varepsilon$ [2507.06557]. This suggests that the precise form of “commutator scaling” for MPF is now an active subject of refinement rather than a closed matter.

## 4. Implementation models and conditioning

The original MPF proposal is naturally implemented through a linear combination of unitaries. One prepares an ancilla superposition encoding the weights, applies controlled product-formula unitaries, and uncomputes the ancilla so that the desired linear combination is obtained on success. In this formulation, the success probability scales as $1/\|a\|_1^2$, which makes conditioning of the coefficient vector central [2207.11268].

A hybrid implementation avoids coherent LCU altogether. Instead of constructing the linear combination unitarily, one runs the distinct product-formula circuits independently, measures the same observable on each circuit, and classically combines the expectation values. The paper introducing this approach states that it has the same approximation bounds as the fully quantum MPFs, but requires no additional qubits, no controlled operations, and is not probabilistic [2207.11268].

Conditioning is quantified by the coefficient 1-norm, $\|a\|_1=\sum_j |a_j|$. Large $\|a\|_1$ amplifies both sampling error and hardware noise in the final post-processed estimate, so well-conditioned MPFs are constructed by searching over integer tuples $(k_1,\dots,k_l)$ and discarding those whose $\|a\|_1$ exceeds a chosen threshold. In the hardware-friendly study, the authors report that for $l\le 15$ one can in all practical cases achieve $\|a\|_1<1.7$ for $S_2$-based MPFs or $\|a\|_1\le 3$ for $S_1$-based MPFs [2207.11268].

That same work integrates MPF with three mitigation layers intended for noisy devices: Pauli Twirling, pulse-efficient transpilation, and zero-noise extrapolation based on scaled cross-resonance pulses. The intended role of MPF in this setting is not only asymptotic improvement, but also circuit-depth reduction under restricted hardware resources.

## 5. Dynamic and dual-channel variants

One line of development treats the MPF coefficients as time-dependent variables rather than fixed extrapolation weights. Dynamic MPF introduces coefficients $c_j(t)$ chosen to minimize a computable proxy for the Frobenius-norm projection error. Because the exact coupling term depends on the unknown target state, it is replaced by an approximate update derived from previous time steps, yielding a discrete-time linear system for the coefficient vector [2306.12569].

To stabilize this procedure in the presence of uncertainty, the same work proposes Minimax MPF. At each time step, the coefficients are obtained from a convex program with an $\ell_2$ penalty,
\[
\hat c(t_j)=\arg\min_c
\Bigl\{
\|\bar M(t_j)c-\bar A(t_j)\hat c(t_{j-1})\|_2+\varepsilon\|c\|_2
\Bigr\},
\qquad
^\top c=1,
\]
where $\varepsilon$ models bounded hardware and sampling noise. The resulting estimator is accompanied by a rigorous error bound and is described as robust to both algorithmic Trotter errors and bounded sampling and hardware noise [2306.12569].

A different recent variant is the dual-channel multi-product formula (DCMPF). It combines a regular product formula $T_p(t)$ with the reversed-sequence formula
\[
\overline T_p(t)=T_p^\dagger(-t),
\]
and defines
\[
U_{\mathrm{dual}}(\Delta t)
=
\frac12\sum_{k=1}^K c_k
\Bigl[
T_p^{\,n_k}(\Delta t/n_k)+\overline T_p^{\,n_k}(\Delta t/n_k)
\Bigr].
\]
Because the forward-plus-reversed combination is symmetric under $t\mapsto -t$, the error expansion contains only even powers of $t$, and the method cancels twice as many orders with the same number of channels $K$ [2602.01713].

For a conventional MPF built from a regular $p$th-order product formula, the leading error scales as $O((\Delta t)^{p+1+K})$; for DCMPF it scales as $O((\Delta t)^{p+1+2K})$. The same analysis states that, for the same algorithmic error, the dual-channel construction uses approximately half the circuit depth of a symmetric-PF MPF, while the sampling error remains essentially unchanged [2602.01713].

## 6. Representative applications and benchmarks

Several studies use MPF to quantify trade-offs among system size, evolution time, precision, circuit depth, and hardware noise. The examples span asymptotic analyses, numerical benchmarks, and small-scale experiments.

| Setting | Reported MPF result | Source |
|---|---|---|
| Plane-wave electronic structure on $n$ spin-orbitals | $\alpha_{\comm,j}=O(n^j)$, $\mu=O(n)$, and MPF cost $O(n^2T\,\poly\log(1/\epsilon))$ | [2403.08922] |
| 1D Heisenberg chain benchmark | Gate counts scaling roughly $n^{1.3}$–$n^{1.5}$, improving over the $n^2$-scaling of Trotter | [2403.08922] |
| Spin-boson resource study | Deepest circuit has $k_l=O\!\bigl(\tfrac{t\log^2(N_q/\epsilon_t)}{\log^2\log(N_q/\epsilon_t)}\bigr)$; reductions up to a factor of 8 for modest systems | [2207.11268] |
| Five-qubit transverse-field Ising model | Final errors $\sim10^{-2}$ with only 2–4 Trotter steps; a single PF with $k=24$ needed for similar accuracy; up to $12\times$ depth reduction | [2207.11268] |
| Heisenberg spin chain with dynamic MPF analysis | MPF can reduce required circuit depth by orders of magnitude in the near-term regime $n\lesssim 100$ | [2306.12569] |

In the plane-wave example, the explicit-commutator analysis compares second-order Trotter, post-Trotter methods, and MPF. The reported scalings are $O(n^{5/2}T^{3/2}\epsilon^{-1/2})$ for second-order Trotter, $O(n^2T\log(1/\epsilon))$ for post-Trotter methods such as qubitization, and $O(n^2T\,\poly\log(1/\epsilon))$ for MPF. The stated conclusion is a polynomial speedup in $T$ and an exponential speedup in $1/\epsilon$ over second-order Trotter, with competitiveness against qubitization when commutator structure is favorable [2403.08922].

For near-term hardware, the five-qubit transverse-field Ising demonstration is significant because it pairs well-conditioned MPF coefficients with explicit error-mitigation procedures. The reported well-conditioned pairs are $[1,2]$, $[1,3]$, and $[2,4]$, with weights $[-1,2]$, $[-\tfrac12,\tfrac32]$, and $[-1,2]$. These combinations produced errors of order $10^{-2}$ using only shallow circuits, whereas a single product formula required substantially more Trotter steps to reach similar accuracy [2207.11268].

## 7. Other mathematical uses of the term

Although current usage on arXiv is dominated by Hamiltonian simulation, “multi-product formula” is not exclusive to quantum algorithms. In stochastic analysis, the term denotes an explicit product expansion for multiple stochastic integrals with respect to the compensated random measure of a Lévy process. Under an integrability condition on contracted kernels, the product
\[
\prod_{j=1}^N I_{m_j}(f^{(j)})
\]
is expanded as a finite sum over contraction data $(\ell,\ell^\circ)\in D_N$:
\[
\prod_{j=1}^N I_{m_j}(f^{(j)})
=
\sum_{(\ell,\ell^\circ)\in D_N}
\frac{1}{m_1!\cdots m_N!\,\ell!\,\ell^\circ!}
I_{|\ell|+|\ell^\circ|}
\Bigl(f^{(1)}\star_{(\ell,\ell^\circ)}\cdots\star_{(\ell,\ell^\circ)}f^{(N)}\Bigr).
\]
That formula is used to derive moments, cumulants, mixed moments, and central limit theorems [2309.11150].

In algebraic combinatorics, the term also appears in the theory of multivariate Rogers–Szegö polynomials. There the product
\[
\tilde H_k(t)\,\tilde H_n(t)
\]
is expanded as
\[
\sum_m
(-1)^{\mathrm{wt}(m)}
\theta_{m,k,n}(q)
\Bigl(\prod_{i=2}^{\ell} e_i^{m_i}\Bigr)
\tilde H_{k+n-\mathrm{wt}(m)}(t),
\]
with recursively defined coefficients $\theta_{m,k,n}(q)$. In that setting, the formula is reinterpreted in symmetric-function terms as a statement about structure constants [1305.2404].

These non-quantum usages are terminologically related but mathematically independent. In contemporary quantum-computing literature, MPF almost always refers to the linear-combination construction for Hamiltonian simulation; in other fields, it names a product-expansion identity specific to the objects under study.

Source: https://www.emergentmind.com/topics/multi-product-formula-mpf