Automated detection of symmetry-protected subspaces in quantum simulations

The analysis of symmetry in quantum systems is of utmost theoretical importance, useful in a variety of applications and experimental settings, and is difficult to accomplish in general. Symmetries imply conservation laws, which partition Hilbert space into invariant subspaces of the time-evolution operator, each of which is demarcated according to its conserved quantity. We show that, starting from a chosen basis, any invariant, symmetry-protected subspaces which are diagonal in that basis are discoverable using transitive closure on graphs representing state-to-state transitions under $k$-local unitary operations. Importantly, the discovery of these subspaces relies neither upon the explicit identification of a symmetry operator or its eigenvalues nor upon the construction of matrices of the full Hilbert space dimension. We introduce two classical algorithms, which efficiently compute and elucidate features of these subspaces. The first algorithm explores the entire symmetry-protected subspace of an initial state in time complexity linear to the size of the subspace by closing local basis state-to-basis state transitions. The second algorithm determines, with bounded error, if a given measurement outcome of a dynamically-generated state is within the symmetry-protected subspace of the state in which the dynamical system is initialized. We demonstrate the applicability of these algorithms by performing post-selection on data generated from emulated noisy quantum simulations of three different dynamical systems: the Heisenberg-XXX model and the $T_6$ and $F_4$ quantum cellular automata. Due to their efficient computability and indifference to identifying the underlying symmetry, these algorithms lend themselves to the post-selection of quantum computer data, optimized classical simulation of quantum systems, and the discovery of previously hidden symmetries in quantum mechanical systems.


I. INTRODUCTION
The analysis of symmetry is a central tool in physics and has enabled some of the most profound discoveries in the field.Noether's Theorem famously connects the symmetries of a system's action with conservation laws to which that system's equations of motion are subject [1].Generally, the analysis of symmetry, or the breaking thereof, allows one to constrain theories [2], solve equations of motion more efficiently [3], and identify phases of matter [4].Applications of symmetry analysis in quantum information include, but are not limited to: quantum error correction [5], error mitigation on quantum hardware [6][7][8][9][10], and quantum machine learning model design [11].
Quantum computing can efficiently simulate quantum dynamics in regimes where classical simulation becomes impossible [12].However, current quantum processors operate in a regime severely constrained by noise, with error rates not yet sufficiently below most error correction thresholds [13].Error mitigation will therefore be critical in the interim before fault-tolerant architectures can be scaled [14][15][16].In recent work, and despite its limitations [17], the technique of post-selection has proven useful to mitigate errors and extract useful results from quantum simulation experiments (see e.g., [18][19][20]).Post-selection works by identifying measured states that could only have come from error processes and excluding them from the statistics used to calculate output quantities.The most obvious example is a conserved quantity such as particle number.In such a case, any measured state that does not preserve the conserved quantity must be the result of errors.Due to the connection established by Noether, a more fundamental way to describe post-selection is with respect to symmetry.According to this description, post-selection works by checking the eigenvalue of the simulation's "fiducial" (i.e., initial) state under the symmetry operator against the corresponding symmetry operator eigenvalues (e.g., value of a conserved quantity) of individual measurement results in the dynamically-generated output state of the simulation.If a particular measurement outcome registers a different eigenvalue under the symmetry operator than the fiducial state does, then the measurement is a result of error and can be discarded.This procedure is restricted in the scope of its application, as symmetries of a quantum system and their corresponding operators (i.e., conserved quantities) are typically either engineered into the dynamics "by hand" or identified by clever theoretical intuition.For a generic quantum system, the relevant symmetry operator(s) may not be obvious a priori and may be difficult to identify.Being able to perform postselection in a manner that does not require explicit identification of a symmetry operator would greatly increase the technique's applicability; this is the subject of this paper.
As a corollary, such an operator-free method for error detection also enables additional applications such as more efficient classical simulation of quantum systems via computational basis state reduction.For example, particle number conservation in hardcore boson models can be used both for post-selection in quantum simulation and for reducing the basis state set size from 2 n to n N , where n is the number of lattice sites and N the number of particles, in classical simulations [21].Interestingly, the identification of symmetry or conserved quantities in some instances can make classical simulation so efficient that it can obviate the need for quantum computation altogether [22].Finally, in certain special cases, one may be able to infer the explicit form of a symmetry operator by inspection of the reduced basis set.
In this paper, we provide algorithms to efficiently make use of symmetry in an operator-free manner.To do this, our methods create the subspace of measurement basis states which would share a conserved quantity of some commuting symmetry operator that is diagonal in that measurement basis, without needing to explicitly create that operator.We call such a space a symmetry-protected subspace (SPS).In the language of linear algebra, these are invariant subspaces of the evolution operator; they are subspaces, determined by the initial state, from which the evolution cannot escape.To reap the benefits of symmetry we only need to find the SPS of the initial state, not a conserved quantity, much less an explicit symmetry.However, naively, to find an SPS we need to actually evolve the system in the full Hilbert space, which is exponentially large in the number of qubits (particles, spins, etc.).Sections III through V describe the formulation and algorithms by which we avoid this exponential scaling, but here we provide a non-technical overview.
First, note that most Hamiltonians and resulting Unitary evolution operators are built from a number of local operators.For example, the Heisenberg-XXX model described below consists only of nearest-neighbor interactions.So, at some level, we have an intuition that the dynamics, thus the SPSs, should be derivable, like the Hamiltonian itself, from a combination of local operations, and that local operations are inherently less computationally expensive to work with.This is indeed the case, as shown below.
Next, note that to say that a wavefunction is in a symmetry-protected or invariant subspace is to say that it is and remains throughout dynamic evolution a linear combination of basis states in that subspace and that subspace alone.And if we care only about finding the subspace, we do not need to keep track of the actual linear combination (i.e., both the basis vectors and their amplitudes) but only the basis vectors.This "binarization" of the evolution is critical, because it allows us to adopt a graph-theoretic framework that is vastly more efficient for finding and searching SPSs.
This also leads to an important restriction in our work; using our methods, we can only automate the discovery of symmetry which is diagonal in a chosen basis.We work in the computational Z basis throughout this paper, though extensions to other bases are of course possible by rewriting the time evolution operators in the new basis and performing the same procedure we describe below.Our automated methods should thus be viewed as a tool to find symmetry-protected subspaces within a given basis (if, of course, they exist), and to potentially improve classical and noisy quantum simulations based on those discovered subspaces.However, they still require an intelligently guessed initial basis as a starting point.For the problems we consider in this work, the computational basis is sufficient to derive novel results, though more complex choices can be required in other cases.
With this caveat in mind, once the basis is chosen we create an undirected and unweighted graph, called the state interaction graph, which describes all possible state-to-state transitions over a single application of a unitary evolution operator.The transitive closure of this graph fragments the Hilbert space of the system, represented in a particular measurement basis, into a cluster graph, whose subgraphs are each a symmetry-protected subspace.
Our main results are two classical algorithms that efficiently construct and work within these subspaces.The first, Algorithm IV.1, uses "transitive closure" on local operations to explore and explicitly construct the full SPS of an initial state, which enables the partition of the Hilbert space into a set of disjoint SPSs.This algorithm scales linearly in both the number of local operators from which the global operator is constructed and the size of the SPS, which is a huge improvement over the exponential scaling of the naive "full evolution" approach.However, because the SPS itself can be exponentially large (albeit with an asymptotically smaller prefactor), the second algorithm, Algorithm V.1, finds a path of local operations through a set of SPS graphs from an initial to final (i.e., measured) state to determine if they lie within the same SPS (and thus the final state is valid).This algorithm scales as the number of local operations raised to a small integer power (that can be tuned for accuracy/performance) times the length of the path, thus completely eliminating any exponential scaling.This paper is structured as follows.In Sec.II, we define symmetry-protected subspaces and the three quantum systems that will serve as our "exemplars" throughout.In Sec.III we outline our novel graph theoretical approach to quantum simulations.Section IV discusses the algorithm for computing an entire symmetry-protected subspace.Section V provides a more efficient algorithm to verify if two states exist within the same symmetryprotected subspace.Finally, we demonstrate the effectiveness of symmetry-protected subspaces to mitigate error in quantum simulations on a classical emulator in Sec.VI.

