UIDP/00618/2020, PTDC/FIS-PAR/31000/2017, CERN/FISPAR /0002/2017, CERN/FIS-
PAR/0014/2019. The work of FE and MM is supported by the BMBF-Project
05H21VKCCA.

## References

  * [1] ATLAS collaboration, G. Aad et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1–29, [1207.7214].
  * [2] CMS collaboration, S. Chatrchyan et al., _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30–61, [1207.7235].
  * [3] F. Egle, M. Mühlleitner, R. Santos and J. a. Viana, _One-loop corrections to the Higgs boson invisible decay in a complex singlet extension of the SM_ , _Phys. Rev. D_ 106 (2022) 095030, [2202.04035].
  * [4] ATLAS collaboration, M. Aaboud et al., _Combination of searches for invisible Higgs boson decays with the ATLAS experiment_ , _Phys. Rev. Lett._ 122 (2019) 231801, [1904.05105].
  * [5] R. Costa, M. Mühlleitner, M. O. P. Sampaio and R. Santos, _Singlet Extensions of the Standard Model at LHC Run 2: Benchmarks and Comparison with the NMSSM_ , _JHEP_ 06 (2016) 034, [1512.05355].
  * [6] A. Djouadi, J. Kalinowski and M. Spira, _Hdecay: a program for higgs boson decays in the standard model and its supersymmetric extension_ , _Computer Physics Communications_ 108 (Jan, 1998) 56–74.
  * [7] A. Djouadi, J. Kalinowski, M. Mühlleitner and M. Spira, _Hdecay: Twenty++ years after_ , _Computer Physics Communications_ 238 (May, 2019) 214–231.
  * [8] sHDECAY collaboration, R. Costa, M. Mühlleitner, M. Sampaio and R. Santos, _Webpage for download of sHDECAY_ , _https://www.itp.kit.edu/ maggie/sHDECAY/_ .
  * [9] M. Krause, R. Lorenz, M. Muhlleitner, R. Santos and H. Ziesche, _Gauge-independent Renormalization of the 2-Higgs-Doublet Model_ , _JHEP_ 09 (2016) 143, [1605.04853].
  * [10] M. Krause, M. Muhlleitner, R. Santos and H. Ziesche, _Higgs-to-Higgs boson decays in a 2HDM at next-to-leading order_ , _Phys. Rev. D_ 95 (2017) 075019, [1609.04185].
  * [11] M. Krause, M. Mühlleitner and M. Spira, _2HDECAY —A program for the calculation of electroweak one-loop corrections to Higgs decays in the Two-Higgs-Doublet Model including state-of-the-art QCD corrections_ , _Comput. Phys. Commun._ 246 (2020) 106852, [1810.00768].
  * [12] M. Krause and M. Mühlleitner, _Impact of Electroweak Corrections on Neutral Higgs Boson Decays in Extended Higgs Sectors_ , _JHEP_ 04 (2020) 083, [1912.03948].
  * [13] M. Krause, D. Lopez-Val, M. Muhlleitner and R. Santos, _Gauge-independent Renormalization of the N2HDM_ , _JHEP_ 12 (2017) 077, [1708.01578].
  * [14] M. Krause and M. Mühlleitner, _ewN2HDECAY - A program for the Calculation of Electroweak One-Loop Corrections to Higgs Decays in the Next-to-Minimal Two-Higgs-Doublet Model Including State-of-the-Art QCD Corrections_ , 1904.02103.
  * [15] A. Denner, L. Jenniches, J.-N. Lang and C. Sturm, _Gauge-independent $\overline{MS}$ renormalization in the 2HDM_, _JHEP_ 09 (2016) 115, [1607.07352].
  * [16] L. Altenkamp, S. Dittmaier and H. Rzehak, _Renormalization schemes for the Two-Higgs-Doublet Model and applications to h $\rightarrow$ WW/ZZ $\rightarrow$ 4 fermions_, _JHEP_ 09 (2017) 134, [1704.02645].
  * [17] L. Altenkamp, S. Dittmaier and H. Rzehak, _Precision calculations for $h\to WW/ZZ\to 4$ fermions in the Two-Higgs-Doublet Model with Prophecy4f_, _JHEP_ 03 (2018) 110, [1710.07598].
  * [18] A. Denner, S. Dittmaier and J.-N. Lang, _Renormalization of mixing angles_ , _Journal of High Energy Physics_ 2018 (Nov, 2018) 104.
  * [19] M. Fox, W. Grimus and M. Löschner, _Renormalization and radiative corrections to masses in a general Yukawa model_ , _Int. J. Mod. Phys. A_ 33 (2018) 1850019, [1705.09589].
  * [20] W. Grimus and M. Löschner, _Renormalization of the multi-Higgs-doublet Standard Model and one-loop lepton mass corrections_ , _JHEP_ 11 (2018) 087, [1807.00725].
  * [21] J. Fleischer and F. Jegerlehner, _Radiative Corrections to Higgs Decays in the Extended Weinberg-Salam Model_ , _Phys. Rev. D_ 23 (1981) 2001–2026.
  * [22] A. Pilaftsis, _Resonant CP violation induced by particle mixing in transition amplitudes_ , _Nucl. Phys. B_ 504 (1997) 61–107, [hep-ph/9702393].
  * [23] S. Kanemura, Y. Okada, E. Senaha and C. P. Yuan, _Higgs coupling constants as a probe of new physics_ , _Phys. Rev. D_ 70 (2004) 115002, [hep-ph/0408364].
  * [24] J. M. Cornwall and J. Papavassiliou, _Gauge invariant three gluon vertex in qcd_ , _Phys. Rev._ D40 (1989) 3474.
  * [25] J. Papavassiliou, _Gauge independent transverse and longitudinal self-energies and vertices via the pinch technique_ , _Phys. Rev._ D50 (1994) 5958–5970, [hep-ph/9406258].
  * [26] D. Azevedo, P. Gabriel, M. Muhlleitner, K. Sakurai and R. Santos, _One-loop corrections to the Higgs boson invisible decay in the dark doublet phase of the N2HDM_ , _JHEP_ 10 (2021) 044, [2104.03184].
  * [27] A. Bredenstein, A. Denner, S. Dittmaier and M. M. Weber, _Precise predictions for the Higgs-boson decay H — $>$ WW/ZZ —$>$ 4 leptons_, _Phys. Rev. D_ 74 (2006) 013004, [hep-ph/0604011].
  * [28] J. Kublbeck, M. Bohm and A. Denner, _Feyn Arts: Computer Algebraic Generation of Feynman Graphs and Amplitudes_ , _Comput. Phys. Commun._ 60 (1990) 165–180.
  * [29] T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput. Phys. Commun._ 140 (2001) 418–431, [hep-ph/0012260].
  * [30] R. Mertig, M. Bohm and A. Denner, _FEYN CALC: Computer algebraic calculation of Feynman amplitudes_ , _Comput. Phys. Commun._ 64 (1991) 345–359.
  * [31] V. Shtabovenko, R. Mertig and F. Orellana, _New Developments in FeynCalc 9.0_ , _Comput. Phys. Commun._ 207 (2016) 432–444, [1601.01167].
  * [32] T. Hahn and M. Perez-Victoria, _Automatized one loop calculations in four-dimensions and D-dimensions_ , _Comput. Phys. Commun._ 118 (1999) 153–165, [hep-ph/9807565].
  * [33] G. J. van Oldenborgh and J. A. M. Vermaseren, _New Algorithms for One Loop Integrals_ , _Z. Phys. C_ 46 (1990) 425–438.
  * [34] F. Staub, _From Superpotential to Model Files for FeynArts and CalcHep/CompHep_ , _Comput. Phys. Commun._ 181 (2010) 1077–1086, [0909.2863].
  * [35] F. Staub, _Automatic Calculation of supersymmetric Renormalization Group Equations and Self Energies_ , _Comput. Phys. Commun._ 182 (2011) 808–833, [1002.0840].
  * [36] F. Staub, _SARAH 3.2: Dirac Gauginos, UFO output, and more_ , _Comput. Phys. Commun._ 184 (2013) 1792–1809, [1207.0906].
  * [37] F. Staub, _SARAH 4 : A tool for (not only SUSY) model builders_ , _Comput. Phys. Commun._ 185 (2014) 1773–1790, [1309.7223].
  * [38] F. Staub, _Exploring new models in all detail with SARAH_ , _Adv. High Energy Phys._ 2015 (2015) 840780, [1503.04200].
  * [39] N. D. Christensen and C. Duhr, _FeynRules - Feynman rules made easy_ , _Comput. Phys. Commun._ 180 (2009) 1614–1641, [0806.4194].
  * [40] C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201–1214, [1108.2040].
  * [41] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250–2300, [1310.1921].
  * [42] R. Coimbra, M. O. P. Sampaio and R. Santos, _ScannerS: Constraining the phase diagram of a complex scalar singlet at the LHC_ , _Eur. Phys. J. C_ 73 (2013) 2428, [1301.2599].
  * [43] M. Mühlleitner, M. O. P. Sampaio, R. Santos and J. Wittbrodt, _ScannerS: Parameter Scans in Extended Scalar Sectors_ , 2007.02985.
  * [44] LHC Higgs Cross Section Working Group, LHC Twiki page.
  * [45] Particle Data Group collaboration, P. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [46] B. W. Lee, C. Quigg and H. B. Thacker, _Weak Interactions at Very High-Energies: The Role of the Higgs Boson Mass_ , _Phys. Rev. D_ 16 (1977) 1519.
  * [47] M. E. Peskin and T. Takeuchi, _Estimation of oblique electroweak corrections_ , _Phys. Rev. D_ 46 (07, 1992) 381–409.
  * [48] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, _$HiggsSignals$ : Confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC_, _Eur. Phys. J. C_ 74 (2014) 2711, [1305.1933].
  * [49] P. Bechtle, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein and J. Wittbrodt, _HiggsSignals-2: Probing new physics with precision Higgs measurements in the LHC 13 TeV era_ , _Eur. Phys. J. C_ 81 (2021) 145, [2012.09197].
  * [50] P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein and K. E. Williams, _HiggsBounds: Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron_ , _Comput. Phys. Commun._ 181 (2010) 138–167, [0811.4169].
  * [51] P. Bechtle, D. Dercks, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein et al., _HiggsBounds-5: Testing Higgs Sectors in the LHC 13 TeV Era_ , 2006.06007.
  * [52] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model_ , _Comput. Phys. Commun._ 176 (2007) 367–382, [hep-ph/0607059].
  * [53] C. Gross, O. Lebedev and T. Toma, _Cancellation Mechanism for Dark-Matter–Nucleon Interaction_ , _Phys. Rev. Lett._ 119 (2017) 191801, [1708.02253].
  * [54] D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki and R. Santos, _Testing scalar versus vector dark matter_ , _Phys. Rev. D_ 99 (2019) 015017, [1808.01598].
  * [55] D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki and R. Santos, _One-loop contribution to dark-matter-nucleon scattering in the pseudo-scalar dark matter model_ , _JHEP_ 01 (2019) 138, [1810.06105].
  * [56] S. Glaus, M. Mühlleitner, J. Müller, S. Patel, T. Römer and R. Santos, _Electroweak Corrections in a Pseudo-Nambu Goldstone Dark Matter Model Revisited_ , _JHEP_ 12 (2020) 034, [2008.12985].
  * [57] XENON collaboration, E. Aprile et al., _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302, [1805.12562].
  * [58] PandaX-4T collaboration, Y. Meng et al., _Dark Matter Search Results from the PandaX-4T Commissioning Run_ , _Phys. Rev. Lett._ 127 (2021) 261802, [2107.13438].
  * [59] J. Aalbers et al., _First dark matter search results from the lux-zeplin (lz) experiment_ , 2022.  10.48550/ARXIV.2207.03764.
  * [60] P. Basler and M. Mühlleitner, _BSMPT (Beyond the Standard Model Phase Transitions): A tool for the electroweak phase transition in extended Higgs sectors_ , _Comput. Phys. Commun._ 237 (2019) 62–85, [1803.02846].
  * [61] P. Basler, M. Mühlleitner and J. Müller, _BSMPT v2 a tool for the electroweak phase transition and the baryon asymmetry of the universe in extended Higgs Sectors_ , _Comput. Phys. Commun._ 269 (2021) 108124, [2007.01725].
  * [62] A. D. Sakharov, _Violation of CP Invariance, c Asymmetry, and Baryon Asymmetry of the Universe_ , _Pisma Zh. Eksp. Teor. Fiz._ 5 (1967) 32–35.
  * [63] H. Abouabid, A. Arhrib, D. Azevedo, J. E. Falaki, P. M. Ferreira, M. Mühlleitner et al., _Benchmarking di-Higgs production in various extended Higgs sector models_ , _JHEP_ 09 (2022) 011, [2112.12515].
  * [64] R. V. Harlander, S. Liebler and H. Mantler, _SusHi: A program for the calculation of Higgs production in gluon fusion and bottom-quark annihilation in the Standard Model and the MSSM_ , _Comput. Phys. Commun._ 184 (2013) 1605–1617, [1212.3249].
  * [65] S. Liebler, _Neutral Higgs production at proton colliders in the CP-conserving NMSSM_ , _Eur. Phys. J. C_ 75 (2015) 210, [1502.07972].
  * [66] R. V. Harlander, S. Liebler and H. Mantler, _SusHi Bento: Beyond NNLO and the heavy-top limit_ , _Comput. Phys. Commun._ 212 (2017) 239–257, [1605.03190].
  * [67] ATLAS collaboration, M. Aaboud et al., _Search for pair production of Higgs bosons in the $b\bar{b}b\bar{b}$ final state using proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _JHEP_ 01 (2019) 030, [1804.06174].
  * [68] ATLAS Collaboration collaboration, _Search for resonant pair production of Higgs bosons in the $b\bar{b}b\bar{b}$ final state using $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector_, , CERN, Geneva, Aug, 2021.
  * [69] CMS collaboration, _Search for resonant Higgs boson pair production in four b quark final state using large-area jets in proton-proton collisions at $\sqrt{s}=13\leavevmode\nobreak\ \mathrm{TeV}$_, .
  * [70] ATLAS collaboration, M. Aaboud et al., _Search for resonant and non-resonant Higgs boson pair production in the ${b\bar{b}\tau^{+}\tau^{-}}$ decay channel in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. Lett._ 121 (2018) 191801, [1808.00336].
  * [71] ATLAS collaboration, G. Aad et al., _Reconstruction and identification of boosted di- $\tau$ systems in a search for Higgs boson pairs using 13 TeV proton-proton collision data in ATLAS_, _JHEP_ 11 (2020) 163, [2007.14811].
  * [72] ATLAS Collaboration collaboration, _Search for resonant and non-resonant Higgs boson pair production in the $b\bar{b}\tau^{+}\tau^{-}$ decay channel using 13 TeV $pp$ collision data from the ATLAS detector_, , CERN, Geneva, Jul, 2021.
  * [73] ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $\gamma\gamma b\bar{b}$ final state with 13 TeV $pp$ collision data collected by the ATLAS experiment_, _JHEP_ 11 (2018) 040, [1807.04873].
  * [74] ATLAS Collaboration collaboration, _Search for Higgs boson pair production in the two bottom quarks plus two photons final state in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, , CERN, Geneva, Mar, 2021.
  * [75] ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $b\bar{b}WW^{*}$ decay mode at $\sqrt{s}=13$ TeV with the ATLAS detector_, _JHEP_ 04 (2019) 092, [1811.04671].
  * [76] CMS collaboration, A. M. Sirunyan et al., _Search for resonant pair production of Higgs bosons in the $bbZZ$ channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Rev. D_ 102 (2020) 032003, [2006.06391].
  * [77] ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $\gamma\gamma WW^{*}$ channel using $pp$ collision data recorded at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Eur. Phys. J. C_ 78 (2018) 1007, [1807.08567].
  * [78] ATLAS collaboration, M. Aaboud et al., _Search for Higgs boson pair production in the $WW^{(*)}WW^{(*)}$ decay channel using ATLAS data recorded at $\sqrt{s}=13$ TeV_, _JHEP_ 05 (2019) 124, [1811.11028].
  * [79] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * [80] G. Passarino and M. J. G. Veltman, _One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model_ , _Nucl. Phys. B_ 160 (1979) 151–207.
  * [81] J. Alison et al., _Higgs boson potential at colliders: Status and perspectives_ , _Rev. Phys._ 5 (2020) 100045, [1910.00012].
  * [82] A. Djouadi, W. Kilian, M. Muhlleitner and P. M. Zerwas, _Testing Higgs selfcouplings at e+ e- linear colliders_ , _Eur. Phys. J. C_ 10 (1999) 27–43, [hep-ph/9903229].
  * [83] A. Djouadi, W. Kilian, M. Mühlleitner and P. M. Zerwas, _Production of neutral Higgs boson pairs at LHC_ , _Eur. Phys. J._ C10 (1999) 45–49, [hep-ph/9904287].
  * [84] J. Baglio, A. Djouadi, R. Gröber, M. M. Mühlleitner, J. Quevillon and M. Spira, _The measurement of the Higgs self-coupling at the LHC: theoretical status_ , _JHEP_ 04 (2013) 151, [1212.5581].
  * [85] LHC Higgs Cross Section Working Group collaboration, D. de Florian et al., _Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector_ , 1610.07922.
  * [86] E. W. N. Glover and J. J. van der Bij, _HIGGS BOSON PAIR PRODUCTION VIA GLUON FUSION_ , _Nucl. Phys._ B309 (1988) 282–294.
  * [87] D. A. Dicus, C. Kao and S. S. D. Willenbrock, _Higgs Boson Pair Production From Gluon Fusion_ , _Phys. Lett. B_ 203 (1988) 457–461.
  * [88] T. Plehn, M. Spira and P. M. Zerwas, _Pair production of neutral Higgs particles in gluon-gluon collisions_ , _Nucl. Phys._ B479 (1996) 46–64, [hep-ph/9603205].
  * [89] ATLAS collaboration, _Constraining the Higgs boson self-coupling from single- and double-Higgs production with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13$ TeV_, 2211.01216.
  * [90] CMS collaboration, A. Tumasyan et al., _A portrait of the Higgs boson by the CMS experiment ten years after the discovery_ , _Nature_ 607 (2022) 60–68, [2207.00043].


