SQWLib: Efficient Szegedy Walk Simulator

• SQWLib is a Python library that implements Szegedy quantum walk simulations and QPE protocols without constructing full N²×N² unitary matrices.
• It employs a matrix-based state representation with structured arrays to achieve O(N²) time for dense graphs and O(|E|) scaling for sparse graphs.
• The library supports multiple formulations—including single-step, double-step, and similarity-transformed walks—and integrates efficient QPE for advanced algorithm simulations.

Searching arXiv for the cited SQWLib/SQUWALS papers to ground the article. SQWLib is a Python library that implements the full simulation framework developed in “Efficient Simulation of Szegedy Quantum Walk Formulations and Algorithms.” It is a classical, NumPy/SciPy–based simulator specialized for Szegedy-type discrete-time quantum walks and for algorithms that couple these walks to quantum phase estimation (QPE). Its defining design goal is to avoid constructing the full unitary evolution operator of size $N^2\times N^2$ for a graph with $N$ nodes, while supporting several different but spectrally equivalent walk formulations. For dense graphs it achieves $O(N^2)$ time and memory, and for sparse graphs its cost scales as $O(|E|)$, where $|E|$ is the number of undirected edges (Ortega et al., 12 Jun 2026).

1. Formal setting and supported walk formulations

SQWLib is built around Szegedy’s quantization of a classical Markov chain. Given a graph with $N$ nodes and a column-stochastic transition matrix $G$, the walk acts on the doubled Hilbert space

$$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$

For each node $i$, one defines

$$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$

The projector onto the span of the $N$0 is

$N$1

and the basic reflection and swap operators are

$N$2

SQWLib supports the principal Szegedy walk formulations used in the literature.

Formulation	Expression	Role
Single-step Szegedy operator	$N$3	Common one-step evolution
Double-step Szegedy operator	$N$4	Szegedy’s original two-reflection form
Similarity-transformed walk	$N$5	Used in QSA and related algorithms

For reversible Markov chains with stationary distribution $N$6, the walk admits the stationary eigenstate

$N$7

Within the dynamical subspace

$N$8

this is the unique eigenstate with eigenvalue $N$9 for $O(N^2)$0, and measuring the first register of $O(N^2)$1 yields a sample from the classical stationary distribution. SQWLib also supports the update-operator formulation based on $O(N^2)$2 and

$O(N^2)$3

for which the stationary state takes the form

$O(N^2)$4

This state is independent of the Markov chain outside the amplitudes $O(N^2)$5, a property used directly in quantum simulated annealing (Ortega et al., 12 Jun 2026).

2. State representation and primitive operators

The library’s efficiency derives from representing walk states as structured arrays rather than as flat vectors. Any state

$O(N^2)$6

is stored as an $O(N^2)$7 matrix $O(N^2)$8 with entries

$O(N^2)$9

Under this convention, the column index corresponds to the first register and the row index to the second register. This layout makes several basic operations immediate: the swap $O(|E|)$0 becomes matrix transposition, $O(|E|)$1, and the simple reflection $O(|E|)$2 is implemented by multiplying all rows except the first by $O(|E|)$3 (Ortega et al., 12 Jun 2026).

The nontrivial primitive is the reflection $O(|E|)$4 about $O(|E|)$5. SQWLib implements it through an expansion–recomposition scheme. For operators of the form

$O(|E|)$6

with

$O(|E|)$7

the expansion step computes coefficients

$O(|E|)$8

and the recomposition step forms $O(|E|)$9. Numerically, expansion is an element-wise multiplication followed by summation along an axis, while recomposition is a broadcasting operation. For $|E|$0, one uses $|E|$1, with

$|E|$2

SQWLib further introduces two concrete unitary implementations of the update operator $|E|$3. The first, $|E|$4, is defined on the update subspace

$|E|$5

and acts as the identity on $|E|$6. The second, $|E|$7, is a block-diagonal Householder construction,

$|E|$8