II. PRELIMINARIES
A. Symmetry-protected subspaces Consider a quantum system undergoing unitary evolution according to the operator U (t) for a time t.The operator U (t) can represent continuous time evolution, but also includes other cases, such as discrete time evolution.We say that U (t) is invariant under the action of an operator S if [S, U (t)] = 0 for all times t.In this instance, S is a symmetry operator.For a basis in which S is diagonal, states can be labeled by their eigenvalues under the action of S: S|s, b = s|s, b , where b is some other (set of) label(s), which could represent, for example, a computational basis state integer encoding or a many-body eigenvalue.Suppose now that we initialize the dynamics in a state of definite s: , and evolve under U (t).Given the commutativity of the symmetry and evolution operators, the action of the symmetry operator on the output state is: , indicating that the eigenvalue s is conserved under the evolution for all time.Our methods identify states that would share an eigenvalue under S without explicitly knowing S. To do this, we use the notion of a symmetryprotected subspace, also known as an invariant subspace, which we define below.Note that while our definition emphasizes the connection to symmetries, an SPS, as defined above, is indeed an invariant subspace according to the usual definition [23], which is simply that where we have used the commutativity of P G and U (t) and the fact that by definition the result of applying P G to anything is in G .Another consequence of Def.In Sec.IV we present an algorithm that prescriptively constructs a subspace, denoted G |ψ0 , of a particular initial state, |ψ 0 .We also show in the corresponding theorem, Thm. 1, that subspaces so constructed share the property described in Lem.0.1, which relies on Def. 1 for its proof.Hence, G |ψ0 constructed according to the procedure in Sec.IV are symmetry-protected subspaces.
Similar to that of the symmetry operator S, the commutation relation involving P G also identifies a dynamical invariance, since for any |g ∈ G, P G [U (t)|g ] = [U (t)|g ].That is, the state U (t)|g is an eigenvector of the projector with eigenvalue 1.Therefore, while the projector P G does not identify conserved symmetry eigenvalues, it does indicate when such an eigenvalue exists.In Sec.III, we will show how discovering P G , rather than S directly, empowers our algorithms to discover underlying invariances.

B. Quantum simulations
The general form for unitary evolution operators we assume is where t is the total duration of evolution, O op denotes some operator ordering, such as time ordering, and p is the number of time-steps used to evolve to t.Each U i (τ j ) acts locally on k-qubits and is parameterized by a real time coordinate τ j .Evolution unitaries of the form in Eq. (2) can evolve over discrete time, where each τ j is finite, or discretized-continuous time, where each τ j is (ideally infinitessimally) small in order to minimize Trotter error.
In the instance where a k-local Hamiltonian, H = i h i , is known, Eq. ( 2) results from the local dynamics governed by the h i via the Trotter-Suzuki formula [24], and p corresponds to the number of Trotter steps.We refer to the operator U as the relevant quantum system and examine the dynamics and associated symmetryprotected subspaces of three exemplary systems.

Heisenberg-XXX
The one-dimensional Heisenberg-XXX model for n spin-1/2 particles with nearest-neighbor interactions is given by the following Hamiltonian: where X i , Y i , and Z i are Pauli operator acting on spin i.The model conserves total spin in the Z basis, repre-sented by the operator as well as the correspondingly-defined operators S x and S y .Quantum simulation of the Heisenberg model on a digital quantum processor can be achieved via exponentiation and Trotterization of Eq. ( 3).In such quantum simulation experiments, one usually picks one or a number of qubit bases in which to measure.The symmetry operators S x,y,z can be used to mitigate errors in the allqubit X, Y, Z measurement bases, respectively, via postselection.In the context of classical simulation of the XXX model, the symmetries can be used to constrain the number of basis states included in the dynamics.

T6 quantum cellular automata
The one-dimensional T 6 quantum cellular automata (QCA) rule has recently come to interest within the context of quantum complexity science as a dynamical smallworld mutual information network generator [20] and a QCA Goldilocks rule [25].Its discrete-time unitary update can be derived from a parent Hamiltonian, but it is more natural to define the system by specifying the simulation unitary for a discrete time t = p, directly: where P (α) i = |α i α i | for α = 0, 1 is the projection operator onto the corresponding state of qubit i, H i is the Hadamard operator, and δ α+β,1 is the Kronecker delta function.At each time step, a Hadamard is applied to a qubit only if exactly one of its neighbors is in the |1 state (i.e., α + β = 1) and does nothing otherwise.It has a known Z basis symmetry related to domain-wall conservation: where it should be understood that indices i = 1, . . ., n refer to dynamical, computational qubits while the indices i = 0 and i = n + 1 refer to non-dynamical qubits fixed to the |0 state.

F4 quantum cellular automata
The one-dimensional F 4 QCA with nearest-and next-nearest-neighbor connectivity is another Goldilocks rule [25].It is also most easily specified by its simulation unitary for discrete-time duration t = p: where if the time step index, j, is even (10) Equation ( 10) applies a Hadamard to a qubit if exactly two out of its neighbors or next-nearest neighbors are in the |1 state.There are no analytically known symmetries for this rule.As shown in Fig. 3 in Sec.IV C, our methods discover previously-unknown, symmetryprotected subspaces, indicating a hitherto hidden symmetry of the system.

III. GRAPH THEORY APPROACH TO QUANTUM SIMULATIONS
In this section, we show how graph theory coupled with transitive closure discovers symmetry-protected subspaces.We describe how to "binarize" interactions between basis states through the dynamical system by discarding amplitude and phases to solely highlight if such an interaction exists, describing these interactions as transitive relations, and finding the states connected by a transitive relation to form the SPSs of the system.
Directly, this method still requires unitary matrix multiplication in Hilbert space to establish basis state interactions through the dynamical system.Thus, to apply our methods, we need an efficient way to describe stateto-state interactions.We do this by creating a structure we call a string edit map L U for a unitary operator U , which relies on the observation that quantum systems are typically structured by local interactions.This string edit map returns the basis states available through any operations of a single local unitary operator on a single basis state in near-constant time complexity, allowing us to inexpensively find basis vectors available in local quantum dynamics.This section will proceed as follows: in Sec.III A, we first define the concept of a state interaction graph, whose edges indicate nonzero amplitudes on transitions  We begin by defining the state interaction graph D U for a unitary operator U over a Hilbert space basis B(H d ) and showing how to construct it.For now we leave the form of U general and will specify particular forms when necessary.

