从点集中估计变换的迭代量子方法

5290从点集中估计变换的迭代量子方法0Natacha Kuete Meli Florian Mannel Jan Lellmann Institute ofMathematics and Image Computing, University of Luebeck0https://www.mic.uni-luebeck.de0摘要0我们提出了一种迭代方法，利用绝热量子计算从点集中估计刚性变换。与现有的量子方法相比，我们的方法依赖于自适应方案来高精度解决问题，并且不会受到不一致的旋转矩阵的影响。在实验中，我们的方法在多个2D和3D数据集上表现出鲁棒性，即使在高异常值比率下也是如此。01. 引言和相关工作0量子计算（QC）的出现提供了一种完全新的计算范式，利用量子机制的原理，承诺在选定的问题上指数级加速，以及数据存储空间等资源的大幅减少[12, 25, 31,36]。直观地说，利用量子系统可以假设混合态，基本上同时存在于几个纯态中，允许计算同时作用于所有这些状态。这种效应被称为量子并行性，使得量子计算机与它们的经典模拟器有所不同，后者只能执行顺序计算[28]。绝热量子计算（AQC）是QC的一个子领域，已经成为解决在经典计算机上通常很难的组合优化问题的一种有前途的方法[21,22]。AQC优化算法通常解决的问题类别是所谓的二次无约束二进制优化（QUBO）问题，其形式为0min q ∈{ 0 , 1 } n  q � Wq  +  c � q, (1)0其中W是耦合矩阵，c是偏置向量。将实际优化问题编码为QUBO问题的形式通常是开发基于AQC的算法的第一步和最关键的一步。这一步骤之后是将问题嵌入到量子硬件中，采用诸如[8, 9]的技术。0（a）2D0（b）3D0图1. 我们的迭代量子方法（IQT）从（a）2D点集（MNIST[26]和Lena [20]）和（b）3D点集（Stanford bunny[30]和completion3D[35]）估计变换的结果。绿色点表示参考点，红色点表示模板点。左侧显示初始对齐，右侧显示注册结果。0问题随后在量子硬件上得到解决，解决方案的比特串q被解嵌入，原问题的解码被解码。已经为一些计算机视觉问题开发了基于QUBO的算法，例如匹配问题[2, 3,5]，这些问题自然处理二进制变量。在这项工作中，我们将重点关注一个相关问题，即对齐点集的问题。根据相应的地标确定两个坐标系之间最佳刚性变换被称为变换估计（TE）[11, 15,16]。它在机器人学和图像处理等计算机视觉领域中有实际应用。例如，在图像配准中，一个准备性的关键步骤是根据相应的地标粗略对齐要注册的图像[23]。从方法解决TE问题的角度来看，有一些期望的属性。它们包括方法准确地建模对齐点集的变换的能力，处理高维和尺寸增加的数据的能力，以及对真实数据挑战（如噪声和异常值）的鲁棒性[11]。虽然计算机视觉中有大量强大的经典算法来解决TE问题，例如ICP[4, 37]̸5300尽管已经有一些方法可以在经典计算机上解决2D变换估计问题，例如ICP [24]和CPD[24]，但在量子计算机上可用的方法仍然非常有限。Golyanik和Theobalt[14]首创了一种基于量子的方法（QA）来解决2D变换估计问题，该方法通过将旋转矩阵近似为基矩阵的线性组合来计算二进制变量，其中线性系数由AQC计算得出。QA的局限性包括其精度固定、不能处理噪声以及由于自由度和所需量子比特数量的增加而难以扩展到3D。此外，它生成的矩阵通常不是正交的，算法也不能确保它接近正交矩阵。在这项工作中，我们在[14]的基础上提出了一种新方法，可以在2D和3D点集上解决变换估计问题，即使在存在噪声的情况下，也可以达到任意所需的精度，并且其矩阵在精度增加时收敛到正交矩阵。本工作的主要贡献如下：0•基于K位二进制表示，我们提出了一种逼近特殊正交群的方案，以便使用绝热量子计算估计旋转矩阵（第3.1节）。0•我们使用逼近方案首先推导出变换估计问题的QUBO公式（第3.2节），然后提出了一种算法，通过适应逼近方案使其在实践中的误差消失（第3.3节）。0•我们使用D-Wave量子退火机进行2D和3D数据集的配准实验证明了该算法（第4节）。其中一些结果如图1所示。02. 数学准备02.1. 从点集中估计变换0设N，D∈N。给定参考集X={xi}Ni=1和模板集Y={yi}Ni=1，变换估计（TE）的目的是确定一个刚性变换f�: RD →RD，解决以下问题：0min f∈SE(D) d(X, f(Y)), (RP)0其中d是适当的相似度度量，SE(D)是D维特殊欧几里德群。后者由旋转和平移的组合形成；其自由度通常计算为P=D(D+1)/2。任意变换f∈SE(D)可以表示为f(y)=Ry+t，其中R∈RD,D是旋转矩阵，t∈RD是平移向量。选择d为0最小二乘距离为0d(X, f(Y)) =0i=1 ∥xi − f(yi)∥220i=1 ∥xi − Ryi −0（2）众所周知[11,16]，通过引入质心居中的质点坐标˜X=X−¯X和˜Y=Y−¯Y，其中¯X和¯Y表示原始点集的质心，刚性配准问题可以简化为寻找最佳旋转矩阵R�，即以下问题的解：0minR∈SO(D)0i=1 ∥˜xi − R˜yi∥22. (RPR)0这将问题的变换空间减小到了D维特殊正交群SO(D)，其自由度为P=D(D−1)/2。最优平移t�由R�通过t�=¯X−R�¯Y推导得出。02.2. 绝热量子计算0绝热量子计算（AQC）是一种基于绝热定理[1,17]的优化范式，该定理表明机械系统的最低能态，即基态，是一个优化问题[18]的解。理论上，n个粒子的量子系统在时间t∈[0, T]内的演化，通常重新缩放为s=t/T∈[0,1]，可以用希尔伯特空间H=C�2n上的哈密顿量H(s)来描述，系统的状态由H中的一个单位向量|ψ(s)�给出。0量子比特。量子比特（qubit）|ψ�是量子信息的基本单位[25, 31]。它是一个单粒子量子系统的状态，根据|ψ� = α|0� +β|1�的超位置表示，其中|α|² + |β|² = 1，α，β ∈ C。状态|0�= (1, 0)�和|1� = (0,1)�是系统的基态，也称为本征态或稳态。之所以这样称呼是因为当量子比特进行测量时，它以概率|α|²或|β|²不可逆地坍缩到|0�或|1�。我们用|+�表示均匀叠加态，即|+� = (|0� + |1�) /√2。02.这个状态矢量的定义可以扩展到n粒子量子系统，通过考虑n个单比特态的张量积|ψ� = �ni=1|ψi�。系统的本征态由{|0�,|1�}�n给出。0横向场IsingHamiltonian。用于优化的绝热量子算法通常使用所谓的随机Hamiltonians[1]。在计算基 {|0�, |1�}�n中，如果Hamiltonian H满足以下条件，则它是随机的：0�q|H|q'� ≤ 0 �q, q' ∈ {0, 1}n with q ≠ q'. (3),(4)5310使用随机Hamiltonians的AQC的主要思想是用与横向场成比例的Hamiltonian H0初始化系统，即0H0 = -A0�n�0k = 1 σxk0�0其中σxk表示作用在系统第k个粒子上的x方向的Pauli矩阵[25]，A是一个正常数。H0的基态是均匀叠加态|ψ(0)� =|+��n。然后，根据Ising能量函数准备另一个HamiltonianH1，其中粒子k和l之间有耦合强度Wkl，并且粒子k上有外部偏置hk。0H1 = -B0�0k,l = 1 Wklσzk�σzl+0k = 1 ckσzk0�0�, (5)0其中σzk再次表示作用在第k个粒子上的z方向的Pauli矩阵，B是一个正常数。如果参数Wkl和ck被选择为QUBO问题的二次和线性参数，如式(1)所示，那么H1的基态|ψ(1)�将编码该问题的解。随着时间s从0到1的演化，初始HamiltonianH0转化为H1，描述了一个时间相关的随机系统Hamiltonian。0H(s) = (1 - s)H0 + sH1. (6)0量子系统演化。由H(s)生成的系统在时间s∈[0,1]上的演化由Schrödinger方程控制：0i�∂|0T∂s = H(s)|ψ(s)�, (7)0其解决方案定义了初始状态矢量|ψ(0)�的时间相关的幺正变换U(s)[1]。如果s满足绝热条件，即s变化足够缓慢，那么绝热定理表明U(s)将以很高的概率将H0的基态映射到H1的基态，这是优化问题的解[17]。0D-Wave量子退火。D-Wave量子计算机[34]旨在使用形式为式(6)的横向场IsingHamiltonians解决形式为式(1)的QUBO优化问题。量子处理器根据问题规模创建一个逻辑比特网络，然后将其嵌入到量子硬件中。网络在s =0处以全局叠加的方式开始，覆盖所有可能的本征态。在称为量子退火的过程中，随着s→1，提供的耦合和偏置被改变为0通过改变状态景观的磁场，强调最有可能成为优化问题解的状态。D-Wave量子退火机通过Leap量子云服务[32]提供，并且可以使用Ocean软件[33]在Python中实现D-Wave量子算法。03. 迭代量子变换估计（IQT）0像许多计算机视觉问题一样，(RPR)本质上不是一个二进制优化问题，因此不能直接在当前的量子退火机上求解。我们展示了通过线性化旋转矩阵并引入二进制系数来近似最小二乘配准问题的方法可以将其转化为QUBO问题。由于2D和3D情况非常相似，我们从一个一般的角度来呈现我们的方法，并且只在必要时切换到逐个讨论的情况。与[14]不同，我们的方法将量子比特数K作为用户可以选择的自由参数。这对于在当前的量子硬件上进行计算是至关重要的。尽管量子比特的数量是固定的，但我们的算法逐步增加了旋转参数的近似精度，从而使其能够以任意精度进行计算。另一个问题是在旋转群SO(D)上进行优化需要对旋转矩阵R施加正交性约束。避免这个约束的一种方法，如[10,19]中已经介绍的，是利用李群SO(D)与其代数so(D)之间的一一对应关系。这允许在线性空间so(D)上优化(RPR)，该空间是维数为D的反对称矩阵M的集合。具体来说，我们将使用指数映射R =exp(M)将M∈so(D)与旋转矩阵R∈SO(D)关联起来。0定点表示。为了以逐渐增加的精度表示解，我们使用K位来近似表示区间[a, b]中的点x，具体如下：0x = a + b - a0s0k=0 qk 2k, (8)0其中二进制变量qk∈{0, 1}，缩放因子s =2K-1。显然，x属于通过将区间[a,b]离散化使用s个bin得到的网格。我们将读者引用到[27,29]中关于量子退火机上固定点和浮点表示的替代方法。1)̸(17)(18)(20)53203.1. 旋转矩阵的近似02D中的旋转矩阵。在2D中，唯一可能的反对称矩阵M∈so(2)是单参数矩阵0M(θ) = θS = θ ∙ 0 -1 1 00�, (9)0其中θ∈R表示旋转角度。相关的旋转矩阵R为0R = eM(θ) = (cosθ)I + (sinθ)S, (10)0任何2D旋转矩阵都可以用这种方式表示[13]。在[14]中，作者提出将cosθ和sinθ视为独立变量（如ω1，ω2），然后对两者进行优化。这样做的缺点是得到的矩阵R = ω1I +ω2S可能不是旋转矩阵，因为这要求ω21 + ω22 =1。为了避免这个问题，我们将R视为θ的非线性函数，并在多步骤过程中确定最优θ。给定当前角度θc，我们在θc附近线性化R：0R(θ) = g(θ)I + h(θ)S ≈ g(θc) + g'(θc)(θ-θc)I0+ h(θc) + h'(θc)(θ-θc)S，0其中 g(θ) = cosθ，h(θ) =sinθ。为了得到一个QUBO问题，我们根据公式(8)将θ量化：对于一个固定的搜索窗口大小∆>0，我们将θ∈[θc-∆, θc+∆]表示为θ = θc-∆+2∆0s ∈ K − 1 k =0 q k 2 k. Inserting this representationinto Eq. ( 11 ) and using ˆ v := 2∆ s ∈ K − 1 k =0 q k 2k, we obtain the approximation0R (ˆ v ) ≈ R c + g ′ ( θ c ) I + h ′ ( θ c ) S ˆ v,(12)0where R c := [ g ( θ c ) − g ′ ( θ c )∆] I + [ h ( θ c ) − h ′ (θ c )∆] S is a term that depends on θ c and ∆, but isindependent of the unknown ˆ v and hence of theoptimization variables q k ∈ { 0 , 1 }, k = 0 , . . . , K − 1.For any of the mass-centered points ˜ y i, i = 1 , . . . , N,as in ( RPR ), we thus find0R ˜ y i ≈ R c ˜ y i + R i ˆ v, (13)0where0R i := g ′ ( θ c ) I + h ′ ( θ c ) S ˜ y i (14)0is independent of ˆ v and hence of the unknowns q k.The approximation ( 13 ) is at the heart of the QUBOproblem that we solve.0Rotation Matrix in 3D. In the 3D case, the completerotation parameters can be recorded in a vector v = ( v 1 ,v 2 , v 3 ) � ∈ R 3 encoding the rotation angle θ = ∥ v ∥ 2and the rotation axis x = v/ ∥ v ∥ 2 = ( x 1 , x 2 , x 3 ) �. Theskew-symmetric matrix M ∈ so (3) has the form0M ( v ) = θ ∙0� 0 − x 3 x 2 x 3 0 x 1 − x 2 x 1 00�0� , (15)0and the exponential map associates to M ( v ) therotation matrix0R = e M ( v ) = I + si0θ M ( v ) + 1 −θ0θ 2 M 2 ( v ) , (16)0where θ = θ ( v ) = ∥ v ∥ 2, see [ 13 ]. For convenience,we introduce the functions g ( v ) = (sin ∥ v ∥ 2 ) / ∥ v∥ 2 and h ( v ) = (1 − cos ∥ v ∥ 2 ) / ∥ v ∥ 2 2 for v ≠ 0and g (0) = 1, h (0) = 1 / 2, g ′ (0) = h ′ (0) = 0. Since M( v ) is linear in v, we have M ′ ( v c )( w ) = M ( w ) forany v c, w ∈ R 3. Hence, linearization of R around thecurrent vector v c yields0R ( v ) = I + g ( v ) M ( v ) + h ( v ) M 2 ( v )0≈ I + g ( v c ) + g ′ ( v c )( v − v c ) M c + g ( v c ) M( v 0+ h ( v c ) + h ′ ( v c )( v − v c ) M 2 c0+ h ( v c ) M ( v − v c ) M c + M c M ( v− v c ),0where we used M c := M ( v c ). Assuming that the components v j, j = 1 , 2 , 3 arein the search window [( v c ) j − ∆ , ( v c ) j + ∆] for some ∆ > 0, we make use ofthe discretization v j = ( v c ) j − ∆ + 2∆0s ∈ K − 1 k =0 q ( j ) k 2 k from Eq. ( 8 ). Introducing qk := ( q (1) k , q (2) k , q (3) k ) � and ˆ v := 2∆ s ∈ K − 1k =0 q k 2 k we can, in slight abuse of notation, rewritethis as v = v c − ∆+ˆ v. Thus, we arrive at theapproximation0R (ˆ v ) ≈ R c + g ( v c ) + h ( v c ) M cM (ˆ v)0+ h ( v c ) M (ˆ v ) M c + M c g ′ ( v c )ˆ v + M 2 c h ′ ( v c )ˆ v,0where again the term R c depends on v c and ∆, butnot on the optimization variables q k ∈ { 0 , 1 } 3 , k =0 , . . . , K − 1. It is straightforward to check that M (ˆ v )w = − M ( w )ˆ v for any w ∈ R 3, so we deduce0R ˜ y i ≈ R c ˜ y i + R i ˆ v, (19)0where the matrices R i are given for i = 1 , . . . , N by0R i = M c ˜ y i g ′ ( v c ) + M 2 c ˜ y i h ′ ( v c )0− h ( v c ) M ( M c ˜ y i ) − g ( v c ) + h ( v c )M c M (˜ y i),(21)where forkP,Pminq∈{0,1}KPminv∈RP v⊤Wv + c⊤v,(25)5330and are independent of ˆ v. As final step, we introduce0q =0�0� 0q 0 q 1... q K −10�0� � � � ∈ { 0 , 1 } KP and0U = 20s ∈ R P,P K0矩阵 D_k ∈ R^P,P03.2. QUBO 表述0根据设计，根据公式 (13) 和 (19)，我们有 R˜y_i ≈ R_c˜y_i+ R_iˆv。关键是，ˆv 和因此 q在这个近似中呈线性关系，因此 ∥˜x_i - R˜y_i∥2是 q的二次函数。通过忽略与 q 无关的项，我们发现 (RPR)可以用 ˆv 的形式近似为0N∙0i = 10∥R_iˆv∥2 + 2�R_c˜y_i - ˜x_i, R_iˆv�. (22)0根据公式 (21) 推导，我们有 ˆ v =Uq。因此，我们提出以下近似 (RPR) 的 QUBO：0min q ∈ {0, 1}^KP q^TWq +c^Tq，(QUBO)0其中，使用符号 ˆ y_i = R_c˜y_i - ˜x_i，0W = U^T∙N∙0i = 1 R^T_iR_i0∙0U 和 c = 2U^TN∙0i = 10(23) 要在量子退火器上解决(QUBO)，我们只需要传递耦合参数 W ∈ R^KP,KP 和偏置c ∈ R^KP。注意，W 和 c 的维度不依赖于点的数量N。在退火过程结束时，测量得到的比特串 q被用于根据公式计算旋转参数0θ = θ_c - ∆ + Uq，或者 v = v_c - ∆ + Uq. (24)03.3. 量子变换估计0我们的迭代量子变换估计策略 (IQT)，描述如算法1，通过解决 (QUBO) 来估计不同的 Taylor 近似点 θ_c(二维0在二维情况下，我们首先对旋转参数进行初始猜测，例如 θc = 0 ∈ R；在三维情况下，我们对 v c进行初始猜测，例如 v c = 0 ∈R^3。对于搜索区间，我们使用初始半径 ∆ = π。组装矩阵W 和 c 后，解决 (QUBO)。解 q 通过公式 (24)提供下一个迭代的 θ 或 v。我们算法的一个重要部分是将 ∆视为自适应参数；我们在可能的情况下减小∆，这会导致对解和相关旋转矩阵的逼近越来越好。当当前旋转参数与上一个旋转参数之间的误差小于某个与分箱大小2∆/s 相关的阈值函数时，我们减小∆，这被视为已达到当前 [v c - ∆, v c + ∆] 区间离散化的 K位的最佳解的指示器。0算法 1：基于点集的变换估计问题的迭代量子方法 (IQT)0数据：˜ X，˜ Y：要注册的点集 maxit：最大迭代次数 θ c 或 vc：初始旋转参数 ∆：搜索区间的半径 κ：∆ 的减小阈值K：量子比特数 j ← 0，初始化 R c，并令 τ = κ∆02K - 1. 当 j < maxit 时0为 i = 1，...，N 构造 R i，U 和 ˆ yi。根据公式 (23) 计算 W 和 c。解决(QUBO) 并测量解 q。根据公式 (24) 计算θ 或 v。如果 |θ - θ c| < τ，或者 ∥v - v0∆ ← ∆/2. τ ← τ/2. end θ_c ← θ，或者 v_c← v. R_c ← R_c(θ_c, ∆)，或者 R_c ←R_c(v_c, ∆)。j ← j + 1. end 返回当前参数θ_c，或者 v_c，以及相关的 R (通过公式(11)，或者公式 (17))。03.4. 经典变换估计0在第4节中，我们将IQT算法与其经典松弛版本进行比较，该版本被称为迭代经典变换估计（ICT），它是通过将QUBO问题替换为一个经典的无约束二次规划（QP）问题得到的。具体而言，根据公式（13）和公式（19），使用连续变量v ∈ R P ，我们构造了QP问题，其中 R ˜ y i  =  R c ˜ y i  + R i v 。2D3D5340QPU/拓扑结构 2000Q Advantage 1.10量子比特数  K 3 5 10 3 5 100Synthetic | θ  −  θ � | 2 . 42 e  −  5 1 . 06 e  −  10 1 . 66e  −  14 2 . 28 e  −  4 1 . 13 e  −  8 2 . 58e  −  12 ∥ R  −  R � ∥ F 3 .42 e  −  5 2 . 26 e  −  10 2 . 24e  −  14 4 . 00 e  −  4 1 . 06 e  −  8 3 . 65e  −  12 平均链断裂 0 0 0 0 0 00MNIST | θ  −  θ � | 1 . 71 e  −  4 1 . 35 e  −  8 2 . 55e  −  15 4 . 70 e  −  5 1 . 00 e  −  8 2 . 65e  −  9 ∥ R  −  R � ∥ F 1 . 28 e−  4 1 . 91 e  −  8 4 . 71e  −  15 6 . 65 e  −  5 1 . 41 e  −  8 3 . 74e  −  9 平均链断裂 0 0 0 0 0 00Lena | θ  −  θ � | 2 . 53 e  −  4 7 . 94 e  −  9 5 . 32e  −  15 2 . 08 e  −  4 1 . 44 e  −  9 6 . 33e  −  12 ∥ R  −  R � ∥ F 3 . 58 e−  4 1 . 12 e  −  8 7 . 25e  −  15 2 . 95 e  −  4 2 . 03 e  −  9 8 . 95e  −  12 平均链断裂 0 0 0 0 0 00Synthetic ∥ v  −  v � ∥ 2 4 . 20 e  −  4 4 . 00e  −  6 8 . 76 e  −  1 9 . 70 e  −  5 6 . 71e  −  7 0 . 80 e  −  2 ∥ R  −  R � ∥ F 3 .80 e  −  4 1 . 63e  −  6 6 . 39 e  −  1 1 . 73 e  −  4 1 . 45e  −  6 1 . 96 e  −  1 平均链断裂 0 0 0 . 01 0 0 0 . 0020Bunny ∥ v  −  v � ∥ 2 1 . 37 e  −  5 3 . 61e  −  7 1 . 25 e  −  1 3 . 11 e  −  4 9 . 52e  −  8 1 . 45 e  −  1 ∥ R  −  R � ∥ F 2 . 30e  −  4 1 . 51e  −  6 1 . 82 e  −  1 3 . 70 e  −  4 1 . 20e  −  7 3 . 12 e  −  1 平均链断裂 0 0 0 . 01 0 0 0 . 0010Completion ∥ v  −  v � ∥ 2 9 . 58 e  −  5 8 . 59e  −  8 5 . 54 e  −  1 7 . 75 e  −  5 6 . 89e  −  5 5 . 24 e  −  1 ∥ R  −  R � ∥ F 1 .28 e  −  4 1 . 12e  −  7 7 . 15 e  −  1 1 . 22 e  −  4 9 . 53e  −  5 4 . 85 e  −  1 平均链断裂 0 0 0 . 01 0 0 . 004 0 . 030表2. IQT的重构旋转参数θ和矩阵R与真值（θ �和v �，R�）之间的距离，对于不同数量的量子比特K在两个D-Wave图拓扑（2000Q和Advantage1.1）上进行了比较。结果是在2D和3D情况下相同问题的15次优化步骤之后报告的。此外，提供了在15次迭代中平均链断裂的比例。在2D情况下，使用K = 10个量子比特观察到最小的误差，在3D情况下使用K =5个量子比特观察到最小的误差。对于3D问题，较大的拓扑结构会导致链断裂，从而降低了优化过程的效果。0其中  W  = � N i =1 R � i R i  和  c  = 2  � N i =1 R � i ˆ y i。这个QP问题的解等价于线性方程组  2 Wv  =  − c 的解。04. 实验结果0在本节中，我们实验性地评估了IQT算法的准确性、鲁棒性和运行时间的混合量子-经典实现。耦合参数W和偏置c在Python 3.9.2上准备，在一台配备16GB RAM的Inteli7-7700KCPU机器上，随后在D-Wave量子退火器上使用默认的退火时间20微秒和每次优化迭代100次读取来解决QUBO问题。每次迭代中选择能量最低的解作为系统的最低能量本征态。实验中使用的数据集列在表1中。除了图1中显示的点集之外，还包括2D中的MNIST [26]和Lena edges[20]点集，以及3D中的Stanford bunny[30]和completion3D[35]点集。点的数量从合成点集的150个变化到Lena的4845个。0关于 K的消融研究。我们方法的设计选择包括用于数值表示的量子比特数  K 。0数据集名称 点的数量02D合成点 150  (可变) MNIST [26] 150 Lena边缘[20] 484503D合成点 150  (可变) Stanford兔子[30] 500 completion3D[35] 20480表1. 实验中使用的数据集及其包含的点数。0未知变量的表示和D-Wave图拓扑的评估。我们对于  K  ∈ {3 ,  5 ,  10 }的2D和3D点集进行了15次迭代，评估了重建的旋转参数和矩阵的精度。请注意，编码QUBO问题需要总共 KP个量子比特，其中在2D中  P  = 1 ，在3D中  P  = 3。我们在使用Chimera拓扑[7]的D-Wave 2000QQPU和使用Pegasus拓扑[6]的D-Wave Advantage 1.1QPU上运行实验。我们将计算得到的旋转参数和矩阵与相应的真值进行比较。结果显示在表2中。我们观察到020406080Outlier Ratio (% template)0.00.20.40.60.8Alignment Error (eA)020406080Outlier Ratio (% template)1021101810151012109106103100103Consistency Error (eR)IQT 2D (iter 5)ICT 2D (iter 5)IQT 2D (iter 10)ICT 2D (iter 10)IQT 2D (iter 15)ICT 2D (iter 15)QA 2D020406080Outlier Ratio (% template)0.00.20.40.60.81.0Alignment Error (eA)020406080Outlier Ratio (% template)101810151012109106103100103Consistency Error (eR)IQT 3D (iter 5)ICT 3D (iter 5)IQT 3D (iter 10)ICT 3D (iter 10)IQT 3D (iter 15)ICT 3D (iter 15)5350更多的量子比特并不一定会增加精度。例如，在3D情况下， K  = 10 时的误差要比  K  = 3  和  K  = 5时的误差大得多。这种差异与平均链断裂比例有关，意味着更多的量子比特会导致系统更加嘈杂，从而增加误计算的数量。在后续的实验中，我们在2D中使用  K  = 10，在3D中使用  K  = 5 ，并在D-Wave 2000QQPU上执行算法，这似乎是完成任务的最佳配置。0对异常值的敏感性。为了验证我们的方法对异常值的鲁棒性，我们随机向一部分模板点添加噪声；该部分被称为异常值比例。为了衡量所提方法的准确性，我们采用了[14]提出的两个度量标准：0并且 e R  := �� I  −  R � R �� F (一致性误差)0并且 e A  := ∥  ˜ X −  R ˜ Y∥ F0∥ ˜ X∥ F (对齐误差) ,0其中  ∥∙∥ F是Frobenius范数。一致性误差衡量了重建的旋转矩阵的正交性，对齐误差表示模板点映射到参考点的准确程度。我们注意到，如果  R � 是通过算法1生成的  v c  的最优解，并且v c  →  v � ，那么由等式（12）右边和等式（18）右边给出的旋转矩阵的近似分别收敛到R �0具有误差阶  O ( ∥ θ c  −  θ � ∥ 2 2 ) ，分别为  O ( ∥ v c−  v � ∥ 2 2 )。在图2中，我们显示了IQT和其经典对应物ICT的配准误差作为异常值比例的函数。我们绘制了2D和3D合成点集上IQT的第5、第10和第15次迭代的误差。在2D中，与QA方法[14]相比，我们将其与完整数据集作为交互点的QA方法的误差进行了比较，即  K  =  N。如预期的那样，IQT和ICT都会收敛到一个正交矩阵，同时减小对齐误差，而一致性误差几乎不受异常值比例的影响。相比之下，基准方法QA[14]在异常值比例增加时更容易产生非正交矩阵。接下来，我们强调了IQT和ICT的收敛性和鲁棒性的相似之处，从而证明了K位二进制离散化的有用性，以及量子计算机通过这种离散化来逼近实数的能力。图3展示了IQT在completion3D数据集[35]上的配准结果。我们以无噪声和50%异常值比例的点对结果进行了说明性展示。我们观察到，在强噪声存在的情况下，IQT能够找到一个良好的变换，其准确性与经典计算（ICT）相当，同时生成一致的（即0图2.IQT方法（绿色）及其经典版本ICT（红色）对抗异常值的鲁棒性。左图显示对于合成2D（上排）和3D（下排）点集的异常值比例，对齐误差eA（左图）和一致性误差eR（右图，注意对数刻度）的变化。所提出的IQT方法在避免了QA方法（蓝色）[14]引入的一致性错误的同时，始终达到与经典ICT方法相同的准确性。0几乎正交）旋转矩阵，而现有的QA方法都不具备这些特性。0计算成本和时间。我们方法中最耗费计算资源的操作与[14]相同，即每次迭代时耦合和偏置的准备工作。准备矩阵W和c需要O(Nξ2K2P4)和O(Nξξ′KP2)的操作，其中ξ和ξ′分别代表计算Ri和Rc˜yi−˜xi的成本。表3比较了IQT的QPU访问时间和IQT解决QP的CPU访问时间在不同参数配置下的总运行时间，平均值是基于10个问题实例和一次迭代。对于IQT和ICT，我们也提供了矩阵准备时间。我们观察到，无论不同的组合，QPU和CPU访问时间都相对稳定。相反，矩阵准备时间随着点的数量N的增加而增长，并且在适度的N下超过了QPU和CPU访问时间。0局限性。IQT应被视为向高效变换估计迈出的第一步，最终目标是点云N1501500200001501500200005360(a) 完美数据。0(b) 数据受到异常值的破坏。0图3.IQT在注册(a)完美点集对和(b)参考点集到模板点集的性能上。点集来自completion3D数据集[35]。绿色点表示参考点，红色点表示模板点。左侧显示初始对齐，右侧显示注册结果。即使在强噪声下，IQT算法也能提供稳健的变换估计。0量子比特数K IQT：矩阵准备/QPU访问时间（毫秒）ICT：矩阵准备/CPU访问时间（毫秒）02D K = 10 2 . 50 / 34 . 67 23 . 20 / 34 . 67 311 . 25 / 34 . 69 6 . 99 / 0 . 10 55 . 56 / 0 . 08 334 . 7 / 0 . 06 K = 20 2 . 51 / 35. 32 23 . 00 / 34 . 79 305 . 19 / 35 . 26 K = 40 2 . 65 / 35 . 42 23 . 14 / 35 . 39 306 . 12 / 35 . 3803D K = 5 3 . 43 / 35 . 04 32 . 07 / 34 . 92 425 . 56 / 34 . 96 6 . 58 / 0 . 01 34 . 84 / 0 . 01 421 . 77 / 0 . 01 K = 10 3 . 50 / 35 .23 32 . 16 / 35 . 23 420 . 18 / 35 . 24 K = 15 3 . 43 / 35 . 39 31 . 68 / 35 . 44 419 . 02 / 35 . 400表3.使用我们的量子公式（IQT，左）和其经典对应物（ICT，右）在合成数据上解决二次问题（QUBO）的运行时间比较。所有值都是基于N个问题实例的平均值。实验结果表明，IQT的运行时间与量子比特数K无关。矩阵准备占据了解决系统的QPU/CPU访问时间的主导地位。0量子硬件上的注册。目前，运行时间主要由量子和经典情况下矩阵的准备所主导。目前量子硬件的有限可用性为转移到量子计算服务提供商引入了额外的时间惩罚。最终，拥有一种不依赖AQC的内在（电路模型）量子公式将是有益的，这是我们留给未来的工作。05. 结论0我们研究了一种用于点集变换估计的绝热量子公式，并提出了一种迭代的量子方法来优