CN114417543A - 用于优化的设备和方法 - Google Patents

Abstract

本申请涉及用于优化的设备和方法。一种优化设备通过运行在施加到表示目标问题的伊辛模型的磁场减小时发生的伊辛模型的状态变化的模拟来求解伊辛模型的基态。在这样做时，所述优化设备将对应于噪声的值添加到模拟中使用的一些系数。然后，所述优化设备执行随着模拟中的时间推进而减小磁场强度的实数时间传播的第一处理和基于虚数时间传播方法减小伊辛模型的能量的第二处理。

Description

用于优化的设备和方法

技术领域

本文讨论的实施方式涉及用于优化的设备和方法。

背景技术

在自然科学和社会科学中经常遇到用于求解评估(也称为能量、函数)的最小值(或产生最小值的评估函数的状态变量的值的组合)的最小化问题(或组合优化问题)。近年来，使用描述相互作用磁自旋行为的伊辛模型(Ising model)来系统阐述这一问题的趋势正在加速。

一种常用的解决方案被称为模拟退火(SA)，其使用马尔科夫蒙特卡洛(MarkovChain Monte Carlo：MCMC)引入随机过程，定义了诸如玻尔兹曼分布的特定分布，然后引入温度参数并人为地降低温度。然而，为了通过模拟退火得到精确的解，温度计划需要与时间的对数成反比。也就是说，快速冷却并不总是得到最优解决方案。

在这种情况下，基于量子力学的计算技术和机器的发展近年来变得活跃。这一趋势的技术背景是伊辛量子计算机的实现。伊辛量子计算机以量子退火(QA)为其逻辑基础。量子计算机有望解决通用的带有电子电路的经典计算机在现实的时间范围内难以解决的问题。在现阶段，并不是所有的问题都能通过量子计算机得到快速的解决，而是部分问题能通过量子计算机得到快速的解决。这类问题的典型示例是素数因式分解问题，这是诸如Rivest-Shamir-Adleman(RSA)密码系统的密码通信的基础。因此，由于量子计算机(尽管是伊辛模型计算机)的实现，探索量子计算机应用范围的趋势加快了。例如，请参见以下文档。

日本特开2018-067200号公报

日本特开2020-64535号公报

日本特开2020-9301号公报

Sanroku Tsukamoto、Motomu Takatsu、Satoshi Matsubara和Hirotaka Tamura，“组合优化问题的加速器结构”，《富士通科技期刊》，第53卷第5期，第8-13页，2017年9月。

S.Kirkpatrick、C.D.Gelatt和M.P.Vecchi，“模拟退火优化”，《科学》，第220卷，第4598期，第671-680页，1983年5月13日。

Koji Hukushima和Koji Nemoto，“交换蒙特卡罗方法及其在自旋玻璃模拟中的应用”，《日本物理学会期刊》，第65卷第6期，第1604-1608页，1996年6月。

E.D.Dahl，“D-Wave编程：地图着色问题”，《D-Wave系统白皮书》，2013年11月。

Tadashi Kadowaki和Hidetoshi Nishimori，“横向伊辛模型中的量子退火”，Phys，修订版，第58卷第5期，第5355-5363页，1998年11月。

Hidetoshi Nishimori和Kabuki Takada，“非随机哈密顿量对量子退火效率的指数增强”，《ICT前沿》，第4卷第2条，2017年2月。

T.Albash和D.A.Lidar，“量子退火炉相对于模拟退火的缩放优势的演示”，arXiv：1705.07452v3，2018年8月6日。

J.R.Abo-Shaeer、C.Raman、J.M.Vogels和W.Ketterle，“玻色-爱因斯坦凝聚体中涡旋晶格的观察”，《科学》，第292卷第5516期，第476-479页，2001年4月20日。

Masaki Mutou、Toru Morishita、Daisuke Kushibe、Shinichi Watanabe和MichioMatsuzawa，“二维玻色-爱因斯坦凝聚体中量子化涡旋的数值模拟”，第五届亚洲国际原子与分子物理研讨会，第202-203页，2002年10月。

D.L.Feder、M.S.Pindzola、L.A.Collins、B.I.Schneider和C.W.Clark，“各向异性陷阱中玻色-爱因斯坦凝聚体的暗孤子态”，Phys，Rev.A，62,053606(2000)，2000年10月19日。

Naomichi Hatano和Masuo Suzuki，“求解高阶指数乘积公式”，量子退火和其它优化方法，第37-68页，2005年11月16日。

Mikio Nakahara、Masamitsu Bando和Shu Tanaka，“绝热量子计算解决旅行商问题”，(Kinki大学)科学技术研究所年度报告29，第1-9页，2017年2月28日。

让我们考虑在由N个二进制变量(N是大于或等于1的整数)定义的最小化问题中，通过执行其原理基于量子力学的计算来求解系统的基态(最佳值)和基态的波函数(产生最佳值的状态变量的组合)的情况。这里的系统是指与周围环境分离的子空间。当转换为伊辛模型时，这个问题可以使用量子退火哈密顿量来解决。用于解决这种最小化问题的设备被称为优化设备。

在使用量子退火的探索中，优化设备设置其中N个自旋受到横向场的作用的初始状态。优化设备还将初始状态设置为实验上可制备的量子态。从初始状态开始，优化设备随着时间推移逐渐减小横向场。最后，横向场完全消失，得到了目标哈密顿量的基态，即，最低能态。

然而，在使用量子退火的探索中，对如何减少横向场施加了很强的限制，以避免产生一种称为一阶量子相变的现象。即，优化设备能够通过缓慢地减小横向场来防止一阶量子相变的发生。然而，为了防止一阶量子相变而对横向场强度的这种缓慢降低导致了对最小化问题解决方案的长时间探索。

发明内容

实施方式的一个方面是提供对伊辛模型的基态的有效确定。

根据一个实施方式，提供了一种优化设备，其包括处理装置，该处理装置用于在通过运行在施加到伊辛模型的磁场减小时发生的伊辛模型的状态变化的模拟来求解伊辛模型的基态时，将与噪声对应的值添加到模拟中使用的一些系数，所述伊辛模型表示目标问题，以及执行随着模拟中的时间推进而减小磁场强度的实数时间传播的第一处理和基于虚数时间传播方法减小所述伊辛模型的能量的第二处理。

附图说明

图1示出了根据第一实施方式的示例性优化方法；

图2示出了优化设备的示例性硬件平台；

图3是使用虚数时间传播(ITP)方法的热耗散的示意图；

图4示出了4城市旅行商问题中的示例性城市位置；

图5示出了根据G₂值的数目的计算结果的差异；

图6是能量值相对于横向场强度的示意图；

图7是涉及噪声添加的QTSA计算的示意图；

图8是示出优化设备的功能的框图；

图9是示出示例性解决方案探索过程的流程图；

图10是示出示例性哈密顿量信息读取过程的流程图；

图11是示出示例性自旋信息初始化过程的流程图；

图12是示出示例性量子试行计算过程的流程图；

图13是示出噪声添加到波函数的示例性过程的流程图；

图14是示出使用虚数时间传播方法的量子退火的示例性过程的流程图；

图15是示出用于获得对应经典状态的示例性过程的流程图；

图16示出了对4城市TSP的200个量子试行的示例性结果；

图17示出了自旋的波函数的示例性时间变化；

图18示出了示例性的8城市TSP；

图19示出了对8城市TSP的量子试行的示例性结果；

图20示出了对8城市TSP的波函数的示例性时间传播；

图21示出了在8城市TSP上的量子试行的计算结果与经典伊辛模型MCMC的计算结果之间的示例性比较；

图22示出了根据试行次数的示例性能量转变；

图23示出了QTSA和经典MCMC在ulysses16上获得的解的示例性比较；

图24示出了量子和经典理论方法之间的解改进率的示例性比较；

图25示出了典型量子试行的时间传播；

图26示出了量子理论方法中的能量耗散；

图27示出了关于ulysses22的排序信息的示例；

图28示出了如何在ulysses22上更新解；

图29示出了对于每个能量分辨率度的量子试行的示例性结果；以及

图30示出了在噪声被添加到横向场的情况下通过QTSA的量子退火的示例性结果。

具体实施方式

下面将参照附图描述几个实施方式。这些实施方式可以彼此组合，除非它们具有矛盾的特征。

(a)第一实施方式

第一实施方式针对一种优化方法，该优化方法将量子力学时间传播算子的时间扩展到复数，以引入热耗散效应来求解伊辛模型的基态，并且基于虚数时间传播(ITP)方法实现热耗散。热耗散机制的引入允许弛豫到基态而不严格遵循绝热跃迁条件。然而，注意，使用ITP方法引入热耗散机制并不总是导致获得基态。鉴于此，根据第一实施方式的优化方法通过向伊辛模型添加干扰波函数对称性的噪声来增加达到基态的可能性。

图1示出了根据第一实施方式的示例性优化方法。图1示出了用于实现伊辛模型1的能量优化方法的优化设备10。优化设备10例如通过运行优化程序来实现优化方法。

优化设备10包括存储单元11和处理单元12。存储单元11例如是设置在优化设备10中的存储器或其它存储设备。处理单元12例如是设置在优化设备10中的处理器或运算电路。

存储单元11在其中存储表示要解决的问题的目标问题信息11a。目标问题例如是可转换为伊辛模型1的最小化问题。作为目标问题信息11a，提供关于表示目标问题的伊辛模型1的信息。伊辛模型1由N个二进制变量定义，并且伊辛模型1的最终基态表示目标问题的最优解。

处理单元12使用量子力学的变分原理来求解伊辛模型1的基态。即，处理单元12模拟在施加到代表目标问题的伊辛模型1的磁场减小时发生的伊辛模型1的状态变化。在求解伊辛模型1的基态时，处理单元12将对应于噪声的值加到模拟中使用的一些系数。然后，处理单元12执行随着模拟中的时间流逝而减小施加到伊辛模型1的磁场的强度的实数时间传播处理(第一处理)。在实数时间传播过程中，施加到伊辛模型1的磁场逐渐减小。然后，当磁场达到零时，伊辛模型1处于基态。

处理单元12在实数时间传播过程期间还执行基于ITP方法降低伊辛模型1的能量的热耗散过程(第二处理)。即，处理单元12将模拟时钟上的时间转换为虚数并使虚数时间前进，以求解伊辛模型1的与在给定时间点施加的磁场相对应的基态。伊辛模型1的与磁场对应的基态是伊辛模型1在稳态(稳定基态)下的基态，在该点处的实数时间被视为固定参数。实数时间传播过程行进，然后磁场强度达到零。伊辛模型1在此时的基态是目标基态。

因此，可以通过将噪声添加到初始状态并运行组合实数时间传播和虚数时间传播的计算来有效地求解伊辛模型1的基态。其原因如下所述。

让我们考虑在例如由N个二进制变量(N是自然数)定义的最小化问题上，通过执行其原理基于量子力学的计算来求解系统的基态(最佳值)和基态的波函数(产生最佳值的状态变量的组合)。在该方法中，可以采用量子退火哈密顿量作为哈密顿量。在量子退火中，初始状态中N个自旋受到横向场的作用。例如，一个实验上可制备的量子态被布置为初始状态。从初始状态开始，横向场随时间推移而减小。最后，当横向场完全消失时，得到了所寻求的哈密顿基态，即，最小能量。对于其数值解，通常使用量子蒙特卡罗(QMC)方法。在量子退火计算中，系统保持在基态，磁场减小得足够慢，不会将系统激发脱离基态。这是对计算的一个非常强的约束。此外，在减小磁场的过程中可能发生被称为一阶量子相变的现象，在该现象下，基态和激发态之间的能隙被指数地抑制。在这种情况下，横向场需要以指数方式缓慢地减小。

至于一阶量子相变，最近报告了多体相互作用加速了一阶量子相变(见前述日本特开2018-067200号公报和“非量子哈密顿量对量子退火效率的指数增强”)。

众所周知，量子系统的动态行为可以通过在虚数时间内求解由量子伊辛哈密顿量和描述自旋翻转的横向场项描述的含时薛定谔方程来模拟。该技术对应于量子化模拟退火，其在这里被称为量子热模拟退火(QTSA)方案。QTSA方案直接模拟量子系统的时间传播，因此不区分非平衡状态和平衡状态。

然而，有时会出现当QTSA方案直接应用于量子退火计算时未获得基态的情况。鉴于此，第一实施方式向优化设备10引入以下两个概念。