Definition 2 (State Interaction Graph). Given a basis
B and a unitary operator U , define a vertex set V ≡ B and an undirected edge set Then, the state interaction graph is defined by the ordered tuple D U ≡ (V, E).In other words, the basis states of B are assigned to vertices (nodes) in the interaction graph and edges are created between vertex states only where the matrix element of the evolution unitary between the two states is nonzero.
Strictly speaking, D U should be a directed graph, where an edge points from |b to |b if b | U |b = 0 and from |b to |b if b| U |b = 0.However, for a symmetry operator S and time-dependent simulation unitary U (t), one can show that [S, U (t)] = 0 ⇐⇒ [S, U † (t)] = 0, meaning that symmetries of the evolution operator, and their associated protected subspaces, are invariant under time-reversal.To remain consistent with this observation, we treat every directed relationship b | U |b = 0 ⇒ (|b → |b ) as an undirected edge (|b ↔ |b ).This treatment is equivalent to assuming that the true, directed state interaction graph corresponding to Def. 2 always has a cycle that leads back to every node, such as is the case in Fig. 1(a).Formally, we assume that if an edge (|b → |b ) ∈ D U , there exists a path {|b → • • • → |b } ⊆ D U .This global cyclicity assumption enables us to use the notion of transitive closure in Sec.III B, and subsequently, in an uncomplicated manner that respects the time-reversal invariance of the resulting subspaces.It is justified on two points: 1.) Most simulation unitaries have a very regular, repetitive structure so that directed, acyclic state interaction graphs are likely only to arise in extremely pathological instances, and 2.) The failure of the global cyclicity assumption will only ever result in the artificial enlargement of a symmetry-protected subspace, and while such a failure leads to underconstrained subspaces, which is bad for the efficacy of e.g., post-selection, it will never result in the corruption of simulation fidelity by overconstraining or throwing out good simulation data.
To construct D U , we use the following steps.First, choose a set of Hilbert space basis vectors B(H d ) = {|b }.Any symmetry-protected subspaces that are discovered must be formed by the basis vectors of this basis.In the context of quantum simulation, B dictates the basis in which a quantum computer will be measured.For example, parallel readout in the computational Z basis will result in bit strings, B comp. = {|0 . . .00 , |0 . . .01 , . . ., |1 . . .11 }, which are Pauli Zstring eigenvectors.In the context of classical simulation, the basis furnishes a representation for the d-dimensional vector of complex amplitudes that stores the evolving many-body wavefunction.When the evolution operator is applied to a basis vector for a single time-step the resulting state |ψ(τ ) has a basis vector decomposition where we suppress the time-ordering subscript j in τ j for simplicity.With this decomposition for each |b ∈ B, we create the state interaction graph D U of the operator U (τ ): for each pair |b , |b such that α b (τ ) = 0 in Eq. ( 11), we add an edge (|b ↔ |b ) to D U .If one has a direct d × d matrix representation of U (τ ) on-hand, then the adjacency matrix for D U can be read off directly as where if an entry in U (τ ) or U † (τ ) becomes non-zero for any value of t, its complex value is replaced by 1 under the operation bit[. ..].This adjacency matrix can be formed, for example, by directly exponentiating a d × d Hamiltonian matrix with a time parameter.In practice, however, constructing or storing an entire evolution unitary in memory is costly, since the size of Hilbert space grows exponentially in the number of qubits, or spins, in the simulation: d = 2 n .Indeed, one of the main advantages of digital quantum simulation is the ability to break global evolution unitaries into sequences of local unitaries, at the expense of introducing error, which are then implemented as quantum gates.Therefore, being able to extract symmetry-protected subspaces from consideration of local operations, rather than from the global unitary they may approximate, is of clear benefit.
Towards this end, Fig. 1(a) shows the state interaction graph D UTrot for one, potentially very bad depending on θ, Trotter approximation to the one-dimensional, fourqubit hopping unitary We take the Trotterization to be where it should be understood that rightmost operators are applied first, and multi-qubit states are ordered as |q 0 q 1 q 2 q 3 .Notice that each iSWAP i,i+1 (θ) is parameterized by an arbitrary θ, and as such we expect each iSWAP operation to be a fractional operation that leaves some residual state behind, i.e., the operator has an identity component.The system is small enough that the state interaction graph can be checked by hand in this case, and the main observation to be made is that it is comprised of four, disjoint subgraphs, each only containing transitions between states of fixed particle number (i.e., number of |1 s), and that all nodes in each subgraph have a path to all other nodes in the subgraph.The zero-and four-particle states are isolated, while the one-, two-, and three-particle states form directed, incomplete, isolated subgraphs.As we will see, this "incompleteness" feature is a pathology of the Trotter approximation which will be rectified in Sec.III B via transitive closure.It is also worth noting that constructing the state interaction graph using the Trotter-approximated unitary is not yet useful, since on a classical computer it currently still requires the storage and evolution of a 2 n -dimensional wavefunction.We will demonstrate the utility of constructing state interaction graphs from component k−local operators in Sec.III C.There are cases where the amplitudes of the states cancel, due to destructive interference, and break this transitive property on the level of individual basis state to basis state interactions.By ignoring the amplitudes of the basis states, we run the risk of including states in the subspace which would be removed via destructive interference.This risk comes with the benefit of efficiently knowing which states are reachable in the quantum simulation, and for our applications it does not add any error to a simulation, as it does not break the underlying commuting subspace projection operator P G ; the subspaces are simply not as restrictive as they could be.See Appendix A for a proof.This also allows one to define the symmetry-protected subspaces to be for any parameterization of the simulation unitary. The for some time/operator exponent t.The transitive closure of an undirected, unweighted graph is a cluster graph, or a set of complete subgraphs; as discussed at the end of Sec.III A, we treat the graph D U as undirected just for this purpose.The transitive closure of the state interaction graph for the Trotterized unitary in Fig. 1(a) can be seen in Fig. 1(b), which in this case turns out to be state interaction graph for the original hopping unitary, or D + UTrot = D UHop .As with D UTrot , within each complete subgraph, total particle number is conserved.
If the edge (|b ↔ |b ) ∈ D + U , then |b and |b share a conserved quantity of the underlying unknown symmetry.This cluster graph structure also makes it apparent that if a wavefunction is initialized as a linear combination of vectors in one subgraph G of D + U : |ψ 0 = b∈G α b |b , it will remain in that subgraph: Therefore, each complete subgraph G within D + U is a symmetry-protected subspace.We will formally prove that transitive closure on the state interaction graph can give an SPS with Thm. 1 in Sec.IV A.
D + U can be represented as a list of disjoint sets of nodes, where each set has implied all-to-all connectivity.This set construction can be seen for the hopping unitary in Fig. 1(c).Here, the index, s, of each subset, G s counts the number of conserved particles.Formally, we can define the associated symmetry-protected subspaces by constructing their projection operators according to Def. 1: P G = b∈G |b b|.

C. Basis state string edit map
We now have a method to identify symmetry-protected subspaces using the language of graph theory.However, actually computing these subspaces still requires the construction and manipulation of vectors and matrices in an exponentially large Hilbert space.Recall, though, that the systems of interest are defined by Hamiltonians composed of local operations U = i U i where each U i is k-local, meaning it only involves k of the n total qubits in the system, and where in general we will have k n.In this section we describe our mechanism for using this fact to build up SPSs efficiently with what we call the basis string edit map.This map enables computation of subspaces using only k-local operations on basis state vectors, so the computational complexity of operation with this map scales with k instead of n.Definition 4 (Basis String Edit Map).Let U i (t) be a unitary operator that acts for a time t non-trivially on k of n qubits, Q k (i) ⊆ {q 0 , . . ., q n−1 }.The basis string edit map L Ui maps a basis state |b to the set of basis states {|b } to which |b can evolve after an arbitrary amount of time under U i (t).Formally, We can apply this construction to any unitary operator, including any k-local Trotter decomposition.Given a unitary operator which is a product of local operators U = i U i , we form the set of local operators used in the Trotter decomposition, {U i }.Then a basis string edit map can be formed for any subset of operators from {U i }, as long as every operator U i ∈ U is included in at least one string edit map.We will use the decomposition which has one string edit map for each local operator in the Trotter decomposition.
Here we have written the operators as acting on the full 2 n -dimensional Hilbert space and highlighted in bold the qubits that are part of the 2 k dimensional subset of this space upon which L iSWAP1,2(θ) acts.This k-local string edit map does not require information about any states besides those at the relevant indices, 1 and 2 in this case.
Add U † i to make the adjacency matrix undirected Algorithm III.1 is computationally trivial to compute for small unitary operators (k = 2 in the iSWAP example), but very expensive for the large unitary operators encountered in quantum simulations.Throughout the rest of this paper, we will use L Ui defined on small k to compute our subspaces in order to keep a small overhead (thus, in our terminology, we will always be talking about "substring edit maps").
The object L lets us create a version of the state interaction graph from Def. 2 which we call the string interaction graph D L ; this construction compactly shows the state-to-state interactions in the algorithms presented in Sec.IV and V.The mapping itself is defined as Notice that, while edges in D U can capture the action of multiple operators at once, each edge in D L is the action of only a single operator; therefore, if D U ≡ (V, E) is the state interaction graph and D L ≡ (V, E L ) is the string interaction graph, then E L ⊆ E. Later, we will show with the proof in Appendix B for Thm. 1 that despite this inequality, the transitive closure of the graph D + L , is equivalent to the transitively closed state interaction graph D + L ≡ D + U .

