added files
This commit is contained in:
Morten Hjorth-Jensen
Binary file not shown.
@@ -0,0 +1,474 @@
|
||||
% Generated by IEEEtranN.bst, version: 1.14 (2015/08/26)
|
||||
\begin{thebibliography}{89}
|
||||
\providecommand{\natexlab}[1]{#1}
|
||||
\providecommand{\url}[1]{#1}
|
||||
\csname url@samestyle\endcsname
|
||||
\providecommand{\newblock}{\relax}
|
||||
\providecommand{\bibinfo}[2]{#2}
|
||||
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
|
||||
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
|
||||
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
|
||||
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
|
||||
\fontdimen4\font\relax}
|
||||
\providecommand{\BIBforeignlanguage}[2]{{%
|
||||
\expandafter\ifx\csname l@#1\endcsname\relax
|
||||
\typeout{** WARNING: IEEEtranN.bst: No hyphenation pattern has been}%
|
||||
\typeout{** loaded for the language `#1'. Using the pattern for}%
|
||||
\typeout{** the default language instead.}%
|
||||
\else
|
||||
\language=\csname l@#1\endcsname
|
||||
\fi
|
||||
#2}}
|
||||
\providecommand{\BIBdecl}{\relax}
|
||||
\BIBdecl
|
||||
|
||||
\bibitem[Ackley et~al.(1985)Ackley, Hinton, and Sejnowski]{HintonBM1985}
|
||||
D.~H. Ackley, G.~E. Hinton, and T.~J. Sejnowski, ``{A Learning Algorithm for
|
||||
Boltzmann Machines},'' \emph{Cognitive Science}, vol.~9, no.~1, 1985.
|
||||
|
||||
\bibitem[Du and Swamy(2019)]{Du2019BM}
|
||||
K.-L. Du and M.~N.~S. Swamy, ``{Boltzmann Machines},'' in \emph{Neural Networks
|
||||
and Statistical Learning}.\hskip 1em plus 0.5em minus 0.4em\relax London:
|
||||
Springer London, 2019.
|
||||
|
||||
\bibitem[Boltzmann(1877)]{Boltzmann1877}
|
||||
L.~Boltzmann, ``{{\"U}ber die Natur der Gasmolek{\"u}le},'' in
|
||||
\emph{Wissenschaftliche Abhandlungen, Vol. I, II, and III}, 1877.
|
||||
|
||||
\bibitem[Gibbs(1902)]{gibbs02}
|
||||
J.~W. Gibbs, \emph{Elementary principles in statistical mechanics}.\hskip 1em
|
||||
plus 0.5em minus 0.4em\relax Cambridge University Press, 1902.
|
||||
|
||||
\bibitem[Liu and Webb(2010)]{Liu2010}
|
||||
B.~Liu and G.~I. Webb, ``{Generative and Discriminative Learning},'' in
|
||||
\emph{Encyclopedia of Machine Learning}, C.~Sammut and G.~I. Webb, Eds.\hskip
|
||||
1em plus 0.5em minus 0.4em\relax Boston, MA: Springer US, 2010.
|
||||
|
||||
\bibitem[Carleo et~al.(2018)Carleo, Nomura, and Imada]{CarleoRBMsQManyBody18}
|
||||
G.~Carleo, Y.~Nomura, and M.~Imada, ``Constructing exact representations of
|
||||
quantum many-body systems with deep neural networks,'' \emph{Nature
|
||||
Communications}, vol.~9, 2018.
|
||||
|
||||
\bibitem[Carleo and Troyer(2017)]{CarleoRBMsQuantumManyBody17}
|
||||
G.~Carleo and M.~Troyer, ``Solving the quantum many-body problem with
|
||||
artificial neural networks,'' \emph{Science}, vol. 355, 2017.
|
||||
|
||||
\bibitem[Nomura et~al.(2017)Nomura, Darmawan, Yamaji, and Imada]{YusukeRBM17}
|
||||
Y.~Nomura, A.~Darmawan, Y.~Yamaji, and M.~Imada,
|
||||
``{Restricted-Boltzmann-Machine Learning for Solving Strongly Correlated
|
||||
Quantum Systems},'' \emph{Physical Review B}, vol.~96, 2017.
|
||||
|
||||
\bibitem[Anshu et~al.(2020)Anshu, Arunachalam, Kuwahara, and
|
||||
Soleimanifar]{AnshuSample-efficientQManyBody}
|
||||
A.~Anshu, S.~Arunachalam, T.~Kuwahara, and M.~Soleimanifar, ``Sample-efficient
|
||||
learning of quantum many-body systems,'' \emph{arXiv preprint -
|
||||
arXiv:2004.07266}, 2020.
|
||||
|
||||
\bibitem[Melko et~al.(2019)Melko, Carleo, Carrasquilla, and Cirac]{Melko2019}
|
||||
R.~G. Melko, G.~Carleo, J.~Carrasquilla, and J.~I. Cirac, ``{Restricted
|
||||
Boltzmann machines in quantum physics},'' \emph{Nature Physics}, vol.~15,
|
||||
no.~9, 2019.
|
||||
|
||||
\bibitem[Hrasko et~al.(2015)Hrasko, Pacheco, and
|
||||
Krohling]{HRASKO2015RBMTimeSeries}
|
||||
R.~Hrasko, A.~G. Pacheco, and R.~A. Krohling, ``{Time Series Prediction Using
|
||||
Restricted Boltzmann Machines and Backpropagation},'' \emph{3rd International
|
||||
Conference on Information Technology and Quantitative Management}, vol.~55,
|
||||
2015.
|
||||
|
||||
\bibitem[Tubiana et~al.(2019)Tubiana, Cocco, and
|
||||
Monasson]{Tubiana19RBMProteins}
|
||||
J.~Tubiana, S.~Cocco, and R.~Monasson, ``{Learning Compositional
|
||||
Representations of Interacting Systems with Restricted Boltzmann Machines:
|
||||
Comparative Study of Lattice Proteins},'' \emph{Neural Computation}, vol.~31,
|
||||
2019.
|
||||
|
||||
\bibitem[Liu et~al.(2013)Liu, Liu, Sun, Liu, and Wang]{LiuRBMsSocialNetworks13}
|
||||
F.~Liu, B.~Liu, C.~Sun, M.~Liu, and X.~Wang, \emph{{Deep Learning Approaches
|
||||
for Link Prediction in Social Network Services}}.\hskip 1em plus 0.5em minus
|
||||
0.4em\relax Springer Berlin Heidelberg, 2013.
|
||||
|
||||
\bibitem[Mohamed and Hinton(2010)]{Mohamed10RBMSignal}
|
||||
A.-r. Mohamed and G.~Hinton, ``{Phone recognition using Restricted Boltzmann
|
||||
Machines},'' \emph{IEEE International Conference on Acoustics, Speech and
|
||||
Signal Processing - Proceedings}, 2010.
|
||||
|
||||
\bibitem[{Assis} et~al.(2018){Assis}, {Pereira}, {Carrano}, {Ramos}, and
|
||||
{Dias}]{Assis18RBMFin}
|
||||
C.~A.~S. {Assis}, A.~C.~M. {Pereira}, E.~G. {Carrano}, R.~{Ramos}, and
|
||||
W.~{Dias}, ``{Restricted Boltzmann Machines for the Prediction of Trends in
|
||||
Financial Time Series},'' in \emph{International Joint Conference on Neural
|
||||
Networks}, 2018.
|
||||
|
||||
\bibitem[Carreira-Perpinan and Hinton(2005)]{Hinton05CD}
|
||||
M.~Carreira-Perpinan and G.~Hinton, ``{On Contrastive Divergence Learning},''
|
||||
\emph{Artificial Intelligence and Statistics}, 2005.
|
||||
|
||||
\bibitem[Murphy(2012)]{MurphyML12}
|
||||
K.~P. Murphy, \emph{Machine learning: a probabilistic perspective}.\hskip 1em
|
||||
plus 0.5em minus 0.4em\relax Cambridge, MA: MIT Press, 2012.
|
||||
|
||||
\bibitem[Hinton(2002)]{Hinton2002TrainingPO}
|
||||
G.~E. Hinton, ``{Training Products of Experts by Minimizing Contrastive
|
||||
Divergence},'' \emph{Neural Computation}, vol.~14, 2002.
|
||||
|
||||
\bibitem[Besag(1975)]{Besag1975}
|
||||
J.~Besag, ``{Statistical Analysis of Non-Lattice Data},'' \emph{Journal of the
|
||||
Royal Statistical Society. Series D (The Statistician)}, vol.~24, no.~3,
|
||||
1975.
|
||||
|
||||
\bibitem[Tieleman(2008)]{Tieleman08}
|
||||
T.~Tieleman, ``{Training Restricted Boltzmann Machines using Approximations to
|
||||
the Likelihood Gradient},'' \emph{Proceedings of the 25th International
|
||||
Conference on Machine Learning}, 2008.
|
||||
|
||||
\bibitem[Sutskever and Tieleman(2010)]{Sutskever10}
|
||||
I.~Sutskever and T.~Tieleman, ``{On the Convergence Properties of Contrastive
|
||||
Divergence},'' \emph{Proceedings of the Thirteenth International Conference
|
||||
on Artificial Intelligence and Statistics}, vol.~9, 2010.
|
||||
|
||||
\bibitem[Amin et~al.(2018)Amin, Andriyash, Rolfe, Kulchytskyy, and
|
||||
Melko]{QBMAmin18}
|
||||
M.~H. Amin, E.~Andriyash, J.~Rolfe, B.~Kulchytskyy, and R.~Melko, ``{Quantum
|
||||
Boltzmann Machine},'' \emph{Phys. Rev. X}, vol.~8, 2018.
|
||||
|
||||
\bibitem[Ansch{\"u}tz and Cao(2019)]{Anschtz2019RealizingQB}
|
||||
E.~R. Ansch{\"u}tz and Y.~Cao, ``{Realizing Quantum Boltzmann Machines Through
|
||||
Eigenstate Thermalization},'' \emph{arXiv preprint - arXiv:1903.01359}, 2019.
|
||||
|
||||
\bibitem[Kieferov\'a and Wiebe(2017)]{QBMWiebe17}
|
||||
M.~Kieferov\'a and N.~Wiebe, ``{Tomography and Generative Training with Quantum
|
||||
Boltzmann Machines},'' \emph{Phys. Rev. A}, vol.~96, 2017.
|
||||
|
||||
\bibitem[Kappen(2018)]{Kappen18QBM}
|
||||
H.~Kappen, ``Learning quantum models from quantum or classical data,''
|
||||
\emph{arXiv preprint - arXiv:1803.11278}, 2018.
|
||||
|
||||
\bibitem[Wiebe and Wossnig(2019)]{Wiebe2019GenerativeTO}
|
||||
N.~Wiebe and L.~Wossnig, ``{Generative training of quantum Boltzmann machines
|
||||
with hidden units},'' \emph{arXiv preprint - arXiv:1905.09902}, 2019.
|
||||
|
||||
\bibitem[Thompson(1965)]{ThompsonInequality1965}
|
||||
C.~J. Thompson, ``{Inequality with Applications in Statistical Mechanics},''
|
||||
\emph{Journal of Mathematical Physics}, vol.~6, no.~11, 1965.
|
||||
|
||||
\bibitem[Golden(1965)]{Golden65Helm}
|
||||
\BIBentryALTinterwordspacing
|
||||
S.~Golden, ``Lower bounds for the helmholtz function,'' \emph{Phys. Rev.}, vol.
|
||||
137, pp. B1127--B1128, Feb 1965. [Online]. Available:
|
||||
\url{https://link.aps.org/doi/10.1103/PhysRev.137.B1127}
|
||||
\BIBentrySTDinterwordspacing
|
||||
|
||||
\bibitem[Troyer and Wiese(2005)]{Troyer05}
|
||||
M.~Troyer and U.-J. Wiese, ``{Computational Complexity and Fundamental
|
||||
Limitations to Fermionic Quantum Monte Carlo Simulations},'' \emph{Phys. Rev.
|
||||
Lett.}, vol.~94, 2005.
|
||||
|
||||
\bibitem[Hangleiter et~al.(2019)Hangleiter, Roth, Nagaj, and
|
||||
Eisert]{Hangleiter2019EasingTM}
|
||||
D.~Hangleiter, I.~Roth, D.~Nagaj, and J.~Eisert, ``{Easing the Monte Carlo sign
|
||||
problem},'' \emph{arXiv preprint - arXiv:1906.02309}, 2019.
|
||||
|
||||
\bibitem[Okunishi and Harada(2014)]{OkunishiSignProblem14}
|
||||
K.~Okunishi and K.~Harada, ``Symmetry-protected topological order and
|
||||
negative-sign problem for $\mathrm{SO(}n)$ bilinear-biquadratic chains,''
|
||||
\emph{Phys. Rev. B}, vol.~89, 2014.
|
||||
|
||||
\bibitem[Li et~al.(2016)Li, Jiang, and Yao]{Li2016SignProblemFreeQMC}
|
||||
Z.-X. Li, Y.-f. Jiang, and H.~Yao, ``{Majorana-Time-Reversal Symmetries: A
|
||||
Fundamental Principle for Sign-Problem-Free Quantum Monte Carlo
|
||||
Simulations},'' \emph{Physical Review Letters}, vol. 117, 2016.
|
||||
|
||||
\bibitem[Alet et~al.(2016)Alet, Damle, and Pujari]{SignProblemAlet16}
|
||||
F.~Alet, K.~Damle, and S.~Pujari, ``{Sign-Problem-Free Monte Carlo Simulation
|
||||
of Certain Frustrated Quantum Magnets},'' \emph{Phys. Rev. Lett.}, vol. 117,
|
||||
2016.
|
||||
|
||||
\bibitem[Li et~al.(2015)Li, Jiang, and Yao]{PhysRevB.91.241117}
|
||||
Z.-X. Li, Y.-F. Jiang, and H.~Yao, ``{Solving the fermion sign problem in
|
||||
quantum Monte Carlo simulations by Majorana representation},'' \emph{Phys.
|
||||
Rev. B}, vol.~91, 2015.
|
||||
|
||||
\bibitem[Ortiz et~al.(2001)Ortiz, Gubernatis, Knill, and
|
||||
Laflamme]{OrtizQAFermionic01}
|
||||
G.~Ortiz, J.~Gubernatis, E.~Knill, and R.~Laflamme, ``Quantum algorithms for
|
||||
fermionic simulations,'' \emph{Physical Review A}, vol.~64, 2001.
|
||||
|
||||
\bibitem[McArdle et~al.(2019)McArdle, Jones, Endo, Li, Benjamin, and
|
||||
Yuan]{VarSITEMcArdle19}
|
||||
S.~McArdle, T.~Jones, S.~Endo, Y.~Li, S.~C. Benjamin, and X.~Yuan,
|
||||
``{Variational ansatz-based quantum simulation of imaginary time
|
||||
evolution},'' \emph{npj Quantum Information}, vol.~5, no.~1, p.~75, 2019.
|
||||
|
||||
\bibitem[Yuan et~al.(2019)Yuan, Endo, Zhao, Benjamin, and
|
||||
Li]{Simon18TheoryVarQSim}
|
||||
X.~Yuan, S.~Endo, Q.~Zhao, S.~Benjamin, and Y.~Li, ``Theory of variational
|
||||
quantum simulation,'' \emph{Quantum}, vol. 3, 191, 2019.
|
||||
|
||||
\bibitem[McLachlan(1964)]{McLachlan64}
|
||||
A.~McLachlan, ``{A variational solution of the time-dependent Schr\"{o}dinger
|
||||
equation},'' \emph{Molecular Physics}, vol.~8, no.~1, 1964.
|
||||
|
||||
\bibitem[Ising(1925)]{Ising1925}
|
||||
E.~Ising, ``{Beitrag zur Theorie des Ferromagnetismus},'' \emph{Zeitschrift
|
||||
f{\"u}r Physik}, vol.~31, no.~1, 1925.
|
||||
|
||||
\bibitem[Peierls(1936)]{peierls_1936}
|
||||
R.~Peierls, ``{On Ising's model of ferromagnetism},'' \emph{Mathematical
|
||||
Proceedings of the Cambridge Philosophical Society}, vol.~32, no.~3, 1936.
|
||||
|
||||
\bibitem[Younes(1996)]{YOUNES1996109}
|
||||
L.~Younes, ``{Synchronous Boltzmann Machines can be universal approximators},''
|
||||
\emph{Applied Mathematics Letters}, vol.~9, no.~3, 1996.
|
||||
|
||||
\bibitem[Fischer and Igel(2012{\natexlab{a}})]{Fischer12RBM}
|
||||
A.~Fischer and C.~Igel, ``{An Introduction to Restricted Boltzmann Machines},''
|
||||
in \emph{Progress in Pattern Recognition, Image Analysis, Computer Vision,
|
||||
and Applications}, L.~Alvarez, M.~Mejail, L.~Gomez, and J.~Jacobo, Eds.\hskip
|
||||
1em plus 0.5em minus 0.4em\relax Springer Berlin Heidelberg, 2012.
|
||||
|
||||
\bibitem[Roux and Bengio(2010)]{Roux10}
|
||||
N.~Roux and Y.~Bengio, ``{Deep Belief Networks Are Compact Universal
|
||||
Approximators},'' \emph{Neural computation}, vol.~22, 03 2010.
|
||||
|
||||
\bibitem[Mont{\'u}far(2018)]{RBM_Montufar_2018}
|
||||
G.~Mont{\'u}far, ``{Restricted Boltzmann Machines: Introduction and Review},''
|
||||
in \emph{{Information Geometry and Its Applications}}, N.~Ay, P.~Gibilisco,
|
||||
and F.~Mat{\'u}{\v{s}}, Eds.\hskip 1em plus 0.5em minus 0.4em\relax Springer
|
||||
International Publishing, 2018.
|
||||
|
||||
\bibitem[Hinton(2012)]{Hinton2012}
|
||||
G.~E. Hinton, ``{A Practical Guide to Training Restricted Boltzmann
|
||||
Machines},'' in \emph{Neural Networks: Tricks of the Trade: Second Edition},
|
||||
G.~Montavon, G.~B. Orr, and K.-R. M{\"u}ller, Eds.\hskip 1em plus 0.5em minus
|
||||
0.4em\relax Berlin, Heidelberg: Springer Berlin Heidelberg, 2012.
|
||||
|
||||
\bibitem[Fischer and Igel(2012{\natexlab{b}})]{Fischer2012}
|
||||
A.~Fischer and C.~Igel, ``{An Introduction to Restricted Boltzmann Machines},''
|
||||
in \emph{{Progress in Pattern Recognition, Image Analysis, Computer Vision,
|
||||
and Applications}}, L.~Alvarez, M.~Mejail, L.~Gomez, and J.~Jacobo,
|
||||
Eds.\hskip 1em plus 0.5em minus 0.4em\relax Berlin, Heidelberg: Springer
|
||||
Berlin Heidelberg, 2012.
|
||||
|
||||
\bibitem[Fischer(2015)]{Fischer2015}
|
||||
A.~Fischer, ``{Training Restricted Boltzmann Machines},'' \emph{KI -
|
||||
K{\"u}nstliche Intelligenz}, vol.~29, no.~4, 2015.
|
||||
|
||||
\bibitem[Magnus(1954)]{MagnusITE54}
|
||||
W.~Magnus, ``On the exponential solution of differential equations for a linear
|
||||
operator,'' \emph{Communications on Pure and Applied Mathematics}, vol.~7,
|
||||
no.~4, 1954.
|
||||
|
||||
\bibitem[Gupta et~al.(2002)Gupta, Roy, and Deb]{GuptaITE02}
|
||||
N.~Gupta, A.~K. Roy, and B.~M. Deb, ``One-dimensional multiple-well
|
||||
oscillators: A time-dependent quantum mechanical approach,'' \emph{Pramana},
|
||||
vol.~59, no.~4, 2002.
|
||||
|
||||
\bibitem[Auer et~al.(2001)Auer, Krotscheck, and Chin]{ITEAuer01}
|
||||
J.~Auer, E.~Krotscheck, and S.~A. Chin, ``{A fourth-order real-space algorithm
|
||||
for solving local Schr\"{o}dinger equations},'' \emph{The Journal of Chemical
|
||||
Physics}, vol. 115, no.~15, 2001.
|
||||
|
||||
\bibitem[Koczor and Benjamin(2019)]{QNGNonUnitary19Simon}
|
||||
B.~Koczor and S.~Benjamin, ``Quantum natural gradient generalised to
|
||||
non-unitary circuits,'' \emph{arXiv preprint - arXiv:1912.08660}, 2019.
|
||||
|
||||
\bibitem[Nielsen and Chuang(2010)]{nielsen10}
|
||||
M.~A. Nielsen and I.~L. Chuang, \emph{Quantum Computation and Quantum
|
||||
Information}.\hskip 1em plus 0.5em minus 0.4em\relax Cambridge University
|
||||
Press, 2010.
|
||||
|
||||
\bibitem[Somma et~al.(2002)Somma, Ortiz, Gubernatis, Knill, and
|
||||
Laflamme]{LaflammeSimulatingPhysPhenom02}
|
||||
R.~Somma, G.~Ortiz, J.~E. Gubernatis, E.~Knill, and R.~Laflamme, ``Simulating
|
||||
physical phenomena by quantum networks,'' \emph{Phys. Rev. A}, vol.~65, 2002.
|
||||
|
||||
\bibitem[Bravyi et~al.(2006)Bravyi, DiVincenzo, Oliveira, and
|
||||
Terhal]{bravyi06LocalHam}
|
||||
S.~Bravyi, D.~DiVincenzo, R.~Oliveira, and B.~Terhal, ``{The Complexity of
|
||||
Stoquastic Local Hamiltonian Problems},'' \emph{Quantum Information and
|
||||
Computation}, vol.~8, 2006.
|
||||
|
||||
\bibitem[Gibbs(2010)]{gibbs_2010}
|
||||
J.~W. Gibbs, \emph{Elementary Principles in Statistical Mechanics: Developed
|
||||
with Especial Reference to the Rational Foundation of Thermodynamics}, ser.
|
||||
Cambridge Library Collection - Mathematics.\hskip 1em plus 0.5em minus
|
||||
0.4em\relax Cambridge University Press, 2010.
|
||||
|
||||
\bibitem[Pauli(1927)]{Pauli1927}
|
||||
W.~Pauli, ``{{\"U}ber Gasentartung und Paramagnetismus},'' \emph{Zeitschrift
|
||||
f{\"u}r Physik}, vol.~41, no.~2, 1927.
|
||||
|
||||
\bibitem[Temme et~al.(2011)Temme, Osborne, Vollbrecht, Poulin, and
|
||||
Verstraete]{Temme2011QuantumMS}
|
||||
K.~Temme, T.~J. Osborne, K.~G.~H. Vollbrecht, D.~Poulin, and F.~Verstraete,
|
||||
``{Quantum Metropolis Sampling},'' \emph{Nature}, vol. 471, 2011.
|
||||
|
||||
\bibitem[Yung and Aspuru-Guzik(2012)]{YungQuantumMetropolis12}
|
||||
M.-H. Yung and A.~Aspuru-Guzik, ``{A quantum{\textendash}quantum Metropolis
|
||||
algorithm},'' \emph{Proceedings of the National Academy of Sciences}, vol.
|
||||
109, no.~3, 2012.
|
||||
|
||||
\bibitem[Poulin and Wocjan(2009)]{PoulinThermalQGibbs09}
|
||||
D.~Poulin and P.~Wocjan, ``{Sampling from the Thermal Quantum Gibbs State and
|
||||
Evaluating Partition Functions with a Quantum Computer},'' \emph{Phys. Rev.
|
||||
Lett.}, vol. 103, 2009.
|
||||
|
||||
\bibitem[Abrams and Lloyd(1999)]{AbramsQPE99}
|
||||
D.~S. Abrams and S.~Lloyd, ``{Quantum Algorithm Providing Exponential Speed
|
||||
Increase for Finding Eigenvalues and Eigenvectors},'' \emph{Phys. Rev.
|
||||
Lett.}, vol.~83, 1999.
|
||||
|
||||
\bibitem[Motta and et~al.(2020)]{MottaQITE20}
|
||||
M.~Motta and et~al., ``Determining eigenstates and thermal states on a quantum
|
||||
computer using quantum imaginary time evolution,'' \emph{Nature Physics},
|
||||
vol.~16, no.~2, 2020.
|
||||
|
||||
\bibitem[Brand{\~a}o and
|
||||
Kastoryano(2019)]{brandaoFiniteCorrLengthEfficientPrep19}
|
||||
F.~G. S.~L. Brand{\~a}o and M.~J. Kastoryano, ``{Finite Correlation Length
|
||||
Implies Efficient Preparation of Quantum Thermal States},''
|
||||
\emph{Communications in Mathematical Physics}, vol. 365, no.~1, 2019.
|
||||
|
||||
\bibitem[Kastoryano and Brand{\~a}o(2016)]{BrandaoGibbsSampler16}
|
||||
M.~J. Kastoryano and F.~G. S.~L. Brand{\~a}o, ``{Quantum Gibbs Samplers: The
|
||||
Commuting Case},'' \emph{Communications in Mathematical Physics}, vol. 344,
|
||||
no.~3, 2016.
|
||||
|
||||
\bibitem[Chowdhury et~al.(2020)Chowdhury, Low, and
|
||||
Wiebe]{WiebeVariationalGibbs2020}
|
||||
A.~Chowdhury, G.~H. Low, and N.~Wiebe, ``{A Variational Quantum Algorithm for
|
||||
Preparing Quantum Gibbs States},'' \emph{arXiv preprint - arXiv:2002.00055},
|
||||
2020.
|
||||
|
||||
\bibitem[Brassard et~al.(2002)Brassard, Hoyer, Mosca, and Tapp]{brassardQAE02}
|
||||
G.~Brassard, P.~Hoyer, M.~Mosca, and A.~Tapp, ``{Quantum Amplitude
|
||||
Amplification and Estimation},'' \emph{Contemporary Mathematics}, vol. 305,
|
||||
2002.
|
||||
|
||||
\bibitem[Noh et~al.(2017)Noh, You, Mun, and Han]{Noh2017RegularizingDN}
|
||||
H.~Noh, T.~You, J.~Mun, and B.~Han, ``Regularizing deep neural networks by
|
||||
noise: Its interpretation and optimization,'' in \emph{NIPS}, 2017.
|
||||
|
||||
\bibitem[Farhi and Neven(2018)]{Farhi2018_gradients}
|
||||
E.~Farhi and H.~Neven, ``{Classification with Quantum Neural Networks on Near
|
||||
Term Processors},'' \emph{arXiv preprint - arXiv:1802.06002}, 2018.
|
||||
|
||||
\bibitem[Mitarai et~al.(2018)Mitarai, Negoro, Kitagawa, and
|
||||
Fujii]{Fujii2018_qcircuitLearn}
|
||||
K.~Mitarai, M.~Negoro, M.~Kitagawa, and K.~Fujii, ``Quantum circuit learning,''
|
||||
\emph{Phys. Rev. A}, vol.~98, 2018.
|
||||
|
||||
\bibitem[Dallaire-Demers and Killoran(2018)]{killoran2018}
|
||||
P.-L. Dallaire-Demers and N.~Killoran, ``Quantum generative adversarial
|
||||
networks,'' \emph{Phys. Rev. A}, vol.~98, 2018.
|
||||
|
||||
\bibitem[Schuld et~al.(2019)Schuld, Bergholm, Gogolin, Izaac, and
|
||||
Killoran]{SchuldQuantumGradients19}
|
||||
M.~Schuld, V.~Bergholm, C.~Gogolin, J.~Izaac, and N.~Killoran, ``Evaluating
|
||||
analytic gradients on quantum hardware,'' \emph{Phys. Rev. A}, vol.~99, 2019.
|
||||
|
||||
\bibitem[Zoufal et~al.(2019)Zoufal, Lucchi, and Woerner]{Zoufal2019}
|
||||
C.~Zoufal, A.~Lucchi, and S.~Woerner, ``Quantum generative adversarial networks
|
||||
for learning and loading random distributions,'' \emph{npj Quantum
|
||||
Information}, vol.~5, no.~1, 2019.
|
||||
|
||||
\bibitem[Dembo and Steihaug(1983)]{TNCDembo1983}
|
||||
R.~S. Dembo and T.~Steihaug, ``{Truncated-Newton algorithms for large-scale
|
||||
unconstrained optimization},'' \emph{Mathematical Programming}, vol.~26,
|
||||
no.~2, pp. 190--212, 1983.
|
||||
|
||||
\bibitem[Kingma and Ba(2014)]{Kingmaadam14}
|
||||
D.~Kingma and J.~Ba, ``{Adam: A Method for Stochastic Optimization},''
|
||||
\emph{International Conference on Learning Representations}, 2014.
|
||||
|
||||
\bibitem[ibm()]{ibmQX}
|
||||
\BIBentryALTinterwordspacing
|
||||
``{IBM Q Experience}.'' [Online]. Available:
|
||||
\url{https://quantumexperience.ng.bluemix.net/qx/experience}
|
||||
\BIBentrySTDinterwordspacing
|
||||
|
||||
\bibitem[Tikhonov et~al.(1995)Tikhonov, Goncharsky, Stepanov, and
|
||||
Yagola]{Tikhonov:1620560}
|
||||
A.~N. Tikhonov, A.~V. Goncharsky, V.~V. Stepanov, and A.~G. Yagola,
|
||||
\emph{{Numerical methods for the solution of ill-posed problems}}, ser.
|
||||
Mathematics and Its Applications.\hskip 1em plus 0.5em minus 0.4em\relax
|
||||
Springer, 1995.
|
||||
|
||||
\bibitem[Tibshirani(2011)]{lassoTibshirani11}
|
||||
R.~Tibshirani, ``Regression shrinkage and selection via the lasso: a
|
||||
retrospective,'' \emph{Journal of the Royal Statistical Society: Series B
|
||||
(Statistical Methodology)}, vol.~73, no.~3, 2011.
|
||||
|
||||
\bibitem[Hansen(2000)]{Hansen00thel-curve}
|
||||
P.~C. Hansen, \emph{{The L-Curve and its Use in the Numerical Treatment of
|
||||
Inverse Problems}}.\hskip 1em plus 0.5em minus 0.4em\relax WIT Press, 2000.
|
||||
|
||||
\bibitem[Aleksandrowicz and et~al.(2019)]{qiskit}
|
||||
G.~Aleksandrowicz and et~al., ``{Qiskit: An Open-source Framework for Quantum
|
||||
Computing},'' 2019.
|
||||
|
||||
\bibitem[Dewes et~al.(2012)Dewes, Ong, Schmitt, Lauro, Boulant, Bertet, Vion,
|
||||
and Esteve]{dewes2012readout}
|
||||
A.~Dewes, F.~R. Ong, V.~Schmitt, R.~Lauro, N.~Boulant, P.~Bertet, D.~Vion, and
|
||||
D.~Esteve, ``{Characterization of a Two-Transmon Processor with Individual
|
||||
Single-Shot Qubit Readout},'' \emph{Phys. Rev. Lett.}, vol. 108, 2012.
|
||||
|
||||
\bibitem[Stamatopoulos et~al.(2019)Stamatopoulos, Egger, Sun, Zoufal, Iten,
|
||||
Shen, and Woerner]{Stamatopoulos2019}
|
||||
N.~Stamatopoulos, D.~J. Egger, Y.~Sun, C.~Zoufal, R.~Iten, N.~Shen, and
|
||||
S.~Woerner, ``Option pricing using quantum computers,''
|
||||
\emph{arXiv:1905.02666}, 2019.
|
||||
|
||||
\bibitem[Sashank~J. et~al.(2018)Sashank~J., Satyen, and Sanjiv]{amsgrad}
|
||||
R.~Sashank~J., K.~Satyen, and K.~Sanjiv, ``{On the Convergence of Adam and
|
||||
Beyond},'' \emph{International Conference on Learning Representations}, 2018.
|
||||
|
||||
\bibitem[Altman(2019)]{Altman2019}
|
||||
E.~Altman, ``Synthesizing credit card transactions,'' \emph{arXiv preprint -
|
||||
arXiv:1910.03033}, 2019.
|
||||
|
||||
\bibitem[Pedregosa et~al.(2011)Pedregosa, Varoquaux, Gramfort, Michel, Thirion,
|
||||
Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos,
|
||||
Cournapeau, Brucher, Perrot, and Duchesnay]{scikit-learn2011}
|
||||
F.~Pedregosa, G.~Varoquaux, A.~Gramfort, V.~Michel, B.~Thirion, O.~Grisel,
|
||||
M.~Blondel, P.~Prettenhofer, R.~Weiss, V.~Dubourg, J.~Vanderplas, A.~Passos,
|
||||
D.~Cournapeau, M.~Brucher, M.~Perrot, and E.~Duchesnay, ``Scikit-learn:
|
||||
Machine learning in {P}ython,'' \emph{Journal of Machine Learning Research},
|
||||
vol.~12, 2011.
|
||||
|
||||
\bibitem[sci()]{scikitClassifierComp}
|
||||
\BIBentryALTinterwordspacing
|
||||
``{Classifier Comparison}.'' [Online]. Available:
|
||||
\url{https://scikit-learn.org/stable/auto_examples/classification/plot_classifier_comparison.html}
|
||||
\BIBentrySTDinterwordspacing
|
||||
|
||||
\bibitem[Gentile et~al.(2020)Gentile, Flynn, Knauer, Wiebe, Paesani, Granade,
|
||||
Rarity, Santagati, and Laing]{LearningModels20}
|
||||
A.~A. Gentile, B.~Flynn, S.~Knauer, N.~Wiebe, S.~Paesani, C.~Granade,
|
||||
J.~Rarity, R.~Santagati, and A.~Laing, ``Learning models of quantum systems
|
||||
from experiments,'' \emph{arXiv preprint - arXiv:2002.06169}, 2020.
|
||||
|
||||
\bibitem[Spieksma(1995)]{Spieksma1995}
|
||||
F.~C.~R. Spieksma, ``{Boltzmann Machines},'' in \emph{Artificial Neural
|
||||
Networks: An Introduction to ANN Theory and Practice}, P.~J. Braspenning,
|
||||
F.~Thuijsman, and A.~J. M.~M. Weijters, Eds.\hskip 1em plus 0.5em minus
|
||||
0.4em\relax Springer Berlin Heidelberg, 1995.
|
||||
|
||||
\bibitem[Farhi et~al.(2014)Farhi, Goldstone, and Gutmann]{FarhiQAOA14}
|
||||
E.~Farhi, J.~Goldstone, and S.~Gutmann, ``{A Quantum Approximate Optimization
|
||||
Algorithm Applied to a Bounded Occurrence Constraint Problem},'' \emph{arXiv
|
||||
preprint - arXiv:1411.4028}, 2014.
|
||||
|
||||
\bibitem[Barkoutsos et~al.(2020)Barkoutsos, Nannicini, Robert, Tavernelli, and
|
||||
Woerner]{Barkoutsos20VQOCVaR}
|
||||
P.~Barkoutsos, G.~Nannicini, A.~Robert, I.~Tavernelli, and S.~Woerner,
|
||||
``{Improving Variational Quantum Optimization using CVaR},'' \emph{Quantum},
|
||||
vol.~4, 2020.
|
||||
|
||||
\bibitem[Kardestuncer(1975)]{Kardestuncer1975}
|
||||
H.~Kardestuncer, ``{Finite Differences},'' in \emph{{Discrete Mechanics A
|
||||
Unified Approach}}.\hskip 1em plus 0.5em minus 0.4em\relax Springer Vienna,
|
||||
1975.
|
||||
|
||||
\end{thebibliography}
|
||||
@@ -0,0 +1,815 @@
|
||||
\documentclass[twocolumn, aps, pra, superscriptaddress, floatfix]{revtex4}
|
||||
\usepackage[T1]{fontenc}
|
||||
\usepackage[latin9]{inputenc}
|
||||
\usepackage{color}
|
||||
\usepackage{bm}% bold math
|
||||
\usepackage{xcolor}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{amssymb}
|
||||
\usepackage{mathtools}
|
||||
\usepackage{graphicx}
|
||||
\usepackage[sort&compress]{natbib}
|
||||
\usepackage{braket}
|
||||
\usepackage{psfrag}
|
||||
\usepackage{tikz}
|
||||
\usepackage[normalem]{ulem}
|
||||
\usepackage{qcircuit}
|
||||
\usepackage{subcaption}
|
||||
\captionsetup{compatibility=false}
|
||||
\usepackage{accents}
|
||||
\usepackage{cleveref}
|
||||
\usepackage{soul}
|
||||
\usepackage{algorithm}
|
||||
\usepackage{algpseudocode}
|
||||
\usepackage{multirow}
|
||||
\usepackage{makecell}
|
||||
|
||||
% Number fields
|
||||
\newcommand{\N}{\mathbb{N}}
|
||||
\newcommand{\Z}{\mathbb{Z}}
|
||||
\newcommand{\R}{\text{Re}}
|
||||
\newcommand{\C}{\mathbb{C}}
|
||||
\DeclareMathOperator{\EX}{\mathbb{E}}
|
||||
|
||||
% Comments
|
||||
\newcommand{\OUF}[1]{\textcolor{orange}{#1}}
|
||||
\newcommand{\WOR}[1]{\textcolor{red}{#1}}
|
||||
\newcommand{\citneeded}{\textcolor{blue}{$^\text{[citation needed]}$}}
|
||||
\newcommand{\varqbm}{VarQBM}
|
||||
|
||||
\begin{document}
|
||||
|
||||
|
||||
\title{Variational Quantum Boltzmann Machines}% Force line breaks
|
||||
|
||||
\author{Christa Zoufal}%
|
||||
\affiliation{IBM Quantum, IBM Research -- Zurich}%, Rueschlikon 8803, Switzerland}
|
||||
\affiliation{ETH Zurich}% 8092, Switzerland}
|
||||
|
||||
\author{Aur\'{e}lien Lucchi}
|
||||
\affiliation{ETH Zurich}% 8092, Switzerland}
|
||||
|
||||
\author{Stefan Woerner}
|
||||
\email{wor@zurich.ibm.com}
|
||||
\affiliation{IBM Quantum, IBM Research -- Zurich}%, Rueschlikon 8803, Switzerland}
|
||||
|
||||
\date{\today}
|
||||
|
||||
\begin{abstract}
|
||||
This work presents a novel realization approach to Quantum Boltzmann Machines (QBMs).
|
||||
The preparation of the required Gibbs states, as well as the evaluation of the loss function's analytic gradient is based on Variational Quantum Imaginary Time Evolution, a technique that is typically used for ground state computation.
|
||||
In contrast to existing methods, this implementation facilitates near-term compatible QBM training with gradients of the actual loss function for arbitrary parameterized Hamiltonians which do not necessarily have to be fully-visible but may also include hidden units.
|
||||
The variational Gibbs state approximation is demonstrated with numerical simulations and experiments run on real quantum hardware provided by IBM Quantum.
|
||||
Furthermore, we illustrate the application of this variational QBM approach to generative and discriminative learning tasks using numerical simulation.
|
||||
\end{abstract}
|
||||
|
||||
\maketitle
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Introduction}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
Boltzmann Machines (BMs) \cite{HintonBM1985, Du2019BM} offer a powerful framework for modelling probability distributions.
|
||||
These types of neural networks use an undirected graph-structure to encode relevant information.
|
||||
More precisely, the respective information is stored in bias coefficients and connection weights of network nodes, which are typically related to binary spin-systems and grouped into those that determine the output, the visible nodes, and those that act as latent variables, the hidden nodes.
|
||||
Furthermore, the network structure is linked to an energy function which facilitates the definition of a probability distribution over the possible node configurations by using a concept from statistical mechanics, i.e., Gibbs states \cite{Boltzmann1877, gibbs02}.
|
||||
The aim of BM training is to learn a set of weights such that the resulting model approximates a target probability distribution which is implicitly given by training data.
|
||||
This setting can be formulated as discriminative as well as generative learning task \cite{Liu2010}.
|
||||
Applications have been studied in a large variety of domains such as the analysis of quantum many-body systems, statistics, biochemistry, social networks, signal processing and finance, see, e.g., \cite{CarleoRBMsQManyBody18, CarleoRBMsQuantumManyBody17, YusukeRBM17, AnshuSample-efficientQManyBody, Melko2019, HRASKO2015RBMTimeSeries, Tubiana19RBMProteins, LiuRBMsSocialNetworks13, Mohamed10RBMSignal, Assis18RBMFin}.
|
||||
However, BMs are complicated to train in practice because the loss function's derivative requires the evaluation of a normalization factor, the partition function, that is generally difficult to compute.
|
||||
Usually, it is approximated using Markov Chain Monte Carlo methods which may require long runtimes until convergence \cite{Hinton05CD, MurphyML12}. Alternatively, the gradients could be estimated approximately using contrastive divergence \cite{Hinton2002TrainingPO} or pseudo-likelihood \cite{Besag1975} potentially leading to inaccurate results \cite{Tieleman08, Sutskever10}.
|
||||
|
||||
%Existing work
|
||||
Quantum Boltzmann Machines (QBMs) \cite{QBMAmin18} are a natural adaption of BMs to the quantum computing framework. Instead of an energy function with nodes being represented by binary spin values, QBMs define the underlying network using a Hermitian operator, a parameterized Hamiltonian
|
||||
\begin{equation*}
|
||||
H_{\theta}=\sum_{i=0}^{p-1}\theta_ih_i,
|
||||
\end{equation*}
|
||||
with $\theta\in\mathbb{R}^p$ and $h_i=\bigotimes_{j=0}^{n-1}\sigma_{j, i}$ for $\sigma_{j, i}\in\set{I, X, Y, Z}$ acting on the $j^{\text{th}}$ qubit.
|
||||
The network nodes are hereby characterized by the Pauli matrices $\sigma_{j, i}$.
|
||||
This Hamiltonian relates to a quantum Gibbs state, $\rho^{\text{Gibbs}} = {e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)}}/{Z}$
|
||||
with $\text{k}_{\text{B}}$ and $\text{T}$ denoting the Boltzmann constant and the system temperature, and $Z=\text{Tr}\left[e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)}\right]$.
|
||||
It should be noted that those qubits which determine the model output are referred to as visible and those which act as latent variables as hidden qubits.
|
||||
The aim of the model is to learn Hamiltonian parameters such that the resulting Gibbs state reflects a given target system.
|
||||
In contrast to BMs, this framework allows the use of quantum structures which are potentially inaccessible classically.
|
||||
Equivalently to the classical model, QBMs are suitable for discriminative as well as generative learning.
|
||||
|
||||
We present here a QBM implementation that circumvents certain issues which emerged in former approaches.
|
||||
The first paper on QBMs \cite{QBMAmin18} and several subsequent works \cite{Anschtz2019RealizingQB, QBMWiebe17, Kappen18QBM, Wiebe2019GenerativeTO} are incompatible with efficient evaluation of the loss function's analytic gradients if the given model has hidden qubits and
|
||||
\begin{equation*}
|
||||
\exists j: \:\left[H_{\theta}, \frac{\partial H_{\theta}}{\partial\theta_j}\right] \neq 0.
|
||||
\end{equation*}
|
||||
Instead, the use of hidden qubits is either avoided, i.e., only fully-visible settings are considered \cite{QBMWiebe17, Kappen18QBM, Wiebe2019GenerativeTO}, or the gradients are computed with respect to an upper bound of the loss \cite{QBMAmin18, Anschtz2019RealizingQB, QBMWiebe17}, which is based on the Golden-Thompson inequality \cite{ThompsonInequality1965, Golden65Helm}.
|
||||
It should be noted that training with an upper bound, renders the use of transverse Hamiltonian components, i.e., off-diagonal Pauli terms, difficult and imposes restrictions on the compatible models.
|
||||
|
||||
Further, we would like to point out that, in general, it is not trivial to evaluate a QBM Hamiltonian with a classical computer, i.e., using exact simulation with Quantum Monte Carlo methods \cite{Troyer05}, because the underlying Hamiltonian can suffer from the so-called \emph{sign-problem} \cite{Hangleiter2019EasingTM, OkunishiSignProblem14, Li2016SignProblemFreeQMC, SignProblemAlet16, PhysRevB.91.241117}. As already discussed in \cite{OrtizQAFermionic01}, evaluations on quantum computers can avoid this problem.
|
||||
|
||||
%Our QBM approach
|
||||
Our QBM implementation works for generic Hamiltonians $H_{\theta}$ with real coefficients $\theta$ and arbitrary Pauli terms $h_i$, and furthermore, is compatible with near-term, gate-based quantum computers.
|
||||
The method exploits \emph{Variational Quantum Imaginary Time Evolution} \cite{VarSITEMcArdle19, Simon18TheoryVarQSim} (VarQITE), which is based on McLachlan's variational principle \cite{McLachlan64}, to not only prepare approximate Gibbs states, $\rho_{\omega}^{\text{Gibbs}}$, but also to train the model with gradients of the actual loss function.
|
||||
During each step of the training, we use VarQITE to generate an approximation to the Gibbs state underlying $H_{\theta}$ and to enable automatic differentiation for computing the gradient of the loss function which is needed to update $\theta$.
|
||||
This \emph{Variational QBM} algorithm (\varqbm) is inherently normalized which implies that the training does not require the explicit evaluation of the partition function.
|
||||
|
||||
We focus on training quantum Gibbs states whose sampling behavior reflects a classical probability distribution. However, the scheme could be easily adapted to an approximate quantum state preparation scheme by using a loss function which is based on the quantum relative entropy \cite{QBMWiebe17, Kappen18QBM, Wiebe2019GenerativeTO}.
|
||||
Hereby, the approximation to $\rho^{\text{Gibbs}}$ is fitted to a given target state $\rho^{\text{data}}$.
|
||||
Notably, this approach is not necessarily suitable for learning classical distributions. More precisely, we do not need to train a quantum state that captures all features of the density matrix $\rho^{\text{data}}$ but only those which determine the sampling probability.
|
||||
It follows that fitting the full density matrix may impede the training.
|
||||
|
||||
The remainder of this paper is structured as follows. Firstly, we review classical BMs and VarQITE in Sec.~\ref{sec:pre}.
|
||||
Then, we outline \varqbm{} in Sec.~\ref{sec:QBM}. Next, we illustrate the feasibility of the Gibbs state preparation and present QBM applications in Sec.~\ref{sec:results}.
|
||||
Finally, a conclusion and an outlook are given in Sec.~\ref{sec:discussion}.
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Preliminaries}
|
||||
\label{sec:pre}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
This section introduces the concepts which form the basis of our \varqbm{} algorithm. First, classical BMs are presented in Sec.~\ref{sec:classBM}. Then, we discuss VarQITE, the algorithm that \varqbm{} uses for approximate Gibbs state preparation, in Sec.~\ref{sec:VarQITE}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Boltzmann Machines}
|
||||
\label{sec:classBM}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
Here, we will briefly review the original concept of classical BMs \cite{HintonBM1985}.
|
||||
A BM represents a network model that stores the learned knowledge in connection weights between network nodes.
|
||||
More explicitly, the connection weights are trained to generate outcomes according to a probability distribution of interest, e.g., to generate samples which are similar to given training samples or to output correct labels depending on input data samples.
|
||||
|
||||
Typically, this type of neural network is related to an Ising-type model \cite{Ising1925, peierls_1936} such that each node $i$ corresponds to a binary variable $z_i \in \set{-1, +1}$.
|
||||
Now, the set of nodes may be split into visible and hidden nodes representing observed and latent variables, respectively.
|
||||
Furthermore, a certain configuration $z=\left\{v,\: h\right\}$ of all nodes -- visible and hidden -- determines an energy, which is given as
|
||||
\begin{equation*}
|
||||
E_{z = \left\{v,\: h\right\}}= -\sum\limits_i\tilde{\theta}_iz_i - \sum\limits_{i, j}\theta_{ij}z_iz_j,
|
||||
\end{equation*}
|
||||
with $\tilde{\theta}_i, \theta_{ij}\in\mathbb{R}$ denoting the weights and $z_i$ representing the value taken by node $i$. It should be noted that the parameters $\theta_{ij}$ correspond to the weights of connections between different nodes. More explicitly, if two nodes are connected in the network, then a respective term appears in the energy function.
|
||||
The probability to observe a configuration $v$ of the visible nodes is defined as \begin{equation}
|
||||
\label{eq:gibbs_dis}
|
||||
p^{BM}_v = \frac{e^{-E_v/\left(\text{k}_{\text{B}}\text{T}\right)}}{Z},
|
||||
\end{equation}
|
||||
where $E_v = \sum_h E_{z = \left\{v,\: h\right\}}$, $k_B$ is the Boltzmann constant, $T$ the system temperature and $Z$ the canonical partition function
|
||||
\begin{equation*}
|
||||
Z=\sum\limits_{z= \left\{v,\: h\right\}}e^{-E_z/\left(\text{k}_{\text{B}}\text{T}\right)}.
|
||||
\end{equation*}
|
||||
We would like to point out that BMs adopt a concept from statistical mechanics.
|
||||
Suppose a closed system that is in thermal equilibrium with a coupled heat bath at constant temperature. The possible configuration space is determined by the canonical ensemble, i.e., the probability for observing a configuration is given by the Gibbs distribution \cite{Boltzmann1877, gibbs02} which corresponds to Eq.~\eqref{eq:gibbs_dis}.
|
||||
|
||||
Now, the goal of a BM is to fit the target probability distribution $p^{\text{data}}$ with $p^{BM}$.
|
||||
Typically, this training objective is achieved by optimizing the cross-entropy
|
||||
\begin{equation}
|
||||
\label{eq:crossEnt}
|
||||
L = -\sum\limits_{v}p_v^{\text{data}}\log{p_v^{BM}}.
|
||||
\end{equation}
|
||||
In theory, fully-connected BMs have interesting representation capabilities \cite{HintonBM1985, YOUNES1996109, Fischer12RBM}, i.e., they are universal approximators \cite{Roux10}.
|
||||
However, in practice they are difficult to train as the optimization easily gets expensive.
|
||||
Thus, it has become common practice to restrict the connectivity between nodes which relates to restricted Boltzmann Machines (RBMs) \cite{RBM_Montufar_2018}.
|
||||
Furthermore, several approximation techniques, such as contrastive divergence \cite{Hinton2002TrainingPO}, have been developed to facilitate BM training.
|
||||
However, these approximation techniques typically still face issues such as long computation time due to a large amount of required Markov chain steps or poor compatibility with multimodal probability distributions \cite{MurphyML12}.
|
||||
For further details, we refer the interested reader to \cite{Hinton2012, Fischer2012, Fischer2015}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Variational Quantum Imaginary Time Evolution}
|
||||
\label{sec:VarQITE}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
Imaginary time evolution (ITE) \cite{MagnusITE54} is an approach that is well known for (classical) ground state computation \cite{VarSITEMcArdle19, GuptaITE02, ITEAuer01}.
|
||||
|
||||
Suppose a starting state $\ket{\psi_0}$ and a time-independent Hamiltonian $H=\sum_{i=0}^{p-1}\theta_ih_i$ with real coefficients $\theta_i$ and Pauli terms $h_i$. Then, the normalized ITE propagates $\ket{\psi_0}$ with respect to $H$ for time $\tau$ according to
|
||||
\begin{equation*}
|
||||
\label{eq:ITE}
|
||||
\ket{\psi_\tau}= C\left(\tau\right)e^{-H\tau}\ket{\psi_0},
|
||||
\end{equation*}
|
||||
where $C\left(\tau\right) = {1}/{\sqrt{\text{Tr}\left[e^{-2H\tau}\ket{\psi_0 }\bra{\psi_0}\right]}}$ is a normalization.
|
||||
The differential equation that describes this evolution is the Wick-rotated Schr\"odinger equation
|
||||
\begin{equation}
|
||||
\label{eq:WickSchroed}
|
||||
\frac{d\ket{\psi_\tau}}{d \tau} =-\left( H - E_{\tau} \right) \ket{\psi_\tau},
|
||||
\end{equation}
|
||||
where $E_{\tau} = \bra{\psi_{\tau}} H\ket{\psi_\tau}$ originates from the normalization of $\ket{\psi_{\tau}}$.
|
||||
The terms in $e^{-H\tau}$, corresponding to small eigenvalues of $H$, decay slower than the ones corresponding to large eigenvalues. Due to the continuous normalization, the smallest eigenvalue dominates for $\tau \rightarrow \infty$. Thus, $\ket{\psi_\tau}$ converges to the ground state of $H$ given that there is some overlap between the ground and starting state.
|
||||
Furthermore, if ITE is only evolved to a finite time, $\tau = {1}/{2\left(\text{k}_{\text{B}}\text{T}\right)}$, then it enables the preparation of Gibbs states, see Sec.~\ref{sec:gibbsVarQITE}.
|
||||
|
||||
As introduced in \cite{VarSITEMcArdle19, Simon18TheoryVarQSim}, an approximate ITE can be implemented on a gate-based quantum computer by using McLachlan's variational principle \cite{McLachlan64}.
|
||||
The basic idea of the method is to introduce a parameterized trial state $\ket{\psi_{\omega}}$ and to project the temporal evolution of $\ket{\psi_\tau}$ to the parameters, i.e., $\omega \coloneqq \omega(\tau)$.
|
||||
We refer to this algorithm as VarQITE and, now, discuss it in more detail.
|
||||
|
||||
First, we define an input state $\ket{\psi_{\text{in}}}$ and a quantum circuit $V\left(\omega\right) = U_q\left(\omega_q\right)\cdots U_1\left(\omega_1\right)$ with parameters $\omega \in \mathbb{R}^{q}$ to generate the parameterized trial state
|
||||
\begin{equation*}
|
||||
\label{eq:varSite}
|
||||
\ket{\psi_{\omega}} \coloneqq V\left(\omega\right)\ket{\psi_{\text{in}}}.
|
||||
\end{equation*}
|
||||
Now, McLachlan's variational principle
|
||||
\begin{equation}
|
||||
\label{eq:McLachlan}
|
||||
\delta \left\lVert \left(d/d\tau + H - E_{\tau}\right) \ket{ \psi_{\omega}} \right\rVert = 0
|
||||
\end{equation}
|
||||
determines the time propagation of the parameters $\omega(\tau)$.
|
||||
This principle aims to minimize the distance between the right hand side of Eq.~\eqref{eq:WickSchroed} and the change $d\ket{ \psi_{\omega}} / d \tau$.
|
||||
Eq.~\eqref{eq:McLachlan} leads to a system of linear equations for $\dot{\omega}= d \omega / d\tau$, i.e.,
|
||||
\begin{align}
|
||||
\label{eq:sle}
|
||||
A\dot{\omega} = C
|
||||
\end{align}
|
||||
with
|
||||
\begin{equation}
|
||||
\begin{split}
|
||||
\label{eq:AC}
|
||||
A_{pq}\left(\tau\right) &= \text{Re}\left(\text{Tr}\left[\frac{\partial V^{\dagger}\left(\omega\left(\tau\right)\right)}{\partial\omega\left(\tau\right)_p}\frac{\partial V\left(\omega\left(\tau\right)\right)}{\partial\omega\left(\tau\right)_q}\rho_{\text{in}}\right]\right) \\
|
||||
C_p\left(\tau\right) &= -\sum\limits_i\theta_i\text{Re}\left(\text{Tr}\left[\frac{\partial V^{\dagger}\left(\omega\left(\tau\right)\right)}{\partial\omega\left(\tau\right)_p}h_iV\left(\omega\left(\tau\right)\right)\rho_{\text{in}}\right]\right),
|
||||
\end{split}
|
||||
\end{equation}
|
||||
where $\text{Re}\left(\cdot\right)$ denotes the real part and $\rho_{\text{in}} = \ket{\psi_{\text{in}}}\bra{\psi_{\text{in}}}$.
|
||||
The vector $C$ describes the derivative of the system energy $ \bra{ \psi_{\omega}}H \ket{ \psi_{\omega}}$ and $A$ is proportional to the classical Fisher information matrix, a metric tensor that reflects the system's information geometry \cite{QNGNonUnitary19Simon}.
|
||||
To evaluate $A$ and $C$, we compute expectation values with respect to quantum circuits of a particular form which is illustrated and discussed in Appendix \ref{app:a_c}.
|
||||
|
||||
This evaluation is compatible with arbitrary parameterized unitaries in $V\left(\omega\right)$ because all unitaries can be written as $U\left(\omega\right) = e^{iM\left(\omega\right)}$, where $M\left(\omega\right)$ denotes a parameterized Hermitian matrix.
|
||||
Further, Hermitian matrices can be decomposed into weighted sums of Pauli terms, i.e., $M\left(\omega\right) = \sum_pm_p\left(\omega\right)h_p$ with $m_p\left(\omega\right)\in\mathbb{R}$ and $h_p=\bigotimes\limits_{j=0}^{n-1}\sigma_{j, p}$ for $\sigma_{j, p}\in\set{I, X, Y, Z}$ \cite{nielsen10} acting on the $j^{\text{th}}$ qubit. Thus, the gradients of
|
||||
$U_k\left(\omega_k\right)$ are given by
|
||||
\begin{equation*}
|
||||
\frac{\partial U_k\left(\omega_k\right)}{\partial\omega_k} = \sum\limits_pi \frac{\partial m_{k,p}\left(\omega_k\right)}{\partial\omega_k}U_k\left(\omega_k\right)h_{k_p}.
|
||||
\end{equation*}
|
||||
This decomposition allows us to compute $A$ and $C$ with the techniques described in \cite{LaflammeSimulatingPhysPhenom02, VarSITEMcArdle19, Simon18TheoryVarQSim}.
|
||||
Furthermore, it should be noted that Eq.~\eqref{eq:sle} is often ill-conditioned and may, thus, require the use of regularized regression methods, see Sec.~\ref{sec:impRel}.
|
||||
|
||||
Now, we can use, e.g., an explicit Euler method to evolve the parameters as
|
||||
\begin{equation*}
|
||||
\omega\left(\tau\right) = \omega\left(0\right) + \sum\limits_{j = 1}^{\tau/\delta\tau} \dot{\omega}\left(\tau\right)\delta\tau.
|
||||
\end{equation*}
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Quantum Boltzmann Machine Algorithm}
|
||||
\label{sec:QBM}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
A QBM is defined by a parameterized Hamiltonian $H_{\theta}=\sum_{i=0}^{p-1}\theta_ih_i $ where $\theta\in\mathbb{R}^p$ and $h_i=\bigotimes_{j=0}^{n-1}\sigma_{j, i}$ for $\sigma_{j, i}\in\set{I, X, Y, Z}$ acting on the $j^{\text{th}}$ qubit.
|
||||
Equivalently to classical BMs, QBMs are typically represented by an Ising model \cite{Ising1925}, i.e., a $2$-local system \cite{bravyi06LocalHam} with nearest-neighbor coupling that is defined with regard to a particular grid. In principle, however, any Hamiltonian compatible with Boltzmann distributions could be used.
|
||||
|
||||
In contrast to BMs, the network nodes, given by the Pauli terms $\sigma_{j, i}$, do not represent the visible and hidden units. These are defined with respect to certain sub-sets of qubits. More explicitly, those qubits which determine the output of the QBM are the visible qubits, whereas the others correspond to the hidden qubits.
|
||||
Now, the probability to measure a configuration $v$ of the visible qubits is defined with respect to a projective measurement $\Lambda_v = \ket{v}\bra{v}\otimes I$ on the quantum Gibbs state
|
||||
\begin{equation*}
|
||||
\label{eq:QGibbs}
|
||||
\rho^{\text{Gibbs}} = \frac{e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)}}{Z}
|
||||
\end{equation*}
|
||||
with $Z=\text{Tr}\left[e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)}\right]$, i.e., the probability to measure $\ket{v}$ is given by
|
||||
\begin{equation*}
|
||||
p_v^{\text{QBM}} = \text{Tr}\left[\Lambda_v\rho^{\text{Gibbs}} \right].
|
||||
\end{equation*}
|
||||
For the remainder of this work, we assume that $\Lambda_v$ refers to projective measurements with respect to the computational basis of the visible qubits. Thus, the configuration $v$ is determined by $v_i\in \set{0, 1}$.
|
||||
It should be noted that this formulation does not require the evaluation of the configuration of the hidden qubits.
|
||||
|
||||
Our goal is to train the Hamiltonian parameters $\theta$ such that the sampling probabilities of the corresponding $\rho^{\text{Gibbs}}$ reflect the probability distribution underlying given classical training data. For this purpose, the same loss function as described in the classical case, see Eq.~\eqref{eq:crossEnt}, can be used
|
||||
\begin{equation}
|
||||
\label{eq:loss_qbm}
|
||||
L = -\sum\limits_{v}p_v^{\text{data}}\log{p_v^{\text{QBM}}},
|
||||
\end{equation}
|
||||
where $p_v^{\text{data}}$ denotes the occurrence probability of item $v$ in the training data set.
|
||||
|
||||
To enable efficient training, we want to evaluate the derivative of $L$ with respect to the Hamiltonian parameters. Unlike existing QBM implementations, \varqbm{} facilitates the use of analytic gradients of the loss function given in Eq.~\eqref{eq:loss_qbm} for generic QBMs.
|
||||
The presented algorithm involves the following steps.
|
||||
First, we use VarQITE to approximate the Gibbs state, see Sec.~\ref{sec:gibbsVarQITE} for further details.
|
||||
Then, we compute the gradient of $L$ to update the parameters $\theta$ with automatic differentiation, as is discussed in Sec.~\ref{sec:implementation}. The parameters are trained with a classical optimization routine where one training step consists of the Gibbs state preparation with respect to the current parameter values and a consecutive parameter update, as illustrated in Fig.~\ref{fig:varqbm}.
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=0.8\linewidth]{varqbm.png}
|
||||
\end{center}
|
||||
\caption{The \varqbm{} training includes the following steps. First, we need to fix the Pauli terms for $H_{\theta}$ and choose initial parameters $\theta$. Then, VarQITE is used to generate $\rho_{\omega}^{\text{Gibbs}}$ and compute $\partial\omega/ \partial\theta$.
|
||||
The quantum state and the derivative are needed to evaluate $p_v^{\text{QBM}}$ and $\partial p_v^{\text{QBM}}/\partial\theta$. Now, we can find $\partial L/\partial\theta$ to update the Hamiltonian parameters with a classical optimizer.}
|
||||
\label{fig:varqbm}
|
||||
\end{figure}
|
||||
|
||||
In the remainder of this section, we discuss Gibbs state preparation with VarQITE in Sec.~\ref{sec:gibbsVarQITE} and VarQBM in more detail in Sec.~\ref{sec:implementation}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Gibbs State Preparation with VarQITE}
|
||||
\label{sec:gibbsVarQITE}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
The Gibbs state $\rho^{\text{Gibbs}}$ describes the probability density operator of the configuration space of a system in thermal equilibrium with a heat bath at constant temperature $T$ \cite{gibbs_2010}. Originally, Gibbs states were studied in the context of statistical mechanics but, as shown in \cite{Pauli1927}, the density operator also facilitates the description of quantum statistics.
|
||||
|
||||
Gibbs state preparation can be approached from different angles. Hereby, different techniques not only have different strengths but also different drawbacks.
|
||||
Some schemes \cite{Temme2011QuantumMS, YungQuantumMetropolis12, PoulinThermalQGibbs09} use Quantum Phase Estimation \cite{AbramsQPE99} as a subroutine, which is likely to require error-corrected quantum computers.
|
||||
Other methods enable the evaluation of quantum thermal averages \cite{MottaQITE20, brandaoFiniteCorrLengthEfficientPrep19, BrandaoGibbsSampler16} for states with finite correlations. However, since QBM-related states may exhibit long-range correlations, these methods are not the first choice for the respective preparation.
|
||||
A thermalization based approach is presented in \cite{Anschtz2019RealizingQB}, where the aim is to prepare a quantum Gibbs state by coupling the state register to a heat bath given in the form of an ancillary quantum register.
|
||||
Correct preparation requires a thorough study of suitable ancillary registers for a generic Hamiltonian as the most useful ancilla system is not a-priori known.
|
||||
Further, a variational Gibbs state preparation method has been presented \cite{WiebeVariationalGibbs2020} which is based on the fact that Gibbs states minimize the free energy of a system at constant temperature.
|
||||
Thus, the goal is to fit a parameterized quantum state such that it minimizes the free energy.
|
||||
The parameter update is hereby conducted with a finite difference method instead of analytic gradients which may impair the training accuracy. Additionally, the method requires the application of Quantum Amplitude Estimation \cite{brassardQAE02}, as well as matrix exponentiation of the input state, and thus, is not well suited for near-term quantum computing applications.
|
||||
|
||||
In contrast to these Gibbs state preparation schemes, VarQITE is compatible with near-term quantum computers, and is neither limited to states with finite correlations nor requires ambiguous ancillary systems.
|
||||
In the following, we discuss how VarQITE can be utilized to generate an approximation of the Gibbs state $\rho^{\text{Gibbs}}$ for a generic $n-$qubit Hamiltonian $H_{\theta}=\sum_{i=0}^{p-1}\theta_ih_i $ with $\theta\in\mathbb{R}^p$ and $h_i=\bigotimes_{j=0}^{n-1}\sigma_{j, i}$ for $\sigma_{j, i}\in\set{I, X, Y, Z}$ acting on the $j_i^{\text{th}}$ qubit.
|
||||
|
||||
First, we need to choose a suitable variational quantum circuit $V\left(\omega\right)$, $\omega \in \mathbb{R}^{q}$, and set of initial parameters $\omega(0)$ such that the initial state is
|
||||
\begin{equation*}
|
||||
\ket{ \psi_{0}} = V\left(\omega\left(0\right)\right)\ket{0}^{\otimes 2n} = \ket{\phi^{+}}^{\otimes n}
|
||||
\end{equation*}
|
||||
where $\ket{\phi^{+}} = \frac{1}{\sqrt{2}}\left(\ket{00}+\ket{11}\right)$ represents a Bell state.
|
||||
We define two $n$-qubit sub-systems $a$ and $b$ such that the first and second qubit of each $\ket{\phi^{+}}$ is in $a$ and $b$, respectively. Accordingly, an effective $2n$-qubit Hamiltonian $H_{\text{eff}} = H_{\theta}^a + I^b$, where $H_{\theta}$ and $I$ act on sub-system $a$ and $b$, is considered.
|
||||
It should be noted that tracing out sub-system $b$ from $\ket{ \psi_{0}}$ results in an $n$-dimensional maximally mixed state
|
||||
\begin{equation*}
|
||||
\text{Tr}_{b}\left[\ket{\phi^{+}}^{\otimes n} \right] = \frac{1}{2^n}I.
|
||||
\end{equation*}
|
||||
Now, the Gibbs state approximation $\rho_{\omega}^{\text{Gibbs}}$ can be generated by propagating the trial state with VarQITE with respect to $H_{\text{eff}}$ for $\tau = 1/2\left(\text{k}_{\text{B}}\text{T}\right)$.
|
||||
The resulting state
|
||||
\begin{equation*}
|
||||
\ket{ \psi_{\omega}} = V\left(\omega\left(\tau\right)\right)\ket{0}^{\otimes 2n}
|
||||
\end{equation*}
|
||||
gives an approximation for the Gibbs state of interest
|
||||
\begin{equation*}
|
||||
\rho_{\omega}^{\text{Gibbs}} = \text{Tr}_{b}\left[\ket{ \psi\left({\omega\left(\tau\right)}\right)}\bra{ \psi\left({\omega\left(\tau\right)}\right)} \right] \approx \frac{e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)}}{Z}
|
||||
\end{equation*}
|
||||
by tracing out the ancillary system $b$.
|
||||
We would like to point out that the VarQITE propagation relates $\omega$ to $\theta$ via the energy derivative $C$ given in Eq.~\eqref{eq:AC}.
|
||||
|
||||
Equivalently to Eq.~\eqref{eq:sle}, Eq.~\eqref{eq:sle_derivative} is also prone to being ill-conditioned. Thus, the use of regularization schemes may be required.
|
||||
|
||||
Notably, this is an approximate state preparation scheme that relies on the representation capabilities of $\ket{ \psi_{\omega}}$. However, since the algorithm is employed in the context of machine learning we do not necessarily require perfect state preparation. The noise may even improve the training, as discussed e.g., in \cite{Noh2017RegularizingDN}.
|
||||
|
||||
McLachlan's variational principle is not only the key component for Gibbs state preparation. It also enables the QBM training with gradients of the actual loss function for generic Pauli terms in $H_{\theta}$, even if some of the qubits are hidden. Further details are given in Sec.~\ref{sec:implementation}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Variational QBM}
|
||||
\label{sec:implementation}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
In the following, \varqbm{} and the respective utilization of McLachlan's variational principle and VarQITE is discussed.
|
||||
We consider training data that takes at most $2^n$ different values and is distributed according to a discrete probability distribution $p^{\text{data}}$.
|
||||
The aim of a QBM is to train the parameters of $H_{\theta}$ such that the sampling probability distribution of the corresponding $\rho_{\omega}^{\text{Gibbs}}= e^{-H_{\theta}/\left(\text{k}_{\text{B}}\text{T}\right)} / Z$ for $\ket{v}, \: v\in{0, \ldots, 2^n-1}$ with
|
||||
\begin{equation*}
|
||||
p_v^{\text{QBM}} = \text{Tr}\left[\Lambda_v\rho_{\omega}^{\text{Gibbs}}\right],
|
||||
\end{equation*}
|
||||
approximates $p^{\text{data}}$.
|
||||
The QBM model is trained to represent $p^{\text{data}}$ by minimizing the loss, given in Eq.~\eqref{eq:loss_qbm}, with respect to the Hamiltonian parameters $\theta$, i.e.,
|
||||
\begin{equation*}
|
||||
\underset{\theta}{\min} \: L = \underset{\theta}{\min} \: \left(-\sum\limits_{v}p_v^{\text{data}}\log{p_v^{\text{QBM}}}\right).
|
||||
\end{equation*}
|
||||
Now, \varqbm{} facilitates gradient-based optimization with the derivative of the actual loss function
|
||||
\begin{equation}
|
||||
\label{eq:loss_derivative}
|
||||
\begin{split}
|
||||
\frac{\partial L}{\partial\theta_i} &= \frac{\partial\left( - \sum\limits_{v}p_v^{\text{data}}\log{p_v^{\text{QBM}}}\right)}{\partial\theta_i} \\ &= - \sum\limits_{v}p_v^{\text{data}}\frac{\partial p_v^{\text{QBM}}/\partial\theta_i}{p_v^{\text{QBM}}}
|
||||
\end{split}
|
||||
\end{equation}
|
||||
by using the chain rule, i.e., automatic differentiation.
|
||||
More precisely, the gradient of $L$ can be computed by using the chain rule for
|
||||
\begin{align}
|
||||
\label{eq:gradGibbsState}
|
||||
\begin{split}
|
||||
\frac{\partial p_v^{\text{QBM}}}{\partial\theta_i} &=
|
||||
\frac{\partial p_v^{\text{QBM}}}{\partial\omega\left(\tau\right)}\frac{\partial\omega\left(\tau\right)}{\partial\theta_i} \\
|
||||
&=\sum\limits_{k=0}^{q-1} \frac{\partial p_v^{\text{QBM}} }{\partial\omega_k\left(\tau\right)}\frac{\partial\omega_k\left(\tau\right)}{\partial\theta_i}.
|
||||
\end{split}
|
||||
\end{align}
|
||||
Firstly, $\partial p_v^{\text{QBM}}/\partial\omega_k\left(\tau\right)=\partial \text{Tr}\left[\Lambda_v\rho_{\omega}^{\text{Gibbs}}\right] / \partial\omega_k\left(\tau\right)$ can be evaluated with quantum gradient methods discussed in \cite{Farhi2018_gradients, Fujii2018_qcircuitLearn,killoran2018, SchuldQuantumGradients19, Zoufal2019} because the term has the following form $\partial \text{Tr}\left[\hat{O}\ket{\phi\left(\alpha\right)}\bra{\phi\left(\alpha\right)}\right] / \partial\alpha$.
|
||||
Secondly, ${\partial\omega_k\left(\tau\right)}/{\partial\theta_i}$ is evaluated by computing the derivative of Eq.~\eqref{eq:sle} with respect to the Hamiltonian parameters
|
||||
\begin{equation*}
|
||||
\frac{\partial A\dot{\omega}\left(\tau\right)}{\partial{\theta_i}} = \frac{\partial C}{\partial{\theta_i}}.
|
||||
\end{equation*}
|
||||
This gives the following system of linear equations
|
||||
\begin{align}
|
||||
\label {eq:sle_derivative}
|
||||
A\left(\frac{\partial\dot{\omega}\left(\tau\right)}{\partial{\theta_i}}\right) = \frac{\partial C}{\partial{\theta_i}} - \left(\frac{\partial A}{\partial{\theta_i}}\right)\dot{\omega}\left(\tau\right).
|
||||
\end{align}
|
||||
Now, solving for $\partial\dot{\omega}\left(\tau\right)/\partial{\theta_i}$ in every time step of the Gibbs state preparation enables the use of, e.g., an explicit Euler method to get
|
||||
\begin{equation}
|
||||
\label{eq:omega_derivative}
|
||||
\begin{split}
|
||||
\frac{\partial \omega_k\left(\tau \right)}{\partial_{\theta_i}} &= \frac{\partial\omega_k\left(\tau-\delta\tau\right)}{\partial_{\theta_i}}+\frac{\partial\dot{\omega}_k\left(\tau-\delta\tau\right)}{\partial_{\theta_i}}\delta\tau \\
|
||||
&=\frac{\partial \omega_k\left(0\right)}{\partial_{\theta_i}} + \sum\limits_{j=1}^{\tau/\delta\tau}\frac{\partial \dot{\omega}_k\left(j\delta\tau\right)}{\partial_{\theta_i}}\delta\tau.
|
||||
\end{split}
|
||||
\end{equation}
|
||||
We discuss the structure of the quantum circuits used to evaluate $\partial_{\theta_i}A$ and $\partial_{\theta_i}C$, in Appendix~\ref{app:a_c}.
|
||||
|
||||
In principle, the gradient of the loss function could also be approximated with a finite difference method. If the number of Hamiltonian parameters is smaller than the number of trial state parameters, this requires less evaluation circuits. However, given a trial state that has less parameters than the respective Hamiltonian, the automatic differentiation scheme presented in this section is favorable in terms of the number of evaluation circuits.
|
||||
A more detailed discussion on this topic can be found in Appendix~\ref{app:complexity}.
|
||||
|
||||
An outline of the Gibbs state preparation and evaluation of ${\partial \omega_k\left(\tau \right)}/{\partial_{\theta_i}}$ with VarQITE is presented in Algorithm \ref{algo:VarQITE}.
|
||||
|
||||
\begin{algorithm}
|
||||
\caption{VarQITE for \varqbm} \label{algo:VarQITE}
|
||||
\begin{algorithmic}[0]
|
||||
\State \textbf{input}
|
||||
\State $H_{\text{eff}} = H_{\theta}^a + I^b$
|
||||
\State$\tau = {1}/{2\left(\text{k}_{\text{B}}\text{T}\right)}$
|
||||
\State $\ket{\psi\left( \omega\left(0\right)\right)} = V\left(\omega\left(0\right)\right)\ket{0}^{\otimes 2n} = \ket{\phi^{+}}^{\otimes n}$
|
||||
\State with $\ket{\phi^{+}} = \left(\ket{00}+\ket{11}\right)/{\sqrt{2}}$
|
||||
\State \textbf{procedure}
|
||||
\For{$t\in\set{\delta\tau, 2\delta\tau, \ldots, \tau}$}
|
||||
\State Evaluate $A\left(t\right)$ and $C\left(t\right)$
|
||||
\State Solve $A\dot{\omega}\left(t\right) = C$
|
||||
\For{$i\in\set{0, \ldots, p-1}$}
|
||||
\State Evaluate $\partial_{\theta_i}C$ and $\partial_{\theta_i}A$
|
||||
\State \label{item:linEqDer}Solve $A\left(\partial_{\theta_i}\dot{\omega}\left(t\right)\right) = \partial_{\theta_i}C - \left(\partial_{\theta_i}A\right)\dot{\omega}\left(t\right)$
|
||||
\State \label{item:grad_omega_tau}Compute $ \: \partial_{\theta_i}{\omega}\left(t\right) = \partial_{\theta_i}{\omega}\left(t-\delta\tau\right)+\partial_{\theta_i}\dot{\omega}\left(t\right)\delta\tau$
|
||||
\EndFor
|
||||
\State \label{item:computeUpdateDer}Compute $\omega\left(t+\delta\tau\right) = \omega\left(t\right) + \dot{\omega}\left(t\right)\delta\tau $
|
||||
\EndFor
|
||||
\State \textbf{return} $\omega\left(\tau\right), \: \partial\omega\left(\tau\right)/\partial\theta$
|
||||
\end{algorithmic}
|
||||
\end{algorithm}
|
||||
|
||||
Now, using a classical optimizer, such as Truncated Newton \cite{TNCDembo1983} or Adam \cite{Kingmaadam14}, allows the parameters $\theta$ to be updated according to ${\partial L}/{\partial\theta}$ from Eq.~\eqref{eq:loss_derivative}.
|
||||
The \varqbm{} training is illustrated in Fig.~\ref{fig:varqbm}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Results} \label{sec:results}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
In this section, the Gibbs state preparation with VarQITE is demonstrated using numerical simulation as well as the quantum hardware provided by IBM Quantum \cite{ibmQX}.
|
||||
Furthermore, we present numerically simulated QBM training results for a generative and a discriminative learning task.
|
||||
First, aspects which are relevant for the practical implementation are discussed in Sec.~\ref{sec:impRel}.
|
||||
Next, experiments of quantum Gibbs state preparation with VarQITE are shown in Sec.~\ref{sec:VarQITEResults}.
|
||||
Then, we illustrate the training of a QBM with the goal to generate a state which exhibits the sampling behavior of a Bell state, see Sec.~\ref{subsec:generative}, and to classify fraudulent credit card transactions, Sec.~\ref{subsec:discriminative}.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Methods}
|
||||
\label{sec:impRel}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
To begin with, we discuss the choice of a suitable parameterized trial state consisting of $V\left(\omega\right)$ and $\ket{\psi_{\text{in}}}$.
|
||||
Most importantly, the initial state $\ket{\psi_{\text{in}}}$ must not be an eigenstate of $V\left(\omega\right)$ as this would imply that the circuit could only act trivially onto the state.
|
||||
Furthermore, the state needs to be able to represent a sufficiently accurate approximation of the target state. If we have to represent, e.g., a non-symmetric Hamiltonian, the chosen trial state needs to be able to generate non-symmetric states.
|
||||
Moreover, $V\left(\omega\right)$ should not exhibit too much symmetry as this may lead to a singular $A$ which in turn causes ill-conditioning of Eq.~\eqref{eq:sle}.
|
||||
Assume, e.g., that all entries of $C$ are zero and, thus, that Eq.~\eqref{eq:sle} is homogeneous. If $A$ is singular, infinitely many solutions exist and it is difficult for the algorithm to estimate which path to choose. If $A$ is non-singular, the solution is $\dot{\omega} = 0$ and the evolution stops although we might have only reached a local extreme point.
|
||||
Another possibility to cope with ill-conditioned systems of linear equations are least-squares methods in combination with regularization schemes. We test Tikhonov regularization \cite{Tikhonov:1620560} and Lasso regularization \cite{lassoTibshirani11} with an automatic parameter evaluation based on L-curve fitting \cite{Hansen00thel-curve}, as well as an $\epsilon$-perturbation of the diagonal, i.e., $A\rightarrow A+\epsilon I$. It turns out that all regularization methods perform similarly well.
|
||||
|
||||
The results discussed in this section employ Tikhonov regularization.
|
||||
Furthermore, we use trial states which are parameterized by Pauli-rotation gates. Therefore, the gradients of the QBM probabilities with respect to the trial state parameters
|
||||
\begin{equation*}
|
||||
\frac{ \partial p_v^{\text{QBM}}}{\partial\omega_k\left(\tau\right)}=\frac{\partial \text{Tr}\left[\Lambda_v\rho_{\omega}^{\text{Gibbs}}\right]} {\partial\omega_k\left(\tau\right)}
|
||||
\end{equation*}
|
||||
can be computed using a $\pi/2-$shift method which is, e.g., described in \cite{Zoufal2019}.
|
||||
All experiments employ an additional qubit $\ket{0}_{\text{add}}$ and parameter $\omega_{\text{add}}$ to circumvent a potential phase mismatch between the target $\ket{\psi_{\tau}}$ and the trained state $\ket{ \psi\left(\omega\left(\tau\right)\right)}$ \cite{VarSITEMcArdle19, Simon18TheoryVarQSim, QNGNonUnitary19Simon} by applying
|
||||
\begin{equation*}
|
||||
R_Z\left(\omega_{\text{add}}\right)\ket{0}_{\text{add}}.
|
||||
\end{equation*}
|
||||
Notably, the additional parameter increases the dimension of $A$ and $C$ by one.
|
||||
The effective temperature, which in principle acts as a scaling factor on the Hamiltonian parameters, is set to $\left(\text{k}_{\text{B}}\text{T}\right) = 1$ in all experiments.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Gibbs State Preparation with VarQITE}
|
||||
\label{sec:VarQITEResults}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=1\linewidth]{ansatz_gibbs.pdf}
|
||||
\end{center}
|
||||
\caption{The depicted circuits illustrate the initial trial state for the Gibbs state preparation of (a) $\rho^{\text{Gibbs}}_1 $ (b) $\rho^{\text{Gibbs}}_2$ using VarQITE.
|
||||
}
|
||||
\label{fig:ansaetze}
|
||||
\end{figure}
|
||||
|
||||
To demonstrate that VarQITE is able to generate suitable approximations to Gibbs states, we illustrate the convergence of the state fidelity with respect to the target state for the following two simple one- and two-qubit Hamiltonians
|
||||
\begin{align*}
|
||||
H_1 &= 1.0 Z, \\
|
||||
H_2 &= 1.0 ZZ - 0.2 ZI - 0.2 IZ + 0.3 XI + 0.3 IX.
|
||||
\end{align*}
|
||||
corresponding to
|
||||
\begin{align*}
|
||||
\rho^{\text{Gibbs}}_1 &= \left( \begin{array}{cc}
|
||||
0.12 & 0. \\
|
||||
0. & 0.88\\\end{array}\right)
|
||||
, \\
|
||||
\rho^{\text{Gibbs}}_2 &= \left( \begin{array}{cccc}
|
||||
0.10 & -0.06 & -0.06 & 0.01 \\
|
||||
-0.06 & 0.43 & 0.02 & -0.05 \\
|
||||
-0.06 & 0.02 & 0.43 & -0.05 \\
|
||||
0.01 & -0.05 & -0.05 & 0.05\\\end{array}\right).
|
||||
\end{align*}
|
||||
The results are computed using the parameterized quantum circuit shown in Fig.~\ref{fig:ansaetze}.
|
||||
|
||||
The algorithm is executed for $10$ time steps on different backends: an ideal simulator and the he \emph{ibmq$\_$johannesburg $20$-qubit} backend.
|
||||
Notably, readout error-mitigation \cite{qiskit, dewes2012readout, Stamatopoulos2019} is used to obtain the final results run on real quantum hardware. Fig.~\ref{fig:VarQITE} depicts the results considering the fidelity between the trained and the target Gibbs state for each time step.
|
||||
It should be noted that the fidelity for the quantum backend evaluations employ state tomography.
|
||||
The plots illustrate that the method approximates the states, we are interested in, reasonably well and that also the real quantum hardware achieves fidelity values over $0.99$ and $0.96$, respectively.
|
||||
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=0.8\linewidth]{VarQTEHW_stateTomo.pdf}
|
||||
\end{center}
|
||||
\caption{Fidelity between trained and target Gibbs state with VarQITE for (a) $\rho^{\text{Gibbs}}_1$ (b) $\rho^{\text{Gibbs}}_2$ trained with an ideal simulator and real quantum hardware, i.e., the \emph{ibmq$\_$johannesburg $20$-qubit} backend. Each simulation used $10$ time steps.
|
||||
}
|
||||
\label{fig:VarQITE}
|
||||
\end{figure}
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Generative Learning}
|
||||
\label{subsec:generative}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
Now, the results from an illustrative example of a generative QBM model are presented.
|
||||
More explicitly, the QBM is trained to mimic the sampling statistics of a Bell state $(\ket{00} + \ket{11})/\sqrt{2}$, which is a state that exhibits non-local correlations.
|
||||
Numerical simulations show that the distribution can be trained with a fully visible QBM which is based on the following Hamiltonian
|
||||
\begin{equation*}
|
||||
H_{\theta}= \theta_0 ZZ + \theta_1 IZ + \theta_2 ZI.
|
||||
\end{equation*}
|
||||
We draw the initial values of the Hamiltonian parameters $\theta$ from a uniform distribution on $\left[-1, 1\right]$.
|
||||
The optimization runs on an ideal simulation of a quantum computer using AMSGrad \cite{amsgrad} with initial learning rate $0.1$, maximum number of iterations $200$, first momentum $0.7$, and second momentum $0.99$ as optimization routine.
|
||||
The Gibbs state preparation uses the initial trial state shown in Fig.~\ref{fig:ansatz_Bell} and $10$ steps per state preparation.
|
||||
|
||||
The training is run $10$ times using different randomly drawn initial parameters.
|
||||
The averaged values of the loss function as well as the distance between the target distribution $p^{\text{data}}=\left[0.5, 0., 0., 0.5\right]$ and the trained distribution $p^{\text{QBM}}$ with respect to the $\ell_1$ norm are illustrated over $50$ optimization iterations in Fig.~\ref{fig:qrbm}.
|
||||
The plot shows that loss and distance converge toward the same values for all sets of initial parameters. Likewise, the trained parameters $\theta$ converge to similar values.
|
||||
Furthermore, Fig.~\ref{fig:bell_prob} illustrates the target probability distribution and for the best and worst of the trained distributions. The plot reveals that the model is able to train the respective distribution very well.
|
||||
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=1\linewidth]{ansatz_bell.png}
|
||||
\end{center}
|
||||
\caption{We train a QBM to mimic the sampling behavior of a Bell state.
|
||||
The underlying Gibbs state preparation with VarQITE uses the illustrated parameterized quantum circuit to prepare the initial trial state. The first two qubits represent the target system and the last two qubits are ancillas needed to generate the maximally-mixed state as starting state for the evolution.}
|
||||
\label{fig:ansatz_Bell}
|
||||
\end{figure}
|
||||
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=0.8\linewidth]{qrbm_plot.png}
|
||||
\end{center}
|
||||
\caption{The figure illustrates the training progress of a fully-visible QBM model which aims to represent the measurement distribution of a Bell state.
|
||||
The green function corresponds to the loss and the pink function represents the distance between the trained and target distribution with respect to the $\ell_1$ norm at each step of the iteration.
|
||||
Both measures are computed for $10$ different random seeds.
|
||||
The points represent the mean and the error bars the standard deviation of the results.
|
||||
}
|
||||
\label{fig:qrbm}
|
||||
\end{figure}
|
||||
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=0.8\linewidth]{bell_prob.png}
|
||||
\end{center}
|
||||
\caption{The figure illustrates the sampling probability of the Bell state (blue), as well as the best (pink) and worst (purple) probability distribution achieved from $10$ different random seeds.
|
||||
}
|
||||
\label{fig:bell_prob}
|
||||
\end{figure}
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\subsection{Discriminative Learning}
|
||||
\label{subsec:discriminative}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
QBMs are not only applicable for generative but also for discriminative learning. We discuss the application to a classification task, the identification of fraudulent credit card transactions.
|
||||
|
||||
To enable discriminative learning with QBMs, we use the input data points $x$ as bias for the Hamiltonian weights. More explicitly, the parameters of the Hamiltonian
|
||||
\begin{equation}
|
||||
H_{\theta}\left(x\right)=\sum\limits_{i}f_i\left(\theta, x\right)h_i
|
||||
\end{equation}
|
||||
are given by a function $f_i\left(\theta, x\right)$ which maps $\theta$ and $x$ to a scalar in $\mathbb{R}$.
|
||||
Now, the respective loss function reads
|
||||
\begin{equation*}
|
||||
\label{eq:loss_sup}
|
||||
\begin{split}
|
||||
L = -\sum\limits_{x}p_x^{\text{data}}\sum\limits_{v}p_{v|x}^{\text{data}}\log{p_{v|x}^{\text{QBM}}}
|
||||
\end{split}
|
||||
\end{equation*}
|
||||
with
|
||||
\begin{equation*}
|
||||
p_{v|x}^{\text{QBM}} = \text{Tr}\left[\Lambda_v\rho\left(x\right)_{\omega}^{\text{Gibbs}}\right],
|
||||
\end{equation*}
|
||||
where $\rho\left(x\right)_{\omega}^{\text{Gibbs}}$ denotes the approximate Gibbs state corresponding to $H_{\theta}\left(x\right)$.
|
||||
The model encodes the class labels in the measured output configuration of the visible qubits $v$ of $\rho\left(x\right)_{\omega}^{\text{Gibbs}}$.
|
||||
Now, the aim of the training is to find Hamiltonian parameters $\theta$ such that, given a data sample $x$, the probability of sampling the correct output label from $\rho\left(x\right)_{\omega}^{\text{Gibbs}}$ is maximized.
|
||||
|
||||
The training is based on $500$ artificially created credit card transactions \cite{Altman2019} with about $15\%$ fraudulent instances. To avoid redundant state preparation, the training is run for all unique item instances in the data set and the results are averaged according to the item's occurrence counts.
|
||||
The dataset includes the following features: location (ZIP code), time, amount, and Merchant Category Code (MCC) of the transactions.
|
||||
To facilitate the training, the features of the given data set are discretized and normalized as follows. Using k-means clustering, each of the first three features are independently discretized to $3$ reasonable bins. Furthermore, we consider MCCs $<10 000$ and group them into $10$ different categories.
|
||||
The discretization is discussed in more detail in Table \ref{tbl:discr_data_preproc}.
|
||||
Furthermore, for each feature, we map the values $x$ to $x' = \frac{x-\mu}{\sigma}$ with $\mu$ denoting the mean and $\sigma$ denoting the standard deviation.
|
||||
\begin{table}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
%\begin{tabular}{|l|l|l|}
|
||||
\begin{tabular}{c|c|c}
|
||||
Feature & Condition & Value \\
|
||||
\hline
|
||||
\multirow{2}{*}{Time} & $0$AM $- 11$AM & 0 \\
|
||||
& $11$AM $- 6$PM & 1 \\
|
||||
& $6$PM - $0$AM & 2 \\
|
||||
\hline
|
||||
\multirow{2}{*}{Amount} & amount < $\$50 $ & 0 \\
|
||||
& amount in $\$50-150 $ & 1 \\
|
||||
& amount > $\$150$ & 2 \\
|
||||
\hline
|
||||
\multirow{3}{*}{ZIP} & east & 0 \\
|
||||
& central & 1 \\
|
||||
& west & 2 \\
|
||||
\end{tabular}
|
||||
\caption{The table discusses the clustering of a transaction fraud data set which is used to train a discriminative QBM model. MCC refers to the merchant category code and ZIP to zone improvement plan.
|
||||
Notably, the given values are approximate.}
|
||||
\label{tbl:discr_data_preproc}
|
||||
\end{table}
|
||||
|
||||
The complexity of this model demands a Hamiltonian that has sufficient representation capabilities. Our choice is the following
|
||||
\begin{equation}
|
||||
\label{eq:H_disc}
|
||||
\begin{split}
|
||||
H_{\theta}\left(x\right) =&\:f_0\left(\theta, x\right) ZZ + f_1\left(\theta, x\right)ZI +\\
|
||||
&\:f_2\left(\theta, x\right) IZ + f_3\left(\theta, x\right)XI + f_4\left(\theta, x\right) IX,
|
||||
\end{split}
|
||||
\end{equation}
|
||||
where $f_i\left(\theta, x\right) = \vec{\theta}_i\cdot\vec{x}$ corresponds to the dot product of the vector corresponding to the data item $\vec{x}$ and a parameter vector $\vec{\theta}_i$ of equal length. Additionally, the first and second qubit correspond to a hidden and visible qubit, respectively.
|
||||
|
||||
Since the numerical simulation of variational Gibbs state preparation for various $H_{\theta}\left( x\right)$, with $x$ corresponding to all unique data items, is computationally expensive, we decided to train the parameters $\theta$ using an exact representation of the quantum Gibbs states.
|
||||
The resulting $\theta$ are then used for Gibbs state preparation with VarQITE.
|
||||
Even though the parameters are not trained with variational Gibbs state preparation, the results discussed in this section demonstrate that we can find suitable parameters $\theta$ such that \varqbm{} corresponds to a well-performing discriminative model.
|
||||
|
||||
The exact training uses a Truncated Netwon optimization routine \cite{TNCDembo1983} with a maximum iteration number of $100$ and the step size for the numerical approximation of the Jacobian being set to $10^{-6}$.
|
||||
The initial values for the Hamiltonian parameters are drawn from a uniform distribution on $\left[-1, 1\right]$.
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=1\linewidth]{ansatz_disc.png}
|
||||
\end{center}
|
||||
\caption{Given a transaction instance, the measurement output of the QBM labels it as being either fraudulent or valid. The underlying Gibbs state preparation with VarQITE uses the illustrated parameterized quantum circuit as initial trial state. The first qubit is the visible node that determines the QBM output, the second qubit represents the hidden unit, and the last two qubits are ancillas needed to generate the maximally-mixed state as starting state for the evolution.}
|
||||
\label{fig:ansatz_Disc}
|
||||
\end{figure}
|
||||
|
||||
Given a test data set consisting of $250$ instances, with about $10\%$ fraudulent transactions, the Gibbs states, corresponding to the unique items of the test data, are approximated using VarQITE with the trained parameters $\theta$ and the trial state shown in Fig.~\ref{fig:ansatz_Disc}.
|
||||
To predict the labels of the data instances, we sample from the states $\rho_{\omega}^{\text{Gibbs}}$ and choose the label with the highest sampling probability.
|
||||
These results are, then, used to evaluate the accuracy, precision, recall and F$_1$ score.
|
||||
It should be noted that we choose a relatively simple quantum circuit to keep the simulation cost small.
|
||||
However, it can be expected that a more complex parameterized quantum circuit would lead to further improvement in the training results.
|
||||
|
||||
The resulting values are compared to a set of standard classifiers defined in a \emph{scikit-learn} \cite{scikit-learn2011} classifier comparison tutorial \cite{scikitClassifierComp}, see Tbl.~\ref{tbl:measuresQBM}. The respective classifiers are used with the hyper parameters defined in this tutorial.
|
||||
Notably, the Linear SVM does not classify any test data item as fraudulent and, thus, the classifier sets precision and recall score to $0$.
|
||||
The comparison reveals that the QBM performs similarly well to the classical classifiers considering accuracy, is competitive regarding precision, and even outperforms them in terms of recall. The best F$_1$ score is achieved with \varqbm.
|
||||
|
||||
\begin{table}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
{\renewcommand{\arraystretch}{1.2}
|
||||
\begin{tabular}{c | c | c | c | c}
|
||||
Model & Accuracy & Recall & Precision & F$_1$\\
|
||||
\hline
|
||||
Nearest Neighbours & $0.94$ & $0.54$ & $0.72$ & $0.31$ \\
|
||||
Linear SVM & $0.90 $& $0$ & $0$ & $0$ \\
|
||||
RBF SVM & $0.94 $& $0.42$ & $0.83$ & $0.28$ \\
|
||||
Gaussian Process & $0.94$ & $0.46$ & $0.85$ & $0.30$\\
|
||||
Gaussian Naive Bayes & $0.91$ &$ 0.42 $& $0.56$ & $0.24$\\
|
||||
Decision Tree & $0.94$ & $0.42$ & $0.83$ & $0.28$\\
|
||||
Random Forrest & $0.93 $& $0.29$ & $1.00$ & $0.22$\\
|
||||
Multi-layer Perceptron & $0.94$ & $0.38$ &$0.9$ & $0.27$\\
|
||||
AdaBoost & $0.94$ & $0.54$ & $0.81$ & $0.32$\\
|
||||
QDA & $0.92$ & $0.46$ & $0.61$ & $0.26$\\
|
||||
\textbf{\varqbm} & $\mathbf{0.95}$ &
|
||||
$\mathbf{0.63}$ & $\mathbf{0.83}$ & $\mathbf{0.36}$\\
|
||||
\end{tabular}
|
||||
}
|
||||
\caption{This table presents performance measures for scikit-learn standard classifiers, as well as the trained QBM. The Nearest Neighbours classifier uses a $3$ nearest neighbours vote. The Linear and RBF Support Vector Machine (SVM) are based on a linear and radial kernel, respectively. The Linear SVM uses a regularization term of $0.25$ and for the RBF SVM the kernel coefficient is set to $2$. The maximum depth of the Decision Tree as well as the Random Forrest is set to 5. Furthermore, the Random Forrest classifier uses $10$ trees and uses $1$ feature to search for the best spit. The Multi-layer Perceptron uses $\ell_2$ regularization with coefficient $1$ and a maximum iteration number of $1000$. QDA refers to Quadratic Discriminant Analysis. It should be noted that the remaining classifier properties are default settings.}
|
||||
\label{tbl:measuresQBM}
|
||||
\end{table}
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%}
|
||||
\section{Conclusion and Outlook}
|
||||
\label{sec:discussion}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
This work presents the application of McLachlan's variational principle to facilitate \varqbm, a variational QBM algorithm, that is compatible with generic Hamiltonians and can be trained using analytic gradients of the actual loss function even if some of the qubits are hidden. Suppose a sufficiently powerful variational trial state, the presented scheme is not only compatible with local but also long-range correlations and for arbitrary system temperatures.
|
||||
|
||||
We outline the practical steps for utilizing VarQITE for Gibbs state preparation and verify that it can train states which are reasonably close to the target using simulation as well as real quantum hardware.
|
||||
Moreover, applications to generative learning and classification are discussed and illustrated with further numerical results.
|
||||
The presented model offers a versatile framework which facilitates the representation of complex structures with quantum circuits.
|
||||
|
||||
An interesting question for future research is the investigation of performance measures that improve our understanding of the model's representation capabilities.
|
||||
Furthermore, QBMs are not limited to the presented applications. They could also be utilized to train models for data from experiments with quantum systems. This is a problem that has recently gained interest, see e.g., \cite{LearningModels20}. Additionally, they might be employed for combinatorial optimization. Classical BMs have been investigated in this context \cite{Spieksma1995} and developing and analyzing quantum algorithms for combinatorial optimization is an active area of research \cite{FarhiQAOA14, Barkoutsos20VQOCVaR}.
|
||||
|
||||
All in all, there are many possible applications which still have to be explored.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Acknowledgments}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
We would like to thank Erik Altman for making the synthetic credit card transaction dataset available to us.
|
||||
Moreover, we are grateful to Pauline Ollitrault, Guglielmo Mazzola and Mario Motta for sharing their knowledge and engaging in helpful discussion.
|
||||
Furthermore, we thank Julien Gacon for his help with the implementation of the algorithm and all of the IBM Quantum team for its constant support.
|
||||
|
||||
Also, we acknowledge the support of the National Centre of Competence in Research \textit{Quantum Science and Technology} (QSIT).
|
||||
|
||||
IBM, the IBM logo, and ibm.com are trademarks of International Business Machines Corp., registered in many jurisdictions worldwide. Other product and service names might be trademarks of IBM or other companies. The current list of IBM trademarks is available at \url{https://www.ibm.com/legal/copytrade}.
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%\section{Author Contributions}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%All authors researched, collated, and wrote this paper.
|
||||
%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%\section{Competing Interests}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%The authors declare that there are no competing interests.
|
||||
%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%\section{Data Availability}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%The data that support the findings of this study are available from the corresponding author upon reasonable request.
|
||||
|
||||
|
||||
\appendix
|
||||
|
||||
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Evaluation of A, C and their gradients}
|
||||
\label{app:a_c}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
The elements of the matrix $A$ and the vector $C$, see Eq.~\eqref{eq:AC} are of the following form
|
||||
\begin{equation}
|
||||
\label{eq:expValueEval}
|
||||
\text{Re}\left(e^{i\alpha}\text{Tr}\left[U^{\dagger}V\rho_{\text{in}}\right]\right)
|
||||
\end{equation}
|
||||
with $\text{Re}\left(\cdot\right)$ denoting the real part and $\rho_{\text{in}} = \ket{\psi_{\text{in}}}\bra{\psi_{\text{in}}}$.
|
||||
As discussed in \cite{VarSITEMcArdle19, Simon18TheoryVarQSim, LaflammeSimulatingPhysPhenom02}, such terms can be computed by sampling the expectation value of an observable $Z$ with respect to the quantum circuit shown in Fig.~\ref{fig:expValueCircuit}.
|
||||
|
||||
\begin{figure}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
\begin{center}
|
||||
\includegraphics[width=1\linewidth]{eval_circ_fig.pdf}
|
||||
\end{center}
|
||||
\caption{Quantum circuit to evaluate $\text{Re}\left(e^{i\alpha}\text{Tr}\left[U^{\dagger}V\rho_{\text{in}}\right]\right)$ with $\rho_{\text{in}} = \ket{\psi}\bra{\psi_{\text{in}}}$.}
|
||||
\label{fig:expValueCircuit}
|
||||
\end{figure}
|
||||
|
||||
Notably, the phase $e^{i\alpha}$ in the first qubit is needed to include phases which may occur from gate derivatives. In our case, $\alpha$ needs to be set to $0$ respectively $\pi/2$ when computing the terms of $A$ or $C$. More precisely, the first qubit is initialized by an $H$ gate for $A$ and $H$ followed by an $S$ gate for $C$.
|
||||
These phases come from the fact that the trial states, used in this work, are constructed via Pauli rotations, i.e., $U\left(\omega\right) = R_{\sigma_l}\left(\omega\right)$ with $\sigma_l \in \set{X, Y, Z}$, which leads to
|
||||
\begin{equation}
|
||||
\frac{\partial U\left(\omega\right)}{\partial\omega} = -\frac{i}{2}{\sigma_l}R_{\sigma_l}\left(\omega\right).
|
||||
\end{equation}
|
||||
Furthermore, this method can be applied for the evaluation of $\partial A / \partial \theta$ and $\partial C / \partial \theta$, i.e., the respective terms can be written in the form of Eq.~\eqref{eq:expValueEval}.
|
||||
More precisely,
|
||||
\begin{equation*}
|
||||
\label{eq:dH_A}
|
||||
\begin{split}
|
||||
&\partial_{\theta_i}A_{p,q}\left(\tau\right) =\\
|
||||
&\sum\limits_s \frac{\partial\omega_s\left(\tau\right)}{\partial_{\theta_i}}\text{Re}\left(\text{Tr}\left[\left(\frac{\partial^2 V^{\dagger}\left({\omega\left(\tau\right)}\right)}{\partial\omega_p\left(\tau\right) \partial\omega_s\left(\tau\right)}\frac{\partial V\left({\omega\left(\tau\right)}\right)}{\partial\omega\left(\tau\right)_q}\right. \right.\right.\\
|
||||
&\left.\left.\left. + \frac{\partial V^{\dagger}\left({\omega\left(\tau\right)}\right)}{\partial\omega\left(\tau\right)_p}\frac{\partial^2 V\left({\omega\left(\tau\right)}\right)}{\partial\omega_q\left(\tau\right)\partial\omega_s\left(\tau\right)}\right)\rho_{\text{in}}\right]\right)
|
||||
\end{split}
|
||||
\end{equation*}
|
||||
and
|
||||
\begin{equation*}
|
||||
\label{eq:dH_C}
|
||||
\begin{split}
|
||||
&\partial_{\theta_j}C_p = \\
|
||||
&-\text{Re}\left(\text{Tr}\left[\frac{\partial V^{\dagger}\left({\omega\left(\tau\right)}\right)}{\partial\omega\left(\tau\right)_p}h_jV\left({\omega\left(\tau\right)}\right)\rho_{\text{in}}\right]\right)\\
|
||||
&\left. - \sum\limits_{i,s}\theta_i\frac{\partial\omega_s\left(\tau\right)}{\partial_{\theta_j}}\text{Re}\left( \text{Tr}\left[\left(\frac{\partial V^{\dagger}\left({\omega\left(\tau\right)}\right)}{\partial\omega_p\left(\tau\right)}h_i\frac{\partial V\left({\omega\left(\tau\right)}\right)}{\partial\omega_s\left(\tau\right)} \right.\right. \right.\right.\\
|
||||
&+ \left.\left. \left.\frac{\partial^2 V^{\dagger}\left({\omega\left(\tau\right)}\right)}{\partial\omega_p\left(\tau\right)\partial\omega_s\left(\tau\right)} h_iV\left({\omega\left(\tau\right)}\right)\right)\rho_{\text{in}} \right] \right).
|
||||
\end{split}
|
||||
\end{equation*}
|
||||
Hereby, $\alpha$ must be set to $\pi/2$ respectively $0$ for the terms in $\partial A/ \partial\theta$ respectively $\partial C/ \partial\theta$. This is achieved with the same gates as mentioned before.
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Complexity Analysis}
|
||||
\label{app:complexity}
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
To compute the gradient $ \partial L / \partial \theta$ of the loss function, given in Eq.~\eqref{eq:loss_qbm}, we could use either a numerical finite differences method \cite{Kardestuncer1975}, or the analytic, automatic differentiation approach that is presented in this paper. In the following, we discuss the number of circuits that have to be evaluated for those gradient implementations for a trial state with $q$ parameters, an $n$-qubit Hamiltonian with $p$ parameters, and VarQITE for Gibbs state preparation using $t$ steps.
|
||||
|
||||
The number of circuits that need to be evaluated for Gibbs state preparation with VarQITE are $\Theta\left(tq^2 \right)$ and $\Theta\left(tqp \right)$ for $A$ and $C$, respectively.
|
||||
Therefore, the overall number of circuits is $\Theta\left(tq(q+p)\right)$.
|
||||
Now, computing the gradient with forward finite differences reads
|
||||
\begin{equation*}
|
||||
\frac{\partial L}{\partial \theta} \approx \frac{L\left( \theta + \epsilon \right) - L\left( \theta \right)}{\epsilon},
|
||||
\end{equation*}
|
||||
for $0 < \epsilon \ll 1$.
|
||||
For this purpose, VarQITE must be run once with $\theta$ and $p$ times with an $\epsilon$-shift which leads to a total number of $\Theta\left(tpq(q+p)\right)$ circuits.
|
||||
|
||||
The automatic differentiation gradient, given in Eq.~\eqref{eq:loss_derivative}, corresponds to
|
||||
\begin{equation*}
|
||||
\begin{split}
|
||||
\frac{\partial L}{\partial\theta} = - \sum\limits_{v}\sum\limits_{k=0}^{q-1}\frac{p_v^{\text{data}}}{\left\langle \Lambda_v\right\rangle} \frac{\partial \left\langle \Lambda_v \right\rangle}{\partial \omega_k} \frac{\partial \omega_k} {\partial \theta}
|
||||
\end{split}
|
||||
\end{equation*}
|
||||
with $\langle \ldots \rangle = \text{Tr}\left[\rho_{\omega}^{Gibbs}\ldots\right]$.
|
||||
VarQITE needs to be run once to prepare $\rho_{\omega}^{\text{Gibbs}}$. Furthermore, the evaluation of $\partial\omega_k / \partial\theta$ requires that $\partial A / \partial\theta$ and $\partial C / \partial\theta$ are computed for every step of the Gibbs state preparation. This leads to $\Theta\left(tq^2(q+p)\right)$ circuits.
|
||||
The resulting overall complexity of the number of circuits is $\Theta\left(tq^2(q+p)\right)$.
|
||||
|
||||
The results are summarized in Tbl.~\ref{tbl:complexity}.
|
||||
Automatic differentiation is more efficient than finite differences if $q < p$. For $q > p$, on the other hand, focusing mainly on computational complexity, one should rather use finite differences.
|
||||
Considering, e.g., a $k$-local Ising model that corresponds to a Hamiltonian with $\mathcal{O}\left(n^k\right)$ parameters. Suppose that we can find a reasonable variational $n$-qubit trial state with $\mathcal{O}\left(n\right)$ layers of parameterized and entangling gates, which results in $q = \mathcal{O}\left(n^2\right)$ parameters, then, automatic differentiation would outperform finite differences for $k>2$.
|
||||
|
||||
\begin{table}[h!]
|
||||
\captionsetup{singlelinecheck = false, format= hang, justification=raggedright, font=footnotesize, labelsep=space}
|
||||
{\renewcommand{\arraystretch}{1.2}
|
||||
\begin{tabular}{ c | c }
|
||||
Method & Number Circuits \\
|
||||
\hline
|
||||
Finite Diff & $\Theta\left(tqp(q+p)\right)$\\
|
||||
Automatic Diff & $\Theta\left(tq^2(q+p)\right)$\\
|
||||
\end{tabular}
|
||||
}
|
||||
\caption{Comparing the number of circuits needed to train a QBM with VarQITE using either finite differences or automatic differentiation. The number of Hamiltonian parameters is $p$, the number of trial state parameters is $q$ and the number of time steps during the Gibbs state preparation is $t$.}
|
||||
\label{tbl:complexity}
|
||||
\end{table}
|
||||
|
||||
\bibliographystyle{IEEEtranN}
|
||||
\bibliography{references}
|
||||
|
||||
|
||||
|
||||
\end{document}
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,694 @@
|
||||
\section{Theory and methods}
|
||||
|
||||
\subsection{Test case}
|
||||
The function to be fitted as verification of the regression methods is the famous Franke's function \(f:\qty[0,1]^2\to\mathbb{R}\) given by
|
||||
\begin{equation}
|
||||
\begin{alignedat}{2}
|
||||
f(x_1,x_2) &= \frac{3}{4}\oldexp(-\frac{(9x_1-2)^2}{4} - \frac{(9x_2-2)^2}{4})
|
||||
+ \frac{3}{4}\oldexp(-\frac{(9x_1+1)^2}{49}- \frac{(9x_2+1)}{10} ) \\
|
||||
&+ \frac{1}{2}\oldexp(-\frac{(9x_1-7)^2}{4} - \frac{(9x_2-3)^2}{4})
|
||||
- \frac{1}{5}\oldexp(-(9x_1-4)^2 - (9x_2-7)^2)
|
||||
\end{alignedat}\label{eq:franke}
|
||||
\end{equation}
|
||||
and shown in~\vref{fig:franke}.
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\includegraphics{franke.pdf}
|
||||
\caption{The Franke's function to be fitted, given by~\vref{eq:franke}.}\label{fig:franke}
|
||||
\end{figure}
|
||||
|
||||
\subsection{Fitting problem statement}
|
||||
The general problem is to fit a given set of data \(D=\qty{(\vec{x}_i,y_i)}_{i=1}^N\), where \(\vec{x}_i\) are inputs of some dimensionality (\(2\) in this project) and \(y_i\) are scalar outputs.
|
||||
A set of basis functions \(\qty{\phi_j}_{j=1}^p\) is chosen, and the goal is to approximate the input data with the linear combination \(\beta_j \phi_j(\vec{x})\).
|
||||
If the data set were generated by such a linear combination, there would exist parameters \(\qty{\beta_j}\) such that \(\beta_j\phi_j(\vec{x}_i)=y_i\), which can be rewritten as the matrix equation
|
||||
\begin{equation}
|
||||
X\vec{\beta} = \vec{y},
|
||||
\end{equation}
|
||||
where \(X\in\mathbb{R}^{N\times p}\) is the matrix with elements \(x_{ij}=\phi_j(\vec{x}_i)\), and \(\vec{\beta}\) and \(\vec{y}\) are the vectors with elements \(\beta_j\) and \(y_i\).
|
||||
|
||||
In general, it will not be possible to find parameters \(\beta_j\) which fulfill this, since the data set will not be generated from the simple basis functions.
|
||||
Consequently, one must find the \(\beta_j\)s which in some sense make \(X\vec{\beta}\) as close to \(\vec{y}\) as possible.
|
||||
The way in which to measure deviation from the perfect solution is called the \emph{cost function}, frequently written \(Q(\vec{\beta};D)\).
|
||||
Regression methods differ in the choice of cost function, while their aim is always to find the parameters \(\beta_j\) which minimise the cost function.
|
||||
|
||||
When a regression method has been used to find \(\vec{\beta}\), predicted values are denoted \(\tilde{y}_i = X_{ij}\beta_j\), while the original, exact values are denoted \(y_i\).
|
||||
|
||||
% _ _ _ _
|
||||
% ___ _ __ __| (_)_ __ __ _ _ __ _ _ | | ___ __ _ ___| |_
|
||||
% / _ \| '__/ _` | | '_ \ / _` | '__| | | | | |/ _ \/ _` / __| __|
|
||||
% | (_) | | | (_| | | | | | (_| | | | |_| | | | __/ (_| \__ \ |_
|
||||
% \___/|_| \__,_|_|_| |_|\__,_|_| \__, | |_|\___|\__,_|___/\__|
|
||||
% ___ __ _ _ _ __ _ _ __ ___ ___|___/
|
||||
% / __|/ _` | | | |/ _` | '__/ _ \/ __|
|
||||
% \__ \ (_| | |_| | (_| | | | __/\__ \
|
||||
% |___/\__, |\__,_|\__,_|_| \___||___/
|
||||
% |_|
|
||||
\subsection{Ordinary Least Squares}
|
||||
The ordinary least squares method is the simplest and most intuitive, since its cost function is simply the square of the error made by the fit,
|
||||
\begin{equation}
|
||||
Q(\beta;D) = \norm{\vec{y}-X\vec{\beta}}_2^2,
|
||||
\end{equation}
|
||||
where \(\norm{\cdot}_p\) denotes the usual \(p\)-norm.
|
||||
|
||||
\subsubsection{Geometric view}
|
||||
The geometric view of the equation \(X\vec{\beta}=\vec{y}\) is to find the linear combination of the columns of \(X\) equal to \(\vec{y}\).
|
||||
This is not necessarily possible if the columns of \(X\) do not span \(\mathbb{R}^N\), which is not possible if \(p<N\).
|
||||
Ordinary least squares seeks to find the linear combination of the columns of \(X\) as close to \(\vec{y}\) as possible.
|
||||
This is achieved when \(X\vec{\beta}\) is equal to the projection of \(\vec{y}\) onto the column space of \(X\), i.e.
|
||||
\begin{equation}
|
||||
X\vec{\beta} = \Proj_{\Col{X}}{\vec{y}}.
|
||||
\end{equation}
|
||||
By construction, this equation will always have a solution, although it may not be unique if the columns of \(X\) are not linearly independent.
|
||||
Since \(X\vec{\beta}\) is the projection of \(\vec{y}\) onto the column space of \(X\), the error, \(\vec{y}-X\vec{\beta}\), is orthogonal to the columns of \(X\), i.e.\ the rows of \(X^T\).
|
||||
Consequently,
|
||||
\begin{equation}
|
||||
X^T\qty(\vec{y}-X\vec{\beta}) = \vec{0}
|
||||
\implies
|
||||
X^T X \vec{\beta} = X^T \vec{y},\label{eq:olsnormal}
|
||||
\end{equation}
|
||||
which is a simple linear set of equations which can be solved for \(\vec{\beta}\).
|
||||
This set of equations is called the \emph{normal equations}.
|
||||
If \(X\) is non-singular, the minimisation problem in ordinary least squares has the closed form solution
|
||||
\begin{equation}
|
||||
\vec{\beta}_{\text{OLS}} = \qty(X^T X)^{-1} X^T \vec{y}.
|
||||
\end{equation}
|
||||
|
||||
\subsubsection{Minimisation view}
|
||||
The cost function for ordinary least squares can be written as
|
||||
\begin{equation}
|
||||
Q(\beta;D) = \norm{\vec{y}-X\vec{\beta}}_2^2
|
||||
= \sum_{i=1}^N \qty(y_i - X_{ij}\beta_j)^2.
|
||||
\end{equation}
|
||||
When this is minimised, \(\nabla Q = \vec{0}\), where the gradient denotes differentiation with respect to the parameters \(\beta_j\).
|
||||
The components of the gradient can straightforwardly be calculated from the cost function,
|
||||
\begin{equation}
|
||||
\nabla_k Q = \pdv{\beta_k}(\qty(y_i - X_{ij}\beta_j)\qty(y_i - X_{ij}\beta_j))
|
||||
= 2\qty(y_i-X_{ij}\beta_j) \pdv{\beta_k}(y_i-X_{ij}\beta_j)
|
||||
= -2X_{ik}\qty(y_i-X_{ij}\beta_j),\label{eq:olsgradcomp}
|
||||
\end{equation}
|
||||
which can be rewritten on vector form as
|
||||
\begin{equation}
|
||||
\nabla Q = -2X^T\qty(\vec{y}-X\vec{\beta}).\label{eq:olsgrad}
|
||||
\end{equation}
|
||||
Setting \(\nabla Q = \vec{0}\) gives~\vref{eq:olsnormal}.
|
||||
|
||||
Ordinary least squares has been implemented by calling \lstinline{dgelss}, which uses a singular value decomposition to prevent numerical instabilities.
|
||||
|
||||
% _ _
|
||||
% _ __(_) __| | __ _ ___
|
||||
%| '__| |/ _` |/ _` |/ _ \
|
||||
%| | | | (_| | (_| | __/
|
||||
%|_| |_|\__,_|\__, |\___|
|
||||
% |___/
|
||||
\subsection{Ridge regression}
|
||||
The cost function used in Ridge regression is
|
||||
\begin{equation}
|
||||
Q_\lambda(\vec{\beta};D) = \norm{\vec{y}-X\vec{\beta}}_2^2 + \lambda \norm{\vec{\beta}}_2^2
|
||||
= \sum_{i=1}^N \qty(y_i - X_{ij}\beta_j)^2 + \lambda\sum_{i=1}^p \beta_i^2,
|
||||
\end{equation}
|
||||
which penalises solution vectors \(\vec{\beta}\) where some coefficients are large.
|
||||
\(\lambda\) is a parameter which should be chosen with great care.
|
||||
The error is as small as possible for \(\lambda=0\), as this will reproduce the solution found by ordinary least squares, but a non-zero value of \(\lambda\) will give \(\beta_j\)s which yield more reasonable predictions for other values of \(\vec{x}\) than those contained in the training data set \(D\)\cite{mehta}.
|
||||
|
||||
Using the derivative of the first term from~\vref{eq:olsgrad}, the gradient of the cost function is
|
||||
\begin{equation}
|
||||
\nabla Q_\lambda = -2 X^T \qty(\vec{y}-X\vec{\beta}) + 2\lambda \vec{\beta}
|
||||
= -2 X^T \vec{y} + 2 X^T X \vec{\beta} + 2\lambda \vec{\beta}
|
||||
= -2 X^T \vec{y} + 2\qty(X^T X + \lambda I)\vec{\beta}.
|
||||
\end{equation}
|
||||
Minimisation requires \(\nabla Q_\lambda = \vec{0}\), which gives
|
||||
\begin{equation}
|
||||
\qty(X^T X + \lambda I)\vec{\beta} = X^T \vec{y}.
|
||||
\end{equation}
|
||||
This can be solved for \(\vec{\beta}\), and the closed-form solution to the minimisation of the Ridge cost function is
|
||||
\begin{equation}
|
||||
\vec{\beta}_\text{Ridge} = \qty(X^T X + \lambda I)^{-1} X^T \vec{y}.
|
||||
\end{equation}
|
||||
Having a non-zero \(\lambda\) clearly reduces problems with linearly dependent columns of \(X\) and singularity of \(X^T X\).
|
||||
|
||||
Ridge regression has been implemented by calculating \(X^T X\) and adding \(\lambda\) on the diagonal, calculating \(X^T y\) and using \lstinline{dposv} to find \(\vec{\beta}\) using a Cholesky-decomposition, since \(X^T X + \lambda I\) is positive definite.
|
||||
|
||||
% _
|
||||
%| | __ _ ___ ___ ___
|
||||
%| |/ _` / __/ __|/ _ \
|
||||
%| | (_| \__ \__ \ (_) |
|
||||
%|_|\__,_|___/___/\___/
|
||||
\tikzexternaldisable
|
||||
\subsection{LASSO regression}
|
||||
The cost function used in LASSO regression is
|
||||
\begin{equation}
|
||||
Q_\lambda(\vec{\beta};D) = \norm{\vec{y}-X\vec{\beta}}_2^2 + \lambda \norm{\vec{\beta}}_1
|
||||
= \sum_{i=1}^N \qty(y_i - X_{ij}\beta_j)^2 + \lambda\sum_{i=1}^p \abs{\beta_i},
|
||||
\end{equation}
|
||||
which also penalises solution vectors \(\vec{\beta}\) where some coefficients are large.
|
||||
LASSO regression has the advantage that certain choices of \(\lambda\) will give a sparse solution, i.e. \(\beta_j=0\) for some values of \(j\)\cite{wieringen}.
|
||||
|
||||
As with the other regression methods, the gradient of the cost function should now be differentiated and set to zero.
|
||||
Mathematicians will now point out that the absolute value is not differentiable at zero and start deriving ingenious minimisation methods to circumvent the problem.
|
||||
This is usually a non-issue in a physics course --- I will now \emph{define} the derivative to be
|
||||
\begin{equation}
|
||||
\dv{\abs{x}}{x} \coloneqq \sgn{x} \coloneqq \begin{cases}1, & x > 0\\ 0, & x = 0 \\ -1, & x < 0 \end{cases}
|
||||
\end{equation}
|
||||
and proceed happily.
|
||||
|
||||
The gradient of the cost function is thus
|
||||
\begin{equation}
|
||||
\nabla Q_\lambda = -2 X^T \qty(\vec{y} - X\vec{\beta}) + \lambda \sgn\qty(\vec{\beta}),
|
||||
\end{equation}
|
||||
and \(\nabla Q_\lambda = \vec{0}\) does unfortunately not have a closed-form solution.
|
||||
A separate minimisation algorithm must therefore be used to find the parameters \(\beta_j\) which minimise the cost function.
|
||||
The calculated gradient can be put straight into the gradient descent algorithm with a suitable step length.
|
||||
Alternatively, Newton's method can be used with the second derivative (Hessian) matrix of the cost function.
|
||||
The second derivative of the absolute value is even more problematic, which is solved by setting it equal to zero.
|
||||
|
||||
The Hessian matrix of the LASSO cost function is thus the same as the second derivative of the ordinary least squares cost function,
|
||||
\begin{equation}
|
||||
H_{kl} = H_{lk} = \pdv{Q_\lambda}{\beta_l}{\beta_k}
|
||||
= \pdv{\beta_l}(\nabla_k Q)
|
||||
\,\husk{gradref}{=}\, \pdv{\beta_l}(-2X_{ik}\qty(y_i-X_{ij}\beta_j))
|
||||
= 2X_{ik}X_{il},
|
||||
\tikz[>=latex,overlay,remember picture,thick]{\draw[<-,gronn] (gradref) .. controls ++(0,0.75) and ++(-1,0) .. ++(1,0.75) node[anchor=west] {\small\ref{eq:olsgradcomp}};}
|
||||
\end{equation}
|
||||
which is the component form of the matrix equation
|
||||
\begin{equation}
|
||||
H = 2 X^T X.
|
||||
\end{equation}
|
||||
The gradient can now be rewritten as
|
||||
\begin{equation}
|
||||
\nabla Q_\lambda = -2 X^T \vec{y} + H\vec{\beta} + \lambda \sgn\qty(\vec{\beta}).
|
||||
\end{equation}
|
||||
\Vref*{eq:newtonstep} in Newton's method can then transformed into
|
||||
\begin{equation}
|
||||
H\qty(\vec{\beta}_{i+1} - \vec{\beta}_i) = 2X^T \vec{y} - H\vec{\beta}_i - \lambda \sgn\qty(\vec{\beta}_i)
|
||||
\implies H\vec{\beta}_{i+1} = 2X^T\vec{y} - \lambda \sgn\qty(\vec{\beta}_i). \label{eq:lassonewton}
|
||||
\end{equation}
|
||||
Since \(H=2X^T X\), this reduces to the normal equations of ordinary least squares (\vref*{eq:olsnormal}) in the case of \(\lambda=0\), as it should.
|
||||
|
||||
Minimisation of the LASSO cost function has the benefit that the second derivative is independent of \(\vec{\beta}\).
|
||||
\(H\) can therefore be Cholesky-decomposed once, an \(\mathcal{O}\qty(p^3)\) operation, and the result can be used to solve the linear set of equations for each iteration, which is \(\mathcal{O}\qty(p^2)\) with a pre-Cholesky-decomposed matrix.
|
||||
The Cholesky-decomposition can be used instead of an LU-decomposition, reducing the number of floating point operations by a factor of two, because \(H=2X^T X\) is positive definite when \(X\) is non-singular.
|
||||
|
||||
Evaluation of the gradient only involves matrix-vector products, which is also an \(\mathcal{O}\qty(p^2)\) operation.
|
||||
The minimisation can, therefore, be done with one initial \(\mathcal{O}\qty(p^3)\) and then only \(\mathcal{O}\qty(p^2)\) operations per iteration, while the more general problem requires both the construction and the Cholesky-decomposition of \(H\), an \(\mathcal{O}\qty(p^3)\) operation, for every single iteration.
|
||||
|
||||
LASSO regression has been implemented by calculating the Cholesky-decomposition of \(H\) using \lstinline{dpptrf}, guessing \(\beta_j=1\) and then repeatedly solving~\vref{eq:lassonewton} using \lstinline{dpptrs} until convergence, as illustrated by the following snippet.
|
||||
\lstinputlisting[linerange={lassostart-lassoend}]{lib/lasso.f90}
|
||||
|
||||
This implementation gives both good performance and accurate results for small values of \(\lambda\).
|
||||
Larger values, however, give rise to issues.
|
||||
The theory suggests that a large \(\lambda\) should force some of the parameters \(\beta_j\) to become zero, but the cost function is not differentiable, and certainly not double differentiable, when \(\beta_j\) is zero.
|
||||
Consequently, Newton's method struggles to converge when some \(\beta_j\)s are zero in the true minimum of the LASSO cost function.
|
||||
Smarter methods such as coordinate descent\cite{friedman} should therefore be used, as in e.g.\ scikit-learn.
|
||||
|
||||
|
||||
% _ _ _ _ _
|
||||
% _ __ ___ (_)_ __ (_)_ __ ___ (_)___ __ _| |_(_) ___ _ __
|
||||
%| '_ ` _ \| | '_ \| | '_ ` _ \| / __|/ _` | __| |/ _ \| '_ \
|
||||
%| | | | | | | | | | | | | | | | \__ \ (_| | |_| | (_) | | | |
|
||||
%|_| |_| |_|_|_| |_|_|_| |_| |_|_|___/\__,_|\__|_|\___/|_| |_|
|
||||
\subsection{Minimisation methods}
|
||||
The important step in a regression method, or indeed in any statistical learning method, is to minimise the cost function.
|
||||
Certain cost functions, such as those of ordinary least squares and Ridge regression, admit analytical closed-form solutions for the parameters \(\beta_j\) which minimise \(Q(\vec{\beta};D)\), while most methods, including LASSO regression and logistic regression, require usage of some general minimisation algorithm.
|
||||
|
||||
\subsubsection{Newton's method}\label{subsubsec:newton}
|
||||
Since the gradient of the cost function, \(\nabla Q\), is a perfectly normal function from \(\mathbb{R}^p\) to \(\mathbb{R}^p\), it can be Taylor expanded around some point \(\vec{\beta}_0\).
|
||||
This Taylor expansion can then be evaluated at a point \(\vec{\beta}_0 + \delta\vec{\beta}\). To first order in \(\delta\vec{\beta}\),
|
||||
\begin{equation}
|
||||
\nabla Q\qty(\vec{\beta}_0 + \delta\vec{\beta}) = \nabla Q\qty(\vec{\beta}_0) + H\qty(\vec{\beta}_0)\delta\vec{\beta},
|
||||
\end{equation}
|
||||
where \(H(\vec{\beta}_0)\) is the derivative of \(\nabla Q\), i.e.\ the Hessian of \(Q\), evaluated at \(\beta_0\).
|
||||
Given a guess \(\beta_0\) somewhere near the minimum of \(Q\), a better estimate can be found by solving for the \(\delta\vec{\beta}\) which makes \(\nabla Q(\vec{\beta}_0 + \delta\vec{\beta}) = \vec{0}\), i.e.
|
||||
\begin{equation}
|
||||
H\qty(\vec{\beta}_0)\delta\vec{\beta} = -\nabla Q\qty(\vec{\beta}_0),
|
||||
\end{equation}
|
||||
and, in general, solving
|
||||
\begin{equation}
|
||||
H\qty(\vec{\beta}_i)\delta\vec{\beta}_i = -\nabla Q\qty(\vec{\beta}_i),\label{eq:newtonstep}
|
||||
\end{equation}
|
||||
letting \(\vec{\beta}_{i+1} = \vec{\beta}_i + \delta\vec{\beta}_i\) and continuing the process until the method has converged sufficiently.
|
||||
The above equation is a simple linear set of equations.
|
||||
|
||||
|
||||
\subsubsection{Gradient descent}
|
||||
Evaluating the Hessian matrix and inverting it can sometimes be infeasible, for example due to poor time complexity or being woefully undefined.
|
||||
In such cases, one can replace the Hessian \(H\) with a number \(1/\alpha\).
|
||||
\Vref{eq:newtonstep} is then rewritten as
|
||||
\begin{equation}
|
||||
\vec{\beta}_{i+1} = \vec{\beta_i} - \alpha \nabla Q\qty(\vec{\beta}_i),
|
||||
\end{equation}
|
||||
with a simple geometric interpretation:
|
||||
By going a small step in the opposite direction of the gradient, the value of the cost function will decrease and \(\vec{\beta}\) will approach the true minimum.
|
||||
A small step will guarantee convergence for a convex function at the cost of slow convergence, while a larger value may give faster convergence or not converge at all.
|
||||
|
||||
% __
|
||||
% _ __ ___ _ __ / _| ___ _ __ _ __ ___ __ _ _ __ ___ ___
|
||||
% | '_ \ / _ \ '__| |_ / _ \| '__| '_ ` _ \ / _` | '_ \ / __/ _ \
|
||||
% | |_) | __/ | | _| (_) | | | | | | | | (_| | | | | (_| __/
|
||||
% | .__/ \___|_| |_| \___/|_| |_| |_| |_|\__,_|_| |_|\___\___|
|
||||
% |_|
|
||||
\subsection{Performance of regression methods}
|
||||
While the cost function is minimised by all regression methods, its value is not particularly meaningful.
|
||||
In particular, the cost functions in Ridge and LASSO regression depend on the parameter \(\lambda\), and a measure independent of \(\lambda\) should be used to determine the optimal \(\lambda\).
|
||||
Other functions are therefore introduced to measure the performance of a given regression method and/or a choice of parameters.
|
||||
|
||||
The simplest measure is the mean square error,
|
||||
\begin{equation}
|
||||
\MSE\qty(\vec{\beta},D) = \frac{1}{N}\norm{\vec{y}-\vec{\tilde{y}}}_2^2 = \frac{1}{N} \sum_{i=1}^N \qty(y_i - \tilde{y}_i)^2,
|
||||
\end{equation}
|
||||
which should be as small as possible. Another measure is the \(R^2\) score, defined as
|
||||
\begin{equation}
|
||||
R^2(\vec{\beta},D) = 1 - \frac{\norm{\vec{y} - \vec{\tilde{y}}}_2^2}{\norm{\vec{y} - \bar{y}}_2^2}
|
||||
= 1 - \frac{\sum_{i=1}^N \qty(y_i - \tilde{y}_i)^2}{\sum_{i=1}^N \qty(y_i - \bar{y})^2},
|
||||
\end{equation}
|
||||
where \(\bar{y}\) is the mean of the measured values. The \(R^2\) score should be as close to \(1\) as possible.
|
||||
|
||||
Prediction and these measures of performance can be applied to two main types of data.
|
||||
Firstly, it can be applied to the data from which \(\vec{\beta}\) was derived, which is called the training data.
|
||||
Ordinary least squares, corresponding to \(\lambda=0\) for the other methods, will, by definition, give the best results for this data set.
|
||||
Secondly, the performance can be measured for values not among the training data, called test data. Ridge and LASSO are expected to outperform ordinary least squares for small, non-zero values of \(\lambda\) for this category of data.
|
||||
|
||||
% _ _
|
||||
% _ __ ___ ___ __ _ _ __ ___ _ __ | (_)_ __ __ _
|
||||
% | '__/ _ \/ __|/ _` | '_ ` _ \| '_ \| | | '_ \ / _` |
|
||||
% | | | __/\__ \ (_| | | | | | | |_) | | | | | | (_| |
|
||||
% |_| \___||___/\__,_|_| |_| |_| .__/|_|_|_| |_|\__, |
|
||||
% |_| |___/
|
||||
\subsection{Resampling methods}
|
||||
Resampling methods are techniques to improve the prediction accuracy and obtain estimates for quantities such as the variance of \(\vec{\beta}\) by repeatedly dividing the data set into training and test data.
|
||||
Two simple examples are \(k\)-fold cross-validation and bootstrapping.
|
||||
The former partitions the data set into \(k\) partitions.
|
||||
One of these subsets is chosen as test data on which the performance is measured, while the rest are used as training data.
|
||||
This process is the repeated \(k\) times, so that each subset is used once as test data.
|
||||
Bootstrapping, on the other hand, repeatedly creates training data sets by randomly selecting \(N\) values from the data set with replacement, while another, separate data set is used as test data each time, as illustrated by the snippet below.
|
||||
\lstinputlisting[linerange={bootstrapstart-bootstrapend}]{lib/bootstrap.f90}
|
||||
|
||||
% _ _ _
|
||||
% | |__ (_) __ _ ___ __ _ _ __ __| |
|
||||
% | '_ \| |/ _` / __| / _` | '_ \ / _` |
|
||||
% | |_) | | (_| \__ \ | (_| | | | | (_| |
|
||||
% |_.__/|_|\__,_|___/ \__,_|_| |_|\__,_|
|
||||
% __ ____ _ _ __(_) __ _ _ __ ___ ___
|
||||
% \ \ / / _` | '__| |/ _` | '_ \ / __/ _ \
|
||||
% \ V / (_| | | | | (_| | | | | (_| __/
|
||||
% \_/ \__,_|_| |_|\__,_|_| |_|\___\___|
|
||||
\subsection{Bias and variance}
|
||||
The choice of the number of basis functions, \(p\), determines the complexity of the model to which the data set is fitted.
|
||||
A higher complexity makes the model a better approximation to the training data, which is said to reduce the \emph{bias} of the model.
|
||||
On the other hand, a higher complexity may decrease the model's ability to predict reasonable values for test data, since a more complex model will be more affected by noise in the model.
|
||||
This is said to increase the model's \emph{variance}.
|
||||
The best prediction performance is therefore achieved for a complexity which balances bias and variance, which is called the bias-variance trade-off\cite{mehta}.
|
||||
|
||||
Following~\cite{mehta}, a mathematical manifestation of the bias-variance trade-off can be derived by assuming that the measured data set \(D=\qty{\qty(\vec{x}_i,y_i)}_{i=1}^N\) is generated by a combination of some model \(f\qty(\vec{x})\) (for example Franke's function) and random noise, specifically
|
||||
\begin{equation}
|
||||
y_i = f\qty(\vec{x}_i) + \varepsilon_i \eqqcolon f_i + \varepsilon_i,
|
||||
\end{equation}
|
||||
where the variables \(\varepsilon_i\) are independent and normally distributed with variance \(\sigma^2\) around zero.
|
||||
Given a data set \(D\), a prediction \(\tilde{y}_i^D\) can be made by any of the regression methods discussed above.
|
||||
An expected mean square error can be found by averaging over all datasets \(D\) and all noises \(\vec{\varepsilon}\),
|
||||
\begin{alignat}{2}
|
||||
E_{D,\varepsilon}\qty[\norm{\vec{y}-\vec{\tilde{y}}_D}_2^2]
|
||||
&= E_{D,\varepsilon}\qty[\norm{\vec{y} - \vec{f} + \vec{f} - \vec{\tilde{y}}_D}_2^2] \label{eq:tmp}
|
||||
= E_{D,\varepsilon}\qty[\norm{\vec{y} - \vec{f}}_2^2 + \norm{\vec{f} - \vec{\tilde{y}}_D}_2^2
|
||||
+ 2 \qty(\vec{y} - \vec{f})\cdot\qty(\vec{f} - \vec{\tilde{y}}_D)],
|
||||
\shortintertext{and using the linearity of the expectation value,}
|
||||
&= E_{D,\varepsilon}\qty[\norm{\vec{y} - \vec{f}}_2^2] + E_{D,\varepsilon}\qty[\norm{\vec{f} - \vec{\tilde{y}}_D}_2^2]
|
||||
+ 2 E_{D,\varepsilon}\qty[\qty(\vec{y} - \vec{f})\cdot\qty(\vec{f} - \vec{\tilde{y}}_D)]\\
|
||||
&= E_{D,\varepsilon}\qty[\norm{\vec{\varepsilon}}_2^2] + E_{D,\varepsilon}\qty[\norm{\vec{f} - \vec{\tilde{y}}_D}_2^2]
|
||||
+ 2 E_{D,\varepsilon}\qty[\vec{\varepsilon}\cdot\qty(\vec{f} - \vec{\tilde{y}}_D)]\\
|
||||
&= \sigma^2 + E_{D,\varepsilon}\qty[\norm{\vec{f} - \vec{\tilde{y}}_D}_2^2]
|
||||
+ 2 \cancel{E_{\varepsilon}\qty[\vec{\varepsilon}]}\cdot E_{D}\qty[\qty(\vec{f} - \vec{\tilde{y}}_D)]\\
|
||||
&= \sigma^2 + E_{D}\qty[\norm{\vec{f} - \vec{\tilde{y}}_D}_2^2].
|
||||
\end{alignat}
|
||||
This expression can be further decomposed by adding and subtracting the average prediction, \(E_D\qty[\vec{\tilde{y}}_D]\),
|
||||
\begin{alignat}{2}
|
||||
E_{D}\qty[\norm{\vec{y}-\vec{\tilde{y}}_D}_2^2]
|
||||
&= \sigma^2 + E_{D}\qty[\norm{\vec{f}- E_D\qty[\vec{\tilde{y}}_D] + E_D\qty[\vec{\tilde{y}}_D] - \vec{\tilde{y}}_D}_2^2]\\
|
||||
&= \sigma^2 + \norm{\vec{f}- E_D\qty[\vec{\tilde{y}}_D]}_2^2 + E_D\qty[\norm{E_D\qty[\vec{\tilde{y}}_D] - \vec{\tilde{y}}_D}_2^2]\\
|
||||
&\phantom{{}={}} + 2E_D\qty[\qty(\vec{f}- E_D\qty[\vec{\tilde{y}}_D])\cdot\cancel{\qty(E_D\qty[\vec{\tilde{y}}_D] - \vec{\tilde{y}}_D)}]\\
|
||||
&= \sigma^2 + \norm{\vec{f}- E_D\qty[\vec{\tilde{y}}_D]}_2^2 + E_D\qty[\norm{E_D\qty[\vec{\tilde{y}}_D] - \vec{\tilde{y}}_D}_2^2].
|
||||
\end{alignat}
|
||||
The first term, \(\sigma^2\), represents the noise in the generated data, which no model can overcome.
|
||||
The second term measures the deviation of the average prediction from the true, noise-free model, which is the (squared) bias.
|
||||
Lastly, the third term measures how much the predictions vary and is called the variance.
|
||||
A complex model (large \(p\)) will minimise the second term, since a complex model will be better suited to fit the true model, \(f\), while the third term will increase with model complexity since the fitting procedure will be more sensitive to noise in the different data sets.
|
||||
|
||||
Unfortunately, the above bias-variance decomposition requires knowledge of the exact model behind the data, \(f\).
|
||||
A more computationally practical expression could have been derived by adding and subtracting \(E_D\qty[\vec{\tilde{y}}_D]\) in~\vref{eq:tmp}, which gives the expression
|
||||
\begin{equation}
|
||||
E_{D}\qty[\norm{\vec{y}-\vec{\tilde{y}}_D}_2^2]
|
||||
= \norm{\vec{y}- E_D\qty[\vec{\tilde{y}}_D]}_2^2 + E_D\qty[\norm{E_D\qty[\vec{\tilde{y}}_D] - \vec{\tilde{y}}_D}_2^2]
|
||||
\end{equation}
|
||||
by manipulations analogous to the ones used above. This is another formulation of the bias-variance decomposition which can be used without any information about the underlying model. Dividing by \(N\) gives the average \(\MSE\) on the left-hand side.
|
||||
|
||||
|
||||
|
||||
% _ _ _ _ _
|
||||
% (_)_ __ ___ _ __ | | ___ _ __ ___ ___ _ __ | |_ __ _| |_(_) ___ _ __
|
||||
% | | '_ ` _ \| '_ \| |/ _ \ '_ ` _ \ / _ \ '_ \| __/ _` | __| |/ _ \| '_ \
|
||||
% | | | | | | | |_) | | __/ | | | | | __/ | | | || (_| | |_| | (_) | | | |
|
||||
% |_|_| |_| |_| .__/|_|\___|_| |_| |_|\___|_| |_|\__\__,_|\__|_|\___/|_| |_|
|
||||
% |_|
|
||||
\subsection{Implementation}
|
||||
The three regression methods and bootstrapping have been implemented in a polymorphic class hierarchy in Fortran.
|
||||
An object-oriented approach ensures ease of reuse for later projects, while also simplifying certain parts of the code and organisation.
|
||||
For instance, the bootstrapping algorithm can take in a regression method object, which is guaranteed to have methods for prediction and fitting, and not care whether ordinary least squares, Ridge or LASSO regression is used. Fortran was chosen because it combines excellent performance with object-oriented capabilities and a nice syntax for numerical work.
|
||||
\tikzexternalenable
|
||||
\tikzsetnextfilename{uml}
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{umlpackage}{Linear regression and bootstrapping}
|
||||
\umlclass[type=abstract]{regressor}{
|
||||
\lstinline{real :: X(N, p), beta(p)}\\
|
||||
\lstinline{class(basis_function) :: basis(p)}
|
||||
}{
|
||||
\lstinline{subroutine create_X(real :: x_values(N,:))}\\
|
||||
\lstinline{subroutine predict(real :: x(N,:), y(N); optional :: y_exact(N), mse, r2)}\\
|
||||
\lstinline{deferred subroutine fit(real, optional :: x_values(N,:); real :: y_values(N))}
|
||||
}
|
||||
\umlclass[below = 1cm of regressor.south, anchor=north]{ridge}{
|
||||
\lstinline{real :: lambda}
|
||||
}{}
|
||||
\umlemptyclass[left=1cm of ridge.north west, anchor=north east]{ols}
|
||||
\umlclass[right = 1cm of ridge.north east, anchor=north west]{lasso}{
|
||||
\lstinline{real :: lambda}\\
|
||||
\lstinline{real :: tolerance}
|
||||
}{}
|
||||
\umlinherit[geometry=|-|,arm1=0.3cm,arm2=0.3cm,anchor1=90,anchor2=240]{ols}{regressor}
|
||||
\umlinherit[geometry=|-|,arm1=0.3cm,arm2=0.3cm,anchor1=90,anchor2=270]{ridge}{regressor}
|
||||
\umlinherit[geometry=|-|,arm1=0.3cm,arm2=0.3cm,anchor1=90,anchor2=300]{lasso}{regressor}
|
||||
|
||||
\umlclass[type=abstract, below=5.5cm of regressor.west,anchor=north west]{basisfunction}{}{
|
||||
\lstinline{deferred real function eval(real :: x(:))}
|
||||
}
|
||||
\umlclass[right=1cm of basisfunction.north east,anchor=north west]{polynomial2d}{
|
||||
\lstinline{integer :: x1_pow, x2_pow}
|
||||
}{}
|
||||
\umlinherit{polynomial2d}{basisfunction}
|
||||
|
||||
\umluniassoc[pos1=2.9, mult1=1..*, geometry=|-|, anchor1=202, anchor2=165,weight=0.35]{regressor}{basisfunction}
|
||||
|
||||
\umlclass[above=1cm of regressor.north, anchor=south]{bootstrapper}{
|
||||
\lstinline{class(regressor) :: fitter}\\
|
||||
\lstinline{real :: y_predictions(N, num_bootstraps), R2s(num_bootstraps), MSEs(num_bootstraps)}\\
|
||||
\lstinline{real :: betas(p, num_bootstraps), mean_beta(p), beta_variance(p)}
|
||||
}{
|
||||
\lstinline{subroutine bootstrap(real :: x(:,:), y(:); integer :: num_bootstraps)}
|
||||
}
|
||||
\umluniassoc[mult1=1,pos1=0.8]{bootstrapper}{regressor}
|
||||
\end{umlpackage}
|
||||
\end{tikzpicture}
|
||||
\caption{Class hierarchy of implementation. All \lstinline{real} variables are of kind \lstinline{real64} from \lstinline{iso_fortran_env}.}
|
||||
\end{figure}
|
||||
\tikzexternaldisable
|
||||
|
||||
% _ _ __
|
||||
% __ _ _ __ __ _| |_ _ ___(_)___ ___ / _|
|
||||
% / _` | '_ \ / _` | | | | / __| / __| / _ \| |_
|
||||
% | (_| | | | | (_| | | |_| \__ \ \__ \ | (_) | _|
|
||||
% \__,_|_| |_|\__,_|_|\__, |___/_|___/ \___/|_|
|
||||
% _ _|___/
|
||||
% _ __ ___ ___| |_| |__ ___ __| |___
|
||||
% | '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|
|
||||
% | | | | | | __/ |_| | | | (_) | (_| \__ \
|
||||
% |_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/
|
||||
\section{Analysis of methods}
|
||||
The accuracy and viability of the methods can be studied on a test case where the underlying function is known.
|
||||
Such an example is Franke's function shown in~\vref{fig:franke}.
|
||||
|
||||
\subsection{Verification of theory and fitting}
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{subfigure}{\textwidth}
|
||||
\centering
|
||||
\includegraphics{figs/verification_OLS.pdf}
|
||||
\caption{Ordinary least squares regression.}
|
||||
\end{subfigure}
|
||||
\begin{subfigure}{\textwidth}
|
||||
\centering
|
||||
\includegraphics{figs/verification_Ridge.pdf}
|
||||
\caption{Ridge regression.}
|
||||
\end{subfigure}
|
||||
% \begin{subfigure}{\textwidth}
|
||||
% \centering
|
||||
% \includegraphics{figs/verification_LASSO.pdf}
|
||||
% \caption{LASSO regression.}
|
||||
% \end{subfigure}
|
||||
\caption{Comparison of Franke's function with noise (red) and a fifth degree polynomial approximation (multi-coloured). The colour of the polynomial approximation shows its deviation from Franke's function.~\Vref{fig:betas} shows that the parameters are significantly different for the different methods, yet the resulting approximations are approximately equal.}
|
||||
\end{figure}
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[thick,
|
||||
height=3in,
|
||||
width=6in,
|
||||
grid,
|
||||
xlabel = {\(i\)},
|
||||
ylabel = {\(\beta_i\)},
|
||||
legend style={draw=none,at={(0.98,0.98)},anchor=north east},
|
||||
legend cell align=left
|
||||
]
|
||||
\foreach \method in {OLS, Ridge}{
|
||||
\addplot+[only marks,mark=o] plot[error bars/y dir=both, error bars/y explicit] table[y error=uncertainty] {data/verification_mean_beta_\method.dat};
|
||||
\addlegendentryexpanded{\method};
|
||||
}
|
||||
\addplot+[only marks,mark=o] plot[error bars/y dir=both, error bars/y explicit] table[y error=uncertainty] {data/verification_mean_beta_sklearn.dat};
|
||||
\addlegendentryexpanded{LASSO (scikit-learn)};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{Comparison of parameters for the different regression methods when applied to Franke's function with noise and using polynomials up to fifth order. The Ridge coefficients are severely suppressed, yet~\vref{table:performance} shows that the fitting quality is virtually the same. LASSO regression (using scikit-learn) gives a sparse solution, i.e.\ some \(\beta_j\)s are exactly zero. Confidence intervals are estimated as \(\beta_j \pm 2\sigma_{\beta_j}\), where \(\sigma_{\beta_j}\) is the standard deviation of \(\beta_j\) calculated from bootstrapping.}\label{fig:betas}
|
||||
\end{figure}
|
||||
|
||||
\begin{table}[H]
|
||||
\centering
|
||||
\caption{Measure of performance on training and test data sampled from Franke's function with noise for the different methods. Ordinary least squares should, by design, perform best on the training data, while Ridge and LASSO are expected to perform better on test data not used in the fitting due to their regularisation parameter. The advantage of regularised methods is amplified when a smaller number of points is used (\(\num{10000}\) points were used here for plotting purposes). See also~\vref{tab:biasvar}, which shows the average error on test data for many bootstrap samples.}\label{table:performance}
|
||||
\pgfplotstabletypeset[zerofill, fixed,precision=4]{./data/verification_mse_r2.dat}
|
||||
\end{table}
|
||||
|
||||
\subsection{Bias and variance}
|
||||
According to the bias-variance decomposition, the mean squared error should be equal to the sum of the bias (squared) and the variance.
|
||||
The bias and variance are calculated from bootstrapping, and the results show good agreement with the theoretical prediction.
|
||||
\begin{table}[H]
|
||||
\centering
|
||||
\caption{Verification of the bias-variance decomposition via bootstrapping for fitting of Franke's function with noise. As expected, the sum of the bias and the variance is equal to the mean squared error, and the regularised methods perform better on test data than ordinary least squares.}\label{tab:biasvar}
|
||||
\pgfplotstabletypeset[sci, sci zerofill, precision=3]{./data/verification_bias_variance.dat}
|
||||
\end{table}
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[thick,
|
||||
height=3in,
|
||||
width=6in,
|
||||
ymode=log,
|
||||
grid,
|
||||
xtick distance=1,
|
||||
xlabel = {Polynomial degree \(d\)},
|
||||
ylabel = {Error},
|
||||
legend style={draw=none,at={(0.02,0.50)},anchor=west},
|
||||
legend cell align=left
|
||||
]
|
||||
\addplot+ table[y index=1] {data/complexity.dat};
|
||||
\addlegendentry{MSE (OLS)};
|
||||
\addplot+ table[y index=2] {data/complexity.dat};
|
||||
\addlegendentry{Bias (OLS)};
|
||||
\addplot+ table[y index=3] {data/complexity.dat};
|
||||
\addlegendentry{Variance (OLS)};
|
||||
\addplot+ table[y index=4] {data/complexity.dat};
|
||||
\addlegendentry{MSE (Ridge)};
|
||||
\addplot+ table[y index=5] {data/complexity.dat};
|
||||
\addlegendentry{Bias (Ridge)};
|
||||
\addplot+ table[y index=6] {data/complexity.dat};
|
||||
\addlegendentry{Variance (Ridge)};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{Visualisation of bias-variance trade-off. The mean squared error for test data is expected to first decrease with increasing complexity as the bias decreases towards \(\sigma^2\), and then increase as the variance becomes dominant. Additionally, the regularisation of the Ridge regressor is intended to keep the variance moderate even with a high polynomial degree.}\label{fig:biasvar}
|
||||
\end{figure}
|
||||
|
||||
\subsection{Effect of noise}
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[
|
||||
thick,
|
||||
width=6in,
|
||||
height=3in,
|
||||
xlabel = {\(\sigma\)},
|
||||
ylabel = {\(R^2\)},
|
||||
legend cell align=left,
|
||||
legend style={draw=none,at={(0.98,0.98)},anchor=north east},
|
||||
grid
|
||||
]
|
||||
\addplot+ table[y index=1] {data/noise.dat};
|
||||
\addlegendentry{Training data};
|
||||
\addplot+ table[y index=2] {data/noise.dat};
|
||||
\addlegendentry{Test data};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{\(R^2\) score of Ridge regression with \(\lambda=0.001\) when normally distributed noise with a standard deviation of \(\sigma\) is added to Franke's function. \(400\) data points are used, and the performance is clearly better for the training data to which the model was fitted than novel test data. A higher number of data points reduces the overfitting, so that the prediction and its performance are less affected by noise.}
|
||||
\end{figure}
|
||||
|
||||
\subsection{Effect of regularisation}
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[
|
||||
thick,
|
||||
width=6in,
|
||||
height=3in,
|
||||
xmode=log,
|
||||
xlabel = {\(\lambda\)},
|
||||
ylabel = {\(R^2\)},
|
||||
legend cell align=left,
|
||||
legend style={draw=none,at={(0.02,0.40)},anchor=west},
|
||||
grid
|
||||
]
|
||||
\addplot+ table[y index=1] {data/r2_lambda.dat};
|
||||
\addlegendentry{Training data};
|
||||
\addplot+ table[y index=2] {data/r2_lambda.dat};
|
||||
\addlegendentry{Test data};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{\(R^2\)-score as a function of the regularisation parameter \(\lambda\) for Ridge regression. \(\lambda=0\) would correspond to ordinary least squares. While ordinary least squares regression by construction gives the best performance on training data, the results show that a smart choice of \(\lambda\) can give better predicting performance since the regularisation reduces overfitting and variance.}
|
||||
\end{figure}
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[
|
||||
thick,
|
||||
width=6in,
|
||||
height=3in,
|
||||
xmode=log,
|
||||
ymode=log,
|
||||
xlabel = {\(\lambda\)},
|
||||
ylabel = {Error},
|
||||
legend cell align=left,
|
||||
legend style={draw=none,at={(0.98,0.98)},anchor=north east},
|
||||
grid
|
||||
]
|
||||
\addplot+ table[y index=1] {data/bivar_lambda.dat};
|
||||
\addlegendentry{MSE};
|
||||
\addplot+ table[y index=2] {data/bivar_lambda.dat};
|
||||
\addlegendentry{Bias};
|
||||
\addplot+ table[y index=3] {data/bivar_lambda.dat};
|
||||
\addlegendentry{Variance};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{Error, bias and variance as a function of the regularisation parameter \(\lambda\) for Ridge regression, calculated from bootstrapping. The variance, which is the error due to overfitting of the noise, gradually decreases as \(\lambda\) increases, while a too large \(\lambda\) causes the bias to increase because the severe suppression of the \(\beta_j\) coefficients makes it impossible for the polynomial expansion to approximate the model function (Franke's function).}
|
||||
\end{figure}
|
||||
|
||||
|
||||
\pgfplotstableread[skip first n=1]{data/beta_lambda.dat}{\betatable}
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[
|
||||
thick,
|
||||
width=6in,
|
||||
height=3in,
|
||||
xmode=log,
|
||||
xlabel = {\(\lambda\)},
|
||||
ylabel = {\(\beta_i\)},
|
||||
legend cell align=left,
|
||||
legend style={draw=none,at={(0.98,0.98)},anchor=north east},
|
||||
grid
|
||||
]
|
||||
\foreach \index in {1,3,...,44}{
|
||||
\addplot+ table[y index=\index, y error index={\index+1}] {\betatable};
|
||||
}
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{Coefficients \(\beta_j\) as a function of the regularisation parameter \(\lambda\) for Ridge regression. Regularisation suppresses the coefficients due to the penalising \(\lambda\|\vec{\beta}\|_2^2\) term in the cost function, in order to achieve better prediction abilities due to lower variance and effect of noise. LASSO regression would set many of the coefficients to exactly zero, ref.~\vref{fig:betas}.}
|
||||
\end{figure}
|
||||
|
||||
% _
|
||||
% __ _ ___ ___ __ _ _ __ __ _ _ __ | |__ _ _
|
||||
% / _` |/ _ \/ _ \ / _` | '__/ _` | '_ \| '_ \| | | |
|
||||
%| (_| | __/ (_) | (_| | | | (_| | |_) | | | | |_| |
|
||||
% \__, |\___|\___/ \__, |_| \__,_| .__/|_| |_|\__, |
|
||||
% |___/ |___/ |_| |___/
|
||||
\section{Geographical data}
|
||||
The previous section verified the implementation of ordinary least squares and Ridge regression (and scikit-learn's LASSO regression), as well as some theoretical predictions.
|
||||
It is now time to apply them to real-world data, such as the geographical data in the figure below.
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\includegraphics{figs/geography.pdf}
|
||||
\caption{Some Norwegian geography, complete with a fjord.}\label{fig:geo}
|
||||
\end{figure}
|
||||
|
||||
To determine the best possible model for the terrain, I have chosen to iterate over different polynomial degrees and compare ordinary least squares with Ridge regression for a selection of \(\lambda\)s for each degree. The best \(\lambda\) is determined for each degree, and the overall best approximation (measured by MSE for test data in bootstrap) is visualised in~\vref{fig:geoapprox}.
|
||||
\begin{table}[H]
|
||||
\centering
|
||||
\caption{Errors for different approximations to~\vref{fig:geo}. Several values of \(\lambda\) for Ridge regression are tried for each degree, and the one resulting in the smallest test MSE during bootstrapping is reported. Ordinary least squares is seen to perform well for polynomials of a small degree, while higher order polynomials lead to severe overfitting. On the other hand, the regularisation of Ridge manages to keep the variance and overfitting in check and give reasonable results also for higher degrees.}\label{tab:geo}
|
||||
\small
|
||||
\pgfplotstabletypeset[columns/lambda/.style={column name={$\lambda$}},columns/$d$/.style={fixed},sci, sci zerofill,precision=1]{data/geography_mse.dat}
|
||||
\end{table}
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\begin{tikzpicture}
|
||||
\begin{axis}[
|
||||
thick,
|
||||
width=6in,
|
||||
height=3in,
|
||||
ymode=log,
|
||||
xlabel = {\(d\)},
|
||||
ylabel = {Error},
|
||||
legend cell align=left,
|
||||
legend style={draw=none,at={(0.98,0.02)},anchor=south east},
|
||||
grid
|
||||
]
|
||||
\addplot+ table[skip first n=1, y index=4] {data/geography_mse.dat};
|
||||
\addlegendentry{MSE};
|
||||
\addplot+ table[skip first n=1,y index=5] {data/geography_mse.dat};
|
||||
\addlegendentry{Bias};
|
||||
\addplot+ table[skip first n=1,y index=6] {data/geography_mse.dat};
|
||||
\addlegendentry{Variance};
|
||||
\end{axis}
|
||||
\end{tikzpicture}
|
||||
\caption{Error, bias and variance as a function of polynomial degree \(d\) for Ridge regression. The parameter \(\lambda\) is estimated for each \(d\) by bootstrapping for different \(\lambda\)s and choosing the one which gives the smallest mean squared error for the test data. Regularisation keeps the variance low even for a high polynomial degree, which is not the case for ordinary least squares according to~\vref{tab:geo}.}
|
||||
\end{figure}
|
||||
Bias and variance still add up to the mean squared error with approximately the same accuracy as when approximating Franke's function with noise. This was not necessarily expected due to the assumption of identically, independently and normally distributed noise in the derivation, which was added artificially. Geographical data is certainly not independent, and the other assumptions may also not hold.
|
||||
|
||||
Ideally, there should be some polynomial degree for which the error is minimised, since overfitting will occur and the variance will increase as the degree of the polynomial grows and the approximation becomes more sensitive to noise. Overfitting requires that the number of basis functions, \(p=(d+1)(d+2)/2\), becomes non-negligible compared to the number of data points. I have here used approximately \(\num{65000}\) data points, meaning that a polynomial must be of a very high degree in order to overfit. Ridge regression scales as \(\mathcal{O}(p^3) = \mathcal{O}(d^6)\) and ordinary least squares possibly even worse, so the runtime increases drastically when the polynomial degree becomes large. The simulation giving the figures and tables in this section took roughly 100 CPU hours.
|
||||
|
||||
\begin{figure}[H]
|
||||
\centering
|
||||
\includegraphics{figs/geography_approx.pdf}
|
||||
\caption{Polynomial approximation to the geographical data in\vref{fig:geo} using the best model from~\vref{tab:geo}. The model captures key elements such as the fjord, while missing the finer details of the mountainous landscape. In particular, the steep climb from the fjord is almost step-like, which no polynomial of finite degree can approximate.}\label{fig:geoapprox}
|
||||
\end{figure}
|
||||
|
||||
% _ _
|
||||
% ___ ___ _ __ ___| |_ _ ___(_) ___ _ __
|
||||
% / __/ _ \| '_ \ / __| | | | / __| |/ _ \| '_ \
|
||||
%| (_| (_) | | | | (__| | |_| \__ \ | (_) | | | |
|
||||
% \___\___/|_| |_|\___|_|\__,_|___/_|\___/|_| |_|
|
||||
\section{Summary and conclusion}
|
||||
This report started with a lengthy discussion of theory, culminating in three different regression methods, some error analysis and a few implementational tidbits.
|
||||
Ordinary least squares is derived from a minimisation of the training error, which has been proven correct.
|
||||
The bias-variance decomposition showed that the predicting error can be decomposed into bias, which measures the approximating function's ability to model the data's underlying function, and variance, which measures how much the approximation varies when applied to different parts of the data set, such as in bootstrapping.
|
||||
Predicting models obtained from ordinary least squares have been inferior to Ridge and LASSO regression when applied to test data not used in the fitting procedure, as these methods have a regularisation term in their cost function which prevents the regressor from being affected by noise and overfitting to the same extent as ordinary least squares.
|
||||
|
||||
Application of the regression techniques to Franke's function showed that the regularised methods perform better than ordinary least squares when a data set with few points and much noise is approximated, while a large number of data points reduces overfitting and makes ordinary least squares on par with Ridge and LASSO regression as long as a model with moderate complexity is used.
|
||||
Newton's method proved ineffective at minimising the LASSO cost function due to its non-differentiability, while the implementations of ordinary least squares, Ridge regression and bootstrapping performed well in terms of both accuracy and performance, the latter being crucial for large data sets with polynomials of a high degree.
|
||||
|
||||
Lastly, geographical data from a Norwegian fjord was fitted. Ordinary least squares proved to be utterly useless unless a moderate order of polynomials is used, while Ridge regression performed well for a variety of polynomial degrees when the optimal regularisation parameter was chosen for each polynomial degree. The main geographical features were reproduced, despite difficult features such as a Heaviside-like incline on the side of the fjord.
|
||||
|
||||
Future work includes using a prediction method which scales better with data and model complexity, such as neural networks. Additionally, a better version of the LASSO technique can be implemented by using a minimisation algorithm such as coordinate descent. A sparse solution with some \(\beta_j\)s exactly equal to zero is then expected due to the use of a soft thresholding operator, and its performance should be evaluated on e.g. Franke's function and geographical data and compared with the other regression methods.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
\clearpage
|
||||
\nocite{*}
|
||||
\printbibliography{}
|
||||
\addcontentsline{toc}{section}{\bibname}
|
||||
\end{document}
|
||||
Reference in New Issue
Block a user