引入的第一个概念是量子试行(quantum trials)。认为量子退火哈密顿量产生时间传播，即，波函数的初始状态随时间演化到其最终状态，这被认为是一种散射问题。优化设备10处理虚数时间中的波函数的散射问题。具体地，优化设备10将取决于系统的初始状态的包括基态的激发态描述为散射问题，并找到最终状态。

要引入的第二个概念是向系统的初始状态提供多样性的技术。该技术在组合优化问题中主动地使用噪声，噪声通常被认为是门模型量子计算机中的问题。一般来说，噪声是有害的，破坏了系统的信息；然而，在QTSA计算中向初始状态添加噪声提供了散射问题初始值的多样性。这允许在求解散射问题时以有限概率找到基态。

至于如何添加噪声，例如，处理单元12将与噪声对应的值(在下文中有时简称为“噪声值”)添加到伊辛模型1的初始状态，并且将添加了噪声值的伊辛模型1的状态设置为模拟的开始状态(添加了噪声的初始状态)。此时，例如，处理单元12将噪声值添加到伊辛模型1中所包括的多个自旋中的一些的波函数中所包括的系数。注意，处理单元12例如从多个自旋中随机地选择要给予噪声值的自旋。

例如，处理单元12对伊辛模型1执行多次噪声值的相加，以产生伊辛模型1的多个状态，每个状态具有不同的噪声模式。然后，处理单元12使用伊辛模型1的每个生成状态作为模拟的开始状态来执行第一处理和第二处理。即，模拟被执行与添加了噪声的初始状态的数目相同的次数。

因此，将各种噪声模式添加到伊辛模型1允许向系统的初始状态提供多样性。添加了噪声的初始状态中的这种多样性增加了与某些添加了噪声的初始状态对应的最终状态成为基态的机会。结果，即使在使用没有添加噪声的初始状态的模拟中没有获得基态作为最终状态的情况下，也可以达到伊辛模型1的基态。此外，噪声的加入破坏了系统的对称性，这避免了量子退火计算中计算费用的急剧增加。即，容易获得基态，并由此提高了解处理的效率。

通常情况下，获得次优解，而不是正确的基态，在次优解中少量的特定自旋各自具有与基态相反的二进制值0或1。鉴于此，处理单元12还可采用使用经典MCMC的计算。根据量子力学，每个自旋的状态被表示为自旋可能状态(例如，0和1)的量子力学叠加。表示为量子力学叠加的自旋状态表示每个状态被观察到的概率。另一方面，在基于经典力学的伊辛模型(经典伊辛模型)中，每个自旋的状态都是可能的状态之一。

例如，处理单元12基于量子力学对其中每个自旋的状态被表示为自旋的可能状态的量子力学叠加的伊辛模型1重复第一处理和第二处理。接下来，基于在伊辛模型1的能量已经收敛之后伊辛模型1的状态，处理单元12获得其中每个自旋的状态已经被确定为可能状态之一的经典伊辛模型的状态。然后，处理单元12执行使用所获得的经典伊辛模型的状态通过基于经典力学的解(例如，MCMC)来求解经典伊辛模型的基态的处理(第三处理)。

因此，结合基于经典力学的求解方法，可以很容易地达到基态。有时情况是，即使仅第一处理和第二处理不能求解正确的基态，基于经典力学的MCMC(以下有时称为经典MCMC)例如在至多几个步骤中也达到基态。在这种情况下，结合经典的MCMC允许获得基态，因此提高了求解基态的效率。

通过破坏系统的对称性来加速计算的另一种方法是向横向场项的强度添加噪声。系统波函数的对称性也通过向横向场的强度添加噪声而被破坏。例如，处理单元12将对应于噪声的值添加到系数，每个系数表示作用在伊辛模型1中的自旋上的磁场的强度。由此，系统的对称性被打破，这避免了量子退火计算中计算费用的急剧增加。

(b)第二实施方式

第二实施方式针对由N个二进制变量定义的最小化问题，通过使用冯·诺依曼计算机来执行其原理基于量子力学的计算，以获得系统的基态(最佳值)和基态的波函数(产生最佳值的状态变量的组合)。

<优化设备的硬件配置>

图2示出了优化设备的示例性硬件平台。优化设备100具有用于控制其整个操作的处理器101。处理器101经由总线109连接到存储器102和其它各种装置和接口。处理器101可以是单个处理装置或包括两个或更多个处理装置的多处理器系统，诸如中央处理单元(CPU)、微处理单元(MPU)和数字信号处理器(DSP)。也可以将处理器101及其程序的处理功能全部或部分地实现到专用集成电路(ASIC)、可编程逻辑器件(PLD)或其它电子电路中，或其任意组合中。

存储器102用作优化设备100中的主存储装置。具体地，存储器102用于临时存储处理器101执行的操作系统(OS)程序和应用程序中的至少一些、以及处理器101用于其处理将使用的各种类型的数据。例如，存储器102可使用随机存取存储器(RAM)或其它易失性半导体存储器装置来实现。

总线109上的其它装置包括存储装置103、图形处理单元(GPU)104、输入装置接口105、光盘驱动器106、外围装置接口107和网络接口108。

存储装置103在其内部存储介质中或内部存储介质上电地或磁地写入和读取数据。存储装置103用作优化设备100中的辅助存储装置，以存储操作系统和应用的程序和数据文件。例如，可以使用硬盘驱动器(HDD)或固态驱动器(SSD)来实现存储装置103。

联接至监视器21的图形处理单元(GPU)104根据来自处理器101的绘图命令产生视频图像并将视频图像显示在监视器21的屏幕上。监视器21例如可以是有机电致发光(OEL)显示器或液晶显示器。

输入装置接口105连接到诸如键盘22和鼠标23之类的输入装置，并且将来自这些装置的信号提供给处理器101。鼠标23是定点装置，其可以用诸如触摸屏、平板电脑、触摸板和轨迹球之类的其它类型的定点装置来代替。

光盘驱动器106通过使用例如激光读出在光盘24上编码的数据或者将数据写入光盘24。光盘24是便携式存储介质，以通过光反射读取的方式将数据记录在便携式存储介质上。光盘24例如可以是数字通用盘(DVD)、DVD-RAM、光盘只读存储器(CD-ROM)、CD-可记录(CD-R)或CD-可重写(CD-RW)。

外围装置接口107是用于将外围装置连接到优化设备100的通信接口。例如，外围装置接口107可以用于连接存储器装置25和存储卡读取器/写入器26。存储器装置25是具有与外围装置接口107通信的能力的数据存储介质。存储卡读取器/写入器26是用于向存储卡27写入数据或从存储卡27读取数据的适配器，存储卡27是卡形式的数据存储介质。

网络接口108连接到网络20，以便与其它计算机或通信装置(未示出)交换数据。

上述硬件平台可以用于实现根据第二实施方式的优化设备100的处理功能。图2的优化设备100的相同硬件配置可以类似地应用于前述第一实施方式的优化设备10。

优化设备100例如通过执行存储在计算机可读存储介质中的计算机程序来提供第二实施方式的各种处理功能。各种存储介质可用于记录要由优化设备100执行的程序。例如，优化设备100可以在其自身的存储装置103中存储程序文件。处理器101从存储装置103读出这些程序的至少一部分，将其加载到存储器102中，并执行加载的程序。用于程序的其它可能的存储位置包括光盘24、存储器装置25、存储卡27和其它便携式存储介质。存储在这种便携式存储介质中的程序在处理器101的控制下被安装在存储装置103中，使得其准备好在请求时被执行。处理器101也可以执行从便携式存储介质读出的程序代码，而无需将其安装在其本地存储装置中。

通过使用具有上述硬件配置的优化设备100，可以对由N个二进制变量定义的最小化问题执行其原理基于量子力学的计算，以获得系统的基态和激发态的波函数。然而，注意使用量子退火延长解探索的一个因素是如上所述的发生称为一阶量子相变的现象。即，解探索需要长时间来避免发生一阶量子相变。下面详细描述当试图控制一阶量子相变的发生时，最小化问题的解决方案的探索花费很长时间的原因。

注意，伊辛模型的基态是由伊辛模型的哈密顿量描述的能量的最低状态。伊辛模型的基态在下文中有时被称为哈密顿量的基态。另外，在量子退火中施加于伊辛模型的磁场称为横向场。

<基态计算技术概述>

在用量子退火求解的过程中，准备了具有横向场的项作为初始状态。在不进行热耗散过程的情况下，从准备好的初始状态开始，被定义为目标N二进制变量问题的优化问题的哈密顿量中的横向场的强度被充分缓慢地减小以最终达到零。此时，系统的波函数从初始状态缓慢地移向目标伊辛模型的基态的波函数，因此磁场减小，使得系统不跳到激发态。这被称为绝热跃迁。由于绝热跃迁，系统最终演化为目标哈密顿量的基态。以这种方式，可以求解目标量子伊辛哈密顿量的基态。

另一方面，量子退火中的一个瓶颈是绝热跃迁耗时较长。在绝热跃迁期间，基态和激发态之间的能隙可能变得太小。当这种情况发生时，计算时间会增加，因此系统的哈密顿量保持在基态，而不是跳到激发态。即，在绝热跃迁过程中，对于中间哈密顿量，到激发态的跃迁源于基态和第一激发态之间的能隙的闭合。在采用量子退火的情况下，如果基态和第一激发态之间的能隙小，则优化设备100通过减小横向场的减量来使哈密顿量演化得足够慢以保持其保留在基态。

这里注意，能隙作为自旋数目的函数而指数地减小的情况被称为一阶量子相变。在一阶量子相变的情况下，满足根据量子退火的绝热跃迁条件的时间传播的时间步的大小取决于基态和激发态之间的能隙。为此，优化设备100指数地减小时间步的大小，这意味着通过量子退火达到最终状态所花费的时间传播运行的次数的指数增加。

另一方面，存在能隙作为自旋数的幂函数而减小的情况。这种情况被称为二阶量子相变。当发生二阶量子相变时，量子退火中绝热跃迁的速度以幂函数的形式减小。在这种情况下，达到最终状态所花费的时间多项式地增加。

然而，由于许多应用上有趣的问题涉及一阶量子相变，因此计算时间呈指数增长这一事实阻碍了其实际应用。鉴于此，已经有研究探索是否可以用二阶量子相变来代替一阶量子相变。一些研究报告说，在称为非随机哈密顿量的特殊哈密顿量系统中，通过设计横向场项，一阶相变被转移为二阶(见前面提到的“非随机哈密顿量对量子退火效率的指数增强”)。因此，在常规量子退火的框架下，如何避免一阶量子相变是进行计算的中心问题。

此外，在量子退火中，需要满足绝热跃迁条件。然而，有报道称，在使用支持量子退火的量子计算机进行计算时，即使违反了绝热跃迁条件，也可能以有限的概率找到基态(见前述“量子退火器相对于模拟退火的缩放优势的演示”)。报告指出，在使用量子退火处理器的计算中存在最佳求解时间(TTS)；然而，存在最佳TTS的原因从物理角度留出了考虑的空间。

从这些上下文判断，认为在量子退火中，由于某种原因可能发生能量耗散(热耗散)，并允许量子退火的机制起作用。因此，考虑到量子理论框架中可能发生热耗散的可能性是实质性理解量子退火过程的关键。

对于热耗散，与外部环境的相互作用是基本的。因此，需要使用热耗散项将外部环境的影响并入哈密顿量中。在核物理领域，提出了一种处理能量耗散的方法，即在虚数时间内求解含时薛定谔方程。通过引入虚数时间将波动方程转化为扩散方程，这是允许处理能量耗散的主要原因。

处理能量耗散的一个具体示例是设计用于处理由稀原子气体中的玻色-爱因斯坦凝聚(BEC)形成的玻色凝聚体的方法。稀原子气体中BEC的理论从多个方面进行了研究，但量子化涡的计算是理论处理中遇到的障碍之一。众所周知，当光学勺子实验性地搅动玻色凝聚体时，根据经验观察到量子化涡旋(见前面提到的“玻色-爱因斯坦凝聚体中涡旋晶格的观察”)。已知具有量子化旋涡的玻色凝聚体的基态可以很好地近似为一个方程，其中角动量项被并入作为平均场近似结果获得的Gross-Pitaevskii(GP)方程。