with $|E|$9 and $N$0. The paper states that $N$1 is more efficient than $N$2, since it requires fewer expansion/recomposition steps and has the self-inverse property (Ortega et al., 12 Jun 2026).

3. Dense and sparse simulation framework

The central algorithmic principle is to simulate the action of primitive operators rather than assemble full $N$3 unitary matrices. In dense settings, all major objects—$N$4, $N$5, and the state matrices $N$6—are $N$7 dense arrays. Each basic operator touches each entry at most a constant number of times, so time and memory per walk step are $N$8. The paper characterizes this as optimal, because the minimal description of a dense Markov chain already requires $N$9 numbers (Ortega et al., 12 Jun 2026).

For sparse graphs, SQWLib exploits the fact that the walk is confined to a reduced subspace

$G$0

where $G$1 is the adjacency matrix of the undirected backbone. Instead of storing an $G$2 matrix, the state is stored as a reduced vector of length $G$3. The implementation uses three auxiliary structures: STARTS, a list of offsets marking the start of each block; SIZES, the block sizes; and PERMUTATION, a permutation array that maps each reduced index to the index corresponding to the swapped edge.

Under this representation, swap $G$4 is implemented by applying PERMUTATION, while $G$5 multiplies all entries by $G$6 and then flips the sign back on the indices corresponding to $G$7. Expansion–recomposition becomes a blockwise computation using numpy.add.reduceat and blockwise broadcasting. Memory consumption is therefore $G$8, and time per operator is also $G$9. For bounded-degree graph families, where $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$0, the simulation becomes linear in the number of nodes. The paper also notes that even when the update operator $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$1 requires adding edges from node $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$2 to all nodes so that $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$3, this adds only $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$4 entries and preserves overall $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$5 scaling (Ortega et al., 12 Jun 2026).

4. Quantum phase estimation layer

A defining extension of SQWLib beyond earlier simulators is its integrated QPE framework. For a unitary $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$6, the library represents a joint QPE state as

$$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$7

where the phase index $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$8 occupies one tensor axis and each column contains an entire Szegedy walk state. This tensorized representation separates the two computational roles: QPE operations act along the phase axis, while walk operators act within each column or batch (Ortega et al., 12 Jun 2026).

Controlled powers of the walk are simulated without constructing controlled-unitary matrices. In the general implementation, one evolves a batch of walk states in parallel: after each step, the leftmost column is removed, so that after $$\mathcal{H}=\mathbb{C}^N\otimes\mathbb{C}^N =\mathrm{span}\{\;|i\rangle_1|j\rangle_2:i,j=0,\ldots,N-1\;\}.$$9 iterations the column with index $i$0 has received exactly $i$1 applications of $i$2. For the common direct-QPE case with input $i$3, SQWLib uses a simpler method: it stores the successive states $i$4 across the phase axis and then applies inverse QFT numerically.

The Hadamard layer $i$5 and QFT/QFT$i$6 are implemented via Walsh–Hadamard transforms and FFTs along the phase axis. The framework also supports multiple phase registers. With $i$7 phase registers of $i$8 qubits each, the total number of controlled-$i$9 applications is $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$0, rather than $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$1. This construction is used in the approximate reflection operator that underlies the walk-based quantum search algorithm simulated in the package (Ortega et al., 12 Jun 2026).

5. Implemented algorithms and representative results

The paper presents SQWLib as a simulator not only for plain walk dynamics but for full algorithmic protocols built from Szegedy walks and QPE.