# Universal Statistical Simulator

Mark Carney  Quantum Threat Group, Cybersecurity Research, Santander Global,
London  Ben Varcoe  Department of Physics and Astronomy, University of Leeds,
Leeds, United Kingdom

###### Abstract

The Quantum Fourier Transform is a famous example in quantum computing for
being the first demonstration of a useful algorithm in which a quantum
computer is exponentially faster than a classical computer. However when
giving an explanation of the speed up, understanding computational complexity
of a classical calculation has to be taken on faith. Moreover, the explanation
also comes with the caveat that the current classical calculations might be
improved. In this paper we present a quantum computer code for a Galton Board
Simulator that is exponentially faster than a classical calculation using an
example that can be intuitively understood without requiring an understanding
of computational complexity. We demonstrate a straight forward implementation
on a quantum computer, using only three types of quantum gate, which
calculates $2^{n}$ trajectories using $\mathcal{O}(n^{2})$ resources. The
circuit presented here also benefits from having a lower depth than previous
Quantum Galton Boards, and in addition, we show that it can be extended to a
universal statistical simulator which is achieved by removing pegs and
altering the left-right ratio for each peg.

##  1 Introduction

The quantum mechanical Galton Board has now been a core concept in quantum
computing for nearly a decade [1] with its roots extending back to the
development of linear optics quantum computing. Apart from the interest in
demonstrating quantum enhanced performance, a Galton Board (GB) is also
interesting for its computational properties. Removing the pegs and altering
the balance of each peg allows for statistical calculations to be performed.
In this way, the GB is a statistical simulator rather than a universal
computer, nevertheless it allows for a number of applications[15] from complex
system simulation[7, 6, 20, 4], modelling random walks across graphs[21, 24,
17, 27, 5], machine learning[18, 9, 8], modelling stock price fluctuations[13,
12, 30, 22], one way functions for cryptography[2, 11, 14, 10], and sampling
and search problems[1, 3, 29, 26, 25].