然而，在涡旋晶格中寻找GP方程的基态是一项非常困难的任务，因为该方程既不能反映涡旋晶格的对称性，也不能反映涡旋的数量。因此，尽管存在基态，但是存在不知道包括对称特性的信息的情况。然后，用于计算在这一点上涡旋晶格基态的方法是ITP方法。当没有可用的包括关于系统基态的对称性的信息时，使用ITP方法允许计算基态(见前述“二维玻色-爱因斯坦凝聚体中量化涡的数值模拟”)。另外，这种方法也可以用来引入一个唯象的热耗散项，并且有一个处理BEC中暗孤子态的研究(见前面提到的“各向异性陷阱中玻色-爱因斯坦凝聚体的暗孤子态”)。ITP方法是一种简单而强大的算法，但是因为在量子退火系统中没有显式地表示动能算符，所以将ITP方法应用于量子退火系统并不容易。

在ITP方法中，通过用虚数表示时间，量子力学时间传播算子从exp(-iEt)变为exp(-Eτ)(i是虚数单位，E是能量，t是实数时间，τ是虚数时间)。因此，使用完整系统|ψ_n>(第n个本征态)将给定波函数展开为方程(1)，当时间由虚数表示时，将其定义为方程(2)。

|Ψ(t)>是表示时间t处自旋状态的波函数，d_n(t)是|ψ_n>的系数，i是虚数单位，En是哈密顿量的第n个能级，E₀是基态的能级，

是狄拉克常数，τ是虚数时间。在方程式(2)中，当取虚数时间τ→∞的极限时，仅保留基态的幅度，这被描述为下面的方程式(3)。

然而，如果保持原样，基态的振幅根据

减小到零。这是因为通过用虚数表示时间而破坏了统一性(概率守恒)。因此，在对每个计算运行重新归一化波函数的同时继续计算，以便强制地保持统一性。以这种方式，可以达到基态。

该ITP方法基于变分原理，并且无论问题如何都是适用的。由于该方法基于变分原理，理论上保证了该方法能达到最优解。采用ITP方法允许以确定性方式处理量子伊辛系统的时间传播。

图3是使用(ITP)方法的热耗散的示意图。图3示出了曲线图30，其中横轴表示λ，纵轴表示由哈密顿量给出的系统的能量。在图3的示例中，实线31表示当λ随着实数时间演化而从0变化到1时，根据常数λ的基态能量，该常数λ对应于横向场的强度。虚线32表示根据λ的第一激发能态。横向场的强度在λ＝0时最大并且在λ＝1时最小。

在一般的量子退火方法中，如果磁场变化太大，则在点P处处于基态的系统被激发到点Q处的第一激发能态。这种情况在一般的量子退火方法中是不允许的，因此需要绝热缓慢地减小磁场。然而，如果可以通过去激励将能级从点Q降低到点R，则不需要严格遵循绝热跃迁条件。优化设备100采用虚数时间传播以通过去激励将能级从点Q降低到点R。

另一方面，即使使用ITP方法，也会发生一阶量子相变。当发生一阶量子相变时，基态和激发态之间的能隙消失，基态在量子退火的其中横向场下降的中间状态变得几乎简并。如果基态和激发态之间的能级是简并的，则基于它们的能级不能区分这两种状态。这使得对于系统是达到基态还是被提升到激发态无法进行控制。即，即使当使用ITP方法时也不总是获得基态，因此，以固定概率来求解。这是一个不仅来自系统波函数的初始状态，而且来自系统的哈密顿量本身的问题。

因此，在量子退火领域，在QTSA框架内达到基态的技术的发展为基础科学和工程提供了价值。第二实施方式是用于系统地计算基态的技术发展的结果。

<伊辛模型与问题哈密顿量>

接下来，对作为优化设备100的求解对象的伊辛最小化问题(伊辛模型问题)进行说明。伊辛模型问题最初是在物理领域中用于研究磁性物质的模型，并且N个自旋的总能量由下面的方程式(1.1)描述。

右侧的第一项是关于N个状态变量的所有组合的相互排斥且共同穷尽的两个状态变量的值(0或1)与权重系数的乘积的总和。x_i是第i个状态变量；x_j是第j个状态变量；W_ij是表示x_i和x_j之间相互作用的强度的权重系数。右侧的第二项是所有状态变量上的偏置系数b_i与每个状态变量的值的乘积的总和。右侧的第三项是常数C。

最小化问题是用于获得由方程式(1.1)给出的哈密顿量的最小值的问题。已知，当使用例如基于蒙特卡罗方法的计算获得由方程式(1.1)给出的哈密顿量的最小值时，计算工作量通常为2^N量级。即，用于获得最优解的计算需要超过多项式的时间并且指数地增加。

为了解决这个问题，引入了量子理论。由方程式(1.1)描述的哈密顿量被转换成由方程式(1.2)、(1.3)和(1.4)给出的量化哈密顿量。

通过量化方程式(1.1)得到左侧为哈密顿量

(下标TH代表目标哈密顿量)的方程式(1.2)。方程式(1.3)和(1.4)中的σ_i ^z(其中σ上附有^)表示Pauli矩阵作用于第i自旋的z分量。

注意，当描述量子力学的动力学时，不可能仅使用上面的方程式(1.2)、(1.3)和(1.4)来描述自旋翻转。因此，如方程式(1.5)所述的具有横向场项的哈密顿量被定义为量子力学哈密顿量。

方程式(1.5)中的T(附有^)称为横向场项，其被引入来描述含时量子力学方程中的自旋翻转。横向场项被定义为下面的方程式(1.6)。

J是描述自旋到自旋强度的系数。在哈密顿量中，G₁是描述纵向场项的强度的常数，而G₂是描述横向场项的强度的常数。在方程式(1.6)中仅取出关于σ_i ^x(其中σ上附有^)的一阶项的情况是上述“横向伊辛模型中的量子退火”中公开的导出量子退火方法。在文献的量子退火方法中，G₁＝λ，并且G₂＝1-λ，其中λ是表示横向场强度的参数，范围从0到1(0≤λ≤1)。然而，通常，G₁和G₂可以彼此无关地独立地取值。

接下来，问题的算子的矩阵和向量符号由以下方程式定义。

方程式(1.7)的定义不同于通常使用的定义。这是为了允许自旋算符具有1和0的特征值，而不是1和-1。|α_i>和|β_i>具有如下公式所定义的彼此正交关系。

接下来，给出了纵向场σ_i ^z(其中σ上附有^)的本征方程。

横向场σ_i ^x(其中σ上附有^)的本征方程如下。

接下来，表示第i个量子比特的单个粒子波函数的ket由方程式(1.16)定义。

C_iα(t)是表示第i自旋的ket态|α_i>的概率振幅和相位的复数。C_iβ(t)是表示第i自旋的ket态|β_i>的概率振幅和相位的复数。单个粒子波函数的bra由下面的方程式(1.17)定义。

C_iα*(t)和C_iβ*(t)分别是C_iα(t)和C_iβ(t)的复合共轭。在量子力学中，波函数是通常由复数定义的函数，因此其系数也由复数定义。这里，方程式(1.18)在单个粒子波函数已被归一化的假设下成立。

接下来，常用物理量定义如下。

N个粒子系统的波函数定义如下。

观察到的第m自旋为α自旋和β自旋的概率分别用P_m，α和P_m，β表示。P_m，α和P_m，β的定义如下。

P_m，α＝|C_mα|² (1.22)

P_m，β＝|C_mβ|² (1.23)

<ITP法的基本方程式>

让我们考虑求解方程式(2.1)，它是由哈密顿方程式(1.5)定义的含时薛定谔方程式。

方程式(2.1)的第n个特征值E_n及其本征态|ψ_n>是下面的方程式(2.2)的解。

利用完备系统|ψ_n>，将给定的波函数|Ψ(t)>展开如下。

这里，使用下面的方程式(2.4)将时间从实数值时间变换为虚数时间τ。

t＝-iτ (2.4)

通过以虚数时间表示，方程式(2.3)被转换为下面的方程式(2.5)。

当取虚数时间τ→∞的极限时，激发态布居数(populations)的衰减速度比基态的衰减速度快。因此，如方程式(2.6)所描述的，仅保留基态的布居数。

然而，当取虚数时间τ→∞的极限时，系数因子

本身变为0。这是因为将实数值时间转换成虚数破坏了时间传播算子的统一性，并且具体地讲，波函数的归一化积分的值低于1。考虑到这一点，在计算进行期间将归一化积分强制设置为1。这个操作称为重整化。利用虚数时间传播和重整化方法，理论上保证了系统达到基态。

其次，针对ITP方法应用于量子伊辛系统的情况，建立了控制方程式。作为时间传播算子，使用由方程式(2.7)给出的线性ITP算子。

在方程式(2.7)中，Δτ的二阶项被包裹在函数O(Δτ²)中。方程式(2.7)被扩展到关于Δt的一阶项，并且被排序以给出下面的方程式(2.8)。

接下来，发展描述波函数的时间传播的方程式。波函数的时间传播由方程式(2.9)给出。

其中只有第m个粒子的波函数被去除的N-1粒子系统的波函数被表示为|Ψ_m(τ)>(m上附有上划线)，以给出下面的方程式(2.10)。

(

表示张量积，且Ψ的下标m附有上划线)应用于方程式(2.9)的两侧。结果，方程式(2.9)的左侧变为方程式(2.11)。

以类似的方式，

(Ψ的下标m附有上划线)应用于方程式(2.9)的两侧。结果，方程式(2.9)的左侧变为方程式(2.12)。

这里，M_i*表示自旋改变其状态的刚性程度，并且由下面的方程式(2.13)给出。

方程式(2.14)定义如下。

Ω_m(τ)被重新排列至Δτ项，以给出下面的方程式(2.15)。

接下来，将

和

(Ψ的下标m附有上划线)应用于方程式(2.9)的两侧以获得右侧，并且得到方程式(2.16)所定义的控制方程式。

方程式(2.16)中的矩阵元素通过下面的方程式(2.17)至方程式(2.19)获得。

H_1，αα＝G₁(E_TH+(1-x_m)h_z，m)+G₂T_αα (2.17)

H_1，ββ＝G₁(E_TH-x_mh_z，m)+G₂T_ββ (2.18)

H_1，αβ＝H_1，βα＝G₂T_αβ (2.19)

下面给出方程式(2.17)至方程式(2.19)中的参数E_TH、h_z，m、T_αα、T_ββ和T_αβ。

E_TH是目标哈密顿量的能量项，并且由方程式(2.20)至方程式(2.22)定义。

E_TH＝E_Q+E_TL+C (2.20)

h_z，m对应于平均场并且由方程式(2.23)和方程式(2.24)给出。

h_m＝v_m+b_m (2.23)

T_αα是横向场项的期望值，它是通过从

(Ψ的下标m附有上划线，并且T附有^)中取C_mα(τ)的系数而获得的。T_αβ是通过从

(Ψ的下标m附有上划线，并且T附有^)中取出C_mβ(τ)的系数而得到的。T_ββ是通过从

(Ψ的下标m附有上划线，并且T附有^)中取出C_mβ(τ)的系数而得到的。

矩阵形式按以下方式进行评估。这里引入方程式(2.25)作为N粒子波函数，其中只有第m个单粒子波函数用

代替。

定义为方程式(2.26)。

|φ_m>＝C′_mα(τ)|α_m>+C′_mβ(τ)|β_m> (2.26)

现在，基于方程式(2.25)，考虑以下积分。

表达式(2.27)的积分被评价为下面的方程式(2.28)。

以上述方式评估的矩阵形式给出了在方程式(2.16)的表示波函数的时间传播的右侧的两行和两列的矩阵元素。T_αα、T_ββ和T_αβ的矩阵元素如方程式(2.29)和(2.30)中给出。

T_αα＝T_ββ＝JX_m (2.29)

T_αβ＝T_βα＝J (2.30)

注意，矩阵元素T_αα、T_αβ、T_βα和T_ββ是根据积分

(T附有^)来计算的。

矩阵元素的具体描述由方程式(2.31)表示。

这里，如方程式(2.32)所示来设置X_m。

接下来描述Ω_m(τ)。α自旋被观测的概率由方程式(1.22)给出，而β自旋被观测的概率由方程式(1.23)给出。因此，通过将观测概率之和关于虚数时间τ微分，得到方程式(2.33)。