IV. ALGORITHM: CREATING THE SYMMETRY-PROTECTED SUBSPACE OF AN INITIAL STATE
We now turn to the construction of the entire symmetry-protected subspace of a given initial state using k-local substring edit maps and transitive closure, a construction consistent with the observation that if a symmetry exists locally everywhere in a quantum circuit, then it will also manifest globally [28].The algorithm by which we do so works by establishing a recurrence relation for computing new states, associated by symmetry-protection to the initial state, in the simulation.This recurrence relation furnishes an efficient way to compute the transitive closure of select subgraphs of the entire Hilbert space, with no extraneous information.For a unitary simulation operator decomposed into enumerated local operations, U (τ j ) = m i=1 U i (τ j ), where each U i (τ j ) is k-local, we will show that a symmetryprotected subspace G |ψ0 of the initial state |ψ 0 costs O((m + 1) × |G |ψ0 |) to compute and O(|G |ψ0 |) to store with a breadth-first search [29].We will give the algorithm for the case where the initial state |ψ 0 is a single measurement basis state (i.e., a product state).If |ψ 0 is a linear combination of measurement basis vectors, the algorithm can be repeated for each basis vector in the sum; this does not impact the asymptotic performance of the algorithm, as it only makes the computed subspace bigger.
We enumerate every state in a symmetry-protected subspace by transitively closing subgraphs created by the local basis substring edit maps established in Sec.III C, which return a set of basis states evolved to by their corresponding unitary operators in O(1) when operating on a single basis state.We recursively build the subspace by checking the set of substring edit maps L ≡ {L Ui : U i ∈ U } on each new state, until none are added.This process can be seen as the transitive closure of a subgraph of the graph D L .For U Trot in Eq. ( 13), Fig. 2(a) shows an example of the subgraph of D L corresponding to the action of each L iSWAPi,i+1(θ) ∈ L Trot starting from the initial state |ψ 0 = |1100 (self-edges are ignored as elsewhere in the manuscript).

A. Recurrence relation
To find the symmetry-protected subspace of an initial state, we begin by computing the string edit map for each k-local unitary, L ≡ {L Ui : U i ∈ U }. Next, we check the set of measurement basis strings generated by operating with each substring edit map on the initial state, notated L(|ψ 0 ) ≡ {L i (|ψ 0 ) : L i ∈ L}.We define the set T 1 |ψ0 ≡ {|ψ 0 } ∪ L(|ψ 0 ).Then, for each new state |φ ∈ L(|ψ 0 ), we check operations under the substring edit maps: This process repeats until no new states are found through the following recurrence relation.
The stop condition activates if no new states are found, that is, when additional operations drawn from L do not unveil any new states.Steps 0-4 in Fig. 2 how the recurrence relation manifests for the input state |ψ 0 = |1100 and the set of substring edit maps generated from the Trotterization in Eq. ( 13).We then define the symmetry-protected subspace G |ψ0 to which the state |ψ 0 belongs via where the definition "≡" in Eq. ( 18) should be taken to mean "all basis states in T i+1 |ψ0 viewed as nodes in a complete graph".We can create the Kleene Closure of the set of substring edit maps, denoted L , which is the set of all finite concatenations of substring edit maps, including the identity.Any arbitrary string of substring edit maps L applied to |ψ 0 will result in a state in G |ψ0 by its definition.Hence, one can write We now state our main result.
Theorem 1.Let U (t) = O op p j=1 m i=1 U i (τ j ) be a quantum simulation unitary of duration t as in Eq. (2), divided into p time-steps, where each U i (τ j ) is k-local and is time-step parameterized by τ j , and some operator ordering (such as time-ordering) is specified.Let L Ui be the string edit map corresponding to any available parameterization of U i (τ j ) and the set of such maps L ≡ {L Ui : i ∈ {1, . . ., m}}.Let B(H 2 n ) ≡ {|b } be the basis in which computations (measurements) are being performed classically (quantumly).Then, given an input state |ψ 0 , expressed in the basis B, if |b / ∈ G |ψ0 , where G |ψ0 is constructed according to Eqs. ( 17)-( 18), then b| U (t) |ψ 0 = 0.
For a proof, see Appendix B. Note that for b f | U (t) |ψ 0 to vanish under these conditions, G |ψ0 must satisfy Def. 1 as demonstrated in Lem.0.1.In other words, for G |ψ0 to be able to exclude particular basis states for arbitrary evolution times, it must be a symmetry-protected subspace.Note that Thm. 1 immediately provides two corollaries.
Corollary 1 (Post-Selection).For simulation on an idealized, noise-free quantum computer, if Hence, if the state |b f is measured in the output of a noisy quantum device, it can be assumed that the state arose as a result of error, and may be discarded.
Corollary 2 (Global Subspace).We assumed a Trotterized form for U (t) in the statement and proof of Theorem 1.However, we can formally recover the corresponding global simulation unitary by taking the limit p → ∞ where τ j = jt/p in the time-ordered case and τ j = t/p ∀j when time-ordering is unnecessary (such as when the Hamiltonian is time-independent).Nothing in the proof of Theorem 1 relies upon the finiteness of p or discreteness of the corresponding time differential t/p.Therefore, our result holds for global simulation unitaries as well.This implies that one can reduce the resource requirements in classical simulations of U (t) by only evolving basis states |b f ∈ G |ψ0 .

B. Pseudocode
With an understanding of the recurrence relation in Eq. ( 17) and how it can compute symmetry-protected subspaces, we present an algorithm which can enumerate these subspaces using a breadth-first search.
Algorithm IV.1 uses the set G, which is eventually the symmetry-protected subspace, to track which states have already been added to Q during the runtime of the algorithm and prevent them from being checked more than once.If this set uses the hash of the basis state's bitstrings, insertion and search will be average-case O(1).Getting the set of single-operator transitions , where m is the number of unitary operators in the system, for a single state |b ; while not all m substring edit maps will provide an edge, as many might act as identity on the state |b , each string edit map must still be checked.This set is computed for each state discovered, and each state discovered is never added to the queue Q more than once, which confirms our original complexity analysis of O((m + 1) This algorithm computes a set equivalent to that described by Eq. (17).See Appendix B for a proof that this set is a symmetry-protected subspace.

C. Usage and limits
As stated in Sec.I, post-selection for quantum simulations is the aim of our methods.In order to perform postselection with the algorithm from this section, the entire symmetry-protected subspace G |ψ0 must be known, and each measurement |b f is verified via Because this requires the computation of the entire subspace, it can still become computationally intractable.
As mentioned in Sec.I, we benchmark this algorithm against the Heisenberg-XXX, T6

