# Fileset

[PhysRevResearch.6.043050.pdf](https://mdr.nims.go.jp/filesets/9c9ef4ff-ebed-4375-bf0f-4edf8e82ba9d/download)

## Creator

Ryosuke Shibukawa, [Ryo Tamura](https://orcid.org/0000-0002-0349-358X), [Koji Tsuda](https://orcid.org/0000-0002-4288-1606)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Boltzmann sampling with quantum annealers via fast Stein correction](https://mdr.nims.go.jp/datasets/a25453bb-0b6f-4365-8fec-3d99a1ce0438)

## Fulltext

Boltzmann sampling with quantum annealers via fast Stein correctionPHYSICAL REVIEW RESEARCH 6, 043050 (2024)Boltzmann sampling with quantum annealers via fast Stein correctionRyosuke Shibukawa,1 Ryo Tamura ,1,2,3,* and Koji Tsuda 1,2,3,†1Graduate School of Frontier Sciences, The University of Tokyo, Chiba 277-8568, Japan2Center for Basic Research on Materials, National Institute for Materials Science, Ibaraki 305-0044, Japan3RIKEN Center for Advanced Intelligence Project, Tokyo 103-0027, Japan(Received 8 September 2023; accepted 30 September 2024; published 21 October 2024)Despite the attempts to apply a quantum annealer to Boltzmann sampling, it is still impossible to performaccurate sampling at arbitrary temperatures. Conventional distribution correction methods such as importancesampling and resampling cannot be applied, because the analytical expression of sampling distribution isunknown for a quantum annealer. Stein correction [Q. Liu and J. Lee, in Proceedings of the 20th InternationalConference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research (PMLR, 2017),Vol. 54, pp. 952–961] can correct the samples by weighting without the knowledge of the sampling distribution,but the naive implementation requires the solution of a large-scale quadratic program, hampering usage in prac-tical problems. In this article, a fast and approximate method based on a random feature map and exponentiatedgradient updates is developed to compute the sample weights and is used to correct the samples generated byD-Wave quantum annealers. In benchmarking problems, it is observed that the residual error of thermal averagecalculations is reduced significantly. If combined with our method, quantum annealers may emerge as a viablealternative to long-established Markov chain Monte Carlo methods.DOI: 10.1103/PhysRevResearch.6.043050I. INTRODUCTIONBoltzmann sampling of the Ising model is central inthe studies of critical phenomena [1,2] and machine learn-ing [3,4]. Although Markov chain Monte Carlo (MCMC)methods can generate samples according to Boltzmann dis-tribution, they often fall short for large models due to slowmixing [5]. To efficiently perform Boltzmann sampling, vari-ous improvement techniques such as exchange Monte Carloand population annealing have been proposed in statisticalphysics [6–8].As an alternative, quantum annealers (QAs) [9] have beenexpected to work as a means to achieve accurate Boltz-mann sampling [10–15]. Theoretically, it has been shown thatthe distribution of quantum annealing samples deviates fromBoltzmann distribution [16]. Nevertheless, scientific discus-sion is still unsettled about whether QA samples can be usedas samples of a Boltzmann distribution in practical terms. Re-cently, Nelson et al. argued that the D-Wave quantum annealerworks as an accurate sampler at certain temperatures [10], butit does not work well at arbitrary temperatures.To accurately obtain the target distribution, conventionaldistribution correction techniques such as importance sam-pling and resampling [5] are the powerful tools when theoriginal distribution is known [17]. However, the sampling*Contact author: tamura.ryo@nims.go.jp†Contact author: tsuda@k.u-tokyo.ac.jpPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.distribution by QAs has not been analytically described, andthese conventional approaches cannot be applied to achieveBoltzmann sampling at arbitrary temperatures. To solve thisproblem, a correction method without knowing the originaldistribution is needed when only the target distribution (inthis case, Boltzmann distribution with arbitrary temperature)is given.Liu and Lee proposed a “black-box” distribution correctionmethod based on Stein statistics, where the analytical form ofthe original distribution is not needed [18]. The original sam-ples are assigned the weights to fit to the target distribution viaquadratic programming. In the original paper, the theoreticalproperties are largely unsolved, but Hodgkinson et al. showedthe convergence of Stein correction for samples generated bya Markov chain [19].This study investigates how well the Stein correction worksin the distribution correction of QA samples to a Boltzmanndistribution with a given temperature. First, we develop a fastapproximate algorithm of Stein correction, because the O(n3)computational cost of quadratic programming for the numberof samples, n, is prohibitive for a large number of samples.Next, in benchmarking studies, we observe that the estimationerror of internal energy, magnetic susceptibility, and Bindercumulant of some Ising models decreases in a large extent byStein correction. This result implies that Stein correction isuseful for improving the sample quality for applications suchas critical phenomena and machine learning.II. METHODA. Boltzmann samplingWe are engaged in sampling from the Boltzmann dis-tribution p(x) ∼ exp[−βHIsing(x)], x ∈ {−1, 1}d , where the2643-1564/2024/6(4)/043050(7) 043050-1 Published by the American Physical Societyhttps://orcid.org/0000-0002-0349-358Xhttps://orcid.org/0000-0002-4288-1606https://ror.org/057zh3y96https://ror.org/026v1ze26https://ror.org/03ckxwf91https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevResearch.6.043050&domain=pdf&date_stamp=2024-10-21https://doi.org/10.1103/PhysRevResearch.6.043050https://creativecommons.org/licenses/by/4.0/SHIBUKAWA, TAMURA, AND TSUDA PHYSICAL REVIEW RESEARCH 6, 043050 (2024)Hamiltonian is described asHIsing(x) = −∑i, j∈EJi jxix j −∑i∈Vhixi. (1)β denotes the inverse temperature, V ⊆ [1, d], and E ⊆[1, d] × [1, d]. Here, we assume that the parameters Ji j andhi are in the range of −1 and 1. The thermal average ofobservable O(x) is defined by〈O(x)〉β = TrO(x) exp[−βHIsing(x)]Tr exp[−βHIsing(x)]. (2)Since the trace calculation is impossible for large models,this trace is approximately replaced to the average of somesamples. When βHIsing(x) is solved by QAs with short an-nealing time, the distribution of samples is more widespreadaround the ground state. If we consider that this distributionis the Boltzmann distribution and n samples (x1, . . . , xn) aregenerated by QAs, the thermal average is calculated as〈O(x)〉QAβ = 1nn∑i=1O(xi ). (3)However, the distribution by QAs deviates from the Boltz-mann distribution, and the correction is needed.B. Kernelized Stein discrepancyStein discrepancy can measure the difference between twodistributions p(x) and q(x). Kernelized Stein discrepancy [20]is defined asS(p, q) = Ex,x′∼q[kp(x, x′)] = Trkp(x, x′)q(x)q(x′), (4)where kp(x, x′) is called the Stein kernel that depends on p(x)and the base kernelk(x, x′) = exp(−∑di=1 I{xi �= x′i}d). (5)The function of I{xi �= x′i} shows 1 for xi �= x′i and 0 for theothers, respectively. Let us define the difference operator as∇x f (x) = [ f (x) − f (¬1x), . . . , f (x) − f (¬d x)], (6)where ¬i denotes the sign flip of the ith variable. In addition,the score function sp(x) ∈ Rd is defined as[sp(x)]i = 1 − p(¬ix)/p(x) (i = 1, . . . , d ).The Stein kernel is then defined askp(x, x′) = sp(x) k(x, x′)sp(x′) − sp(x) ∇x′k(x, x′)− sp(x′) ∇xk(x, x′) + Tr[∇x,x′k(x, x′)]. (7)Notably, kp(x, x′) depends on p(x) only through the scorefunctions. Here, S(p, q) is always non-negative becausekp(x, x′) is a provably positive definite kernel if k(x, x′) is apositive definite kernel. In addition, S(p, q) = 0 holds if andonly if p = q [18].To address the relationship between the probability distri-butions and their score functions, we consider the three-bitIsing model with random interactions. Figure 1 shows theBoltzmann distributions of the same Ising model with dif-ferent temperatures as β = 0.5, 0.4, and 0.05 and their scorefunctions. When the distributions are similar, their score func-tions are close, while when the distributions are dissimilar,their score functions diverge significantly. The Stein discrep-ancies are also shown in Fig. 1, and we have confirmed thatthe relationship between distributions can be correctly evalu-ated by the Stein discrepancy. Thus, we can understand thatminimizing the difference in the score function minimizes thedifference between distributions. Stein correction follows thisidea of minimizing the difference in score functions.C. Stein correctionIn this study, we use the kernelized Stein discrepancy tocorrect the contributions of each sample obtained by QAs. Letp(x) be the Boltzmann distribution p(x) ∼ exp[−βHIsing(x)]and q(x) be an unknown distribution by QAs. Suppose thatn samples are obtained from QAs. However, since q(x) isunknown, the weight of each sample cannot be obtained.Thus, we assume that the weight of each sample is wi insteadof q(xi ), and the Stein discrepancy defined by Eq. (4) isapproximated asS(p, q) n∑i, j=1wiw jkp(xi, x j ) = w Kpw, (8)where Kp is the n × n matrix where elements are kp(xi, x j ),which is called the Stein kernel matrix. In addition, w isthe weight vector for n samples. In Stein correction [18],the weights are adjusted such that the Stein discrepancy isminimized,ŵ = argminw{w Kpw s.t. wi � 0,n∑i=1wi = 1}. (9)Then, the samples weighted by ŵ approximate a Boltzmanndistribution p(x). Using the obtained weights, the thermalaverage of observables O(x) is approximately obtained by〈O(x)〉SCβ =n∑i=1ŵiO(xi ). (10)D. Fast Stein correctionA naive implementation of Stein correction requires O(n2)space and O(n3) time. We reduce the complexity by in-troducing a random feature map [21] and an exponentiatedgradient descent [22]. Using the random feature map, thebase kernel defined by Eq. (5) is approximated as the innerproduct k(x, x′) ≈ φ(x) φ(x′), where the feature map φ(x) :{−1, 1}d → R� is computed as follows. Let us draw � samplesω1, . . . ,ω� fromh(ω) =d∏i=11π(1 + ω2i) , (11)where ω is the d-dimensional vector with each component ωi.Also, b1, . . . , b� are sampled from the uniform distributionover [0, 2π ]. Then, the feature map is defined asφ(x) =[zω1,b1(x + 12d), . . . , zω�,b�(x + 12d)] , (12)043050-2BOLTZMANN SAMPLING WITH QUANTUM ANNEALERS VIA … PHYSICAL REVIEW RESEARCH 6, 043050 (2024)FIG. 1. Relationship between Boltzmann distributions and their score functions in three-bit Ising model with different temperatures asβ = 0.5, 0.4, and 0.05. The Stein discrepancies between distributions are also shown.where zω,b(x) =√2�cos(ω x + b) and 1 is the d-dimensionalvector of (1, . . . , 1). Since the Stein kernel is a linear functionof the base kernel, it can also be approximated as kp(x, x′) ≈φp(x) φp(x), where φp(x) is the concatenation of the follow-ing vectors:θk (x) = p(¬kx)p(x)φ(x) − φ(¬kx) (k = 1, . . . , d ). (13)Using the random feature map, the Stein discrepancy can beapproximated asS(p, q) ∥∥∥∥∥n∑i=1wiφp(xi )∥∥∥∥∥2, (14)which is strictly positive. The optimization problem definedby Eq. (9) is a convex optimization problem with non-negativity and normalization constraints. When the standardgradient descent algorithm is applied, the constraints are vi-olated every time the parameters are updated. In this case,exponentiated gradient descent is known to work well [22],because constraint violation never happens. The update isdescribed aswt+1,i = wt,i exp(−η[∇ f (wt )]i)Zt, (15)wheref (wt ) =n∑i=1wt,iφp(xi ), (16)Zt =n∑i=1wt,i exp{−η[∇ f (wt )]i}, (17)and η is the learning rate. The modification shown above re-duces the space requirement to O(n�). Each update takes onlyO(n) time, enabling us to deal with a large number of samples.The implementation of a fast Stein correction can be found onGitHub (see Ref. [23]). We have confirmed that approximatelythe same accuracy is obtained with the exact calculation ofthe Stein correction by our approximate algorithm (see theSupplemental Material, Fig. S1 [24]). Furthermore, we havemeasured the computation time as a function of the number ofbits in the problem, d , with the number of samples fixed at 2and the number of random feature maps � = 5000. The resultsare shown in Fig. S2 in Ref. [24], and the O(d2) computationtime is required for the problem size.III. RESULTSIn our experiments, we employ 16-bit Ising Hamiltoniansproposed by Nelson et al. [10] for benchmarking. They are043050-3SHIBUKAWA, TAMURA, AND TSUDA PHYSICAL REVIEW RESEARCH 6, 043050 (2024)FIG. 2. Results for Stein correction on GSD_8. (a) Computational time of exact and fast Stein corrections depending on the number ofsamples n when � = 5000. The inset is the log scale figure. (b) Residual error of the Stein kernel matrix Kp against the number of randomfeatures � when n = 1000. Five independent runs were performed, and the means and standard deviations are plotted as lines and error bars,respectively. The target Stein kernel matrix is common and random features are different towards five runs.called GSD_X and the number X indicates the number ofdegenerated ground states. The Hamiltonians are designed tofit the chimera topology of D-Wave 2000Q systems, but herethe samples are generated by Advantage System 6.2, because2000Q is already out of service in their cloud platform. Hamil-tonians are embedded in the Pegasus topology of Advantageusing the MINORMINER PYTHON package [25].First, we evaluate the efficiency and approximation errorof fast Stein correction using GSD_8. Throughout the paper,the learning rate is η = 10−5 and the number of updates is3000 in the fast Stein correction. In addition, an exact cor-rection was performed by CVXOPT [26]. Figure 2(a) showsthe computational time of Stein correction depending on thenumber of samples when � = 5000. The fast correction wasorders of magnitude faster than the exact one against thenumber of samples. Figure 2(b) shows the residual error of theStein kernel matrix Kp, which is defined as ‖Kp − K∗p‖/‖K∗p‖,depending on the number of features �. Here, K∗p is the exactmatrix. The residual error by the random feature map de-creases as the number of features is increased, showing thatthe more random features are preferred for accurate approxi-mation. In the following analysis, calculations are performedwith � = 5000, indicating a sufficiently better residual errorand taking computational time into account. Note that residualerrors depending on the system size d are shown in Fig. S3in Ref. [24], and when the number of features is larger than20 000, the residual error becomes efficiently small up to 1024bits.Next, we observe how fast Stein correction improves theaccuracy of estimating observables on GSD_8, GSD_38, andGSD_F_6. For GSD_F_6, the finite fields hi are imposed. Thecalculated thermal averages of observables are the internalenergy E (β ), the magnetic susceptibility χ (β ), and the Bindercumulant U4(β ), defined byE (β ) = 〈HIsing(x)〉β, (18)χ (β ) = β〈(∑di=1 xi)2〉β, (19)U4(β ) = 1 − 〈(∑di=1 xi)4〉β/3〈( ∑di=1 xi)2〉2β, (20)where d is the system size. For each observable, we computethe residual error as ‖y − y∗‖/y∗, where y is the thermal aver-age of observables calculated by Eqs. (3) and (10) and y∗ is theexact value computed via brute-force enumeration. For eachβ, 10 000 samples are generated by a D-Wave QA by settingnum_reads = 10 000 with an annealing time of 5 µs. For thissetting, when β is increased, the number of unique samplesis decreased in 10 000 (see Fig. S4 in Ref. [24]). In addition,the Metropolis method with single-spin flip, a basic MCMCmethod, has also been applied. In our implementation, theinitial state is randomly generated at each temperature, and theMCMC calculation is quenched at the target temperature. Inthe MCMC calculation, we restart the Markov chain when thelength of the chain reaches 2000. This corresponds to the factthat a single MCMC run consists of 5 independent Markovchains. To match the number of QA samples, the first 1500samples are used for burn-in and the remaining 500 samplesare used for estimation. Computing time comparisons for eachinverse temperature between Stein-corrected QAs and MCMCcalculation are shown in Fig. S5 in Ref. [24].The results over a range of inverse temperatures are shownin Fig. 3. For each case, five independent runs are conductedand the mean values of residual error are evaluated. Theaccuracy of fast Stein correction outperforms the originalQA samples in all cases, showing the effectiveness of ourapproach. The temperature dependence for internal energy be-haves nonmonotonically, and there is a temperature at whichthe residual error reaches a minimum. If appropriate sam-ples are obtained from D-Wave QAs for each temperature,the results by Stein correction will be more accurate. Thus,we consider that the temperature with the minimum error isclose to the effective temperature of the D-Wave QA, whichis considered to give the best sampling. On the other hand,the accuracies by the Stein correction and QA samples arelow at lower inverse temperatures (high temperatures) in theinternal energy. Even if high temperatures are set, the low-energy states are obtained by QAs and these are not importantsamples in the thermal averages. Thus, if the temperature is farfrom the effective temperature of D-Wave QAs, the accuracycan be expected to improve by adding appropriate non-QA043050-4BOLTZMANN SAMPLING WITH QUANTUM ANNEALERS VIA … PHYSICAL REVIEW RESEARCH 6, 043050 (2024)FIG. 3. Residual errors of internal energy E (β ), magnetic susceptibility χ (β ), and Binder cumulant U4(β ). The thermal averages werecalculated by MCMC, a naive average by quantum annealers (QA) and Stein correction (SC), respectively. The number of samples was fixedas n = 10 000, and the results depending on n are shown in Fig. S6 in Ref. [24] when β = 0.4. Five independent runs were performed, and themeans and standard deviations are plotted as lines and shaded areas, respectively.samples. In addition, the residual errors depending on thenumber of samples when β = 0.4 are summarized in Fig. S6in Ref. [24], resulting in the case that the error does not alwaysdecrease with more samples. One reason is that the numberof unique samples is smaller and smaller as the temperatureis decreased. In fact, many samples are the same in the QAsamples (see Fig. S4 in Ref. [24]). In particular, errors onmagnetization and the Binder cumulant are not decreased asthe number of samples is increased. This suggests that expec-tation calculation of magnetization and the Binder cumulant isessentially more difficult and sensitive to samples compared tothe internal energy.Furthermore, the error by fast Stein correction was com-parable to that by MCMC except for the internal energy athigh-temperature regions. These results show that fast Steincorrection has the potential to expand the applicability rangeof quantum annealers significantly and may replace MCMCin diverse tasks of discrete sampling.Since the 16-bit system is used, the exact distribution de-pending on the temperature can be obtained by brute-force043050-5SHIBUKAWA, TAMURA, AND TSUDA PHYSICAL REVIEW RESEARCH 6, 043050 (2024)FIG. 4. Total variation distance between exact distributions and Stein-corrected distributions for energies of states (∑di=1 xi )2 and(∑di=1 xi )4, respectively. The results by QA samples are also shown.technique. Thus, we can compare the exact distributions andthe Stein-corrected distribution directly, and the total varia-tion distance depending on the temperature is summarized inFig. 4. The energies of states (∑di=1 xi )2 and (∑di=1 xi )4 arethe target observables. By comparing with the results of QAs,we confirm that the Stein correction can fit the distributions tomultiple observables simultaneously.IV. CONCLUSIONWe demonstrated that fast Stein correction is a helpfulcompanion of quantum annealers and fundamentally enhancestheir usability. The proposed method will play a role asan efficient postprocessing for quantum annealing [11] toachieve temperature-dependent Boltzmann sampling. Sinceour approach is a general method that can obtain the desireddistribution given arbitrary samples, it can be applied to othersoftware- and hardware-based sampling methods. The advan-tage of QA samples is that they are not locally concentrated,whereas MCMC samples have difficulty covering the wholespace. We believe that our proposed method combining QAsand Stein corrections is useful in many kinds of applicationswhere MCMC does not work efficiently and high-quality sam-ples cannot be obtained by random sampling. For example,our method is effective for constrained optimization prob-lems. When dealing with problems with equality constraints,random walk does not work efficiently because some statessatisfying constraints are not reachable from each other by asingle spin flip. Generating random samples by MCMC is thusdifficult in constrained problems. On the other hand, QAs canconsider equality constraints by using appropriate procedureswhen designing the Ising model [27]. Therefore, high-qualitysamples can be obtained even when such constraints areincluded.Although Stein correction cannot bring the distributionalerror to zero, it would be particularly useful to samplefrom highly constrained spaces [28], where global mixingby MCMC is extremely hard. Our future work involves theapplication of our method to machine learning and statisticalphysics and other highly scalable Ising machines such ascoherent Ising machines [29] and GPU-based algorithms [30].It can also be applied to general quantum computers includingNISQ, where distributional error is unavoidable due to envi-ronmental noise [31].ACKNOWLEDGMENTSR.S. thanks participants at AQC2023 in Albuquerque, NewMexico, for fruitful discussions. This work is supported byAIP Kasoku under Grant No. JPMJCR21U2, JST CRESTunder Grant No. JPMJCR21O2, JST ERATO under Grant No.JPMJER1903, KAKENHI under Grant No. 19H05819, andMEXT under Grant No. JPMXP1122712807.[1] K. Binder and A. P. Young, Spin glasses: Experimental facts,theoretical concepts, and open questions, Rev. Mod. Phys. 58,801 (1986).[2] K. Binder and E. Luijten, Monte Carlo tests of renormalization-group predictions for critical phenomena in Ising models, Phys.Rep. 344, 179 (2001).[3] N. Zhang, S. Ding, J. Zhang, and Y. Xue, An overview onrestricted Boltzmann machines, Neurocomputing 275, 1186(2018).[4] R. G. Melko, G. Carleo, J. Carrasquilla, and J. I. Cirac, Re-stricted Boltzmann machines in quantum physics, Nat. Phys.15, 887 (2019).043050-6https://doi.org/10.1103/RevModPhys.58.801https://doi.org/10.1016/S0370-1573(00)00127-7https://doi.org/10.1016/j.neucom.2017.09.065https://doi.org/10.1038/s41567-019-0545-1BOLTZMANN SAMPLING WITH QUANTUM ANNEALERS VIA … PHYSICAL REVIEW RESEARCH 6, 043050 (2024)[5] A. Doucet, N. De Freitas, N. J. Gordon et al., Sequential MonteCarlo methods in practice (Springer, New York, 2001).[6] K. Hukushima and K. Nemoto, Exchange Monte Carlo methodand application to spin glass simulations, J. Phys. Soc. Jpn. 65,1604 (1996).[7] Y. Iba, Extended ensemble Monte Carlo, Int. J. Mod. Phys. C12, 623 (2001).[8] K. Hukushima and Y. Iba, Population annealing and its applica-tion to a spin glass, AIP Conf. Proc. 690, 200 (2003).[9] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F.Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P.Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi,E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich,M. C. Thom et al., Quantum annealing with manufacturedspins, Nature (London) 473, 194 (2011).[10] J. Nelson, M. Vuffray, A. Y. Lokhov, T. Albash, and C. Coffrin,High-Quality thermal Gibbs sampling with quantum annealinghardware, Phys. Rev. Appl. 17, 044046 (2022).[11] J. Raymond, S. Yarkoni, and E. Andriyash, Global warm-ing: Temperature estimation in annealers, Front. ICT 3, 23(2016).[12] J. Marshall, D. Venturelli, I. Hen, and E. G. Rieffel, Power ofpausing: Advancing understanding of thermalization in experi-mental quantum annealers, Phys. Rev. Appl. 11, 044083 (2019).[13] J. Marshall, E. G. Rieffel, and I. Hen, Thermalization, freeze-out, and noise: Deciphering experimental quantum annealers,Phys. Rev. Appl. 8, 064025 (2017).[14] R. Sandt and R. Spatschek, Efficient low temperature MonteCarlo sampling using quantum annealing, Sci. Rep. 13, 6754(2023).[15] M. Vuffray, C. Coffrin, Y. A. Kharkov, and A. Y. Lokhov, Pro-grammable quantum annealers as noisy Gibbs samplers, PRXQuantum 3, 020317 (2022).[16] Y. Matsuda, H. Nishimori, and H. G. Katzgraber, Ground-statestatistics from annealing algorithms: quantum versus classicalapproaches, New J. Phys. 11, 073021 (2009).[17] S. T. Tokdar and R. E. Kass, Importance sampling: A review,Wiley Interdiscip. Rev. Comput. Stat. 2, 54 (2010).[18] Q. Liu and J. Lee, Black-box Importance Sampling, in Pro-ceedings of the 20th International Conference on ArtificialIntelligence and Statistics, Proceedings of Machine LearningResearch (PMLR, 2017), Vol. 54, pp. 952–961.[19] L. Hodgkinson, R. Salomone, and F. Roosta, The reproduc-ing Stein kernel approach for post-hoc corrected sampling,arXiv:2001.09266.[20] J. Yang, Q. Liu, V. Rao, and J. Neville, Goodness-of-fit testingfor discrete distributions via Stein discrepancy, in Proceedingsof the 35th International Conference on Machine Learning,Proceedings of Machine Learning Research (PMLR, 2018),Vol. 80, pp. 5561–5570.[21] A. Rahimi and B. Recht, Random features for large-scale kernelmachines, Proceedings of the 20th International Conferenceon Neural Information Processing Systems (Curran Associates,Inc., 2007), p. 1177.[22] J. Kivinen and M. K. Warmuth, Exponentiated gradient ver-sus gradient descent for linear predictors, Inf. Comput. 132, 1(1997).[23] GitHub, https://github.com/tsudalab/fast-stein-correction.[24] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.6.043050 for Fig. S1: Residual er-rors by exact and fast calculations of Stein correction; Fig. S2:Computational time depending on the problem size; Fig. S3:Residual error depending on the number of random features;Fig. S4: Number of unique samples depending on the targettemperature; Fig. S5: Elapsed time between Stein-correctionand Metropolis algorithm; Fig. S6: Residual errors of physicalquantities depending on the number of samples.[25] Dwavesystems/minorminer, https://github.com/dwavesystems/minorminer, accessed: 2023-08-03.[26] M. S. Andersen, J. Dahl, L. Vandenberghe et al., CVXOPT: APython package for convex optimization, https://cvxopt.org/ 54(2013).[27] I. Hen, Equation planting: A tool for benchmarking Ising ma-chines, Phys. Rev. Appl. 12, 011003(R) (2019).[28] S. Aoki, H. Hara, and A. Takemura, Markov Bases in AlgebraicStatistics (Springer, Berlin, 2012).[29] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T.Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu,O. Tadanaga, H. Takenouchi, K. Aihara, K.-I. Kawarabayashi,K. Inoue, S. Utsunomiya, and H. Takesue, A coherent Isingmachine for 2000-node optimization problems, Science 354,603 (2016).[30] Z. Mao, Y. Matsuda, R. Tamura, and K. Tsuda, Chemical designwith GPU-based Ising machines, Digital Discovery 2, 1098(2023).[31] E. Pelofske, J. Golden, A. Bärtschi, D. O’Malley, and S.Eidenbenz, Sampling on NISQ devices: “Who’s the fairest oneof all?”, in 2021 IEEE International Conference on QuantumComputing and Engineering (QCE) (IEEE, New York, 2021),pp. 207–217.043050-7https://doi.org/10.1143/JPSJ.65.1604https://doi.org/10.1142/S0129183101001912https://doi.org/10.1063/1.1632130https://doi.org/10.1038/nature10012https://doi.org/10.1103/PhysRevApplied.17.044046https://doi.org/10.3389/fict.2016.00023https://doi.org/10.1103/PhysRevApplied.11.044083https://doi.org/10.1103/PhysRevApplied.8.064025https://doi.org/10.1038/s41598-023-33828-2https://doi.org/10.1103/PRXQuantum.3.020317https://doi.org/10.1088/1367-2630/11/7/073021https://doi.org/10.1002/wics.56https://arxiv.org/abs/2001.09266https://doi.org/10.1006/inco.1996.2612https://github.com/tsudalab/fast-stein-correctionhttp://link.aps.org/supplemental/10.1103/PhysRevResearch.6.043050https://github.com/dwavesystems/minorminerhttps://cvxopt.org/https://doi.org/10.1103/PhysRevApplied.12.011003https://doi.org/10.1126/science.aah4243https://doi.org/10.1039/D3DD00047H