In this paper we present an intuitive approach to generating a Quantum Galton
Board (QGB) that mimics the physical action of the classical GB in a quantum
circuit. The aim of this paper is not necessarily to present a new quantum
algorithm, but rather to provide a clear and transparent demonstration of the
exponential speed up that can be provided with a quantum processor. In
comparison with previous research which replicated the output of the classical
GB, we create a quantum circuit that creates a superposition of all possible
trajectories. In a ”sampling” mode this possesses the statistics of a GB, but
it also allows for individual trajectories to be tailored to meet specific
requirements. Moreover, in contrast with previous work the entire circuit is
performed using combinations of three Clifford operations and ancilla qubits.
To perform an n-level QGB we require $n(2n-1)$ gates and $2n$ qubits ($2n-1$
working qubits, and a recycled control qubit).

The requirement of a large number of ancilla qubits poses a large overhead.
however, as the aim is to presenting an approach that generates the actual
trajectories in a classical GB this is a relatively small price to pay.

##  2 Preliminaries

###  2.1 The classical Galton Board

Figure 1: A standard pyramidal Galton board, from Wolfram MathWorld [28].

The standard construction of a Galton board is presented in figure 1. In the
regular formulation, each column from $1\ldots n+1$ of the output of an $n$
level Galton board has, for $q=1-p=0.5$ the value

$\displaystyle\binom{n}{k}p^{k}q^{n-k}=\frac{1}{2^{n}}\binom{n}{k}$ (1)

Although it is omitted here, the proof that a Galton board produces the normal
distribution can be obtained via the De Moivre-Laplace theorem.

##  3 The Quantum Galton Board

The basis for our approach to the construction of the QGB is to model the
action of an individual peg on a physical Galton board. In the classical GB an
ideal ball hitting the ideal peg produces a 50% probability of going left or
right, then hitting another peg on the board.

To mimic this, we take a set of $\ket{0}$ initialized qubits, and invert the
middlemost qubit using an $X$ gate. We then construct a series of superposed
$SWAP$s that this ‘ball’ can effectively ‘fall through’ mimicking the action
of a real ball through an array of pegs on a Galton board.

We can use this model to get a normal distribution by relying on equation 1 in
the same way that you do on a physical, macro-scale Galton board.

###  3.1 A Quantum Peg

$B$$P$Go rightGo left Figure 2: Basic model of a ball $B$ interacting with a
physical peg $P$ that forms the basis of our ‘quantum peg’ circuit module in
fig 3. $\Qcircuit{&q_{0}\lstick{\ket{0}}\gate{H}\ctrl{1}\targ\ctrl{2}\qw\qw\\\
q_{1}\lstick{\ket{0}}\qw\qswap\qw\qw\meter\cw\\\
q_{2}\lstick{\ket{0}}\gate{X}\qswap\qwx\ctrl{-2}\qswap\qw\qw\\\
q_{3}\lstick{\ket{0}}\qw\qw\qw\qswap\qwx\meter\cw}$ Figure 3: A Quantum
circuit analogue of the physical peg in figure 2. Qubit $q_{2}$ is the input
channel or ‘ball’ and meters (measurements) indicate the output channels of
the ‘quantum peg’.

We first construct a ‘quantum peg’ module (figure 3) that we can then
replicate to build a complete Galton board. The basic peg uses three working
qubits ($q_{1}$, $q_{2}$, and $q_{3}$) and one control qubit ($q_{0}$). We
initialize all qubits to zero, then place the control qubit $q_{0}$ into
superposition with the Hadamard gate and place a ‘ball’ on the middle working
qubit, $q_{2}$ using an $X$-gate.

We then perform a controlled-$SWAP$ operation on $q_{1}$ and $q_{2}$ giving us
an eventual state for $\ket{q_{2}q_{1}q_{0}}$ as

$\frac{1}{\sqrt{2}}(\ket{011}+\ket{100})$

which by an application of an inverted $CNOT$ on $q_{2}$ to $q_{0}$ become

$\frac{1}{\sqrt{2}}(\ket{011}+\ket{101})$

Given our control qubit is now a stable $\ket{1}$ we can apply a swap between
$q_{2}$ and $q_{3}$ to give the desired final state of

$\ket{q_{3}q_{2}q_{1}q_{0}}=\frac{1}{\sqrt{2}}(\ket{0011}+\ket{1001})$

This module can, with some care around the control qubit, be replicated to
reproduce each level of the Galton board successively.

###  3.2 A 3-peg QGB

$\Qcircuit<EMAIL_ADDRESS>q_{1}\lstick{\ket{0}}\qw\qw\qw\qw\qw\qw\qswap\qwx\qw\qw\qw\qw\qw\qw\meter\cw\\\
q_{2}\lstick{\ket{0}}\qw\qswap\qw\qw\qw\qw\qswap\qwx\ctrl{-2}\qswap\qwx\qw\qw\qw\qw\qw\qw\\\
q_{3}\lstick{\ket{0}}\gate{X}\qswap\qwx\ctrl{-3}\qswap\qw\qw\qw\qw\qswap\qwx\ctrl{-3}\qswap\qw\qw\meter\cw\\\
q_{4}\lstick{\ket{0}}\qw\qw\qw\qswap\qwx\qw\qw\qw\qw\qw\qw\qswap\qwx\ctrl{-4}\qswap\qw\qw\\\
q_{5}\lstick{\ket{0}}\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qswap\qwx\meter\cw}$
Figure 4: A Quantum circuit analogue of the physical peg in figure 2.

We detail in figure 4 a 3-peg (2 level) version of a QGB implemented using our
method. Note, for 6 qubits, we get 3 classical bits out, which is half what is
achieved by other methods (see [23]). However, we have what we consider to be
close to a minimal circuit depth for our algorithm.

The circuit in fig 4 can be thought of in three sections. The first is
identical to the individual peg in fig 3. We then reset the control qubit
$q_{0}$ and proceed to take the ‘output’ of the first peg on qubits $q_{2}$
and $q_{4}$, and apply the same quantum peg logic, adding an additional
inverted $CNOT$ (the 10th gate column in our diagram fig 4). This $CNOT$ now
returns the control qubit probability to 50%, which then allows us to proceed
onto implementing a quantum peg for the other ‘half’ of the circuit. Finally,
we have three measurements on $q_{1}$, $q_{3}$, and $q_{5}$, with state of
measurement given by

$\ket{q_{5}q_{3}q_{1}}=\frac{1}{\sqrt{4}}(\ket{001}+2\ket{010}+\ket{100})$

with the full state begin given by

$\ket{q_{5}q_{4}q_{3}q_{2}q_{1}q_{0}}=\frac{1}{2}(\ket{10000}+2\ket{00100}+\ket{00001})\sqrt{3}\ket{1}$

###  3.3 Scaling the QGB

The additional $CNOT$ after our peg module circuit we described previously is
needed to re-balance the control qubit for the successive gates. The
requirement for a mid-circuit reset on the control qubit is also required for
every successive level of the circuit. Thus, we can expect for every ‘peg’ we
want to model, we require at most 4 gates, with optional local optimizations
of only needing 3-gates.

###  3.4 QGB Features

We first note that our outputs will differ from other approaches, as each
output measurement will have precisely a single ‘1’ output. Thus, some post-
processing should be required to determine the final output.

We require $n$-many ancilla qubits for each $n$ bits of desired output, which
is a limitation given the current scale of NISQ devices that are currently
available.

In terms of gate count, the following give the upper bound for $n$ desired
output bits, subject of course to local optimisations:

  * •

We require $n$ Hadamard gates, and $n$ reset ($\ket{0}$) gates for a Galton
board of $n$-many levels, plus one $X$ gate for the ‘ball’.

  * •

We require 4 gates per ‘peg’, the total number of which is governed by
triangle numbers.

  * •

We require $n+1$ measurement gates.

This gives us an upper limit on our gates of

$n+n+4\Big{(}\frac{n^{2}+n}{2}\Big{)}+n+1+1=2n^{2}+5n+2$

We note that for $n+1$-many output bits, this equates to $n$ levels of a
Galton board.

Our method has circuit depth less than half that demonstrated in [23], for
example, who present 167 gates on 5-qubit output, compared to our upper limit
of 76. However, ours requires post-processing to convert it from an array of
$n$ bits with a single 1 to a full 5-bit normally distributed output. The
significant decrease on circuit depth, however, significantly reduces the
amount of error our circuits are prone to experiencing, meaning in part that
an increase in the number of runs does not necessitate an increase in the
level of error.

##  4 Experimental Results for QGB

We performed local and remote simulations using Qiskit, IBM’s python-based
quantum circuit library.

###  4.1 Remote Simulations

Figure 5: Output results from an IBM simulation of the 4-QGB circuit.

We produced a 4-level Quantum Galton Board (4-QGB) circuit and implemented
this on the IBM-QX quantum simulator, as access to quantum computers with more
than 5 real qubits is currently restricted. The IBM-QX simulator, however,
uses quantum random numbers for its simulations, so the fidelity to actual
experiment is likely as close as can be achieved.

We constructed the circuit in OpenQASM - listed in Appendix A. The IBM
transpiler made no alterations to our code, and simulated the circuit, see
figure 13 in Appendix D.

The results are listed in the histogram in figure 5. As you can see, the
expected normal distribution can be found, albeit on 5 output qubits.

###  4.2 Real Simulations of Quantum Peg

Figure 6: Results histogram from a quantum computer hardware run of the quantum peg circuit (fig 3) State | Num. | % (from 8192 shots)
---|---|---
$\ket{000}$ | 975 | $\approx 11.90$
$\ket{001}$ | 2362 | $\approx 28.83$
$\ket{010}$ | 837 | $\approx 10.22$
$\ket{011}$ | 497 | $\approx 6.07$
$\ket{100}$ | 2075 | $\approx 25.33$
$\ket{101}$ | 380 | $\approx 4.64$
$\ket{110}$ | 611 | $\approx 7.46$
$\ket{111}$ | 455 | $\approx 5.55$
Table 1: Results from running the QGB peg circuit on quantum computer
hardware. The desired outcomes are highlighted in grey.

Although simulating a full QGB is not immediately possible, it was possible to
simulate a quantum peg. The circuit was uploaded to IBM-QX and the transpiler
created the circuit depicted in figure 14.

This transpiled circuit, however, clearly has a lot of noise owing to the fact
that 5 active quantum gates have now become 64 gates, which necessarily leads
to a significant increase in error. The results are presented in figure 6 and
table 1.