V. ALGORITHM: VERIFICATION OF A SHARED SYMMETRY-PROTECTED SUBSPACE
Because the worst-case size of a symmetry-protected subspace is still exponential in the number of qubits, Alg.IV.1 is only applicable at relevant system sizes for simulations where G |ψ0 is not exponentially large.Thus, to generalize the usability of symmetry-protected subspaces to any simulation, we present an alternative algorithm in this section that uses an efficient, but greedy, heuristic.Instead of computing the entire exponentially large SPS, this method performs an efficient search for a path of substring edit maps to connect the initial state |ψ 0 and a measured state |b f in the graph D L .The heuristic nature of this algorithm means that it produces only approximate results (albeit with a degree of approximation that can be continuously improved at the cost of more time complexity), and since its runtime scales favorably, it can be used to check the output of quantum simulations well beyond the scale where other classical methods become intractable.
Naive search algorithms in a graph traditionally use a breadth-first search from an initial vertex [30], which is what we described in Alg.IV.1 to enumerate the symmetry-protected subspace through transitive closure; as stated, this is too computationally expensive in many cases.
Our greedy algorithm works as follows: Following all edges (application of the set of substring edit maps L) from a state generates a set of possible states accessible from that state.Working both forward from |ψ 0 and backward from |b f , we have found a path when steps from both directions lead to a shared element.To detect whether this has happened, we use a simple observation about ordered sets: two sets are the same set if they have the same minimal (or maximal) element.A natural order for sets of quantum computational basis states (i.e., binary strings) is just the integer they encode.So, our greedy algorithm works by building two sets (one starting from the initial state and one from the measured state) and comparing them; these sets are subsets of the symmetry-protected subspace corresponding to the initial and final states.If these sets have any common elements (which we check in constant time by looking at their minimal elements), then we have found a path connecting the initial and measured state, and they occupy the same SPS.
In our methods, post-selection's accept/reject decision is made by checking for a measurement result in the symmetry-protected subspace of the initial state.Using the rationale outlined above, we can declare with certainty that two states share a SPS when their paths collide.On the other hand, if the two paths terminate in dissimilar minima we assume the initial and measured states inhabit disjoint subspaces.We say "assume" because the accuracy of our heuristic differs depending on the simulation; when the locally-minimal choice at each step from either state does not build a path to the true minimum element, the two states may build distinct paths while occupying the same SPS.

A. Searching in the string interaction graph
This algorithm builds a single branch, following local minima in depth-limited breadth-first searches, from a single starting vertex in the D L graph, where each vertex |b 's edgeset is given by L(|b ).
Throughout this section we use the integers encoded by the bitstrings b of basis states |b ∈ B(H 2 n ), which provides a natural ordering to states.Let the function min(A) on a set of states A return the state with the smallest encoded integer in that set: To find the minimal element of the symmetry protected subspace, we start with an element |b 0 = |ψ 0 or |b f and build a set towards the minimum of G |b0 , notated min(G |b0 ), with locally optimal decisions.At each step |b curr in the search, we compute T µ |bcurr from Eq. ( 17), which is a breadth-first search to depth µ, or every state in L µ (|b curr ).Then, the starting point for the next step, |b next , is the state with the smallest binary encoded integer in T µ |bcurr , notated min(T µ |bcurr ).Repeat until the set T µ |bcurr does not offer a state smaller than |b curr .Let this process be represented by the recursive function χ: To reiterate, because each step is given by applications of the substring edit maps, χ(|b 0 , µ) and |b 0 must share a symmetry-protected subspace, i.e., χ(|b 0 , µ) ∈ G |b0 .Therefore, when χ(|ψ 0 , µ) = χ(|b f , µ) there must be a sequence of substring edit maps between |ψ 0 and |b f ; i.e., a there is a path When both searches conclude in the true minimal element of their corresponding symmetry-protected subspaces, χ(|ψ 0 , µ) = min(G |ψ0 ) and χ(|b f , µ) = min(G |b f ), we can conclude with certainty that the states do or do not inhabit the same subspace if χ(|ψ 0 , µ) = χ(|b f , µ) or χ(|ψ 0 , µ) = χ(|b f , µ).However, if either search does not conclude in their targeted minimal state, our conclusions can be wrong.Suppose When postselecting with our procedure under these conditions, state |b f would be wrongly rejected because the underlying assumption for this heuristic, that χ(|b 0 , µ) = min(G |b0 ), is wrong.In practice, the true minimal element of the symmetry-protected subspace is unknown; thus, the assertion that χ finds the minimum is always an assumption.
Consequently, if the result of the two searches is a collision, χ(|ψ 0 , µ) = χ(|b f , µ), we know the two states must share a symmetry-protected subspace, even if the searches are at a false minima.If the two searches do not find a common element, we assume the initial and measured states occupy separate symmetry-protected subspaces.As such, measurement outcomes that lie within the intial state's symmetry-protected subspace can be rejected.Thus, each measurement |b f is verified with Formally, we can state that the protected subspace formed by χ around an initial state |ψ 0 is which is approximately equal to the true symmetryprotected subspace, P G ≈ P G .The confidence of this assumption depends on the quantum system being studied, and the depth µ of the local breadth-first searches.Thus, we present arguments to support it for our three exemplar systems: Heisenberg-XXX, T 6 QCA, and F 4 QCA as outlined in Sec.II B. For the Heisenberg-XXX model this assumption is always correct at µ = 1: systems whose substring edit maps L are isomorphic to nearest-neighbor SWAP substring edit maps L SWAP , such as the Heisenberg-XXX model, provably find min(G |b0 ) using χ(|b 0 , µ = 1), as shown in Appendix C 1. We also numerically find that χ(|b 0 , µ = 2) is exact for the T 6 QCA model.Out of our three studied systems the F 4 QCA demonstrates the worst accuracy; χ(|b 0 , µ = 9) still returns false minima, which potentially causes a false rejection.In Sec.V C 1 we analyze these properties of the T 6 and F 4 QCAs more closely.In Sec.VI C we show that post-selection is still effective when χ causes a false rejection to occur.
The computational performance of this method is also model-dependent.Let |Depth| be the depth of our search, or the number of intermediate states traversed to (the second line in Eq. ( 22)) before finding a locally optimal vertex.This means the sequence of substring edit maps between |b 0 and χ(|b 0 , µ) is length O(µ × |Depth|), and the computational complexity of Eq. ( 22) is O(m µ × |Depth|).
Because the reliability and computational complexity of this post-selection method varies, before applying the algorithm to a given quantum simulation, benchmarks should be performed on a small version of the system to understand the µ required to obtain reliable results.

B. Pseudocode
Here we describe the algorithm, Alg.V.1, to find the minimum binary encoded state in the symmetryprotected subspace of a given basis state.This algorithm starts at a state |b 0 and follows the local minima of limited breadth-first searches in an attempt to find the set's minimal element, resulting in χ(|b 0 , µ).
Notice that the innermost while loop is borrowed from Alg. IV.1 to compute each T µ |bcurr , where |b curr is the current best minimum state.This time, the queue Q uses the ordered tuple (|b , η) as its elements, where |b is the current, un-checked, state in the breadth-first search and η is the number of edges (i.e., depth) that |b is away from the current optimal solution |b curr .Each node checked in this loop has O(m) new edges, and we check a maximum depth of µ, giving the inner loop computational complexity O(m µ ).The set T µ |bcurr is checked for every node traversed to in the search, and our search traverses to |Depth| elements, which is how we arrive at an overall complexity of O(m µ × |Depth|) to compute χ(|b 0 , µ).

Algorithm V.1 Greedy pathfinding to the minimum binary encoded state
Require: Path starting state |b0 , substring edit maps L, To post-select with Alg.V.1, run it once for |ψ 0 to get χ(|ψ 0 , µ), again for each measurement |b f to get χ(|b f , µ), and check the results in Eq. (23).
As an aside, if all states computed in the search, T µ |bcurr , are cached with the result of their search, subsequent executions of the algorithm can be preempted with their previously calculated result if a state is found inside the cache.
C. Benchmarking data for Algorithm V.1 This subsection presents benchmark data for Alg.V.1 that demonstrates its practical applicability for our example systems.These benchmarks include heuristic reliability and search depth.We supplement this data with proofs or explicit analytical forms when possible.When this is not possible, we instead rely on extrapolation from the data to inform the asymptotics of these quantities.We find that for the Heisenberg-XXX model, the heuristic is provably exact (Appendix C 1) and the search depth has an explicit equation (Appendix C 2), while the T 6 and F 4 QCA models rely on data-driven intuition for these quantities.