对于观测概率之和，P_m,α+P_m,β＝1被定义为总是真。因此，M_m*是没有实部的纯虚数。

在此考虑对于时间传播控制方程式(2.16)中的波函数的初始值C_mα(τ＝0)和C_mβ(τ＝0)二者选择实数的情况。方程式(2.16)右侧的矩阵元素都是实数。因此，当C_mα和C_mβ是实数时，C_mα(τ+Δτ)Ω_m(τ)也是实数。另一方面，Ω_m(τ)根据M_m*计算，并且在时间τ处的M_m*根据方程式(2.13)也是实数，方程式(2.13)是M_m*的定义方程式。然而，根据方程式(2.13)，M_m*是纯虚数，因此在这种情况下M_m*＝0。因此，Ω_m(τ)＝1并且C_mα(τ+Δτ)也是实数。对于C_mβ(τ+Δτ)也适于完全相同的理论，因此C_mβ(τ+Δτ)也是实数。

接下来，考虑为波函数的初始值C_mα(τ＝0)和C_mβ(τ＝0)选择复数的情况。在这种情况下，M_m*被给出为复数。在这种情况下，通常Ω_m(τ)≠1，并且是一个复数。M_m*是在第m个自旋上确定的值。因此，Ω_m(τ)是基于除第m个感兴趣自旋之外的所有自旋的信息来确定的。因此，C_mα(τ+Δτ)受满足i≠m的所有其它自旋的紧前波函数C_iα(τ)和C_iβ(τ)的影响。即，在相位信息被合并到时间τ＝0的波函数的信息中以使得波函数取复数的情况下，经由虚数分量交换信息。对各自被定义为绝对值的平方的观测概率没有影响的相位分量在自旋之间的相互作用中起着关键作用。然后，自旋通过相互作用以及相位分量耦合到所有其它自旋。

由此，导出了以确定性方式处理量子退火的控制方程式。注意，总能量由方程式(2.34)和(2.35)给出。

E(G₁,C₂)＝G₁E_TH+G₂E_0L (2.34)

至于计算成本，用于虚数时间传播的控制方程式(2.16)的复杂度的阶数是O(2^N)。方程式(2.16)中的矩阵元素的计算涉及目标哈密顿量的总能量的计算，并且总能量计算的复杂度的阶数为O(N²)。因此，假定直到ITP方法收敛为止所需的计算运行的数目为L，那么计算成本为O(LN²)。然而，能量被估计为以exp(ΔE₀₁τ)衰减，然后在虚数时间传播中收敛，从而以指数方式快速下降。

<量子试行的定义>

下文中使用的术语量子试行定义如下。t＝0时的N粒子波函数|Ψ₀>的虚数时间传播算子定义为方程式(3.1)。

此时，获得方程式(3.2)，其中|Ψ(τ)>是作为虚数时间传播的结果而获得的最终状态。

此时，不同的|Ψ₀>通常产生不同的最终状态。术语量子试行是指从初始状态获得最终状态的过程。因此，可以将量子退火视为一种散射过程。

<通过向波函数添加k粒子激发噪声产生初始状态>

存在多种制备初始状态的方法。构成N粒子波函数的单个粒子波函数的初始值如下制备。

这里假设使用方程式(4.1)中的这样的初始值，其保持所有单个粒子波函数的系数相同。即，对于给定的两个自旋i和j(i≠j)，准备确保C_iα(t＝0)＝C_jα(t＝0)且C_iβ(t＝0)＝C_jβ(t＝0)的初始状态。此时，k粒子激励噪声定义如下。在N个自旋中，选择k个自旋，对于每个自旋，将z％的噪声添加到其概率幅度。对于所选择的第一自旋，添加噪声之后的概率幅度由方程式(4.2)和(4.3)定义。

C'_iα(t＝0)＝(1+z)C_ia(t＝0) (4.2)

当在方程式(4.2)中使用C_iα＝0的初始状态时，不添加噪声。因此，在方程式(4.2)和(4.3)的情况下，噪声被添加到C_iβ。

然而，应当注意，由方程式(4.2)和(4.3)表示的用于定义噪声添加的方程式仅仅是示例。重要的是，通过添加噪声，只有某些自旋具有与其它量子自旋不同的单粒子量子态。另外，确保z通常对于每个自旋取不同值是合适的。这是因为，尽管噪声被添加到k个自旋，但是相同的噪声被添加到所有k个自旋是不合逻辑的。因此，优化设备100每次使用例如随机数来确定噪声幅度z。

当根据方程式(4.2)和(4.3)将噪声添加到k个自旋时，概率幅度增大或减小，这可被视为不同状态的部分引入。这可以解释为激发的部分结合。因此，添加的噪声在这里被称为k粒子激励噪声。注意，可以从噪声添加之前的最开始准备具有破坏对称性的初始状态，其中C_iα(t＝0)≠C_jα(t＝0)且C_iβ(t＝0)≠C_jβ(t＝0)。例如，可以将一个自旋指定为C_iα(t＝0)＝1，而将剩余自旋指定为C_jα(t＝0)＝0(i≠j)。在这种情况下，通过噪声添加获得不同的最终状态。

<使用量子试行和噪声添加来解决的挑战>

到目前为止的讨论还没有为量子试行和波函数的噪声添加提供依据。接下来描述的是一种机制，其中这两种方案的引入解决了仅通过ITP方法难以克服的问题，所述问题为即使在ITP方法的框架内允许量子系统随时间演化，ITP方法也不总是能够找到基态。

让我们提出使用ITP方法实际解决4城市旅行商问题(TSP)时发生的具体问题。

图4示出了4城市TSP中的示例性城市位置。在图4的示例中，第一个城市的位置是(0,0)；第二个城市是(3,0)；第三个城市是(3,1)；第四个城市是(1,1)。

对于图4所示的问题，有三个候选的解。候选路径长度依次为7.414、7.812和10.398。

让我们考虑将TSP转换为伊辛模型并在外部磁场中通过使用ITP方法进行量子退火计算的情况，如上述“横向伊辛模型中的量子退火”中所讨论的那样。设G₁＝1–G₂，其中0≤G₂≤1，并且在虚数时间τ＝0时G₂＝1。G₂被给出为离散点，并且以规则的间隔提供划分点。这里，用于计算的分割点n的数目从64增加到1024。第n个G₂值由G₂(n)表示。中间状态由方程式(5.1)定义。

G₂(n)＝1-0.01n (5.1)

对于由方程式(5.1)定义的每个中间状态，执行虚数时间传播过程并且继续计算直到能量已经收敛。设计算收敛的阈值为10^-10。波函数的初始值在这里被设置为C_mα(τ＝0)＝1且C_mβ(τ＝0)＝0，尽管它们将取任何给定值。还设G₂(0)＝1且G₁＝1.0。一旦n＝0的计算结束，n的值增加1，并且重复虚数时间传播过程以找到在具有G₂(1)且G₁(1)的下一个中间点处的基态。这里，G₂(1)的值已经从先前点的值减小。注意，作为新中间点处的波函数的初始状态，使用作为前一点的最终结果而获得的波函数，其中能量以G₂(0)收敛。重复该计算，并且最终当G₂(n＝100)＝0的计算结束时，获得目标问题哈密顿量的基态(或激发态)。这一系列的计算是一次量子试行。

图5示出了根据G₂值的数目的计算结果的差异。图5示出了曲线图33，横轴表示G₂，纵轴表示能量。曲线图33中的三条实线表示与G₂相关的能量值，每条实线对应于具有不同数目的划分点的量子试行。如曲线图33所示，根据划分点的数目，求解具有能量7.414的基态或者求解具有能量7.812或10.398的激发态。

在4城市TSP的情况下，通过减少截断误差和增加划分点的数量来确定基态。然而，例如ulysses16不是这种情况，ulysses16是TSPLIB(TSP工作簿)中包括的更实际的问题。在这样的问题上，即使当以双精度浮点表示能量截断误差并且提供相当大数量的划分点时，也不能达到基态。

因此，尽管存在未达到基态的可能性，但是获得了要在4城市TSP中找到的所有三个候选。该计算不一定遵循绝热跃迁条件，而是支持以固定概率求解基态的上述实验量子退火结果。因此，解释了在实际处理器的量子退火计算中，发生由于热耗散而产生的自然能量耗散。

接下来描述引起这种现象的机理。

图6是能量值相对于横向场的强度的示意图。在图6的曲线图34中，水平轴表示横向场的强度，垂直轴表示能量。曲线图34描绘了典型的一阶量子相变。一阶量子相变是在量子退火期间发生的现象，其中随着横向场的减小，目标哈密顿量的影响变得更强，然后多个量子态实际上变得简并。在定性意义上，能隙像ΔE≈exp(-bN)那样指数地闭合，其中N是自旋数，b是系数。

量子退火是用于在保持由曲线图34中的最底部电位曲线表示的量子态的同时通过减小横向场来求解目标基态的超启发式方法。然而，随着量子比特(自旋)数量的增加，能隙呈指数级关闭。结果，发生能量耗散，并且即使存在与外部环境的相互作用，能隙实际上也可以达到0。双精度浮点数通常用于数值计算，但它们只允许16位数字的精度。因此，在第17位或随后的数字中出现差异的情况下，双精度浮点数的使用不足以表示差异。也就是说，所有的差异都是数值退化的。

第二实施方式从不同角度观察计算分辨率内的能量退化，并认为所有退化状态都是可能的。即，我们认为以有限概率找到基态而不是以100％的概率找到基态是可接受的。在这种思路下，也有可能求解激发态。然而，因为相同的计算过程将产生相同的微观状态，所以故意引入小噪声以在初始状态中产生变化。作为QTSA控制方程式的方程式(2.16)是一种平均场近似方程式，其包括其在势项中在当前时间的自身波函数，因此表现为典型的复杂系统。因此，添加的非常小的噪声在时间传播期间被放大，并引起系统的意外行为。系统随机地达到基态和激发态。根据第二实施方式的优化设备100引入量子试行和噪声，使得系统随机地达到目标基态和激发态。

图7是涉及噪声添加的QTSA计算的示意图。例如，优化设备100首先将预定初始值设置为表示自旋状态的自旋信息40。在自旋信息40中，每个自旋的取向被设置为例如0或1。优化设备100向自旋信息40添加噪声以生成多个自旋信息块41、42、43等。在每个自旋信息块41、42、43等中，一些随机选择的自旋已经改变为除0和1之外的状态(基于量子力学的概率表示)。

优化设备100使用自旋信息块41、42、43等中的每一个作为初始状态来执行QTSA计算。从对每个自旋信息块41、42、43等进行的QTSA计算获得的结果将是基态或激发态。因此，通过生成足够数量的自旋信息块41、42、43等作为初始状态，期望任何计算结果来指示基态。

量子试行和噪声添加的使用也有望缩短计算时间。其原因将在下面描述。

基于方程式(2.5)，在虚数时间传播过程中，激发态的衰减率需要满足表达式(5.2)。

这里，对于哈密顿方程式(1.5)定义能量E_n，该哈密顿方程包括横向场项并且定义个体G₁和G₂。这由E_n(G₁,G₂)表示。注意，具有最慢衰减速率的项限制了计算速率。这是第一激发态。因此，该计算大约花费方程式(5.3)所示的时间量。

这里，从E₀到E₁的变化ΔE₁由方程式(5.4)给出。

ΔE₁＝E₁(G₁,G₂)-E₀(G₁,G₂) (5.4)

如果这里发生一阶量子相变，则存在近似等于由方程式(5.5)指示的能量的能隙。

ΔE₁＝O(exp(-αN)) (5.5)

因此，计算花费由方程式(5.6)指示的时间量。

这表明计算时间实际上呈指数增长。这种计算时间由不确定性原理确定，因此难以避免时间的指数增长。

然而，一旦接受实际上正在发生退化，计算复杂性的增长可以被控制在能量截断误差的范围内。代替地，能量截断误差范围内的所有量子态都变成候选态，与整个解空间的大小相比，候选态仍然很小。这是因为实际上并不是所有的量子态都是简并的。

如上所述，优化设备100通过向每个初始自旋状态添加噪声并使用QTSA计算最小化问题的解来有效地求解。

<在优化设备上实现具有噪声添加的QTSA>

接下来，提供关于优化设备100的功能的具体描述，优化设备100向Q自旋的初始状态添加噪声并使用QTSA方案解决最小化问题。