As is apparent, the output is very noisy, as the desired results represent
only 54.16% of the total output, with the rest being noise. That said, the two
peaks are on the desired outputs we defined earlier, $\ket{100}$ and
$\ket{001}$.111 _NB_ \- the measurement operators were not applied to qubit
$q_{0}$ during the runs, and IBM-QX defaults to a ‘0’ readout in that case.

###  4.3 Local Simulation and Rescaling of QGB Output

Figure 7: The re-scaled output of the 4-QGB locally simulated circuit. The
plot here shows the output of our processing on the simulated QC output
against the normal distribution with parameters $var=5$ (for visibility) and
$\mu=16$.

To demonstrate the post-processing we mentioned in section 3.4, we loaded the
circuit onto a local simulator and generated 20,000 shots, storing each
result. For the outputs we reassigned them with an integer from 0 to 4
(‘00001’ became ‘0’, ‘00010’ became ‘1’, and so on). Under this substitution,
the mean average was 1.9977, the standard deviation was 1.001521, and the
variance was 1.002995.

With this reassignment done, we took blocks of 8 outputs and summed each
block. This scaled our output from 0 to 32, as the maximum sum is $8\times 4$.
We then plotted how many of each sum we had, and compared this to the normal
distribution. The result can be seen in figure 7.

As can be seen, we have successfully scaled our output of individual bits
under normal distribution to a range - this method can be repeated for any
multiple of the highest-value substitution.

##  5 A Biased QGB

Of course, regular normal distributions are only part of the story when it
comes to generating useful random numbers. We present a construction of a
‘biased QGB‘ (B-QGB) where we can execute control over the final distirbution
with significant accuracy.

###  5.1 The Biased Quantum Peg

$\Qcircuit{&q_{0}\text{
}\lstick{\ket{\psi}}\gate{\ket{0}}\gate{R_{x}(\theta)}\ctrl{1}\targ\ctrl{2}\qw\\\
q_{1}\lstick{\ket{0}}\qw\qw\qswap\qw\qw\qw\\\
q_{2}\lstick{\ket{0}}\qw\gate{X}\qswap\qwx\ctrl{-2}\qswap\qw\\\
q_{3}\lstick{\ket{0}}\qw\qw\qw\qw\qswap\qwx\qw}$ Figure 8: A biased quantum
peg, based on figure 2 with the $H$ gate replaced with an $x$-axis rotation
$R_{x}(\theta)$

As before, we follow a modular approach and construct a quantum peg circuit
that is biased, depicted in figure 8.

This circuit follows naturally from the observation that our circuit in figure
3 can have its $H$ gate replaced with a rotation about the $x$-axis, a gate
usually denoted $R_{x}(\theta)$ which is defined as

$R_{x}(\theta)=\begin{pmatrix}\cos{\frac{\theta}{2}}&-i\sin{\frac{\theta}{2}}\\\
-i\sin{\frac{\theta}{2}}&\cos{\frac{\theta}{2}}\end{pmatrix}$

As such, if we replace the $H$ gate with $R_{x}(\frac{\pi}{2})$, then the
resultant circuit operates precisely the same - the $\cos{\frac{\pi}{4}}$
terms replace the $(\frac{1}{\sqrt{2}}-\frac{1}{\sqrt{2}})$ null factors on
the diagonal. Likewise, the $-\frac{i}{\sqrt{2}}$ and $\frac{1}{\sqrt{2}}$
terms are equivalent under the regular norm used in measurement.

With this, we can then see that by introducing a free rotation about the
$x$-axis in our circuit we can gain more control over the nature of the
underlying probabilities affecting the ‘quantum peg’.

However, to manage possibly unique rotations per peg, we also require an
additional prefix of a $RESET$ gate, which is theoretically no problem but
might present some issues on certain quantum computers. However, the payoff of
distribution control, and the fact that our circuits are significantly shorter
still, means that this payoff is worthwhile.

The overall equation for this biased quantum peg is then given as follows.
First note that:

$R_{x}(\theta)\ket{0}=\begin{pmatrix}\cos{\frac{\theta}{2}}\\\
-i\sin{\frac{\theta}{2}}\end{pmatrix}$

and the complimentary rotation

$\overline{R_{x}}(\theta)\ket{0}=R_{x}(\pi-\theta)\ket{0}=\begin{pmatrix}\cos{\frac{\pi-\theta}{2}}\\\
-i\sin{\frac{\pi-\theta}{2}}\end{pmatrix}$

We shall denote these states as $\ket{R_{\theta}}$ and
$\ket{\overline{R_{\theta}}}$ respectively. Thus the state of the biased
quantum peg in figure 8 is given by

$\ket{q_{3}q_{2}q_{1}q_{0}}=(\ket{001R_{\theta}}+\ket{100\overline{R_{\theta}}})$
(2)

The way that these states feed forwards into successive rows on our QGB means
that we can have a given rotational bias on each row that will behave as we
expect under successive compositions.

###  5.2 Biased Quantum Peg Example

We provide the following worked example for a biased peg that has the
probability of measuring a 1 on the upper qubit of 75% and 25% on the lower
qubit - qubits $q_{1}$ and $q_{3}$ respectively in figure 8.

To do this we set our $\theta=\frac{2\pi}{3}$ and note that the final state as
given in equation 2 becomes

$\ket{q_{3}q_{2}q_{1}q_{0}}=\Big{(}\sqrt{\frac{3}{4}}\ket{0011}+\frac{1}{\sqrt{4}}\ket{1001}\Big{)}$

When we apply measurement, these terms will give us the desired output of 75%
chance of a $\ket{0011}$ state, and a 25% chance of a $\ket{1001}$ state.
Should we wish/require to, we can flip this order by letting
$\theta=\frac{\pi}{3}$, as we might expect.

###  5.3 Biased QGB Features

The overall gate count does not increase too significantly. If we suppose that
we want per peg row bias to be set, which is the maximal case for control of
the circuit in our current setup, then we have:

  * •

5 gates per peg -

    * –

one $RESET$,

    * –

one $R_{x}(\theta)$,

    * –

two controlled-$SWAP$s,

    * –

(up to) two $CNOT$.

  * •

$n+1$ measurement gates.

  * •

One $X$ gate for the ‘ball’.

Again, we use the triangle numbers to get a maximum gate count for an
$n$-level QGB given by:

$3(n^{2}+n)+n+2$ (3)

In section 11 we take this idea further and demonstrate fine-grained control
of each individual peg without adding extra ancilla qubits.

##  6 Experimental Results for Biased-QGB

###  6.1 Biased-QGB Peg on QC Hardware

State | Num. | % (from 8192 shots)
---|---|---
$\ket{000}$ | 1108 | $\approx 13.52$
$\ket{001}$ | 2991 | $\approx 36.51$
$\ket{010}$ | 851 | $\approx 10.39$
$\ket{011}$ | 490 | $\approx 5.98$
$\ket{100}$ | 1464 | $\approx 17.87$
$\ket{101}$ | 492 | $\approx 6.00$
$\ket{110}$ | 352 | $\approx 4.30$
$\ket{111}$ | 444 | $\approx 5.42$
Table 2: Results from running a biased QGB peg on quantum computer hardware.
The desired outputs are highlighted in grey. Figure 9: Results histogram from
a real QC run of the biased quantum peg circuit (fig 8)

We ran the biased quantum peg circuit from fig 8 on the `ibmq-manila` quantum
computer. The transpiled circuit was of similar depth to the one generated
previously, also shown in figure 14 \- as we can see, it is broadly the same.
The results are presented in the graph in figure 9.

The number of shots in this particular run was 8192, the maximum, and the
outcome results are listed in table 2. Clearly the two states we are most
interested in, $\ket{001}$ and $\ket{100}$, were the most common outcomes,
with 36.51% and 17.87% of the outcomes respectively. However neighbouring
states $\ket{000}$ and $\ket{010}$ were also very prominent with $>10\%$ of
the overall outcomes each.

Thus we can claim that our circuit has worked on real hardware, but that the
overall noise floor omnipresent in NISQ-era quantum computing poses a real
challenge for extended use of multiple peg QGB’s constructed in the manner we
present.

###  6.2 Locally Simulated Biased-QGB

Figure 10: The re-scaled output of the biased 4-QGB locally simulated circuit.

We extended our methods in section 4.3 to run a biased QGB by means of
swapping out the $H$ gates for $R_{x}(\theta)$ gates, and setting
$\theta=\frac{2\pi}{3}$ as described in section 5.2. As before, we present our
plotted results in figure 10.

Clearly our scaling produces a desired skewed distribution, which validates
our approach in simulations at least. In this run, the mean value produced was
2.66, with a standard deviation of 1.176 and a variance of 1.383.

As before, our OpenQASM is included, found in appendix B.

##  7 Fine-Grained Biased QGB Circuits

$\Qcircuit@C=1em@R=1em{&q_{0}\lstick{\ket{0}}\gate{R_{x}(\theta_{1})}\ctrl{2}\targ\ctrl{3}\gate{\ket{0}}\gate{R_{x}(\theta_{2})}\ctrl{1}\targ\ctrl{2}\gate{\ket{0}}\gate{R_{x}(\theta_{3})}\ctrl{3}\targ\ctrl{4}\qw\qw\qw\\\
q_{1}\lstick{\ket{0}}\qw\qw\qw\qw\qw\qw\qswap\qwx\qw\qw\qw\qw\qw\qw\qw\qw\meter\cw\\\
q_{2}\lstick{\ket{0}}\qw\qswap\qw\qw\qw\qw\qswap\qwx\ctrl{-2}\qswap\qwx\qw\qw\qw\qw\qw\qw\qw\qw\\\
q_{3}\lstick{\ket{0}}\gate{X}\qswap\qwx\ctrl{-3}\qswap\qw\qw\qw\qw\qswap\qwx\qw\qw\qswap\qw\qw\targ\meter\cw\\\
q_{4}\lstick{\ket{0}}\qw\qw\qw\qswap\qwx\qw\qw\qw\qw\qw\qw\qw\qswap\qwx\ctrl{-4}\qswap\ctrl{-1}\gate{\ket{0}}\qw\\\
q_{5}\lstick{\ket{0}}\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qswap\qwx\qw\meter\cw}$
Figure 11: A Quantum circuit for the fine-grained per-peg control of bias on a
QGB. This circuit has the additional $CNOT$ at the end to adjust the outputs.
The final reset is not used here but would be important for subsequent rows.