Search reliability
If searching for the minimal element from |b 0 = |ψ 0 or |b f fails, because χ(|b 0 , µ) = min(G |b0 ), then it is possible to falsely assume that two states occupy disjoint subspaces.As such, we examine the reliability of the searching function χ on each exemplar system.
In the Heisenberg-XXX model, the base case for χ, as shown in Eq. (22), only activates at a state which is always the minimum binary-encoded state in the symmetry-protected subspace.This is because the set of substring edit maps for the Heisenberg-XXX system is isomorphic to the set of substring edit maps for the nearest-neighbor SWAP network: Because the T 6 QCA and F 4 QCA do not have SWAPisomorphic substring edit maps, we do not have any available proofs for their exactness, and instead rely on simulating χ for these models.The results of this are in Fig. 4, where for the T 6 and F 4 QCA models we plot the propor- First, notice Fig. 4(a), which shows search failures for the T 6 QCA model; while µ = 1 causes most searches to fail as n increases, µ = 2 causes every search to succeed for every n in the domain.Thus, it we can assume that the heuristic for Alg.V.1 is accurate for this model and that Next, notice Fig. 4(b), which shows search failures for the F 4 QCA model; we use µ = 1, 3, 5, 7, 9, and find that χ can fail even with high µ.We see that µ = 1, 3 sees majority failures over the domain, that µ = 5, 7, 9 has bounded failure ≈ 0.05, and that no µ completely removes failures.
This metric only determines the reliability of χ, not the reliability of post-selecting with χ.This is because, as discussed previously in Sec.V A, two searches can fail and obtain the same false minima, χ(|ψ 0 , µ) = χ(|b f , µ) = min(G |ψ0 ) = min(G |b f ), i.e., two failed searches can still successfully identify that the states share a conserved quantity.Therefore, the search reliability shown in Fig. 4 should be seen as an upper-bound on the error caused by the heuristic.
Sections VI B and VI C explore this reliability further in the context of post-selection on noisy quantum data.In the context of post-selection, we find that even small µ = 3 is acceptable for the F 4 QCA, as the most likely measurement outcomes also share false minima with the initial state.We discuss this observation more with the data that supports it in Sec.VI C.

Search depth
The time performance of Alg.V.1 is heavily dependent on the length of the path taken.The asymptotic worst-case of the search depth is system-dependent, being identified as n 2 /4 for the Heisenberg-XXX model in Ap- For the F 4 model, shown in Fig. 5(b), the search depth becomes shorter as µ increases; this is because higher µ means more edges are traversed at each step in the search, and thus less depth is needed in the search overall.These path lengths, combined with the µ found in Sec.V C 1, results in the computational complexity of Alg.V.1 on a single state being: Heisenberg-XXX is O(mn 2 ), T 6 QCA is O(m 2 n ln(n)), and F 4 QCA is O(m µ n log 10 (n)).The µ exponent is left in the F 4 QCA computational complexity because, as shown in Sec.V C 1, we find no value of µ which guarantees that a search finds the minimal element, and thus expect µ to vary depending on the symmetry-protected subspace occupied for the simulation.As we can see in Fig. 4(b), µ ≥ 5 is mostly sufficient.

VI. POST-SELECTION IN EMULATED NOISY SIMULATIONS
To demonstrate the power of our algorithm, we perform automated post-selection on simulated quantum computations with noise.We perform discretized time evolutions with our exemplar systems, the Heisenberg-XXX, T 6 QCA, and F 4 QCA models, and show that our post-selection methods restructure noisy measurement probability distributions to be closer to the one found in an ideal simulation.We first run the simulation without errors via U (t) |ψ 0 , constructing the measurement distribution at each time step.Next, we repeated the simulations with depolarizing noise injected after each layer of gates with independent probability /3 = 0.005, 0.01, 0.02, 0.05 of X i , Y i , Z i on each qubit i. See Appendix D 1 for more noise model details.We use Kullback-Liebler divergence [31] to quantify the distance between the ideal and noisy measurement distributions constructed at each measurement layer, with and without post-selection using Alg.V.1 and the symmetry-check in Eq. (23).
We observe a significant increase in accuracy of the data when using our post-selection methods without adding circuit runs, even in simulation subspaces which have imperfect pathfinding.

A. Methods
We simulate discrete time evolution, and for each discrete time step a sequence of measurements is used to construct a probability distribution.Let P (p) be the measurement distribution of the wavefunction after p Trotter steps, such that P (p) ≈ || |ψ(p) || 2 up to shot noise.We compare an ideal simulation measurement distribution P ideal (p), a noisy simulation without post-selection P raw (p), and a noisy simulation with postselection using Alg.V.1 at a given µ, P ps,µ (p).Let P (p) be a normalized probability distribution constructed by a sequence of M measurements after p Trotter steps, and P (p, b) ≈ || b|ψ(p) || 2 , which is the probability amplitude of state |b at that layer p, up to shot noise.We use Kullback-Liebler divergence [31], which measures the distance between two probability distributions, and use the following equation from [32] where P IRN is the incoherent random noise probability distribution which we expect to measure once noise has proliferated in the computation.Equation (26) will return 1 if P sim = P ideal , and 0 if P sim = P IRN .We compute Eq. ( 26) at each measurement layer in the computation, with P sim being either the raw data P raw or the post-selected data P ps .We post-select by removing measurement results that fail under Alg.V.1 and Eq. ( 23) and re-normalizing the probability distribution.If P raw (p) is the measurement distribution with noise at step p of M measurements, P ps (p) is the post-selected measurement distribution, given by checking each state in P raw (p) in Alg.V.1, at the same Trotter-step p with measurements M kept ≤ M .

B. Post-selection with perfect symmetry-protected subspaces
We run a 15 qubit simulation of each of the three exemplar models, using initial conditions seen in the liter-ature, the anti-ferromagnetic state |ψ 0 = |1010 . . .101 for the Heisenberg-XXX model, a single bit-flip on the middle qubit |ψ 0 = |0 . . . 1 . . .0 for the T 6 QCA [20,25], and two bit-flips neighboring the center qubit |ψ 0 = |0 . . .101 . . .0 for the F 4 QCA [25].No false minima are seen past µ = 2 in any of the simulations for these initial conditions, meaning post-selection functions effectively.The data uses M = 30, 000 measurements at each measurement layer and does not add any measurements to replace error-victim circuit runs.
Results can be seen in Fig. 6.We observe that postselected data approaches incoherent noise in the SPS of that initial state P ISPS , defined as Because our simulation's raw data converges to random noise, which is a uniform probability distribution over the Hilbert space, and post-selection only removes measurement outcomes, we should expect that the result of post-selection is still a uniform probability distribution, just over the symmetry-protected subspace instead.In other words, we should observe that lim p→∞ F(P ps (p), P ideal (p)) ≈ F(P ISPS (p), P ideal (p)) (29) given enough measurements.