图8是示出优化设备的功能的框图。优化设备100包括存储单元110和处理单元120。存储单元110使用例如在存储器102或存储装置103中确保的存储区域的一部分来实现。处理单元120例如通过执行预定程序的处理器101来实现。

存储单元110在其中存储能量信息、自旋信息、横向场信息、问题设置信息和哈密顿量信息。

能量信息包括计算的初始能量值、以及在迄今为止计算的能量值中的从最小值开始以升序的预定数量的能量值。能量信息还可以包括与预定数量的能量值中的每一个相对应的各个状态变量的值的组合。自旋信息包括各个状态变量的值和关于自旋初始化方法的规范的信息。横向场信息包括关于作用在系统上的横向场的强度的信息。问题设置信息例如包括表示解决方案目标问题的信息、使用的优化方法(诸如SA和QTSA方案)、用于执行优化方法的各种参数、以及计算结束条件。哈密顿量信息包括例如由方程式(1.1)给出的能量函数的权重系数(W_ij)、偏置系数(b_i)和常数(C)。

处理单元120包括控制单元121、设置读取单元122、自旋初始化单元123、哈密顿量生成单元124、横向场计算单元125、ITP控制单元126、噪声添加单元127和经典状态渲染单元128。

控制单元121对处理单元120的每个单元实施控制，以执行最小值求解处理。

设置读取单元122将上述各种信息从存储单元110读取为控制单元121能够理解的形式。

自旋初始化单元123将包括在伊辛模型中的多个自旋的状态初始化。自旋初始化单元123在存储单元110中存储自旋信息，该自旋信息包括表示初始化的自旋状态的状态变量的值。

哈密顿量生成单元124基于问题设置信息生成哈密顿量信息，该哈密顿量信息表示用于解决目标问题的伊辛模型的哈密顿量。哈密顿量生成单元124将生成的哈密顿量信息存储在存储单元110中。

横向场计算单元125计算随着实数时间的进展而更新的横向场。

ITP控制单元126控制虚数时间传播。例如，每次在固定横向场的强度的情况下计算自旋状态时，ITP控制单元126将虚数时间向前移动每个虚数时间步宽度Δτ。

噪声添加单元127将噪声添加到由自旋初始化单元123初始化的自旋状态。例如，噪声添加单元127从多个自旋中随机选择预定数目的自旋。噪声添加单元127将噪声添加到每个选择的自旋的概率幅度。噪声添加单元127例如使用随机数来确定每个自旋的噪声量。

经典状态渲染单元128从通过量子试行获得的一组解中找到对应于经典伊辛模型的状态，并将该状态设置为经典MCMC过程的初始状态。然后，经典状态渲染单元128通过经典MCMC处理获得基态。

注意，图8的实线互连功能块表示它们的一些通信路径。本领域技术人员将理解，在实际实现中可以存在其它通信路径。图8所示的每个功能块可以被实现为程序模块，使得计算机执行程序模块以提供其编码功能。

接下来描述由优化设备100运行以使用QTSA方案搜索最小化问题的解的过程。

图9是示出示例性解探索过程的流程图。图9中的过程在下面以步骤编号的顺序描述。

[步骤S101]控制单元121执行从存储单元110读取哈密顿量信息的处理。该处理将在后面详细描述(参见图10)。

[步骤S102]自旋初始化单元123使自旋信息初始化。例如，自旋初始化单元123从存储单元110读取指示由用户指定的自旋初始化方法的指定信息。存在选择自旋初始化方法的多个模式，包括：指定模式、键合模式、反键合模式和随机模式。指定模式是用于按用户指定的设置自旋状态的初始化方法。键合模式是一种为每个自旋的状态|1>和状态|0>的系数分配相同的符号的初始化方法。反键合模式是一种为每个自旋的状态|1>和状态|0>的系数分配相反的符号的初始化方法。如果用户在指定模式、键合模式和反键合模式中指定了自旋初始化方法，则自旋初始化单元123通过指定的初始化方法来初始化自旋信息。随机模式是一种随机设置每个自旋以指向上或指向下的初始化方法。在随机模式中，通过将向上和向下的概率分别设置为50％，自旋初始化单元123准备相等数量的向上和向下自旋。稍后将详细描述自旋信息初始化处理(参见图11)。

[步骤S103]控制单元121重复步骤S104和S105，直到满足计算结束条件为止。例如，控制单元121将迭代计算运行的预定次数设置为tmax。然后，控制单元121每执行一次步骤S104和S105，将表示迭代次数的变量tstep从1开始增加计数。当tstep的值小于或等于tmax时，控制单元121重复迭代处理步骤S104和S105。

[步骤S104]横向场计算单元125确定横向场项的系数G₂的值。例如，横向场计算单元125将G₂的初始值设置为1，并且在步骤S104和105的每次迭代之后将G₂的值减小预定量。当tstep的值已经达到tmax时，横向场计算单元125允许G₂的值为0。

[步骤S105]ITP控制单元126、噪声添加单元127和经典状态渲染单元128相互协作执行量子试行计算。通过量子试行计算，得到了电流横向场作用下哈密顿量的基态。关于量子试行计算的详细说明将在后面描述(参见图12)。

[步骤S106]当满足迭代处理的计算结束条件时，控制单元121移动到步骤S107。例如，如果在步骤S105结束时tstep的值等于tmax，则控制单元121确定满足计算结束条件。

[步骤S107]控制单元121输出最终状态。例如，控制单元121输出当横向场的系数G₂的值已经达到0时获得的哈密顿量的基态(最佳值)以及基态中的状态变量的组合。

通过随着虚数时间传播过程的重复而逐渐减小横向场的强度，因此可以找到产生系统最优值的状态变量的组合。得到的状态变量的组合是最小化问题的解。

接下来描述读取哈密顿量信息的过程的进一步细节。

图10是示出示例性哈密顿量信息读取过程的流程图。图10中的过程在下面以步骤编号的顺序描述。

[步骤S201]控制单元121从存储单元110读取目标哈密顿量的系数W_ji和b_i以及常数C。此时，控制单元121还从存储单元110读取在QTSA计算中使用的各种参数。例如，如果预先指定了虚数时间步宽度Δτ，则控制单元121读取该虚数时间步宽度Δτ。注意，可以采用可变时间步宽度作为虚数时间步宽度Δτ。在这种情况下，通过预定的计算获得Δτ，并且因此省略对虚数时间步宽度Δτ的读取。

[步骤S202]自旋初始化单元123读取关于自旋初始化方法的指定信息。

[步骤S203]自旋初始化单元123确定自旋初始化方法是否是指定模式。如果是指定模式，则自旋初始化单元123移动到步骤S204，否则，结束哈密顿量信息读取处理。

[步骤S204]自旋初始化单元123从存储单元110读取由用户指定的自旋值。随后，自旋初始化单元123结束哈密顿量信息读取处理。

接下来描述初始化自旋信息的过程的进一步细节。

图11是示出示例性自旋信息初始化过程的流程图。图11中的过程在下面以步骤编号的顺序描述。

[步骤S301]自旋初始化单元123确定自旋初始化方法是否是指定模式。如果是指定模式，则自旋初始化单元123移动到步骤S302，否则，移动到步骤S303。

[步骤S302]自旋初始化单元123将自旋状态初始化为用户指定的值。例如，自旋初始化单元123基于从存储单元110读取的关于自旋状态的信息来初始化自旋状态。随后，自旋初始化单元123结束自旋信息初始化处理。

[步骤S303]自旋初始化单元123确定自旋初始化方法是否为键合模式。如果是键合模式，则自旋初始化单元123移动到步骤S304，否则，移动到步骤S305。

[步骤S304]自旋初始化单元123初始化自旋状态以进入键合模式。例如，自旋初始化单元123将所有自旋的状态初始化为

随后，自旋初始化单元123结束自旋信息初始化处理。

[步骤S305]自旋初始化单元123确定自旋初始化方法是否为反键合模式。如果是反键合模式，则自旋初始化单元123移动到步骤S306，否则，移动到步骤S307。

[步骤S306]自旋初始化单元123初始化自旋状态以进入反键合模式。例如，自旋初始化单元123将所有自旋的状态初始化为

随后，自旋初始化单元123结束自旋信息初始化处理。

[步骤S307]自旋初始化单元123初始化自旋状态以进入随机模式。例如，自旋初始化单元123将从所有自旋中以50％的概率选择的自旋的状态初始化为

而将其余自旋的状态初始化为

随后，自旋初始化单元123结束自旋信息初始化处理。

以上述方式，自旋信息被初始化。在量子试行计算的早期阶段，基于初始化的自旋信息，将噪声添加到每个量子比特(自旋)集的波函数的一部分。接下来详细描述涉及噪声添加的量子试行计算的过程。

图12是示出示例性量子试行计算过程的流程图。图12中的过程在下面以步骤编号的顺序描述。

[步骤S401]ITP控制单元126设置包括量子试行计数N_q的计算参数。例如，ITP控制单元126将预先指定的计算参数的值设置为对应的变量。

[步骤S402]ITP控制单元126重复步骤S403至S409 N_q次。

[步骤S403]ITP控制单元126为每个自旋设置波函数的初始状态。例如，ITP控制单元126针对每个自旋生成与由自旋信息初始化处理初始化的值相对应的方程式(4.1)的波函数。

[步骤S404]噪声添加单元127对波函数执行噪声添加处理。稍后将详细描述噪声添加处理(参见图13)。

[步骤S405]ITP控制单元126利用虚数时间传播执行量子退火。后面将详细描述使用虚数时间传播的量子退火处理(参见图14)。

[步骤S406]经典状态渲染单元128计算对应于当前自旋状态的经典状态。在计算经典状态的过程中，如果不存在经典状态，则输出“亚稳态”作为处理结果。计算经典状态的过程将在后面详细描述(参见图15)。

[步骤S407]经典状态渲染单元128确定是否存在对应于当前自旋状态的经典状态。例如，如果获得“亚稳态”作为处理结果，则经典状态渲染单元128确定不存在经典状态，否则，确定存在经典状态。如果存在对应的经典状态，则经典状态渲染单元128移动到步骤S408，否则，移动到步骤S409。

[步骤S408]经典状态渲染单元128使用在步骤S406中获得的经典状态作为初始状态，通过经典MCMC方法获得基态。例如，经典状态渲染单元128通过使用随机数随机地确定伊辛模型的经典状态中的哪个自旋要被翻转以及是否接受自旋的翻转。如果已经确定接受自旋翻转，则经典状态渲染单元128将伊辛模型的当前状态改变为自旋翻转的状态。通过重复这种状态转换，经典状态渲染单元128获得具有最小能量的状态。然后，经典状态渲染单元128保存由经典MCMC获得的能量值和自旋状态作为当前处理循环的处理结果。随后，经典状态渲染单元128移动到步骤S410。注意，使用经典MCMC的解决方案可以与上述“组合优化问题的加速器结构”中描述的设备合作进行。

[步骤S409]经典状态渲染单元128输出表明不存在对应的经典状态的处理结果。

[步骤S410]在步骤S403至S409的N_q次迭代完成时，ITP控制单元126结束量子试行计算处理。

[步骤S411]ITP控制单元126输出在步骤S408中保存的所有能量值中的最小能量值以及产生最小能量值的自旋状态作为量子试行计算的结果。

以上述方式，多次执行采用其中添加了噪声的波函数被用作初始状态的虚数时间传播的量子退火，然后获得产生最小能量值的自旋状态。利用随机数随机地实现了对波函数的噪声添加。

图13是示出噪声添加到波函数的示例性过程的流程图。图13中的过程在下面以步骤编号的顺序描述。

[步骤S501]噪声添加单元127设置添加噪声的自旋数M。例如，噪声添加单元127将预先指定的值设置为指示用于噪声添加的自旋数M的变量。

[步骤S502]噪声添加单元127设置噪声的最大幅度z_max[％]。例如，噪声添加单元127将预先指定的值设置为指示噪声z的幅度的变量。

[步骤S503]噪声添加单元127重复步骤S504至S507 M次。

[步骤S504]噪声添加单元127通过使用随机数产生1与N之间的一个整数j。

[步骤S505]噪声添加单元127确定在步骤S504中当前生成的随机数是否与之前生成的随机数重叠。如果存在重复随机数，则噪声添加单元127移动到步骤S504，否则，移动到步骤S506。

[步骤S506]噪声添加单元127通过使用随机数产生0到z_max[％]之间的实数值。然后，噪声添加单元127将产生的实数值设置为z[％]。