###  7.1 Iterating Pegs for Fine-Grained Control on a QGB

In figure 11 we show how we can have per-peg control of the bias in a QGB.
This circuit takes three parameters for the $x$-axis rotations $\theta_{1}$,
$\theta_{2}$, and $\theta_{3}$.

The construction is based on the circuit in figure 4 but with the 10th $CNOT$
gate replaced with a $RESET$ and $R_{x}$ gate preparing the final peg instead.
To adjust for the loss of this correcting $CNOT$ gate we have to use an
additional $CNOT$ gate at the end to add the ‘leftover’ rotation to the
correct output on $q_{3}$, followed by an additional $RESET$. Although we do
not need this here, it would be required if we were to add additional rows of
pegs.

###  7.2 Features and Simulation of a Fine-Grained Biased QGB

####  7.2.1 Features

It is important that the pegs proceed in linear order so that the corrective
$CNOT$ gates can function properly. The correction needs to take place for
each peg except the first. As such the gate count for an $n$-row fine-grained
QGB is roughly:

  * •

One $X$ gate for the circuit.

  * •

Per peg:

    * –

One $RESET$ gate.

    * –

One $R_{x}(\theta)$ gate.

    * –

Two controlled-$SWAP$ gates.

    * –

One $CNOT$ gate

  * •

For row $i$, $(i-1)$-many $CNOT$ and $RESET$ gates for the end of that row’s
pegs, for each row.

  * •

$n$ measurement gates.

Thus overall we have

$1+5\Big{(}\frac{n^{2}+n}{2}\Big{)}+\Big{(}\frac{(n-1)^{2}+(n-1)}{2}\Big{)}+n=3n^{2}+3n+1$

####  7.2.2 Simulations

Figure 12: Results histogram for a fine-grained biased QGB simulation. The
circuit is depicted in figure 13.

The third (lower) circuit in figure 13 is our example of a 4-QGB that has
control over each peg. Although we could only simulate this circuit, we have
also included the OpenQASM in Appendix C.

The results are in figure 12. Although we are only simulating a uniform
75%:25% across all pegs in the QGB the results match our other simulations.

This worked example helps to demonstrate the function of the corrective $CNOT$
and $RESET$ gates we mentioned previously. The left-to-right top-to-bottom
arrangement of the quantum peg circuit modules means that our corrective
additions and resets occurs always on the lower parts of the circuit.

##  8 Conclusions

We have presented a more intuitive and efficient way of generating quantum
Galton boards. We have also shown that these circuits can be simulated
effectively, and are reasonably optimal, if of some limited output power. We
have also shown how these QGB implementations can be biased to produce other
distributions in a reliable way. We have also provided methodologies and code
to achieve these results locally, with a view to expanding research on these
quantum algorithms.

Clearly noise is a limiting factor for circuits using the Fredkin
(controlled-$SWAP$) gates on real hardware owing to the extensive transpiler
work that has to take place, significantly increasing the gate count that in
turn lifts the noise floor. Until such gates are supported natively or error-
corrected qubits are in play, our approach will have limitations.

What we have shown is that the true power of a quantum computer in simulating
this kind of behaviour lies in the randomness that allows the outputs to be
reliable given the fairly minimal circuit construction in play. Thus, we
consider there is value in utilising these QGB circuits as they provide
results based on high quality randomness.

Whilst other approaches give $2^{n}$ possible outputs for $n$ qubit
measurements on $\mathcal{O}(n^{2})$ resources in play, we note that it is
straightforward to transform this normally distributed output to another by
linearly increasing numbers of repeated runs of the circuit.

It is a natural fact that, given the current state of NISQ-era quantum
computers available to us, the most often thought of route for generating
normally-distributed outputs with quantum technology is to get quantum random
bits and feed them into a well known sampling algorithm such as the Ziggurat
[19] or others that are in common currency, see [16, Vol II.].

Though this route still seems to be the best possible option in practical
terms, our aim is that presenting this approach can lend more variety to the
construction of quantum circuits that solve the Galton Board simulation and
Universal Statistical simulation problems.

## References

  * [1] Aaronson, S., and Arkhipov, A. The Computational Complexity of Linear Optics.  In Research in Optical Sciences (2014), Optical Society of America, p. QTh1A.2.
  * [2] Abd El-Latif, A. A., Abd-El-Atty, B., Mehmood, I., Muhammad, K., Venegas-Andraca, S. E., and Peng, J. Quantum-Inspired Blockchain-Based Cybersecurity: Securing Smart Edge Utilities in IoT-Based Smart Cities.  Information Processing & Management 58, 4 (2021), 102549.
  * [3] Avrachenkov, K., Borkar, V. S., Kadavankandy, A., and Sreedharan, J. K. Revisiting random walk based sampling in networks: evasion of burn-in period and frequent regenerations.  Computational Social Networks 5, 1 (2018), 4.
  * [4] Bagrets, D., Kim, K. W., Barkhofen, S., De, S., Sperling, J., Silberhorn, C., Altland, A., and Micklitz, T. Probing the topological Anderson transition with quantum walks.  Physical Review Research 3, 2 (jun 2021), 23183.
  * [5] Blank, C., Park, D. K., and Petruccione, F. Quantum-enhanced analysis of discrete stochastic processes.  npj Quantum Information 7, 1 (2021), 126.
  * [6] Böttcher, L., and Porter, M. A. Classical and Quantum Random-Walk Centrality Measures in Multilayer Networks.  SIAM Journal on Applied Mathematics 81, 6 (jan 2021), 2704–2724.
  * [7] Bouwmeester, D., Marzoli, I., Karman, G. P., Schleich, W., and Woerdman, J. P. Optical Galton board.  Physical Review A 61, 1 (dec 1999), 13410.
  * [8] de Souza, L. S., de Carvalho, J. H. A., and Ferreira, T. A. E. Quantum Walk to Train a Classical Artificial Neural Network.  In 2019 8th Brazilian Conference on Intelligent Systems (BRACIS) (2019), pp. 836–841.
  * [9] Dernbach, S., Mohseni Kabir, A., Pal, S., Gepner, M., and Towsley, D. Quantum walk neural networks with feature dependent coins.  Applied Network Science 4 (sep 2019), 76.
  * [10] El-Latif, A. A. A., Abd-El-Atty, B., Venegas-Andraca, S. E., Elwahsh, H., Piran, M. J., Bashir, A. K., Song, O. Y., and Mazurczyk, W. Providing End-to-End Security Using Quantum Walks in IoT Networks.  IEEE Access 8 (2020), 92687–92696.
  * [11] El-Latif, A. A. A., Iliyasu, A. M., and Abd-El-Atty, B. An Efficient Visually Meaningful Quantum Walks-Based Encryption Scheme for Secure Data Transmission on IoT and Smart Applications, 2021.
  * [12] Emmanoulopoulos, D., and Dimoska, S. Quantum Machine Learning in Finance: Time Series Forecasting, 2022.
  * [13] Fama, E. F. Random Walks in Stock Market Prices.  Financial Analysts Journal 21, 5 (feb 1965), 55–59.
  * [14] Janani, T., and Brindha, M. A secure medical image transmission scheme aided by quantum representation.  Journal of Information Security and Applications 59 (2021), 102832\.
  * [15] Kadian, K., Garhwal, S., and Kumar, A. Quantum walk and its application domains: A systematic review.  Computer Science Review 41 (2021), 100419.
  * [16] Knuth, D. E. The Art of Computer Programming, Vol. 1: Fundamental Algorithms, third ed.  Addison-Wesley, Reading, Mass., 1997.
  * [17] Lahini, Y., Steinbrecher, G. R., Bookatz, A. D., and Englund, D. Quantum logic using correlated one-dimensional quantum walks.  npj Quantum Information 4, 1 (2018), 2.
  * [18] Li-Hua Lu, Y.-Q. L. Quantum Approach to Fast Protein-Folding Time.
  * [19] Marsaglia, G., and Tsang, W. W. The ziggurat method for generating random variables.  Journal of Statistical Software 5, 8 (2000).
  * [20] Max, E., Robert, K., J., M. L., Christoph, D., Matthias, H., and Alexander, S. Exploring complex graphs using three-dimensional quantum walks of correlated photons.  Science Advances 7, 9 (feb 2022), eabc5266.
  * [21] Nguyen, V., Orbell, S. B., Lennon, D. T., Moon, H., Vigneau, F., Camenzind, L. C., Yu, L., Zumbühl, D. M., Briggs, G. A. D., Osborne, M. A., Sejdinovic, D., and Ares, N. Deep reinforcement learning for efficient measurement of quantum devices.  npj Quantum Information 7, 1 (2021), 100.
  * [22] Pistoia, M., Ahmad, S. F., Ajagekar, A., Buts, A., Chakrabarti, S., Herman, D., Hu, S., Jena, A., Minssen, P., Niroula, P., Rattew, A., Sun, Y., and Yalovetzky, R. Quantum Machine Learning for Finance, 2021.
  * [23] Rattew, A. G., Sun, Y., Minssen, P., and Pistoia, M. The Efficient Preparation of Normal Distributions in Quantum Registers.  Quantum 5 (Dec. 2021), 609.
  * [24] Roehsner, M.-C., Kettlewell, J. A., Fitzsimons, J., and Walther, P. Probabilistic one-time programs using quantum entanglement.  npj Quantum Information 7, 1 (2021), 98.
  * [25] Segawa, E., and Yoshie, Y. Quantum search of matching on signed graphs.  Quantum Information Processing 20, 5 (2021), 182.
  * [26] Shenvi, N., Kempe, J., and Whaley, K. B. Quantum random-walk search algorithm.  Physical Review A 67, 5 (may 2003), 52307.
  * [27] Su, Q.-P., Zhang, Y., Yu, L., Zhou, J.-Q., Jin, J.-S., Xu, X.-Q., Xiong, S.-J., Xu, Q., Sun, Z., Chen, K., Nori, F., and Yang, C.-P. Experimental demonstration of quantum walks with initial superposition states.  npj Quantum Information 5, 1 (2019), 40.
  * [28] Weisstein, E. W. Galton board. From MathWorld—A Wolfram Web Resource.  Last visited on 10/1/2022.
  * [29] Xia, F., Liu, J., Nie, H., Fu, Y., Wan, L., and Kong, X. Random Walks: A Review of Algorithms and Applications.  IEEE Transactions on Emerging Topics in Computational Intelligence 4, 2 (2020), 95–107.
  * [30] Yi, B., and Bose, S. Information Flow as Causation Quantifier in Quantum Dynamics, 2022.