C. Post-selection with incorrect symmetry-protected subspaces
We show that even in subspaces of the F 4 QCA model where Alg.V.1 has failures at high µ, which is equivalent to false-positive mislabeling some basis vectors as outside the symmetry-protected subspace, our methods mitigate more errors than they introduce and thus still show promise for NISQ quantum simulations.
We found 18 symmetry-protected subspaces of the F 4 QCA that had false minima at n = 15 with µ ≥ 9 (henceforth called incomplete subspaces) and ran simulations to p = 29 with M = 10000 measurements at each cycle, initialized in one state from each rocky subspace.See Appendix D 4 b for a list of these initial states.This way, Alg.V.1 is highly likely to meet false minima and reject measurements that were in the symmetryprotected subspace of the initial state.We calculate F(P ps,µ (p), P ideal (p)) for each simulation, and run our post-selection method with µ = 5, 7, 9. See Fig. 7(a) for the set of fidelities {F} at /3 = 0.02, compared to postselection with the "Full SPS", which is post-selection using the full symmetry-protected subspace, instead of the heuristic; in the average case our post-selection keeps the simulation data above incoherent noise, even with lower µ.The worst-case fidelity (in red) drops to the fidelity of incoherent noise, which indicates that our methods only offer improvement on particular simulations.However, the mean (orange) and standard deviation (orange bars) remain reliably higher than incoherent noise, with the maximum fidelity keeping close to the ideal simulation.
The other important observation is that our postselection, which is using an inaccurate symmetryprotected subspace, obtains a higher mean and maximum fidelity than post-selection with a well-defined symmetryprotected subspace, with fidelity lowering as µ increases.We theorize the reason for this is as follows: if two states share a false minima with χ, they are likely close to each other in the string interaction graph D L , and thus have a higher probability of transitions between each other in the simulation.On the other hand, if a measured state does not share a false minima with the initial state, it is less likely to be observed in a simulation, but a noisy simulation will artificially amplify the probability of that state.The result is that post-selection with the true SPS will accept these artificially amplified results, where our false SPS will coincidentally reject states that can possibly have this happen.
We also examine our post-selection in a simulation absent any error; this is to see how much a false rejection can degrade a perfect computation.See Fig. 7(b) for post-selection with an error-free simulation.It can be seen that most fidelities remain above ≈ 0.999, meaning most of the time this error-prone post-selection method will have minimal impact even on an ideal simulation.The worst-case degraded fidelity is ≈ 0.992.

VII. DISCUSSION AND CONCLUSION
We introduce a graph theory interpretation of quantum time evolution, which provides a theoretical framework through which symmetry-protected subspaces can be constructed via transitive closure.We identify that these invariant subspaces are an operator-free method for characterizing symmetry in the system, by indicating a conserved quantity of an initial state as it is manipulated by the dynamical system without needing to explicitly identify it.Along these lines, our approach complements recent work in open quantum systems where symmetryprotected density matrices were constructed without the knowledge of explicit symmetry operators [33].
We observe that a symmetry-protected subspace can be used to provide a smaller computational space, postselect noisy quantum simulation data, or be analyzed to deduce a symmetry operator.We identify post-selection as a pertinent application and introduce two main classical algorithms that elucidate the features of a quantum system's symmetry-protected subspaces.These algorithms employ a basis string edit map, which is an efficient construction to provide the local dynamics of an operator by focusing on the presence or absence of basis vectors through its action.
The first algorithm uses transitive closure, calculated with breadth-first search of these basis string edit maps, to enumerate every state within the symmetry-protected subspace of the initial state.Because post-selection with this subspace requires constructing the entire SPS, which can still be exponentially difficult, we provide a second, polynomial-scaling algorithm.The second algorithm attempts to find the smallest basis state (by binary-encoded integer) in an SPS by following a locally-optimal heuristic, and establishing that a path of basis string edit maps exists between an initial and final state if their respective searches collide at any elements.
We conclude by demonstrating post-selection using our second algorithm, which shows that even when the raw simulation data degrades to incoherent noise, our methods effectively recover a probability distribution closer to an ideal computation.Our methods are compatible with any other error-mitigation technique compatible with post-selection, such as zero noise extrapolation, and thus present a further addition to a growing array of techniques to improve noisy quantum computation [34].
We identify a few obvious extensions to this project.First, if Alg.IV.1 concludes with certain classical resources, this implies that the wavefunction can be stored with constant overhead on those same resources; this could shrink the computational memory requirements from the quantum regime to the classical regime.Second, we speculate that more reliable algorithms than Alg.V.1 may fulfill the same function; a string-matching algorithm similar to Needleman-Wunsch [35] with the basis string edit maps, instead of insert/delete/shift edits, may fulfill this function.
This work also presents interesting results on post-selection that are worth further exploration.The first is the convergence to incoherent noise within a symmetryprotected subspace; this shows that, as expected, postselection will not converge to an ideal computation, regardless of the additional circuit runs to replace errorvictim runs detected, as error will still have proliferated within the the SPS.However, such incoherent distributions still have significant complexity and nontrivial structure (see e.g., [20]), making them particularly interesting in cases where the underlying symmetry is not analytically known.The second is the observation that misidentifying measurements within the SPS paradoxically results in a more accurate computation than what results from full knowledge of the SPS.The exact cause of this is up to speculation, and is likely a problem-dependent effect, but is interesting in its own right.Finally, since the algorithms presented herein are specified with respect to a particular simulation or measurement basis, our methods are currently constrained to discovering subspaces that are protected either by a single symmetry generator in that basis or by a collection of Abelian generators, which are diagonal (or have diagonal representations) in the chosen basis.However, constraining many-body dynamics by sets of non-Abelian generators often yields rich physics as has been noted with respect to non-Abelian thermal states [36], bipartite entanglement entropy growth [37], eigenstate thermalization [38], and thermalization in finite-size quantum simulators [39,40].As such, extending our subspace detection framework to subspaces that are protected by multiple non-Abelian generators would be a fruitful direction for future research., where a gate X( /3) is an X gate applied with probability /3 and identity 1 applied with probability 1− /3.In practice, we simulate this with density matrices through the following steps: if we have the wavefunction's density matrix ρ from the last application of the Trotterized circuit layer U (τ j ), we apply the symmetric depolarizing noise with the equation on each qubit i.As stated in Sec.VI, we use /3 = 0.005, 0.01, 0.02, 0.05.

Heisenberg-XXX
Here we give the simulation details for the Heisenberg-XXX model, with Hamiltonian given in Sec.II B 1 and simulation results in Sec.VI B.

a. Circuit implementation
The k-local unitary used (k = 2) is U i,i+1 = iSWAP i,i+1 (θ)Z i Z i+1 (θ), or , where i is the complex coefficient when multiplied and the qubit when indexed, using θ = 0.1.We then do an even-and odd-layered Trotterization for one discrete layer of the simulation.This Trotterization on n = 5 is Previously, no symmetry was known for the F 4 QCA.Below, we outline the non-identity behavior of each local string edit map L F4i : Any state |b which does not appear above will have the behavior L F4i (|b ) = {|b }.It should also be understood that i is the middle qubit in the ordering above.

b. Initial conditions for simulations
We use a variety of initial conditions to test the imperfect path finding of Alg.V.1.Our initial condition |ψ 0 = |000000101000000 , tested in Sec.VI B, has no false rejections (in other words perfect path finding) with µ = 2.

Definition 1 (
Symmetry-Protected Subspace).Let H d be a Hilbert space of dimension d spanned by an orthonormal set of basis vectors, B(H d ) = {|b } .A subspace G ⊆ H d , which is spanned by a subset of B(H d ), is a symmetry-protected subspace of unitary operator U (t) if and only if a projection onto G, P G = b∈G |b b|, obeys the commutation relation [P G , U (t)] = 0.

FIG. 1 .
FIG.1.Example construction of the state interaction graph DU Trot and corresponding symmetry-protected subspaces for the one-dimensional, four-site hopping unitary in Eq.(12).Multi-qubit states are ordered as |q0q1q2q3 .One Trotter approximation for the unitary, UTrot, is given in Eq.(13), where it should be understood that right-most operators act first.For graphical clarity we omit loops, with the understanding they are always implied.(a) State interaction graph DU Trot .Nodes represent the different 4-bit strings, and edges occur where matrix elements of UTrot are nonzero.Note the presence of five disconnected subgraphs.In practice, we treat each of these edges as undirected.(b) Transitive closure of DU Trot resulting in the closed state interaction graph, which in this case corresponds exactly to the state interaction graph for the hopping unitary: DU Hop = D + U Trot .Each disconnected subgraph in DU Trot has become a disconnected complete graph in D + U Trot .(c) Each complete graph in DU Hop = D +U Trot corresponds to a symmetry-protected subspace Gs.In this example, the subspace indices correspond to eigenvalues under particle number, S = 3 i=0 a † i ai, conservation.