Algorithm	Demonstration setup	Reported outcome
Marked-node detection	2D lattice, $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$2, $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$3 or $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$4, $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$5	For $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$6, phase outcome $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$7 has probability $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$8; for $$|\psi_i\rangle=|i\rangle_1\otimes|\omega_i\rangle_2, \qquad |\omega_i\rangle_2=\sum_{k=0}^{N-1}\sqrt{G_{ki}}\;|k\rangle_2.$$9, outcomes concentrate on non-zero phases
Quantum simulated annealing	1D Ising spin glass, $N$00, $N$01, $N$02, 5 temperature steps from $N$03 to $N$04	Marginals match the classical stationary distribution; phase-0 probability at each step is $N$05; joint probability over all 5 steps is $N$06
Quantum search on graphs	Same Ising system, $N$07, $N$08, varying number of phase registers $N$09	Larger $N$10 makes the success-probability curve approach the theoretical amplitude-amplification behavior

For marked-node detection, the underlying transition matrix is modified so that marked vertices become absorbing sinks,

$N$11

The algorithm prepares a uniform superposition over unmarked $N$12 states. If there are no marked nodes, that state is an eigenstate with eigenvalue $N$13, so QPE returns phase $N$14. If marked nodes are present, the state has no support on the eigenspace with eigenvalue $N$15, and QPE yields non-zero phase with high probability given sufficient precision.

For quantum simulated annealing, SQWLib uses the similarity-transformed walk $N$16 built from a Metropolis–Hastings chain. The acceptance probabilities are

$N$17

and the transition matrix is

$N$18

The library simulates the standard QPE-with-postselection protocol across a temperature schedule. In the demonstrated 1D Ising spin-glass instance,

$N$19

the measured marginals agree closely with the classical Boltzmann distributions at each intermediate temperature.

For quantum search on graphs, SQWLib implements approximate amplitude amplification using an oracle

$N$20

together with an approximate reflection around the stationary state $N$21 obtained through QPE. The paper reports that as the number of phase registers increases, the success-probability curve approaches the theoretical one for perfect amplitude amplification; at the first maximum, the measured distribution over marked nodes matches the stationary distribution restricted to the marked subset (Ortega et al., 12 Jun 2026).

6. Relation to SQUWALS, implementation dependence, and limitations

SQWLib is the direct successor to the earlier simulator “SQUWALS: A Szegedy QUantum WALks Simulator,” which already provided a memory-saving classical simulator for the original Szegedy walk and showed $N$22 scaling in both time and memory for dense graphs (Ortega et al., 2023). SQUWALS supported the reflection $N$23, swap $N$24, simple oracles, mixed-state evolution, the semiclassical Szegedy walk, and high-level applications such as quantum PageRank. Its central contribution was the matrix-state representation that replaced $N$25 operator storage with $N$26 array operations (Ortega et al., 2023).

SQWLib extends that earlier framework in four explicit directions. First, it adds operator-level generality by introducing efficient simulation for the primitive operators $N$27, $N$28, and $N$29, which allows simulation of a broader class of Szegedy formulations. Second, it adds a NumPy-based sparse-state representation in the reduced subspace $N$30, bringing time and memory down to $N$31. Third, it integrates a QPE framework, including controlled-$N$32 simulation, multiple phase registers, and approximate reflections. Fourth, it analyzes implementation dependence by providing both $N$33 and $N$34: the reported QSA results are independent of whether $N$35 or $N$36 is used when measurement is included at each step, whereas some heuristic multi-chain behaviors without measurement are implementation-dependent (Ortega et al., 12 Jun 2026).

The package is therefore not a framework for controlling quantum hardware; it is a classical simulator for numerical exploration of Szegedy-walk-based algorithms. The paper also states several limitations. Practical system size remains bounded by memory: $N$37 scaling for dense graphs limits $N$38 to a few thousands depending on the machine, and sparse simulations are limited by $N$39 and by the overhead of reduced-subspace tensors. The implementation assumes column-stochastic transition matrices and is tailored to Markov-chain-derived Szegedy walks. Alternative formulations that strongly exploit the second register may behave differently depending on how $N$40 is defined. Within those constraints, SQWLib is presented as a practical tool for studying Szegedy-walk-based algorithms numerically beyond purely analytical treatments; the code and scripts are available at https://github.com/qDNA-yonsei/SQWLib (Ortega et al., 12 Jun 2026).