##  Appendix A OpenQASM for 4-QGB

We present the OpenQASM code for a 4-level QGB, as simulated on IBM-QX.

    OPENQASM 2.0;
    include "qelib1.inc";

    qreg q[10];
    creg c[10];

    reset q[0];
    x q[5];
    h q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    h q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[6],q[7];
    cx q[6],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    h q[0];
    cswap q[0],q[2],q[3];
    cx q[3],q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[7],q[8];
    cx q[7],q[0];
    cswap q[0],q[6],q[7];
    cx q[6],q[0];
    cswap q[0],q[5],q[6];
    cx q[5],q[0];
    cswap q[0],q[4],q[5];
    reset q[0];
    h q[0];
    cswap q[0],q[1],q[2];
    cx q[2],q[0];
    cswap q[0],q[2],q[3];
    cx q[3],q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    cx q[6],q[0];
    cswap q[0],q[6],q[7];
    cx q[7],q[0];
    cswap q[0],q[7],q[8];
    cx q[8],q[0];
    cswap q[0],q[8],q[9];
    measure q[1] -> c[1];
    measure q[2] -> c[2];
    measure q[3] -> c[3];
    measure q[4] -> c[4];
    measure q[5] -> c[5];
    measure q[6] -> c[6];
    measure q[7] -> c[7];
    measure q[8] -> c[8];
    measure q[9] -> c[9];

##  Appendix B OpenQASM for a Biased 4-QGB

We present the OpenQASM code for a 4-level Biased QGB using the circuit change
example in section 5.2.

    OPENQASM 2.0;
    include "qelib1.inc";

    qreg q[10];
    creg c[10];

    reset q[0];
    x q[5];
    rx(2*pi/3) q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[6],q[7];
    cx q[6],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[2],q[3];
    cx q[3],q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[7],q[8];
    cx q[7],q[0];
    cswap q[0],q[6],q[7];
    cx q[6],q[0];
    cswap q[0],q[5],q[6];
    cx q[5],q[0];
    cswap q[0],q[4],q[5];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[1],q[2];
    cx q[2],q[0];
    cswap q[0],q[2],q[3];
    cx q[3],q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    cx q[6],q[0];
    cswap q[0],q[6],q[7];
    cx q[7],q[0];
    cswap q[0],q[7],q[8];
    cx q[8],q[0];
    cswap q[0],q[8],q[9];
    measure q[1] -> c[1];
    measure q[2] -> c[2];
    measure q[3] -> c[3];
    measure q[4] -> c[4];
    measure q[5] -> c[5];
    measure q[6] -> c[6];
    measure q[7] -> c[7];
    measure q[8] -> c[8];
    measure q[9] -> c[9];

##  Appendix C OpenQASM for Fine-Grained Biased QGB

    OPENQASM 2.0;
    include "qelib1.inc";

    qreg q[10];
    creg c[10];
    reset q[0];
    x q[5];
    rx(2*pi/3) q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[5],q[6];
    cx q[6],q[0];
    cswap q[0],q[6],q[7];
    barrier q[0],q[1],q[2],q[3],q[4];
    barrier q[5],q[6],q[7],q[8],q[9];
    reset q[0];
    cx q[6],q[5];
    rx(2*pi/3) q[0];
    reset q[6];
    cswap q[0],q[2],q[3];
    cx q[3],q[0];
    cswap q[0],q[3],q[4];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[4],q[5];
    cx q[5],q[0];
    cswap q[0],q[5],q[6];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[6],q[7];
    cx q[7],q[0];
    cswap q[0],q[7],q[8];
    barrier q[0],q[1],q[2],q[3],q[4];
    barrier q[5],q[6],q[7],q[8],q[9];
    reset q[0];
    cx q[5],q[4];
    cx q[7],q[6];
    rx(2*pi/3) q[0];
    reset q[5];
    reset q[7];
    cswap q[0],q[1],q[2];
    cx q[2],q[0];
    cswap q[0],q[2],q[3];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[3],q[4];
    cx q[4],q[0];
    cswap q[0],q[4],q[5];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[5],q[6];
    cx q[6],q[0];
    cswap q[0],q[6],q[7];
    reset q[0];
    rx(2*pi/3) q[0];
    cswap q[0],q[7],q[8];
    cx q[8],q[0];
    cswap q[0],q[8],q[9];
    barrier q[0],q[1],q[2],q[3],q[4];
    barrier q[5],q[6],q[7],q[8],q[9];
    cx q[4],q[3];
    cx q[6],q[5];
    cx q[8],q[7];
    reset q[4];
    reset q[6];
    reset q[8];
    measure q[1] -> c[1];
    measure q[2] -> c[2];
    measure q[3] -> c[3];
    measure q[4] -> c[4];
    measure q[5] -> c[5];
    measure q[6] -> c[6];
    measure q[7] -> c[7];
    measure q[8] -> c[8];
    measure q[9] -> c[9];

##  Appendix D IBM-QX Circuit Diagrams

To aid understanding, figure 13 depicts the circuit used for simulation,
derived from the OpenQASM code above in Appendices A, B, and C.

Figure 13: The 4-QGB circuit (upper) and the biased 4-QGB circuit (middle) and
fine-grained biased QGB (lower).

Figure 14: Transpiled Quantum Peg circuits for the ibmq-manila quantum
computer. The unbiased (upper) and biased (lower) transpiled circuits are
presented.


11institutetext:  1 Department of Computer Science,
Johns Hopkins University, Baltimore, MD 21218, USA
11email<EMAIL_ADDRESS>
2 Department of Biomedical Engineering,
Johns Hopkins University School of Medicine, Baltimore, MD 21205, USA
3 Department of Electrical and Computer Engineering,
Johns Hopkins University, Baltimore, MD 21218, USA
4 Center for Neuroscience and Regenerative Medicine, Henry M. Jackson
Foundation, Bethesda, MD 20817, USA

# Deep filter bank regression for super-resolution of anisotropic MR brain
images

Samuel W. Remedios 11 Shuo Han 22 Yuan Xue 33 Aaron Carass 33
Trac D. Tran 33 Dzung L. Pham 44 and Jerry L. Prince 33

###### Abstract

In 2D multi-slice magnetic resonance (MR) acquisition, the through-plane
signals are typically of lower resolution than the in-plane signals. While
contemporary super-resolution (SR) methods aim to recover the underlying high-
resolution volume, the estimated high-frequency information is implicit via
end-to-end data-driven training rather than being explicitly stated and
sought. To address this, we reframe the SR problem statement in terms of
perfect reconstruction filter banks, enabling us to identify and directly
estimate the missing information. In this work, we propose a two-stage
approach to approximate the completion of a perfect reconstruction filter bank
corresponding to the anisotropic acquisition of a particular scan. In stage 1,
we estimate the missing filters using gradient descent and in stage 2, we use
deep networks to learn the mapping from coarse coefficients to detail
coefficients. In addition, the proposed formulation does not rely on external
training data, circumventing the need for domain shift correction. Under our
approach, SR performance is improved particularly in “slice gap” scenarios,
likely due to the constrained solution space imposed by the framework.

###### Keywords:

super-resolution filter bank MRI

##  1 Introduction

Anisotropic magnetic resonance (MR) images are those acquired with high in-
plane resolution and low through-plane resolution. It is common practice to
acquire anisotropic volumes in clinics as it reduces scan time and motion
artifacts while preserving signal-to-noise ratio. To improve through-plane
resolution, data-driven super-resolution (SR) methods have been developed on
MR volumes [12, 7, 1, 2]. The application of SR methods to estimate the
underlying isotropic volume has been shown to improve performance on
downstream tasks [13].

For 2D multi-slice protocols, the through-plane point-spread function (PSF) is
known as the slice profile. When the sampling step is an integer, the through-
plane signals of an acquired MR image can be modeled as a strided 1D
convolution between the slice profile and the object to be imaged [3, 8, 10].
Commonly, the separation between slices is equivalent to the full-width-at-
half-max (FWHM) of the slice profile, but volumes can also be acquired where
the slice separation is less than or greater than the slice profile FWHM,
corresponding to “slice overlap” and “slice gap” scenarios, respectively.

Data-driven SR methods usually simulate low-resolution (LR) data from high-
resolution (HR) data using an assumed slice profile [12, 7, 1, 2], or an
estimated slice profile according to the image data or acquisition [3]. In
either case, SR methods are generally formulated as a classical inverse
problem:

$y=Ax,$ (1)

where $y$ is the LR observation, $A$ is the degradation matrix, and $x$ is the
underlying HR image. Commonly, this is precisely how paired training data is
created for supervised machine learning methods; HR data is degraded by $A$ to
obtain the LR $y$ and weights $\theta$ of a parameterized function $\phi$
(e.g., a neural network) are learned such that $\phi_{\theta}(y)\approx x$.
However, under this framework there is no specification of information lost by
application of $A$; contemporary SR models train end-to-end and are directed
only by the dataset.

In our work, we propose an entirely novel SR framework based on perfect
reconstruction (PR) filter banks. From filter bank theory, PR of a signal $x$
is possible through an $M$-channel filter bank with a correct design of an
analysis bank $H$ and synthesis bank $F$ [11]. Under this formulation, we do
not change Eq. 1 but explicitly recognize our observation $y$ as the “coarse
approximation” filter bank coefficients and the missing information necessary
to recover $x$ as the “detail” coefficients (see Fig. 1). For reference, in
machine learning jargon, the analysis bank is an encoder, the synthesis bank
is a decoder, and the coarse approximation and detail coefficients are
analogous to a “latent space.”

Figure 1: The filter bank observation model. Both $y$ and $H_{0}$ (green) are
given and fixed. In stage 1, filters $H_{1},\ldots,H_{M-1}$ and
$F_{0},F_{1},\ldots,F_{M-1}$ (purple) are learned; in stage 2, a mapping from
$y$ to $d_{1},\ldots,d_{M-1}$ (gold) is learned.