[步骤S507]噪声添加单元127将z[％]的噪声添加到第j个自旋的波函数的概率幅度。

[步骤S508]当完成步骤S504至S507的M次迭代时，噪声添加单元127结束向波函数添加噪声的处理。

以上述方式，随机选择的噪声幅度被添加到M个自旋的波函数，然后使用虚数时间传播过程进行量子退火。注意，相同幅度的噪声可以应用于所有M个自旋。在该情况下，例如，噪声添加单元127将z_max[％]设置为z[％]，而不在步骤S506中基于随机数产生实数值。

图14是示出使用虚数时间传播的量子退火的过程的流程图。图14中的过程在下面以步骤编号的顺序描述。在下面的描述中，用符号k代替表示在上述方程中使用的自旋数的符号i。

[步骤S601]ITP控制单元126基于方程式(1.19)计算所有自旋中的每一个的纵向场算子的期望值x_k(k＝1、2、...、N)。

[步骤S602]ITP控制单元126基于方程式(1.20)计算所有自旋中的每一个的横向场算子的期望值y_k(k＝1、2、...、N)。

[步骤S603]ITP控制单元126基于方程式(2.23)和(2.24)计算所有自旋中的每一个的纵向平均场h_k(k＝1、2、...、N)。

[步骤S604]ITP控制单元126基于方程式(2.35)计算横向场能E_OL。

[步骤S605]ITP控制单元126重复步骤S606至S613，直到满足结束条件为止。例如，ITP控制单元126将迭代计算的预定次数设置为nitl。然后，控制单元121每运行一次步骤S606至S609，将表示虚数时间传播的迭代次数的变量i从1开始依次向上加1。当i的值小于或等于nitl时，控制单元121重复迭代处理步骤S606至S609。

[步骤S606]ITP控制单元126基于方程式(2.16)计算下一个时间步(即，通过将虚数时间向前移动Δτ获得的时间点)的C_kα(τ₀+Δτ)和C_kβ(τ₀+Δτ)。

[步骤S607]ITP控制单元126基于下面的方程式(6.1)执行重整化，以满足归一化条件。

γ_m是归一化常数。将重整化之前的波函数的系数C_mα(t+Δt)和C_mβ(t+Δt)分别乘以归一化常数γ_m，并将这些乘积用作重整化之后的波函数的系数。

[步骤S608]ITP控制单元126基于使用重整化系数的方程式(1.19)计算所有自旋中的每一个的纵向场算子的期望值x_k(k＝1、2、...、N)。

[步骤S609]ITP控制单元126基于使用重整化系数的方程式(1.20)计算所有自旋中的每一个的横向场算子的期望值y_k(k＝1、2、...、N)。

[步骤S610]ITP控制单元126基于方程式(2.23)和(2.24)计算所有自旋中的每一个的纵向平均场h_k(k＝1、2、...、N)。

[步骤S611]ITP控制单元126根据方程式(2.35)计算横向场能E_OL。

[步骤S612]ITP控制单元126计算基于方程式(2.34)的总能量E(G₁,G₂)与前一时间步处的总能量之间的差ΔE。

[步骤S613]ITP控制单元126判断当前时间步与前一时间步之间的总能量差ΔE是否小于或等于计算的截断精度ε。如果ΔE小于或等于ε，则ITP控制单元126结束虚数时间传播处理，否则，移动到步骤S614。

[步骤S614]如果满足迭代过程的计算结束条件，则ITP控制单元126结束虚数时间传播过程。例如，当在步骤S613结束时i的值等于nitl时，ITP控制单元126确定已经满足计算结束条件。

在利用虚数时间传播的量子退火完成之后，执行获得相应经典状态的过程。

图15是示出用于获得对应经典状态的示例性过程的流程图。图15中的过程在下面以步骤编号的顺序描述。

[步骤S701]经典状态渲染单元128设置阈值ε。例如，经典状态渲染单元128将预先指定的值设置为表示阈值ε的变量。

[步骤S702]经典状态渲染单元128在从1开始对变量i进行递增计数的同时，重复步骤S703至S709 N次(N＝自旋数)。

[步骤S703]经典状态渲染单元128将表示概率幅度(参见方程式(4.1))的C_iα和C_iβ的绝对值之间的幅度关系进行比较。如果|C_iα|>|C_iβ|，则经典状态渲染单元128移动到步骤S704，否则移动到步骤S707。

[步骤S704]经典状态渲染单元128确定|C_iα|是否大于1-ε。如果|C_iα|>1–ε，则经典状态渲染单元128移动到步骤S706，否则，移动到步骤S705。

[步骤S705]经典状态渲染单元128输出“亚稳态”作为处理结果，并结束获得相应经典状态的过程。

[步骤S706]经典状态渲染单元128将第i个自旋带入状态1。随后，经典状态渲染单元128移动到步骤S710。

[步骤S707]经典状态渲染单元128确定|C_iβ|是否大于1–ε。如果|C_iβ|>1–ε，则经典状态渲染单元128移动到步骤S709，否则，移动到步骤S708。

[步骤S708]经典状态渲染单元128输出“亚稳态”作为处理结果，并结束获得相应经典状态的过程。

[步骤S709]经典状态渲染单元128将第i个自旋带入状态0。

[步骤S710]当完成步骤S703至S709的N次迭代时，经典状态渲染单元128结束获得对应经典状态的过程。

通过获得经典状态，即使仅通过量子试行没有找到基态，也能够找到波动方程式的基态。这减少了由于使用虚数时间传播过程而不能获得基态的机会。

<组合优化问题的计算示例>

接下来描述实际组合优化问题的实际计算示例。注意，在下面的计算示例中，使用独立的冯·诺依曼计算机。该计算是在实际的时间帧内实现的，既没有并行化的实现，也没有使用诸如图形处理单元的加速器。

<<添加了噪声的初始状态上的量子试行引起的坍缩现象的实例>>

让我们来观察通过优化设备100运行200个量子试行而获得的上述4城市TSP的计算结果。赋予优化设备100的计算条件如下。

优化设备100随机选择k个自旋(k＝1、2等)用于噪声添加。另外，优化设备100解决了G₁固定为1.0而G₂从1改变为0的问题。G₂的划分点的数目为n₀＝128。G₂(n)由下面的方程式(7.1)给出。

优化设备100将QTSA(虚数时间传播)的截断误差δE设置为10^-8(δE＝10^-8)。另外，对于给定的G₂(n)，优化设备100继续QTSA计算，直到计算达到或低于能量收敛过程的截断误差。一旦能量已经收敛，优化设备100将n的值增加1，然后再次重复QTSA计算。最后，得到在G₂＝0时产生的值。

对于波函数的初始值，优化设备100随机地选择16个量子比特(自旋)中的任意两个，并且将高达1％的噪声随机地添加到每个所选比特的概率幅度C_mα(0)。图16中示出了在重复这样的量子试行约200次之后获得的结果。

图16示出了对4城市TSP的200个量子试行的示例性结果。在图16中，左侧的曲线图51表示与横向场的强度相关的能量转变。在曲线图51中，水平轴表示横向场的强度，而垂直轴表示能量。右侧的曲线图52表示沿着虚数时间路径的能量转变。在曲线图52中，水平轴表示虚数时间，而垂直轴表示能量。

在曲线图51和52中，重叠的实线对应于单独的200个量子试行。根据曲线图51，所有试行中的每一个达到作为4城市TSP的解的能量状态(即，7.414、7.812和10.398)中的一个。结果表明，基态不是100％的时间获得的，而是以有限的概率获得的。还可以看出，尽管每个试行都经历了不同的过程来达到其解，但许多试行都经历了高能态。

这种高能量子态在曲线图52中是明显的。一个值得注意的趋势是，能量在有限的时间长度内保持在较高水平，然后坍缩成TSP解水平。这些较高的能级被称为亚稳态。

当虚数时间传播继续时，一些试行直接达到其最终状态，而另一些试行则经历亚稳态。亚稳态并不总是通常只使用1和0来描述的经典状态。

图17示出了自旋的波函数的示例性时间变化。左侧的曲线图53表示当能量到达基态而不经过亚稳态时16个自旋的波函数的时间变化。右侧的曲线图54表示当能量在经过亚稳态之后到达基态时16个自旋的波函数的时间变化。在曲线图53和54中，水平轴表示虚数时间，垂直轴表示波函数绝对值的平方。

在不经历亚稳态的情况下，所有的自旋都在虚数时间传播的早期翻转。另一方面，在经历亚稳态的情况下，16个自旋中有12个在虚数时间传播的早期翻转；然而，剩下的4个自旋是无序的并且保持在中间状态。不存在与这些中间状态相对应的经典状态，这些中间状态是量子力学所独有的。然后，在虚数时间τ≈6处，在由剩余4个自旋组成的两对中，一个自旋翻转到状态0，另一个自旋翻转到状态1，它们分别达到其基态。这是一种称为分叉(bifurcation)的现象。

注意，随着系统中量子比特数的增加，亚稳态的持续时间往往会延长。这表明存在深的势垒。由于这种亚稳态持续时间较长且稳定，因此在能量塌陷到较低能级之前，即使有热耗散效应，计算也可能结束。在这种情况下，不存在对应于亚稳态的经典状态，因此不可能获得对应的经典伊辛模型的解。因此，在提取经典伊辛模型的解时，优化设备100从解候选中剔除其最终结果指示亚稳态的解。然而，请注意，这并不是说亚稳态没有意义。这只是一个无法将这种亚稳态映射到经典对应状态的问题。这个问题源于量子理论方法本质上是经典的扩展。

<<噪声添加破坏对称性的实例>>

接下来给出的是一个实际情况表明可以通过向波函数的初始值添加噪声来故意破坏自旋波函数的对称性。对于64个量子比特的系统大小，我们使用8城市TSP作为目标问题。

图18示出了示例性的8城市TSP。在图18中，将每个城市的位置(x坐标和y坐标)映射到该城市的识别号。问题的最优解之一指示要访问的城市的序列为城市1、2、3、4、5、6、7和8。此时的最佳值为2.5023。

图19示出了对8城市TSP的量子试行的示例性结果。在图19的曲线图55中，水平轴表示横向场的强度，垂直轴表示能量。曲线图55中的线表示与量子试行中的横向场的强度相关的能量转变，每个试行在具有不同噪声模式的多个初始状态之一上运行。

根据方程式(3.2)，初始状态通过散射演化为其最终状态；然而，由于添加了不同的噪声模式，这里的初始状态都彼此不同，因此它们演化为不同的最终状态。结果，获得了包括基态在内的各种量子态。可以看出，即使添加到初始状态的噪声很小，作为量子退火的结果，也不总是获得基态。

如果寻求获得具有100％概率的基态，则这是不期望的。然而，利用噪声的影响，可以通过故意添加噪声来计算包括基态在内的量子态。换句话说，基态不是以100％的概率而是以有限的概率计算的。

接下来描述如何通过噪声添加来破坏波函数的对称性。

图20示出了8城市TSP中的波函数的示例性时间传播。图20描述了对于基态和第一激发态的64个量子比特的波函数的时间传播。左侧曲线图56表示基态的时间传播，右侧曲线图57表示第一激发态的时间传播。在曲线图56和57中，64个的|1>的时间传播用实线表示，64个的|0>的时间传播用虚线表示。曲线图56和57中的水平轴和垂直轴分别表示虚数时间和每个波函数的绝对值的平方。

通过比较基态和第一激发态可以看出，噪声添加在系统的时间传播的早期阶段没有产生太大的差异。然而，随着时间的推移，噪声的影响越来越大。然后，最终放大的噪声影响整个系统的行为，即，某个旋转翻转。因此，将轻微噪声添加到波函数的初始值以确定性方式避免系统演变为基态。

接下来描述的是自旋如何翻转。在时间传播的早期阶段，所有的自旋都开始以协调的方式与其它自旋保持同步地翻转。然而，当只有一个自旋开始以比其它自旋更快的速度翻转时，整个系统的对称性就被打破了。然后，自旋的翻转触发其它自旋的翻转，以这种方式，整个系统的自旋状态就被确定了。

