Subspace diagonalization methods differ primarily in the choice of subspace. In classical computing, one of the common approaches is to construct the subspace by generating correlation via local operators such as the hopping terms for fermions, as in multi-reference configuration interaction45<\/a>. Alternatively, one can use global operators. For instance, the classical Lanczos method employs the power series of the Hamiltonian to construct the subspace as \\({{{\\mathcal{K}}}}_{P}={{\\rm{Span}}}\\{{H}^{\\,j}\\left\\vert {\\psi }_{0}\\right\\rangle \\}\\), which is also referred to as the power or polynomial Krylov space. The main advantage of such a construction is that the accuracy of the solution improves exponentially with the subspace size D46<\/a>,47<\/a>,48<\/a>. The limiting factor in classical Lanczos and related methods is that they inevitably suffer from memory consumption that grows exponentially with the system size, owing to the need to represent entangled quantum states.<\/p>\n While various adaptations of this scheme to quantum computers have been proposed13<\/a>,14<\/a>,15<\/a>,16<\/a>,17<\/a>,18<\/a>,21<\/a>,22<\/a>,23<\/a>,24<\/a>,25<\/a>,26<\/a>,27<\/a>,28<\/a>,29<\/a>,30<\/a>,31<\/a>,32<\/a>,35<\/a>,36<\/a>,37<\/a>,38<\/a>,39<\/a>,40<\/a>,41<\/a>,42<\/a>, the most appropriate for near-term quantum computers is to use real-time evolutions as the global operators to generate the Krylov space:<\/p>\n $${{{\\mathcal{K}}}}_{U}={{\\rm{Span}}}\\{{U}^{\\, j}\\left\\vert {\\psi }_{0}\\right\\rangle \\},\\quad j=0,1,\\ldots,D-1,$$<\/p>\n \n (3)\n <\/p>\n where U = e\u2212iH\u2009dt is the time evolution operator for some timestep dt13<\/a>,14<\/a>,15<\/a>,16<\/a>,17<\/a>,18<\/a>,21<\/a>,22<\/a>,23<\/a>,24<\/a>,25<\/a>,26<\/a>,27<\/a>,28<\/a>. The advantage of this is two-fold: first, time evolutions can be approximated by circuits of short enough depth to be implemented on existing quantum devices. Second, one can show that even in the presence of noise, the error due to projection into this unitary Krylov space converges exponentially quickly with the Krylov dimension, just as in classical Krylov algorithms. The noise simply contributes an additional error term as long as it is not so large that it completely overwhelms the signal27<\/a>,28<\/a>. This means that it is possible to reach convergence of the approximate ground-state energy with a Krylov space of limited dimension.<\/p>\n While evaluation of the Krylov matrices on the quantum computer resolves the issue of memory, which is the main obstacle to scaling on the classical side, the main obstacle on the quantum side is noise. Two major contributions are statistical noise due to finite shot sampling and hardware noise in the device. Algorithmic error from the approximation of time evolutions also enters, but below we show numerically that its effects are below the level of the hardware errors. On the other hand, suppressing and mitigating those hardware errors proves to be crucial in order to scale the size of the simulation: we apply experimental techniques for this purpose (see Supplementary Note\u00a04<\/a> for details) as well as keeping the quantum circuit as shallow as possible while maintaining global coupling structure of the Krylov space.<\/p>\n To simplify our circuits, we exploit the U(1) symmetry possessed by many condensed matter models, including the Heisenberg model we focus on. As a qubit operator, U(1) symmetry can be expressed as conservation of Hamming weight; in terms of spin-1\/2 operators, it corresponds to conservation of the z component of total spin. Equivalently, we can think of the symmetry subspaces as k-particle subspaces, treating \u2191(\u2193) spins as absence (presence) of a particle.<\/p>\n Figure\u00a01<\/a> shows a sequence of circuits that could, in principle, be used to calculate the matrix elements (1<\/a>). Panel (a) shows the standard Hadamard test, which would be the default tool for such a calculation. Panel (b) illustrates how we use spin conservation to avoid implementing the controlled time evolutions present in the conventional Hadamard test: instead, we implement controlled initializations of the reference state \\(\\left\\vert {\\psi }_{0}\\right\\rangle\\), and then rely on the fact that the time evolutions preserve the \u201cvacuum state\u201d \\(\\left\\vert 00\\ldots 0\\right\\rangle\\) up to a classically calculable phase.<\/p>\n Fig. 1: Schematic of Krylov quantum diagonalization.<\/b> a<\/b> Hadamard circuit for computing matrix elements of the form \u3008\u03c8i\u2223P\u2223\u03c8j\u3009, which relies on controlled unitary implementation of Krylov basis states. b<\/b> Simplification of the circuit by exploiting a symmetry such as particle-number conservation. c<\/b> The construction employed in this work. Only one time evolution circuit is required, and the second controlled preparation circuit is absorbed into the basis of the measurement. d<\/b> Classical postprocessing to construct matrices \\(\\tilde{H}\\) and \\(\\tilde{S}\\), which yield a generalized eigenvalue problem. The matrices are Hermitian for the circuits shown in a, b<\/b>, and Toeplitz Hermitian for c<\/b>. Note that the diagonal elements, enclosed by black lines, can be computed classically.<\/p>\n As a second simplification, we note that for the exact time evolutions, \\(\\langle {\\psi }_{0}| {U}_{j}^{{\\dagger} }H{U}_{k}| {\\psi }_{0}\\rangle=\\langle {\\psi }_{0}| H{U}_{k-j}| {\\psi }_{0}\\rangle\\), which gives us two formally equivalent ways to measure the same matrix element, with the second yielding a simpler circuit since it only involves one time evolution. However, once the time evolutions are approximated by Trotterization, these two expressions are no longer equal. In Fig.\u00a01<\/a>c, we show the circuit corresponding to the latter version.<\/p>\n It is not a priori clear whether one should prefer the circuits shown in panels b or c in Fig.\u00a01<\/a>, purely from a Trotter error perspective. One advantage of Fig.\u00a01<\/a>b is that it still corresponds to variational optimization in a subspace, since each matrix element still has the form (1<\/a>). However, even this ceases to be true in the presence of finite sample and device noise28<\/a>. Figure\u00a01<\/a>c, the version in which Toeplitz structure is explicitly enforced, is preferable from the perspective of circuit depth for two reasons: it only requires one time evolution, and as a result, the second controlled initialization can be applied as a Clifford transformation to the Pauli observables in the Hamiltonian rather than explicitly implemented in the circuit. In practice, we do not see dramatic violations of variationality with this method, thanks to the regularization technique used to avoid ill-conditioning of the eigenvalue problem (2<\/a>) (see Supplementary Note\u00a05<\/a> for details). As an example, we compare exact classical simulation results for circuits in Fig.\u00a01<\/a>b, c, which are to be compiled in experiments as shown in Fig.\u00a02<\/a> (see the next section for details). Energy curves for a 20-qubit system shown in Fig.\u00a03<\/a> indicate that the variationality restores quickly with the number of Trotter steps. These findings motivated using the version of the circuits shown in Fig.\u00a01<\/a>c.<\/p>\n Fig. 2: Quantum circuits for Krylov quantum diagonalization.<\/b> a<\/b> Each circuit performs the controlled preparation of an initial state within the target particle sector, followed by a Trotterized time evolution. b<\/b> The controlled preparation prepares a computational basis state in which the Hamming weight corresponds to the number of particles for the given experiment, controlled on the auxiliary qubit. Since the heavy-hex lattice can be edgewise three-colored (colors given in the figure by red, green, and blue), both the controlled preparation and the Trotterized time evolution can be implemented using sequences of three unique two-qubit gate layers interleaved with single-qubit rotations. See the main text for details. c<\/b> Each layer of two-qubit gates is Pauli twirled in order to tailor the noise to a sparse Pauli-Lindblad noise model \u039b43<\/a>,44<\/a>, preceded by its amplification \u039bG for PEA. Note that adjacent layers of single-qubit gates, originating from either the source circuit, the twirling, or the noise amplification layer, are always combined in a single layer; they are left unmerged in the figure for clarity. d<\/b> (12\u2009+\u20091)-qubit example of the CZ layers.<\/p>\n Fig. 3: Numerical investigations of algorithmic errors.<\/b> a<\/b> (20\u2009+\u20091)-qubit layout of the Heisenberg model used for numerical simulations, with the green and red circles indicating the control and excited qubits. b<\/b> Energy versus Krylov space dimension. The dotted and solid lines indicate the results from the circuits in Fig.\u00a01<\/a>b, c, respectively. c<\/b> Heat map of the ground-state energy error \u0394E for k\u2009=\u20095-particle sector with various dt and D, using 4\u2009second-order Trotter steps. The white arrow indicates the value of \u03c0\/||H||.<\/p>\n Large-scale experimental demonstrations<\/p>\n For our experiments, we studied the spin-1\/2 antiferromagnetic Heisenberg model, which is defined for a set of edges E as<\/p>\n $$H={\\sum}_{(i,j)\\in E}{J}_{ij}({X}_{i}{X}_{j}+{Y}_{i}{Y}_{j}+{Z}_{i}{Z}_{j})$$<\/p>\n \n (4)\n <\/p>\n with uniform couplings Jij\u2009=\u20091, where Xi,\u00a0Yi,\u00a0Zi denote the Pauli matrices on the ith site. The set of interactions E is a subset of the heavy-hex lattice (see Fig.\u00a04<\/a>). Note that, while the heavy-hex lattice is bipartite and hence the ground state in the entire Hilbert space can be simulated efficiently using the path-integral Monte Carlo method49<\/a>, the sign problem is present for excited states in general. Among the excited states, we focus on the lowest-energy eigenstates within several k-particle subspaces. The dimension of the k-particle subspace scales as O(Nk). Note that the circuit construction relies on the U(1) symmetry but not on SU(2) symmetry, and hence our method is directly applicable to XXZ model as well.<\/p>\n Fig. 4: Experimental diagonalization of many-body Hamiltonians.<\/b> a<\/b> The energy per site of Heisenberg model for particle numbers k\u2009=\u20091,\u00a03, and 5 in system sizes of N\u2009=\u200956,\u00a044, and 42, respectively. The error bars indicate standard deviations estimated by bootstrapping. The dashed curves indicate the energies from noiseless classical simulations, and solid black lines show the exact lowest energy in the given k-particle subspace. b<\/b>\u2013d<\/b> Qubit layout graphs. The green and red circles indicate the control and initial locations of particles, respectively. e<\/b>\u2013g<\/b> Energy curves for individual particle numbers k\u2009=\u20091,\u00a03, and 5. h<\/b>\u2013j<\/b> Error matrices \\(\\Delta \\tilde{H}\/N :=| {\\tilde{H}}_{\\exp }-{\\tilde{H}}_{{{\\rm{num}}}}| \/N\\) and \\(\\Delta \\tilde{S} :=| {\\tilde{S}}_{\\exp }-{\\tilde{S}}_{{{\\rm{num}}}}|\\), where the subscripts \u201cexp” and \u201cnum” denote data from experiments and numerical calculations, respectively.<\/p>\n We ran experiments in three different k-particle sectors: k\u2009=\u20091,\u00a03,\u00a05. The initial states in all three cases were computational basis states with numbers of \\(\\left\\vert 1\\right\\rangle\\) s given by k: for example, in the single-particle case, \\(\\left\\vert {\\psi }_{0}\\right\\rangle=\\left\\vert 00\\ldots 1\\ldots 0\\right\\rangle\\). The circuit implementations for the different values of k therefore differ in the controlled preparation (see Figs.\u00a01<\/a> and 2<\/a>). The k\u2009=\u20091 case corresponds to generating only one particle in the initial state, which can easily be implemented with a CX gate between the control qubit (the ancilla in the Hadamard test) and an adjacent qubit. For k\u2009>\u20091, we chose locations for the particles that were distributed approximately uniformly over the qubit graph.<\/p>\n The heavy-hex lattice permits a three-coloring of its edges, in which each color corresponds to a layer of two-qubit gates that can be implemented simultaneously (see Fig.\u00a02<\/a>). Since each distinct two-qubit layer requires its own noise learning for probabilistic error amplification (PEA\u2014see below), it is advantageous to minimize the number of distinct layers in the circuits. The controlled preparation circuits can be implemented using a set of two-qubit layers corresponding to the three-coloring of edges in the heavy-hex, with only a constant overhead compared to arbitrary layers (see Supplementary Note\u00a02<\/a> for details). For our Trotterized time evolutions, we partitioned the Hamiltonian terms into the same set of layers. Therefore, we only had to learn the noise models of three unique layers in total for each experiment.<\/p>\n The depth of the controlled-initialization part of the circuit is proportional to the distance between the two furthest apart initial particles in the qubit graph. We used two second-order Trotter steps to approximate the time evolutions in all of our experiments. r second-order Trotter steps with three commuting groups of Hamiltonian terms require 4r\u00a0+1 two-qubit layers (see panel b in Fig.\u00a02<\/a>), yielding 9 layers in our case for the time evolution part of the circuit.<\/p>\n To measure the observables corresponding to real or imaginary parts of the matrix elements in \\(\\tilde{H}\\) and \\(\\tilde{S}\\), we partitioned the observables into as few locally-commuting sets (measurement bases) as possible, since such sets are co-measurable7<\/a>. The shortened circuits, as in the third row of Fig.\u00a01<\/a>, require conjugating the Hamiltonian terms (4<\/a>) by the second controlled-initialization circuit, since it is not physically implemented. This yields the same number of Pauli observables since the controlled-initialization is a Clifford circuit, and one can prove that these observables can be partitioned into 2(k\u2009+\u20092) measurement bases; see Supplementary Note\u00a02<\/a>.<\/p>\n We performed experiments on the Heron R1 processor IBM_montecarlo. This is a 133-qubit device with fixed-frequency transmon qubits connected to each other via tunable couplers. Heron processors have faster two-qubit gates (similar in duration to the single-qubit gates) and lower cross-talk compared to the fixed-coupling devices of earlier generations. To further improve the measured observables (see Fig.\u00a01<\/a>), we used probabilistic error amplification (PEA)44<\/a> and twirled readout error extinction (TREX)50<\/a>, which mitigates SPAM errors, to approximate noise-free expectation values. We additionally employed error suppression, in particular Pauli twirling and dynamical decoupling. Details of the error mitigation and suppression are given in Supplementary Note\u00a04<\/a>.<\/p>\n In our experiments, for each measurement basis, a certain number of twirled instances were generated, and each instance was then repeatedly measured for different values of the noise amplification factor. For the single-particle (k\u2009=\u20091) experiment, we used 300 twirled instances with 500 shots each, at noise amplification factors of 1,\u00a01.5,\u00a03. For k\u2009=\u20093,\u00a05, we used 100 twirled instances with 500 shots each, at noise factors 1,\u00a01.3,\u00a01.6. The reduction in twirled instances for the larger experiments was introduced in order to reduce the total runtime, since the number of measurement bases as well as the circuit sizes increase with k. The adjustment of the noise amplification factors was due to the increased noise rates in the deeper circuits. The controlled-initialization part of the circuit involves creating a maximally entangled state of the control qubit and the initial particle locations. With an increase in the number of particles, this translates to a larger maximally entangled state prepared at the beginning of the circuit, which in turn makes the results more susceptible to noise.<\/p>\n The size of the Krylov space was fixed to D\u2009=\u200910 across all experiments in order to achieve a total runtime of the algorithm within the timescale of device recalibration processes (24\u2009h). For a fixed value of k the experiment was run on a specific qubit subset, chosen according to the current status of the device by using a heuristic routine for optimal qubit mapping51<\/a>. The k\u2009=\u20091 experiment was executed on a 57-qubit subset, the k\u2009=\u20093 experiment on a 45-qubit subset, and the k\u2009=\u20095 experiment on a 43-qubit subset (the layouts are shown in Fig.\u00a04<\/a>). The latter two were partially chosen by hand in order to have five complete heavy hexes in each case.<\/p>\n Although the time step dt theoretically has an optimal value of \u03c0\/||H||27<\/a>,28<\/a>, the restriction to low-particle-number subspaces alters this. Consequently, we chose the time steps heuristically, with values 0.5, 0.022, and 0.1 for k\u2009=\u20091, 3, and 5, respectively.<\/p>\n Results are shown in Fig.\u00a04<\/a>. Panel a summarizes the results on a normalized energy scale, while e, f, and g show the convergence curves for each separate experiment. The corresponding qubit graphs are shown in panels b, c, and d, respectively. These convergence curves are a useful diagnostic tool for assessing the results of noisy KQD experiments. We know from the theoretical analysis that if error rates are low enough to resolve the signal, i.e., to distinguish the lowest energy state in the Krylov space from pure noise, then we should see an exponential decay of the energy towards a value offset from the true ground-state energy by a constant depending on the error rate27<\/a>,28<\/a>. Our results show this behavior up to some fluctuations, which is expected since the theoretical results only provide for an exponentially-decaying upper bound. If the noise had completely dominated the signal, however, the rate of convergence with subspace dimension would have been exponentially slow with respect to system size, and we would not have seen the initial fast convergence in Fig.\u00a04<\/a>. See Supplementary Notes\u00a04<\/a> and 5<\/a> for further details.<\/p>\n In our experimental results, noise and algorithmic error (due to the Trotter approximation as well as the limited Krylov dimension) are still significant limiting factors, as evidenced by the differences between the most accurate estimated energies (at D\u2009=\u200910) and the true values. We estimated standard deviations for our experimental energies using bootstrapping, since the post-processing of solving the regularized, generalized eigenvalue problem (2<\/a>) makes direct error propagation difficult. This yielded the error bars in Fig.\u00a04<\/a>; for further details, see Supplementary Note\u00a05<\/a>. Figure\u00a04<\/a> also shows the energy convergence curves for ideal classical simulations of our circuits, which are tractable by representing vectors and operators only in the restricted particle-number subspaces. While the error bars are large due to the noisy experimental results, our estimated energies for the two larger values of k are consistent with the ideal simulation curves up to these standard deviations at nearly all points.<\/p>\n In the k\u2009=\u20091 experiment, the results deviate below the true lowest energy, indicating that noise has created an effective leakage out of the k\u2009=\u20091 subspace. This illustrates a risk of relying on symmetry conservation to remain in a particular subspace, although studying the global ground state would not be subject to this concern.<\/p>\n Exact diagonalization can also be carried out in the sectors of the Hilbert space studied in the present experiments, though not in the full Hilbert space. However, the experiments did not depend on those particular particle number sectors in any way except for the reduced circuit depth of the controlled initialization, so there are not qualitative or structural obstacles to scaling, only effects of noise. In the specific case we focused on\u2014the ground states of the Heisenberg model on a 2D heavy-hexagonal lattice\u2014it is also still possible to compute precise approximations using classical techniques such as tensor networks.<\/p>\n One may ask why KQD was employed rather than one of the various other algorithms that have been recently developed for ground-state energy estimation in a near-term or early-fault-tolerant setting, e.g.,3<\/a>,4<\/a>,5<\/a>,52<\/a>. One primary reason, in addition to KQD\u2019s relatively well-understood noise tolerance21<\/a>,25<\/a>,27<\/a>,28<\/a>, those alternative methods all extract eigenenergies from the time evolution rather than directly from a projection of the Hamiltonian itself. This is a problem in a setting such as ours where the Trotter circuit is held fixed as the number of timesteps is increased (which was necessary to minimize circuit depth), because the spectra of the Trotter circuits diverges from the spectra of the ideal time evolutions, indeed becoming periodic with a period depending on the\u00a0timestep and the fixed Trotter circuit. Hence, with this constraint\u00a0algorithms depending only on the evolution will at some point cease to converge as the number of timesteps is increased.<\/p>\n","protected":false},"excerpt":{"rendered":"Theory of Krylov quantum diagonalization KQD consists of two main steps. The first is a quantum subroutine to…\n","protected":false},"author":3,"featured_media":12869,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[25],"tags":[10046,10047,492,13632,836,159,67,132,68],"class_list":{"0":"post-12868","1":"post","2":"type-post","3":"status-publish","4":"format-standard","5":"has-post-thumbnail","7":"category-physics","8":"tag-humanities-and-social-sciences","9":"tag-multidisciplinary","10":"tag-physics","11":"tag-quantum-information","12":"tag-quantum-physics","13":"tag-science","14":"tag-united-states","15":"tag-unitedstates","16":"tag-us"},"share_on_mastodon":{"url":"https:\/\/pubeurope.com\/@us\/114742721554975869","error":""},"_links":{"self":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/12868","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/comments?post=12868"}],"version-history":[{"count":0,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/12868\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media\/12869"}],"wp:attachment":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media?parent=12868"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/categories?post=12868"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/tags?post=12868"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}![\"figure](\"https:\/\/www.europesays.com\/us\/wp-content\/uploads\/2025\/06\/41467_2025_59716_Fig1_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/us\/wp-content\/uploads\/2025\/06\/41467_2025_59716_Fig2_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/us\/wp-content\/uploads\/2025\/06\/41467_2025_59716_Fig3_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/us\/wp-content\/uploads\/2025\/06\/41467_2025_59716_Fig4_HTML.png\")
<\/a><\/p>\n