The primary contribution of this work is to reformulate SR to isotropy of
2D-acquired MR volumes as a filter bank regression framework. The proposed
framework has several benefits. First, the observed low-frequency information
is untouched in the reconstruction; thus, our method explicitly synthesizes
the missing high frequencies and does not need to learn to preserve acquired
low frequency information. Second, the downsampling factor $M$ specifies the
number of channels in the $M$-channel filter bank, constraining the solution
space in tougher scenarios such as “slice gap” acquisition recovery. Third,
the analysis filters of PR filter banks necessarily introduce aliasing which
is canceled via the synthesis filters; therefore, we do not need to directly
handle the anti-aliasing of the observed image. Fourth, our architecture has a
dynamic capacity for lower-resolution images. The rationale behind the dynamic
capacity is intuitive: when fewer measurements are taken, more estimates must
be done in recovery and a more robust model is required. Fifth, our method
exploits the nature of anisotropic volumetric data; in-plane slices are HR
while through-plane slices are LR. Thus, we do not rely on external training
data and only need the in-plane HR data to perform internal supervision. In
the remainder of the paper, we describe this framework in detail, provide
implementation details, and evaluate against a state-of-the-art internally
supervised SR technique. We demonstrate the feasibility of formulating SR as
filter bank coefficient regression and believe it lays the foundation for
future theoretical and experimental work in SR of MR images.

##  2 Methods

The analysis bank $H$ and synthesis bank $F$ each consist of $M$ 1D filters
represented in the $z$-domain as $H_{k}$ and $F_{k}$, respectively, with
corresponding spatial domain representations $h_{k}$ and $f_{k}$. As
illustrated in Fig. 1, input signal $X(z)=\mathcal{Z}(x)$111$\mathcal{Z}(x)$
is the $\mathcal{Z}$-transform of $x$[11]. is filtered by $H_{k}$, then
decimated with $\downarrow M$ (keeping every $M^{\text{th}}$ entry) to produce
the corresponding coefficients. These coefficients exhibit aliasing and
distortion which are corrected by the synthesis filters [11]. Reconstruction
from coefficients comes from zero-insertion upsampling with $\uparrow M$,
passing through filters $F_{k}$, then summing across the $M$ channels.

Figure 2: This network architecture, used in the second stage of our
algorithm, has the same structure for both the generator and discriminator but
with different hyperparameters. All convolutional layers used a $3\times 3$
kernel. The generator and discriminator used $16$ and $2$ residual blocks,
respectively. The generator had $128\times M$ features per convolutional layer
while the discriminator had $64\times M$ features per convolutional layer. The
final convolution outputs $M-1$ channels corresponding to the missing filter
bank detail coefficients. The internal structure of the residual block is
encapsulated in green.

Traditional design of $M$-channel PR filter banks involves a deliberate choice
of a prototype low-pass filter $H_{0}$ such that modulations and alternations
of the prototype produce the remaining filters for both the analysis and
synthesis filter banks [11]. $M$ is also chosen based on the restrictions of
the problem at hand. However, for anisotropic 2D-acquired MRI, the slice
profile is the low-pass filter and as such we have a fixed, given $H_{0}$. The
separation between slices is equal to the FWHM of $h_{0}$ plus any further gap
between slices. We denote the slice separation as $M$, corresponding to the
number of channels in the PR filter bank. We use $A\|B$, read “A skip B”, to
denote a FWHM of $A$ mm and slice gap of $B$ mm and note that $M=A+B$. For
this preliminary work, we assume $A,B$, and $M$ are all integer and, without
loss of generality, assume that the in-plane resolution is $1\|0$.

Our goal is to estimate filters $H_{1},\ldots,H_{M-1}$ and
$F_{0},\ldots,F_{M-1}$ and the detail coefficients $d_{1},\ldots,d_{M-1}$
which lead to PR of $x$. We approach this problem in two stages. In stage 1,
we approximate the missing analysis and synthesis filters, assuming there
exists a set of filters to complete the $M$-channel PR filter bank given that
$H_{0}$ and $M$ are fixed and known ahead of time. These must be learned first
to establish the approximate PR filter bank conditions on the coefficient
space. Then, in stage 2, we perform a regression on the missing coefficients.
Both of these stages are optimized in a data-driven end-to-end fashion with
gradient descent. After training, our method is applied by regressing
$d_{1},\ldots,d_{M-1}$ from $y$ and feeding all coefficients through the
synthesis bank to produce $\hat{x}$, our estimate of the HR signal. The stage
2 coefficient regression occurs in 2D, so we construct our estimate of the 3D
volume by averaging stacked 2D predictions from the synthesis bank from both
cardinal planes containing the through-plane axis.

Stage 1: Filter Optimization Previous works assumed the slice profile is
Gaussian with FWHM equal to the slice separation [12, 7]; instead, we estimate
the slice profile, $H_{0}$, directly with
ESPRESO222https://github.com/shuohan/espreso2 [3]. We next aim to estimate the
filters $H_{1},\ldots,H_{M-1}$ and $F_{0},\ldots,F_{M-1}$. To achieve this, we
learn the spatial representations $h_{1},\ldots,h_{M-1}$ and
$f_{0},\ldots,f_{M-1}$ from 1D rows and columns drawn from the high resolution
in-plane slices of $y$, denoted $\mathcal{D}_{1}=\\{x_{i}\\}_{i=1}^{N}$. We
initialize these filters according to a cosine modulation [11] of $h_{0}$,
which is defined as

$f_{k}[n]=h_{k}[n]=h_{0}[n]\sqrt{\frac{2}{M}}\cos{\left[\left(k+\frac{1}{2}\right)\left(n+\frac{M+1}{2}\right)\frac{\pi}{M}\right]},$

for $k\in\\{1,2,\ldots,M-1\\}$. Accordingly, we initialize $f_{0}$ to $h_{0}$.
We estimate $\hat{x}_{i}$ by passing $x_{i}$ through the analysis and
synthesis banks, then (since the entire operation is differentiable) step all
filters except $h_{0}$ through gradient descent. The reconstruction error is
measured with mean squared error loss and the filters are updated based on the
AdamW [6] optimizer with a learning rate of $0.1$, the one-cycle learning rate
scheduler [9], and a batch size of $32$ for $100,000$ steps.

Stage 2: Coefficient Regression From stage 1, we have the analysis and
synthesis banks and now want to estimate the missing detail coefficients given
only the LR observation $y$. With the correct coefficients and synthesis
filters, PR of $x$ is possible. For this stage, we use 2D patches, in spite of
the 1D SR problem, as a type of “neighborhood regularization”. Let
$\mathcal{D}_{2}=\\{x_{i}\\}_{i=1}^{N}$, $x_{i}\in\mathbb{R}^{p\times pM}$;
i.e., the training set for stage 2 consists of 2D $p\times pM$ patches drawn
from the in-plane slices of $y$. The second dimension will be decimated by $M$
after passing through the analysis banks, resulting in
$y,d_{1},\ldots,d_{M-1}\in\mathbb{R}^{p\times p}$. We use the analysis bank
(learned in stage 1) to create training pairs
$\\{(y_{i},(d_{1},d_{2},\ldots,d_{M-1})_{i})\\}_{i=1}^{N}$ and fit a
convolutional neural network (CNN) $G:\mathbb{R}^{p\times
p}\rightarrow\mathbb{R}^{{p\times p}^{M-1}}$ to map $y_{i}$ to
$(d_{1},\ldots,d_{M-1})_{i}$. In this work, we set $p=32$. Since this is an
image-to-image translation task, we adopt the widely used approach proposed in
Pix2Pix [4] including the adversarial patch discriminator.

Empirically, we found more learnable parameters are needed with greater $M$.
Thus, our generator $G$ is a CNN illustrated in Fig. 2 with $16$ residual
blocks and $128\times M$ kernels of size $3\times 3$ per convolutional layer.
The discriminator $D$ has the same architecture but with only $2$ residual
blocks and $64\times M$ kernels per convolutional layer. Our final loss
function for stage 2 is identical to the loss proposed in [4] and is
calculated on the error in $(d_{1},\ldots,d_{M-1})_{i}$. We use the AdamW
optimizer [6] with a learning rate of $10^{-4}$ and the one-cycle learning
rate scheduler [9] for $500,000$ steps at a batch size of $32$.

##  3 Experiments and Results

Table 1: Mean $\pm$ std. dev. of volumetric PSNR values for stage 1 reconstruction of the $30$ subjects. “Self” indicates a reconstruction of the input low-resolution volume on which the filter bank was optimized, while “GT” indicates reconstruction of the isotropic ground truth volume. (L-R), (A-P), and (S-I) are the left-to-right, anterior-to-posterior, and superior-to-inferior directions, respectively. | Self (L-R) | Self (A-P) | GT (L-R) | GT (A-P) | GT (S-I)
---|---|---|---|---|---
$2\|0$ |  $62.24\pm 0.97$ |  $60.19\pm 3.74$ |  $60.63\pm 0.56$ |  $59.59\pm 2.54$ |  $55.47\pm 4.69$
$2\|1$ | $63.01\pm 4.91$ | $62.25\pm 5.09$ | $64.32\pm 0.63$ | $59.49\pm 5.52$ | $53.81\pm 6.50$
$2\|2$ | $62.57\pm 1.59$ | $57.93\pm 5.32$ | $60.62\pm 1.34$ | $59.31\pm 3.65$ | $52.09\pm 4.34$
$4\|0$ | $55.47\pm 3.81$ | $52.36\pm 5.32$ | $48.91\pm 4.65$ | $48.77\pm 4.68$ | $44.08\pm 4.78$
$4\|1$ | $53.03\pm 1.54$ | $50.31\pm 3.41$ | $44.19\pm 1.57$ | $45.65\pm 1.63$ | $44.28\pm 2.14$
$4\|2$ | $54.71\pm 2.61$ | $51.08\pm 4.51$ | $46.75\pm 2.83$ | $46.39\pm 3.27$ | $43.27\pm 2.80$
$6\|0$ | $49.97\pm 1.07$ | $40.18\pm 4.77$ | $40.14\pm 1.35$ | $41.04\pm 1.40$ | $35.76\pm 3.19$
$6\|1$ | $52.35\pm 0.55$ | $45.69\pm 5.24$ | $42.11\pm 0.84$ | $42.74\pm 1.25$ | $39.76\pm 3.47$
$6\|2$ | $53.17\pm 3.17$ | $49.11\pm 3.41$ | $43.66\pm 4.12$ | $44.87\pm 3.99$ | $41.50\pm 2.29$
Figure 3: Estimated PR filters from stage 1 for a single subject at $4\|1$ resolution ($M=5$) in the frequency domain. Note the amplitudes for analysis and synthesis banks are on different scales, DC is centered, and $h_{0}$ is estimated by ESPRESO [3]. Figure 4: Mid-sagittal slice for a representative subject at different resolutions and gaps for each method. The low resolution column is digitally upsampled with k-space zero-filling. $A\|B$ signifies a slice thickness of $A$ mm and a gap of $B$ mm. Fourier magnitude is displayed in dB on every other row. The top two rows correspond to $2\|0$ ($M=2$) for the MR slice and Fourier space, the second two rows are for $4\|1$ ($M=5$), and the bottom two rows are for $6\|2$ ($M=8$).  |  |
---|---|---
(a) |  | (b)
Figure 5: Quantitative metrics PSNR in (a) and SSIM in (b), computed over the
$30$ image volumes. Significance tests are performed between SMORE and our
proposed method with the Wilcoxon signed rank test; $\ast$ denotes $p$-values
$<0.05$; “ns” stands for “not significant”.