在对称性被破坏之后再次添加噪声是没有用的。一旦翻转开始，除非采取措施允许自旋完全翻转，否则系统对于添加一些噪声是稳定的。从物理角度来看，这是很重要的，因为波函数的对称性会影响量子退火的计算稳定性和成败。就在自旋信息初始化之后，所有自旋的波函数都处于相同的状态，因此高度对称，但很容易被破坏。另一方面，一旦对称性被破坏，一些噪声的添加不足以破坏波函数的对称性而使系统处于具有较低能量的状态。因此，如果期望在量子退火中100％的时间内找到基态，则从初始状态到最终状态不应添加噪声直到选择基态的对称性。因此，对系统的转变施加了严格的约束，诸如绝热演化的实现。

另一方面，如果在发生热耗散的前提下主动地采用噪声，则可以获得还包括激发态的状态，尽管到达基态的概率不是100％。

<<量子试行和经典MCMC试行之间的效率比较>>

接下来，将通过噪声添加的量子试行获得的结果与通过处理与经典伊辛模型相同的问题而获得的具有副本交换的经典MCMC方法的结果进行比较。该比较的目的是检验基于量子试行的采样计算与经典的基于MCMC的采样计算的效率。赋予优化设备100的计算条件如下。

经典伊辛模型的数值解是基于MCMC的。采用两种方法进行求解：单自旋翻转法和对自旋翻转法。单自旋翻转方法在N个自旋中随机选择一个自旋并翻转所选择的自旋。成对自旋翻转方法同时进行在N个自旋中随机选择状态1的自旋并将其翻转到状态0，并且随机选择状态0的自旋并将其翻转到状态1，以此保持状态1的自旋数目相同。

因为在TSP的情况下状态1的自旋的数目被确定为城市的数目，所以采用成对翻转方法。成对翻转方法保持答案的数量相同，因此与单翻转方法相比显著减小了搜索空间的大小。Metropolis方法被用来决定是接受还是拒绝通过翻转建立的新状态。另外，因为简单的MCMC太慢，所以经典的MCMC计算中使用了副本交换采样。副本的数量固定为10个。对于每个副本，执行十亿次试行，并且副本交换的频率被设置为每100次一次。手动优化温度。

图21示出了对8城市TSP的量子试行的计算结果与经典伊辛模型MCMC的计算结果之间的示例性比较。在图21的曲线图58中，圆圈表示对8城市TSP的1000个量子试行当中的按能量升序排列的前100个试行的结果。对于8城市TSP上的经典伊辛模型MCMC，能量升序的前100个结果用正方形表示。在曲线图58中，水平轴表示排序(具有较低能量的结果排序较高)，而垂直轴表示能量。

可以理解，基态是通过量子理论和经典理论方法二者获得的。对于激发态，从底部到第11位(在量子理论方法中排至第100位)的激发态是彼此一致的。对于基态，八个基态是简并的。获得八个或更多个基态源自于计算的完全相同的量子态。由此可以看出，量子理论方法能够获得基态和激发态。

图22示出了根据试行次数的示例性能量转变。图22中的曲线图59表示作为试行次数的函数，能量值如何被更新。曲线图60是曲线图59的仅在从0到5的能量范围附近的放大图。在曲线图59和60中，水平轴表示试行次数，垂直轴表示能量。

在图22的示例中，在第75次量子试行中获得基态能量的值，如带圆圈的折线所指示的。然而，请注意，当达到基态时，取决于如何给出随机数而不同。如果随机种子改变，则存在在第三次试行中达到基态的情况。另一方面，在通过基于MCMC的采样来求解经典伊辛模型的情况下，在第47492次试行中获得解，如由具有三角形的折线所指示的。考虑到副本的数量，实际上在474920次试行中达到了基态。

然而，应该注意的是，474920次试行的绝对值没有太大的意义，因为基于MCMC的副本交换可以通过优化而加速10倍到100倍。在基于MCMC的采样中总是看到，最小值在初始老化(初始松弛)结束之后按照使得几乎与时间的对数成反比的慢的速度更新。因此，缩放特性非常重要。

在8城市TSP的情况下，缩放特性的影响并不显著；但是，随着城市数量的增加，缩放特性的差异清晰可见。这一点稍后通过从随着城市数量增加的TSP获得的计算结果来说明。

在量子理论方法的结果中看到的一个重要特征是噪声和热耗散效应引起的去激励。该去激励可被认为是通过隧穿效应从高能态到低能态的转变。因此，获得的最终状态相对接近于基态。

<<ulysses16的情况的缩放评估>>

接下来描述在解更新速率的缩放方面与经典的基于MCMC的采样的比较。这里使用TSPLIB中包括的ulysses16，其最优解已知为6859。与8城市TSP不同，ulysses16涉及大量路径，因此足够复杂以突出经典理论方法和量子理论方法之间的差异。

对于使用经典伊辛模型解的示例性计算，在副本交换中使用10个副本，并且对于每个副本执行10亿次试行，如在8城市TSP的情况中那样。在计算的示例中，副本交换的频率是每100次一次。对于QTSA方案，根据方程式(6.1)，G₁被设置为0.001，而G₂从1减小到0。量子试行的次数是1000次。在单粒子激励的基础上添加噪声。

图23示出了QTSA和经典MCMC在ulysses16上获得的解的示例性比较。左侧曲线图61示出了经典MCMC和QTSA之间的以能量升序的前100个试行中的比较。在曲线图61中，水平轴表示排序，而垂直轴表示能量。右侧的曲线图62描绘了作为横向场的强度G₂的函数的能量(解目标)的变化。在曲线图62中，水平轴表示横向场的强度，而垂直轴表示能量。

参照图23，在基于量子理论的计算中达到基态(曲线图61中的圆圈)。与经典的MCMC(曲线图61中的正方形和三角形)的结果相比，量子理论方法产生更好的结果。由于参数集用于为每个副本运行10亿次试行，经典的MCMC无法到达解。注意，在曲线图61中，SSF和PSF分别代表单自旋翻转和对自旋翻转。SSF的最小值为7442，PSF的最小值为7181。通过量子理论方法获得的最小值是6859，其与官方最小值匹配。

图24示出了量子理论方法和经典理论方法之间的解改进速率的示例性比较。图24表示作为试行次数的函数，解如何被更新。左侧曲线图63描绘了从第一试行开始的解更新的转变。右侧曲线图64描绘了在老化期结束之后解更新的转变。在曲线图63和64中，水平轴表示试行次数，垂直轴表示能量。

曲线图63中的前1000次试行被认为是老化期。老化期间的更新并不是解更新速率的重要部分。在经典的MCMC采样(具有正方形和三角形的折线)的情况下，如曲线图63和64所示，在初始老化期结束之后，更新速率急剧地减慢，并且解更新随着试行次数的增加而对数地进行。另一方面，在量子理论方法(带圆圈的折线)的情况下，测量从一开始就在8000左右开始，然后在第597次试行中达到解。

接下来描述分叉加速计算的机制。这种现象是通过观察自旋微观状态的时间变化来解释的。

图25示出了典型量子试行的时间传播。图25描述了在QTSA期间，在一个典型的量子试行中，自旋微观状态的时间变化。在图25中的曲线图65和66中，水平轴表示虚数时间，垂直轴表示状态|1>下自旋的观察概率(波函数绝对值的平方)。曲线图66是曲线图65在从0到0.2的观察概率的范围附近的放大图。

从曲线图65可以理解，自旋翻转不是同时发生的，而是特定的自旋翻转，并且被设置在该状态中，然后发生连续的自旋翻转。在曲线图65中，在虚数时间τ≈15附近观察到上述分叉。在经典伊辛模型中，当已经发生这种现象时，2^N个候选解被减少到四分之一，即2^{N -2}。这是因为两个自旋的经典状态已经被确定了。另一种现象也可能发生，即，设置一个自旋的状态，然后驱动许多量子比特束一起设置在经典状态。这在曲线图66中可见。将一个量子比特的状态设置为1，然后几乎同时将捆绑在一起的许多量子比特的状态设置为0。束中的量子比特的数量是变化的，有时是2个，有时是10个或更多。重要的是，一旦这种现象发生，给定束中的量子比特数为q，解搜索空间立即变为2^N-q。解空间每次指数地减小。例如，在曲线图66的情况下，一束约20个量子比特一起一次变为状态0。假设一束包括设置为状态|0>的大约20个量子比特，每当这种现象发生时，解空间的大小就会减少到百万分之一的一小部分，因为2²⁰大约相当于百万分之一的一小部分。因为该现象在短时间内连续发生，所以解空间的大小非常快地减小。

这相当于比二叉树更快地搜索大小为2^N的解空间。即，指数大的搜索空间以指数减小。

接下来描述能量如何通过热耗散而损失。

图26示出了量子理论方法中的能量耗散。图26示出了在达到ulysses16的最优解(即，6859)的过程中通过虚数时间传播的热耗散。在图26的曲线图67和68中，水平轴表示虚数时间，垂直轴表示能量。曲线图68是曲线图67在从6840到7000的能量范围附近的放大图。

图26中描绘的现象总是发生在每个量子试行中。首先，当能量接近横向场的基态时，能量逐渐减小。然而，随着虚数时间传播的进行，在一定时段期间，存在以非常短的时间的能量骤降。此时，系统的对称性几乎被确定，这意味着此时确定哪些自旋处于状态1或0。随后，存在能量实际上保持恒定的时间段。如图25的自旋翻转图所示，并不是所有的自旋都在同一时间翻转。有些自旋还没有设置，因此，在这里发生了非常缓慢的能量弛豫。然而，在某一点上，一些自旋被确定为最终处于状态1，然后其它自旋随后处于状态0。然后，能量再次开始指数地减小，并且达到6859的最终状态。这很容易理解，因为这是一种正在发生的基底选择。

如果量子比特的数量增加，这种趋势就更明显了。因此，接下来提供包括比ulysses16更多数量的城市的ulysses22的情况作为示例。ulysses22也是公开的，其最优解已知为7013。当ulysses22被映射到伊辛模型中时，量子比特的数量是484。

图27示出了关于ulysses22的排序信息的示例。图27中的曲线图69描述了在基于QTSA的200个量子试行中，以能量升序排列的前50个试行的结果。给予优化设备100的计算条件如下。

经典的基于MCMC的采样中的副本数量是20个。以与ulysses16中相同的方式实现副本交换，并且在温度空间中执行副本的随机游走。另外，模拟温度被设置为足够高的值。对于每个复本，执行10亿次试行。采用SSF(曲线图69中的三角形)和PSF(正方形)两种方案。惩罚值被设置为4000。

曲线图69中的实心三角形表示满足约束条件的解。在SSF的情况下，尽管能量降低，但由于惩罚值设置为4000，所以不排除违反约束的解。得到的最小能量值为7861，然而，这是不可行的解。在曲线图69中，由实心三角形描绘在所获得的解的前50个中的仅可行解(即满足约束条件的解)。满足约束的解的最小值为8365，其能量比最优解高1300以上。在PSF的情况下，因为自旋翻转使得满足约束条件，所以所有解都满足约束。解的最佳值为7852。

另一方面，通过量子理论方法获得的最小值在200次试行后为7019，其高于7013，但仍充分接近最优解。

图28示出了如何在ulysses22上更新解。在曲线图70和71中，水平轴表示试行次数，垂直轴表示能量。量子试行的结果用圆圈表示，而经典的MCMC的结果用正方形表示。曲线图71是曲线图70在7000至10000之间的范围附近的放大图。

老化期(初始松弛)的步数趋向于随着自由度更高而增加。此外，初始松弛后的松弛进一步减慢，因此经典的MCMC计算面临着严重的困难，其随着自由度的增加而以加速的速度增长。量子试行的情况并非如此，量子试行很快就得出了几乎最佳的解。虽然这两种方法都没有达到最优解，但是使用由量子理论方法获得的解作为经典MCMC计算的输入允许通过交换两个城市来获得最优解。注意，因为试行次数增加10倍导致解的改进很小，所以经典的基于MCMC的计算被终止。在曲线图71中可以看出，在10亿次试行之后通过经典MCMC计算获得的解仍然高于由量子试行方法产生的初始值。

<<量子-经典混合计算>>

对于通过使用量子试行方法求解ulysses22而获得的解7019，假设该问题实际上得到了解决。这可以通过检查自旋的微观状态来理解。在ulysses22的情况下，最优解中的城市序列是：{1,14,13,12,7,6,15,5,11,9,10,19,20,21,16,3,2,17,22,4,18,8}。所得解7013的微观状态具有按顺序交换的城市14和13，即：{1,13,14,12,7,6,15,5,11,9,10,19,20,21,16,3,2,17,22,4,18,8}。