B
. Defining symmetry-protected subspaces with transitive closure Suppose we have states |b and |b , such that b | U |b = 0 and b | U U |b = 0.This requires a transitive relation: b | U U |b = b ∈B b | U |b b | U |b , because |b must first transition to an intermediate state |b to reach its final destination at |b .Therefore, in our state interaction graph D U , as defined in Def. 2, there are edges (|b → |b ), (|b → |b ) ∈ D U , which will have the same transitive relation encoded in the path {|b → |b → |b } ⊆ D U .We use this duality to make the assumption that if the edges (|b → |b ), (|b → |b ) ∈ D U , the transition amplitude b | U U |b = 0.
Algorithm III.1 creates L Ui by taking transitive closure of the adjacency matrix A Ui of D Ui via boolean matrix multiplication [27].This algorithm requires O(2 3k ) time to compute L Ui , O(2 k ) space to store it, and O(1) time to use L Ui .Algorithm III.1 Create a string edit map L Ui Require: Local unitary Ui, orthonormal basis Bi = {|b } Ensure: Dim(Bi) = Dim(Ui) = 2 k and Col(Ui) = Col(Bi)

Definition 5 (
String Interaction Graph).Let D L be a graph defined by an ordered tuple D L ≡ (V L , E L ) and L ≡ {L Ui : U i ∈ U } for some Trotterized unitary U = U i .The vertex set is then given by V L ≡ B(H d ) and the edgeset is given by E L ≡ {(|b ↔ |b ) : |b , |b ∈ B and |b ∈ L(|b )}.

FIG. 2 .
FIG. 2.This figure shows how we can incrementally build the symmetry-protected subspace for U = 2 i=0 iSWAPi,i+1(θ).(a) A single un-closed symmetry-protected subspace, or a subgraph of D L for an iSWAP network.All operations LTrot = {L iSWAP 0,1 (θ) , L iSWAP 1,2 (θ) , L iSWAP 2,3 (θ) } which are non-identity on each node are shown as an edge.Notice the similarities between this graph and the s = 2 subgraph of Fig.1(a): both have the same vertex set, but the edgeset of D L is a subset of the edgeset in DU .(b) Follow the recursion relation in Eq. (17) to iteratively build the symmetry-protected subspace G |1100 , starting from |ψ0 = |1100 .Even though the graph in (a) is not equivalent to the graph in Fig.1(a), their transitively closed graphs are equivalent, as can be seen by the vertices covered by the red line.
Breadthfirst search to enumerate an entire graph (V, E) has computational complexity O(|V | + |E|).In our implementation, there are |G |ψ0 | vertices and we check for m edges at each vertex, giving O(m × |G |ψ0 |) edges in the entire graph.Thus, our breadth-first search to enumerate the symmetry-protected subspace is O(|G |ψ0 | + m × |G |ψ0 |) = O((m + 1) × |G |ψ0 |) Algorithm IV.1 Enumerate symmetry-protected subspace G |ψ0 with a breadth-first search Require: String edit maps L of the simulation operator, initial state |ψ0

FIG. 3 .
FIG. 3. The number of states in symmetry-protected subspaces for each model, on a log 2 scale.Each blue dot is the size of an individual subspace.(a) Heisenberg-XXX.(b) T6 QCA.(c) F4 QCA, which has no known symmetry operator.
Either T pm |ψ0 = T |ψ0 ≡ G |ψ0 , or else T pm |ψ0 ⊂ T |ψ0 .In either case, |b f / ∈ G |ψ0 =⇒ |b f / ∈ T pm |ψ0 , which means that b f |b p,m = 0 ∀ |b p,m ∈ T pm|ψ0 and the full transition amplitude must vanish.with noise, we also added an "error layer" E after each U (τ j ), such that U → U E. It is detailed below in Appendix D 1.

1 .
Simulation noise modelWe use symmetric single-qubit depolarizing noise with probability /3 after each Trotter layer in the circuit.The circuit diagram for this process at an arbitrary Trotter step j is U (τ j ) → U (τ j )

U 2 , 3 2 U 3 U 4 .
(θ) U 3,4 (θ)for a single discrete time-step.We use the initial condition |ψ 0 = |101010101010101 .b. String edit mapsThe string edit maps L HeisXXX for this simulation are SWAP-isomorphic, with each local string edit map L HeisXXXi,i+1 having the behavior:L HeisXXXi,i+1 (|00 ) = {|00 } L HeisXXXi,i+1 (|01 ) = {|01 , |10 } L HeisXXXi,i+1 (|10 ) = {|01 , |10 } L HeisXXXi,i+1 (|11 ) = {|11 } (D2)Each string edit (row in Eq. (D2)) obviously conserves S = i+1 j=i Z j .3.T6 quantum cellular automataHere we give the simulation details of the T 6 quantum cellular automata, with model details in Sec.II B 2 andU (τ j=odd ) =U Notice that gates in the same product in Eq. (D5) and (D6), for example, U 4 and U 7 , commute because only their controls overlap.To get to, for example t = 4, our circuit would look like U (t = 4) = U odd U even U odd U even a.String edit maps (13)ulti-qubit states are ordered as |q0q1q2q3 .One Trotter approximation for the unitary, UTrot, is given in Eq.(13), where it should be understood that right-most operators act first.For graphical clarity we omit loops, with the understanding they are always implied.(a)State interaction graph DU Trot .Nodes represent the different 4-bit strings, and edges occur where matrix elements of UTrot are nonzero.Note the presence of five disconnected subgraphs.In practice, we treat each of these edges as undirected.(b)Transitive closure of DU Trot resulting in the closed state interaction graph, which in this case corresponds exactly to the state interaction graph for the hopping unitary: DU Hop = D + U Trot .Each disconnected subgraph in DU Trot has become a disconnected complete graph in D + U Trot .(c)Each complete graph in DU Hop = D +U Trot corresponds to a symmetry-protected subspace Gs.In this example, the subspace indices correspond to eigenvalues under particle number, S = 3 i=0 a † i ai, conservation.
A. State interaction graph transitive property exists for every state |b in D U , so we can take the transitive closure of the state interaction graph to create a closed state interaction graph D + U .The transitive closure of an edgeset E is a transitively closed edgeset E + , where every pair of states |b , |b ∈ E which can be associated by any transitive relation, in other words can be connected by a path {|b → • • • → |b } ⊆ E, has an edge (|b → |b ) ∈ E + [26].Definition 3 (Transitively Closed State Interaction Graph).Let D U = (V, E) be a potentially non-closed state interaction graph for unitary U and basis B. Define V + ≡ V to be the closed state interaction graph vertex (node) set and E + to be the state interaction graph edge set.An edge, (|b ↔ |b ) ∈ E + , exists in this edge set if and only if there is a path between |b and |b in D U , {|b ↔ • • • ↔ |b } ⊆ E. The transitively closed state interaction graph is then defined as D + U ≡ (V + , E + ).Because the original interaction graph represents single-operator state-to-state transitions, any two basis states which can discover each other through the quantum evolution have an edge in D + U ; in other words, piction of the size of the subspaces, up to 17 qubits.Figure3(c) is especially significant because, as alluded to in Sec.II B 3, this model previously had no known conservation laws, but this data shows the partitioning of Hilbert space into symmetry-protected subspaces.By examining Fig.3, we can see that the worst case of each subspace size, max({|G|}), is log 2 (max({|G|})) ≈ log 2 (|H d |) − k when the model is comprised of k-local operations.