Experiments We performed two experiments to evaluate the efficacy of each
stage in our approach. We randomly selected $30$ T1-weighted MR brain volumes
from the OASIS-3 dataset [5] to validate both stages and simulated LR
acquisition via convolution with a Gaussian kernel with FWHM $\in\\{2,4,6\\}$
and slice gap $\in\\{0,1,2\\}$, yielding nine combinations of FHWM and slice
gap in total. Since $M=A+B$ for a scan of resolution $A\|B$,
$M\in\\{2,3,4,5,6,7,8\\}$. For these experiments, the HR plane was axial while
the cardinal LR planes were sagittal and coronal. We note that both stage 1
and stage 2 are trained for each LR volume separately as our proposed method
does not use external training data, but instead relies on the inherent
anisotropy in the multi-slice volume (i.e., HR in-plane and LR through-plane
data).

Stage 1 Results We trained stage 1 using both cardinal 1D directions from in-
plane data; that is, left-to-right (L-R) and anterior-to-posterior (A-P)
directions. We then performed 1D reconstruction along these cardinal
directions and collated all reconstructions into 3D volumes. In other words,
this is an evaluation of self-auto-encoding. The mean volumetric
reconstruction PSNR $\pm$ std. dev. across the $30$ subjects is shown in Table
1. In addition to applying the learned filters to the LR image itself, we
would also like to test the extent of signal recovery for the HR counterpart
that is the ground truth (GT) of the LR volume. Indeed, the coefficients
generated by our learned analysis bank are what we will regress in stage 2, so
a reconstruction of the GT is also shown in the right three columns of Table
1. This serves as a sort of “upper bound” on our super-resolution estimate and
also answers the question of how well internal training generalizes to
reconstruction of an isotropic volume.

We note that if we had attained PR filters, the PSNR would be $\infty$; our
estimates fall short of this. Notably, reconstruction performance drops in the
(S-I) direction; this is likely due to the fact that signals along this
direction were not included in the training data. Additionally, an example of
learned filters in the frequency domain for one resolution, $4\|1$ ($M=5$), is
shown in Fig. 3. Recall that the fixed filter $h_{0}$ is the slice selection
profile. We observe that our optimization approximated bandpass filters.

Stage 2 Results To evaluate stage 2, we compared our method to two approaches
which also do not rely on external training data: cubic b-spline interpolation
and SMORE [12], a state-of-the-art self-super-resolution technique for
anisotropic MR volumes. For a fair comparison and improving SMORE results,
SMORE was trained with the same slice profile that we use (the ESPRESO
estimate [3]) instead of a Gaussian slice profile used in the original paper.

Qualitative results are displayed in Fig. 4 of a mid-sagittal slice for a
representative subject at $2\|0$, $4\|1$, and $6\|2$. This subject is near the
median PSNR value for that resolution across the $30$ subjects evaluated in
our experiments and for which SMORE outperforms our method at $2\|0$, is on
par with our method at $4\|1$, and is outperformed by our method at $6\|2$.
Also shown in Fig. 4 is the corresponding Fourier space, and we see that our
proposed method includes more high frequencies than the other methods. For
quantitative results, PSNR and SSIM were calculated on entire volumes, as
illustrated in box plots in Fig. 5.

##  4 Discussion and conclusions

In this paper, we have presented a novel filter bank formulation for SR of
2D-acquired anisotropic MR volumes as the regression of filter-specified
missing detail coefficients in an $M$-channel PR filter bank that does not
change the low-frequency sub-bands of the acquired image. We would emphasize
that our approach establishes a new theoretic basis for SR. In theory, these
coefficients exist and give exact recovery of the underlying HR signal.
However, it is unknown whether a mapping of
$y\rightarrow(d_{1},\ldots,d_{M-1})$ exists, and whether it is possible to
find filters to complete the analysis and synthesis banks to guarantee PR. In
practice, we estimate these in two stages: stage 1 estimates the missing
analysis and synthesis filters towards PR and stage 2 trains a CNN to regress
the missing detail coefficients given the coarse approximation $y$. According
to our experiments, as the resolution worsens and slice gap increases our
proposed method better handles the SR task than the competitive approach,
validating the usefulness of our method for super resolving anisotropic MR
images with large slice gaps. Future work will include: 1) deeper
investigation into the limits of the training set in learning the regression;
2) the degree to which the mapping $G$ is valid; 3) more analysis of the
super-resolved frequency space; and 4) develop methods to exactly achieve or
better approximate PR. True PR filter banks should greatly improve the method,
as Table 1 serves as a type of “upper bound” for our method; regardless of the
quality of coefficient regression, even given the ideal ground truth
coefficients, reconstruction accuracy would be limited. Furthermore, our work
suffers two major shortcomings. First, our current assumptions are integer
slice thickness and slice separation, which is not always true in reality. To
address this, the use of fractional sampling rates with filter banks [11] may
be a promising research direction. Second, our model in stage 2 scales the
number of convolution kernels per layer by $M$. This induces a longer training
and testing time when the image is of lower resolution. For reference, SMORE
produced the SR volume in about 86 minutes on a single NVIDIA V100 regardless
of the input resolution, but our proposed method produced the SR volume in 27
minutes for $2\|0$, 85 minutes for $4\|1$, and 127 minutes for $6\|2$.
Additionally, further investigation into improved regression is needed—a model
which can better capture the necessary aliasing in the coefficient domain is
vital for PR.

##  5 Acknowledgements

This material is supported by the National Science Foundation Graduate
Research Fellowship under Grant No. DGE-1746891. Theoretical development is
partially supported by NIH ORIP grant R21 OD030163 and the Congressionally
Directed Medical Research Programs (CDMRP) grant MS190131. This work also
received support from National Multiple Sclerosis Society RG-1907-34570, CDMRP
W81XWH2010912, and the Department of Defense in the Center for Neuroscience
and Regenerative Medicine.

## References

  * [1] Chen, Y., Shi, F., Christodoulou, A.G., Xie, Y., Zhou, Z., Li, D.: Efficient and accurate MRI super-resolution using a generative adversarial network and 3D multi-level densely connected network. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2018. pp. 91–99. Springer International Publishing, Cham (2018)
  * [2] Du, J., He, Z., Wang, L., Gholipour, A., Zhou, Z., Chen, D., Jia, Y.: Super-resolution reconstruction of single anisotropic 3d MR images using residual convolutional neural network. Neurocomputing 392, 209–220 (2020)
  * [3] Han, S., Remedios, S., Carass, A., Schär, M., Prince, J.L.: MR slice profile estimation by learning to match internal patch distributions. In: International Conference on Information Processing in Medical Imaging. pp. 108–119. Springer (2021)
  * [4] Isola, P., Zhu, J.Y., Zhou, T., Efros, A.A.: Image-to-image translation with conditional adversarial networks. CVPR (2017)
  * [5] LaMontagne, P.J., Benzinger, T.L., Morris, J.C., Keefe, S., Hornbeck, R., Xiong, C., Grant, E., Hassenstab, J., Moulder, K., Vlassenko, A.G., Raichle, M.E., Cruchaga, C., Marcus, D.: OASIS-3: Longitudinal neuroimaging, clinical, and cognitive dataset for normal aging and alzheimer disease. medRxiv (2019)
  * [6] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (2017)
  * [7] Oktay, O., Bai, W., Lee, M., Guerrero, R., Kamnitsas, K., Caballero, J., de Marvao, A., Cook, S., O’Regan, D., Rueckert, D.: Multi-input cardiac image super-resolution using convolutional neural networks. In: Medical Image Computing and Computer-Assisted Intervention - MICCAI 2016. pp. 246–254. Springer International Publishing, Cham (2016)
  * [8] Prince, J.L., Links, J.M.: Medical Imaging Signals and Systems. Pearson (2015)
  * [9] Smith, L.N., Topin, N.: Super-convergence: Very fast training of neural networks using large learning rates. In: Artificial intelligence and machine learning for multi-domain operations applications. vol. 11006, p. 1100612. International Society for Optics and Photonics (2019)
  * [10] Sønderby, C.K., Caballero, J., Theis, L., Shi, W., Huszár, F.: Amortised map inference for image super-resolution. arXiv preprint arXiv:1610.04490 (2016)
  * [11] Strang, G., Nguyen, T.: Wavelets and filter banks. Wellesley - Cambridge Press (1997)
  * [12] Zhao, C., Dewey, B.E., Pham, D.L., Calabresi, P.A., Reich, D.S., Prince, J.L.: SMORE: a self-supervised anti-aliasing and super-resolution algorithm for MRI using deep learning. IEEE transactions on medical imaging 40(3), 805–817 (2020)
  * [13] Zhao, C., Shao, M., Carass, A., Li, H., Dewey, B.E., Ellingsen, L.M., Woo, J., Guttman, M.A., Blitz, A.M., Stone, M., Calabresi, P.A., Halperin, H., Prince, J.L.: Applications of a deep learning method for anti-aliasing and super-resolution in MRI. Magnetic Resonance Imaging 64, 132–141 (2019)