这意味着，通过使用解7013作为不同解方法的初始值，在一次试行的最小值之后获得真实解。QTSA方案通过与MCMC的机制不同的机制来提供解候选。因此，当仅通过量子理论方法不能获得解时，将量子理论方法与经典理论方法相结合的QTSA可以能够获得真实解。上面描述的ulysses22的计算是表明量子理论方法和经典理论方法的组合很容易找到真实的解的示例。

根据QTSA方案，随机地获得基态。QTSA方案完全不同于MCMC方案，并且通过一种基于量子力学的隧道效应来求解。一旦经典的基于MCMC的方法陷入局部最小值，就不容易摆脱它。在大的汉明距离下也很难达到状态。然而，这些问题在量子理论方法中是不相关的，因为能量损失是由热耗散引起的，并且就像通过隧穿效应穿过势垒一样，系统能量本身迅速减少并达到接近基态能量的值。

<<设置能量分辨率的效果>>

接下来，说明基于方程式(5.3)预先决定ITP方法的能量分辨率的效果。能量分辨率用作测试虚数时间传播中的能量收敛的阈值。

为了证明决定能量分辨率的效果而进行的数值实验如下。使用图18中的8城市TSP。所检查的能量分辨率为ΔE＝10^-5、10^-6、10^-7、10^-8、10^-9、10^-10、10^-11和10^-12。对于每个能量分辨率，进行1000次量子试行。给定波函数初始状态的噪声在单粒子激发基础上随机添加。噪声的幅度也是随机确定的。因此，即使当噪声被添加到相同的量子比特时，也会达到不同的最终状态。这是因为方程式(2.16)是非线性方程式并且取决于噪声的幅度。

下面描述数值实验的结果的概要。在能量分辨率大于ΔE＝10^-6的情况下，能级没有达到对应于经典能级的状态。这对应于不确定性原理方程式(5.3)的观察时间缩短的情况，并且还意味着，在量子退火中，磁场快速减小导致失败。

在ΔE＝10^-7的情况下，在1000次量子试行中的858次量子试行中，能级收敛到与经典伊辛模型相对应的状态；然而，在剩余的142次试行中，能级未能收敛到与经典状态相对应的状态。在ΔE≤10^-8的情况下，在所有1000次量子试行中，能级成功地收敛到对应于经典状态的状态。

这表明在量子退火中存在要遵守的磁场的最大减小速率。如果减小磁场的速度太高，则能级甚至不收敛到与经典状态相对应的状态，更不用说基态。

接下来描述增加能量分辨率的效果。

图29示出了对于每个能量分辨率的量子试行的示例性结果。图29中的曲线图72是在ΔE＝10^-7和10^-10的不同分辨率的1000次量子试行中获得的能量值的直方图。在曲线图72中，水平轴表示能量，垂直轴表示频率(观察到相应能量值的试行次数)。直方图仅包括所获得的能量状态对应于经典状态的试行结果。在ΔE＝10^-7的情况下，采样失败的概率是固定的，这导致能量状态不能映射到经典状态。不能被映射的能量状态值的数量级为10²。

直方图描述了通过量子试行获得的能量状态值的概率分布。从图29中可以看出，在较高能量分辨率下，更频繁地观察到具有较低能量值的状态，并且基态也更可能被采样。

当能量分辨率低(粗分辨率)时，分布为钟形。因为随着能量分辨率的程度变得更高，成功采样基态的概率增加，所以分布的整体形状向右倾斜，在较低端有更多的观测结果。

虽然随着能量分辨率的程度被设置得越高，在较低端上采样能量状态的成功概率增加，但是计算花费的时间更长。在虚数时间传播中，能量倾向于指数衰减，因此，如果能量分辨率增加到十分之一或百分之一，计算时间最多保持在几倍以内。然而，能量分辨率越高，计算时间越长。

这表明，在采样中，存在用于量子退火的磁场的最佳减小速率。通过更快地减小横向场来减少计算时间；然而，这使得热耗散机制很难工作，这反过来又降低了获得基态的概率。相反，横向场的足够慢的退火允许热耗散，但反过来增加了计算时间。

在实验方面，将热耗散效应作为非酉项(non-unitary term)纳入计算框架，并执行重整化，以避免波函数中途变为0。如果希望使实验系统更接近实际系统，则将热耗散项乘以系数以调节热耗散的幅度。此外，在不进行重整化的情况下，实验系统向实际系统移动。在这种情况下，当计算时间增加得更多时，波函数的振幅消失，并且只有噪声保留在末端。因此，在实际系统中不建议将能量分辨率增加得太高。存在随问题而变化的最佳观察时间。即，存在要遵守的最小退火速率。

执行波函数重整化的框架处理了这样一个理想条件，即，热耗散损失的概率总是从外部放大。

因此，结果是，在即使发生一阶量子相变也允许随机获得系统的基态的情况下，能量分辨率(退火速率)不必符合绝热跃迁条件。

(c)第三实施方式

根据第二实施方式，优化设备100向初始状态下的自旋的波函数添加噪声。也就是说，通过采用这种打破波函数对称性的措施，定义了量子试行。然而，QTSA控制方程是平均场近似。为了确保QTSA的解随机收敛到基态，对称性需要随着时间的推移而被打破。为此目的，优化设备100可以向磁场的系数添加噪声。

例如，优化设备100可以将如下面的方程式(8.1)所定义的噪声添加到第i个自旋的横向场J_i的系数。

这里，J是总横向场的强度。为了建立一个公式，假设作用在每个自旋上的横向场的强度是相同的。现在，优化设备100随机地将扰动作为噪声引入横向场的强度。以这种方式，添加非常小的噪声足以破坏系统的对称性。因此，这产生与向波函数添加噪声的情况类似的结果。

图30示出了在噪声被添加到横向场的情况下通过QTSA的量子退火的示例性结果。左侧的曲线图81表示作为横向场强度的的函数的能量。在曲线图81中，水平轴表示横向场的强度，而垂直轴表示能量。曲线图81中的实心圆表示当没有噪声加到横向场时获得的QTSA的示例性计算结果。空心圆和三角形都表示具有噪声添加的QTSA的示例性计算结果。

右侧曲线图82示出了作为横向场收敛的函数的在虚数时间传播之前所需的计算步骤的数目。在曲线图82中，水平轴表示横向场的强度，垂直轴表示在虚数时间传播收敛之前所需的计算步骤的数目。实心圆描绘了当没有噪声加到横向场时获得的QTSA的示例性计算结果。空心圆描绘了添加了噪声的QTSA的示例性计算结果。

图30的计算示例是在典型的量子退火协议下能级不降低并且系统被捕获在亚稳态的情况。如曲线图81所示，当在计算中向磁场引入不对称性时，横向场减小，然后能量迅速减小。如果横向场在负方向上增加，则系统试图返回到其原始状态。也就是说，当系统的对称性以某种方式被破坏时，系统的能量降低。

同样如曲线图82所示，在不引入不对称性(实心圆)的情况下，虚数时间传播以几乎恒定的试行次数收敛。可以看出，试行次数在G₂≈0的范围内增加。这是因为横向场的强度降低，因此不太可能发生自旋翻转。另一方面，在磁场中引入不对称性突然延长了G₂的多个值的寿命，这代表了许多亚稳态的出现。因此，事实证明，轻微的不对称会对系统的行为带来实质性的改变，并且系统塌陷到较低的能量状态。

以上述方式，将不对称性引入磁场实现了与向波函数添加噪声类似的效果。注意，将噪声添加到上述磁场的方式仅是示例。如在向波函数添加噪声的情况下，可以选择特定的k个自旋，然后向k个自旋添加均匀或非均匀噪声。

(d)其它实施方式

具有噪声添加的QTSA可以允许并行处理。根据量子试行的定义，噪声添加产生不同的结果，两个不同的量子试行彼此独立地运行。因此，哈密顿量信息被给予单独的L个经典计算机器一次(L是大于或等于1的整数)，然后仅需要给出关于要添加的噪声的信息。计算后，然后简单地收集量子比特信息的结果。计算只涉及向量和矩阵运算。因此，期望该计算与基于超标量的并行计算机和诸如GPU的加速器兼容。

当被给定由伊辛模型表示的优化问题的系数时，优化设备100能够通过QTSA方案来计算最优解。即，如果优化问题由伊辛模型描述而不管该问题的区域和行业如何，则优化设备100能够计算最优解。例如，优化设备100能够在包括金融、分销、材料工程、计算化学和计算生物学的领域中找到最小化问题的解决方案。另外，优化设备100能够在包括金融、分销和医疗保健(制药)的行业中找到最小化问题的解决方案。

根据上述第二和第三实施方式，优化设备100由包括处理器101和存储器102的冯·诺依曼计算机结构来实现，然而，优化设备可以例如使用量子计算机来实现。在这种情况下，基于量子计算机的优化设备中的热耗散可以在现象上结合到冯·诺依曼架构中的优化设备100中。这允许基于冯·诺依曼计算机的优化设备100也用作用于设计基于量子计算机的优化设备的硬件的模拟器。

根据一个方面，可以有效地找到伊辛模型的基态。

Claims (9)

1.一种优化设备，该优化设备包括：

处理装置，该处理装置用于：

在通过运行在施加到伊辛模型的磁场减小时发生的所述伊辛模型的状态变化的模拟来求解所述伊辛模型的基态时，将与噪声对应的值添加到所述模拟中使用的系数中的一些，所述伊辛模型表示目标问题，以及

执行随着模拟中的时间推进而减小所述磁场的强度的实数时间传播的第一处理和基于虚数时间传播方法减小所述伊辛模型的能量的第二处理。

2.根据权利要求1所述的优化设备，其中，

所述处理装置将与噪声对应的所述值添加到所述伊辛模型的初始状态，并且使用所述伊辛模型的添加了与噪声对应的所述值的结果状态作为所述模拟的开始状态。

3.根据权利要求2所述的优化设备，其中，

所述处理装置将与噪声对应的所述值加到包括在所述伊辛模型中的多个自旋当中的一些自旋的波函数中所包含的系数。

4.根据权利要求3所述的优化设备，其中，

所述处理装置在所述多个自旋当中随机地选择将被赋予与噪声对应的所述值的一些自旋。

5.根据权利要求2至4中的一项所述的优化设备，其中，

所述处理装置多次对所述伊辛模型执行与噪声对应的所述值的添加，以生成所述伊辛模型的各自具有不同噪声模式的多个状态，并且通过使用所述伊辛模型的所述多个状态中的每个状态作为所述模拟的开始状态来执行所述第一处理和所述第二处理。

6.根据权利要求1所述的优化设备，其中，

所述处理装置将与噪声对应的所述值添加到各自指示作用在所述伊辛模型中包括的自旋上的所述磁场的强度的系数。

7.根据权利要求1所述的优化设备，其中，

所述处理装置：

确定基于量子力学的所述伊辛模型的每个自旋的状态以取基于经典力学的可能值之一，所述伊辛模型的每个自旋的状态通过重复所述第一处理和所述第二处理来获得，并且

执行使用每个自旋具有所确定的值的经典伊辛模型的状态作为初始状态，通过基于所述经典力学的解来获得所述伊辛模型的基态的第三处理。

8.一种非暂时性计算机可读记录介质，在该非暂时性计算机可读记录介质中存储使计算机执行处理的计算机程序，所述处理包括：

在通过运行在施加到伊辛模型的磁场减小时发生的所述伊辛模型的状态变化的模拟来求解所述伊辛模型的基态时，将与噪声对应的值添加到所述模拟中使用的系数中的一些，所述伊辛模型表示目标问题；以及

执行随着所述模拟中的时间推进而减小所述磁场的强度的实数时间传播的第一处理和基于虚数时间传播方法减小所述伊辛模型的能量的第二处理。

9.一种由计算机执行的优化方法，该优化方法包括以下步骤：

在通过运行在施加到伊辛模型的磁场减小时发生的所述伊辛模型的状态变化的模拟来求解所述伊辛模型的基态时，将与噪声对应的值添加到所述模拟中使用的系数中的一些，所述伊辛模型表示目标问题；以及

执行随着所述模拟中的时间推进而减小所述磁场的强度的实数时间传播的第一处理和基于虚数时间传播方法减小所述伊辛模型的能量的第二处理。

Images