基于奇异值分解的矩阵低秩近似量子算法*

王芙蓉 杨帆 张亚 李世中 王鹤峰

1) (中北大学机电工程学院, 太原 030051)

2) (清华大学物理系, 北京 100084)

3) (西安交通大学理学院, 西安 710049)

在大数据时代, 高效的数据处理至关重要, 量子计算具有平行计算能力, 为方便处理数据提供了新的解决途径.本文提出了一个基于奇异值分解的矩阵低秩近似量子算法, 复杂度为 O [log(pq)].在核磁共振量子计算系统完成了算法的原理演示, 选择一个 8 ×8 维的图像矩阵, 实现共振跃迁算法的哈密顿量 H 的时间演化,用量子态层析法分别读出密度矩阵的不同成分, 对密度矩阵进行重构, 保真度为99.84%, 在误差范围内验证了本文提出的矩阵低秩近似量子算法的正确性.而通过奇异值分解计算低秩矩阵的经典算法的复杂度是O[poly(pq)], 量子算法与经典算法相比, 实现了指数加速.

1 引 言

随着互联网的技术不断提高, 大数据时代正在到来.矩阵作为处理大数据的基本工具之一, 常被用于解决许多实际问题.大多数情况下, 矩阵表示在空间和时间复杂度上会随着数据的规模呈二次方增长, 这导致数据处理困难.因此, 去除冗杂信息, 只保留有用的信息, 构造低秩矩阵, 快速高效地近似一个目标矩阵, 可提高机器学习和数据管理的效果.目前, 低秩矩阵已被广泛用于模式识别[1]、图像压缩和检测[2,3]、图像处理[4,5]、人脸识别[6]等应用领域.

矩阵的低秩近似是一种稀疏表示形式, 可以在保持原矩阵的诸多性质同时, 减少冗余和噪声, 降低了数据的存储空间和计算复杂度.常用的矩阵低秩近似方法包括: 低秩矩阵恢复[7], 主成分分析法(PCA)[8]、非负矩阵分解(NMF)[9]、奇异值分解(SVD)、核映射中的核主成分分析(KPCA)和拉普拉斯特征映射(LE)等.

2014年, Tipping和Bishop[10]采用PCA方法降秩, 他们将图像分割成8 × 8个不重叠的块, 进行均匀量化将比特平均分配, 用变换后的前4个主成分对整个图像进行重构, 模型似然度最大化后,将数据分配给重构误差最小的分量进行图像编码,最终的比特率为0.5比特/像素, 误差为 6.2×10-2.与PCA法相比, KPCA法既可以保持全局属性又可以获取非线性信息.2017年, 徐梦珂[11]提出基于二维核主成分分析法的拉普拉斯特征映射(2DKPCA+LE)算法进行人脸识别.在此方法中, 用2D-PCA算法训练样本去相关性, 之后利用核主成分分析法提取非线性特征, 最后用拉普拉斯特征映射再一次降秩, 其算法识别率高, 计算复杂度低.2001年, Tsuge S[12]等提出将NMF[13]应用于逐项文档矩阵中文档向量的降秩.NMF将非负矩阵分解为两个非负矩阵, 将分解后的矩阵之一视为基本向量.通过将文档向量投影到由这些基本向量形成的较低维空间上来实现降秩.大多矩阵的低秩近似方法都避免不了求解特征值或奇异值分解, 因此研究基于奇异值分解的低秩近似方法是提高大数据处理效率的一个核心问题.

量子计算基于量子力学原理, 利用量子叠加性、量子纠缠等特性进行计算, 在处理特定问题时相比经典计算能指数加速.利用叠加态原理, 将某些类型的数据存入量子系统所需要的物理资源也远小于经典系统.近年来, 许多求解线性代数问题的量子算法被提出, 例如量子线性方程组算法[14]、量子支持向量机[15]、量子主成分分析[16]、量子奇异值分解算法[17]等, 也提出了许多量子数据存储、读取的方案.在某些情况下, 这些量子算法和方案在时间复杂度或空间复杂度上, 相比经典计算机有着指数级优势.目前只拥有中等规模、带噪声的量子计算机, 大多数基于通用量子计算模型的量子算法还无法实现.利用有限的物理资源对量子算法进行验证, 是目前量子计算的重要研究方向.

本文提出了一种量子低秩近似算法, 该方法保留了图像数据的主要成分而去掉那些相对不重要或者由噪声引入的成分, 利用主成分分析的原理,实现了对图像数据矩阵“奇异值分解”、“奇异值过滤”、“矩阵重构”的步骤, 相比于经典计算, 本量子算法复杂度有指数加速.首次利用该量子算法替代了以往的相位估计算法, 大量地减少了对辅助比特的需求, 成功地实现了对图像矩阵的低秩近似, 验证了算法的正确性.

2 理 论

2.1 矩阵的低秩近似与奇异值分解

设 A 是一个秩为r的任意矩阵, 若r远小于A的行数和列数, 则称 A 是低秩矩阵.给定一个秩为r的矩阵 A , 欲求其最近似矩阵该问题可形式化为

任何实矩阵 A ∈Rm×n都可以分解为

其中, U ∈Rm×m和 V ∈Rn×n都是幺正矩阵, Σ 是对角元为从大到小依次排列的奇异值 σi的对角矩阵.对矩阵 A 进行奇异值分解后, 将矩阵 Σ 中的r-k 个最小的奇异值置零获得矩阵 Σk, 仅保留最大的k个奇异值, 则

(3)式就是(1)式的最优解(Eckart-Young-Mirsky定理), 其中 Uk和 Vk分别是(2)式中前k列组成的矩阵[18].

2.2 量子低秩近似算法

假设对于矩阵 A =[aij]∈Rp×q, 有1个量子操作黑箱(oracle) OA:

即给定i和j, OA能给出矩阵A对应的元素值.考虑到一般图像在量子态上的编码方式, 使用振幅编码制备如下量子态:

以上是将图像制备成量子态的过程, 也就是初始化的过程.该步可利用量子随机存储器(QRAM)[19,20],采用 O [log(pq)] 个步骤实现, 当然也可以采用其他初始化方法[21].

对于矩阵 A 的奇异值分解:

若只选取由大到小前 r′个大于阈值 τ 的奇异值重构出低秩近似矩阵 A′, 则有

其中, F作用于矩阵 A 时, 会将第 r′+1 到r个成分奇异值取反, 而其他不发生改变.初始化完成后的量子态:

则F可以利用相位估计算法[22]和阈值 τ 对应的反相算符 Oτ高效实现.相位估计算法可以总结为以下算符:

其中寄存器B和C分别用来存储厄米矩阵A的本征值的估计值和量子态 | ψ〉 ,是量子傅里叶变换算符的逆算符, H⊗t是t比特Hadamard门.Oτ可以认为是一个量子过滤器算符,

即对寄存器C中所有量子态二进制表示小于 τ 的态相位取反, 而其他则不变.

(7)式需要计算两个矩阵的差, 早期的酉算子乘积量子计算模型无法直接计算[23], 利用对偶量子计算(DQC)[24-26], 可以在量子计算机上实现(7)式, 对偶量子计算在包括计算空间和辅助比特的大空间中是酉的, 酉演化算符可以写成:

之后对辅助量子比特进行测量, 当测量结果为0时, 就得到了非幺正算符 ( U0+U1)/2 作用在量子态上的结果.

其中 ρAA†是与厄米矩阵 A A†对应的密度矩阵, 二者矩阵元相同, 仅相差一个常系数.

综上, 量子图像低秩近似算法可以概括如下.

输入图像矩阵 A 对应的量子态 | ψA〉 、幺正演化算符、阈值 τ 对应的反相算符 Oτ.

输出A 矩阵的低秩近似矩阵 A′对应的量子态

算法步骤

1) 在3个量子寄存器上准备初始态:

2) 对初始态运行量子相位估计算法, 其中幺正演化算符为.通过求偏迹制备出密度算符 ρAA†, 再根据Lloyd等[16]在2014年提出的方法, 用 ρAA†有效地实现算符, 之后得到量子态:

3) 利用DQC实现(7)式.令 U0=I , U1=Oτ2,得到:

其中前 r′个奇异值 σk大于等于设定的阈值 τ.

4) 运行步骤2)的逆运算, 并测量寄存器R.当测量结果为0时, 忽略寄存器R和C, 剩余部分就是所需要的低秩近似矩阵 A′对应的量子态 | ψA′〉 :

现在分析这一算法的复杂度.利用QRAM制备初始态 | ψ0〉 , 需要的时间是 O [log(pq)].相位估计算法中, 一般有 t0=O(κ/ε).考虑到目标是奇异值大于 τ 的部分, 则 O (κ)=O(σmax/τ).对于低秩近似问题, 通常要求近似后的矩阵和初始矩阵差别很小, 即通常不随矩阵维度呈多项式增长.相位估计算法的复杂度一般是而该算法中相位估计算法中 ε 的大小对于远离阈值 τ 的奇异值几乎不产生影响, 只会对接近 τ 的部分产生误差, 因此常数误差不会对大多数被过滤的成分产生影响, 不影响算法输出的低秩性.综上所述, 该量子低秩近似算法的复杂度为 O [log(pq)].作为对比, 经典计算机解决基于奇异值分解的低秩近似问题需要的时间是 O [poly(pq)].本量子算法的复杂度有指数提高.

2.3 代替相位估计算法的共振跃迁量子算法

本文的量子低秩近似算法, 在实验上实现有以下困难.首先是制备初始态 | ψA〉.在有QRAM的情况下, 可以高效的制备初始态, 但目前为止还没有真正实现QRAM功能的量子器件, 因此要选择其他方式制备.其次是相位估计算法需要一个额外的寄存器存放本征值的估计值, 需要大量的量子比特, 这也是所有基于相位估计的算法的共同问题,使得目前在硬件上实现该类算法非常困难; 除此,实现密度算符 ρAA†和演化算符也需要大量的辅助比特.在设计算法时, 通常假设拥有足够多的比特资源, 只考虑时间复杂度的优势, 但目前量子计算机仍在研发当中, 很难提供足够多的比特资源.

针对以上问题, 分别提出了对应的解决方案:1)利用梯度下降优化脉冲(GRAPE)[28,29], 将量子态从赝纯态演化到初始态 | ψA〉 , 算法直接从初始态|ψA〉开始, 而制备该初始态的方式不会影响到本文算法的正确性.2)为了减少量子比特数, 选择量子共振跃迁算法来替代相位估计算法以减少辅助比特的数量, 只用一个辅助比特即可实现这一步骤.

相位估计在该问题中的作用是先区分目标奇异向量(需要保留的)和非目标奇异向量(不需要保留的), 以便之后通过测量进行筛选, 量子共振跃迁算法也可以实现对特定的量子态的选择.这里,对量子共振跃迁算法简单介绍(图1).

图1 共振跃迁.一个本征值未知的系统与另一个二能级辅助系统存在相互作用.H s 是左边系统的哈密顿量, 假设E0=0.当辅助系统的激发态能级 ε 接近任意特征值Ei时, 两个系统之间会发生共振跃迁Fig.1.Resonant transition.A system with unknown eigenvalues interacts with another two-level auxiliary system.Hs is the Hamiltonian of the system on the left.Assuming E0=0 , when the excited state energy level ε of the auxiliary system is close to any eigenvalue E i , a resonance transition will occur between the two systems.

在两个原本独立的量子系统中, 产生了比较微弱(相比系统本身的能量)的相互作用, 例如两个分子距离逐渐靠近等, 若两个系统的一些能级能量接近, 则部分能量会在两个系统的这些能级上跃迁.基于这一物理现象, Wang[30]设计了一种算法,并在核磁量子计算平台上成功计算了2比特低能等效水分子哈密顿量的基态[31].该算法的步骤如下.

1) 将量子态初始化到 | 0〉⊗|Φ〉.

2) 模拟哈密顿量 H 并进行时间演化, H 具体形式由下式给出:

其中, σ 是泡利算符, c是两个系统的相互作用强度, B 是相互算作算符, 需要求解的系统哈密顿量 为 Hs, 在 HT中 引 入 参 考 点,HT=ε0|0〉〈0|⊗|Φ〉〈Φ|+|1〉〈1|⊗Hs, 其中 | Φ〉 是初始态.

3) 演化一段时间后测量第1个辅助比特, 若发生共振, 则测量结果可能为1或0.当结果为1时得到目标态.结果为0时重复上一步, 再测量,直到测量结果为1.

理论上共振跃迁算法可以替代相位估计算法,利用较少的辅助比特实现目标.以图2(a)的校徽为例, 利用共振跃迁算法进行数值模拟, 结果如图2(b)—(d)所示.

图2 共振跃迁量子算法数值模拟 (a) 3 00×300 的原图像; (b), (c), (d)分别是前10, 30, 50个奇异值恢复的图像.50个奇异值已经可以很好地恢复图像, 可以辨认校徽的细节、文字Fig.2.Numerical simulation of resonance transition by quantum algorithm.Panel (a) is the original image of 300×300.Panels (b), (c) and (d) are the images recovered by the first 10, 30 and 50 singular values respectively.50 singular values can be a good restoration of the image, the details and words of the school emblem are legible.

从数值模拟的结果来看, 共振跃迁量子算法可以很好地实现预期目标.

3 实验演示及结果分析

3.1 实验参数和步骤

为了减少量子比特数, 首先选择1个 8 ×8 维的图像矩阵, 并预先进行厄米化, 表示该矩阵需要的量子比特数是3.之后利用共振跃迁算法, 只需要1个辅助比特, 共计需要4量子比特, 可以在核磁共振量子计算平台进行演示.原始矩阵 M 是由图3所示数字1的灰度图像在Matlab中对应的8×8维实矩阵:

图3 数字1的灰度图像, 共 8 ×8 个像素Fig.3.A grayscale image of the number 1, with a total of 8×8pixels.

考虑到辅助比特只有1个, 令阈值τ=(σ1+σ2)/2并对 M 厄米化后得到 A.

本次实验平台是核磁共振量子计算平台, 实验中使用13C 标记的巴豆酸作为4比特样品, 溶解在d6-丙酮中, 整个过程中1H 解耦.该分子的结构和参数如图4所示.C1—C4表示4个量子位, 选择C1作为辅助量子位.弱耦合近似下的内部哈密顿量为图中 vi为化学位移, Jij为第i和第j个核之间的耦合强度.实验在室温(T = 296.5 K)下的Bruker DRX 400-MHz谱仪上进行.实验过程如下.

图4 1 3C 标记巴豆酸的分子结构和参数.C1是辅助量子比特, C2—C4是目标量子比特.在整个实验中, 1 H 是解耦的.对角元素和非对角元素是化学位移和J耦合(单位:Hz).最下面是相干弛豫时间 T 2 (单位: s).Fig.4.Molecular structure and parameters of crotonic acid are labeled 1 3C.C1 is the auxiliary qubit and C2-C4 is the target qubit.1 H is decoupled throughout the experiment.The diagonal and off-diagonal elements are chemical shifts and J coupling (in Hertz).At the bottom is the coherent relaxation time T 2 (in seconds).

第1步制备赝纯态.这里使用的是空间平均法[32-34], 即利用梯度磁场和多个幺正演化算符实现, 当然也可以采用其他方法[35,36].完成后密度算符有如下形式:

其中单位矩阵I不产生有效信号, 极化度 ϵ ≈10-5.制备赝纯态后, 利用GRAPE将量子态从 | 0000〉 演化到需要的态 | 0〉⊗|ψA〉 , 完成初始态的制备.

第2步实现共振跃迁算法的哈密顿量 H 的时间演化.利用形状脉冲能够实现时间演化算符e-iHt0, 其中 H 由(17)式给出, 考虑到这里的目标态是初始态的一个很好的近似, 令 c = 0.01 ,B=I⊗3, Hs=A , ε =1 并且 ε0=σ1-1.将演化算符多次作用在初始态上, 也就是让 | 0〉⊗ | ψA〉 在H下演化一段时间.

最后, 利用量子态层析方法对密度矩阵进行读出[37-41], 需要用到多个读出脉冲分别读出密度矩阵的不同成分, 最后对密度矩阵重构.

本实验的简要流程如图5所示.

图5 简要量子电路图.制备好赝纯态后, 将经过GRAPE优化好的形状脉冲作用于第2个寄存器, 完成初始化.之后将时间演化算符 e -iHt0 作用于初始态N次, 再测量第1个比特, 结果为 | 1〉 时 得到目标态|ψA′〉Fig.5.A brief quantum circuit diagram.After the pseudomorphic state is prepared, the shape pulse optimized by GRAPE is applied to the second register to complete the initialization.Then, the time evolution operator e -iHt0 is applied to the initial state for N times, and the first bit is measured.When the result is | 1〉 , the target state | ψA′〉 is obtained.

3.2 实验结果及分析

对实验读出数据进行处理并重构密度矩阵后得到的结果如图6所示, 另一个奇异向量可以通过完全相同的方法获得, 只需要交换厄米化时矩阵M 与 M†的顺序.通过第1个奇异值对应的奇异向量对图像进行重构, 从图6能看出实验结果和理论值符合得很好.

图6 (a)通过实验得到的密度矩阵; (b)理论计算的结果; (c)最大的奇异值恢复的矩阵对应图像; 由于辅助比特是 | 1〉 时才给出有意义的结果, 所以这里只考虑其为 | 1〉 的子空间, 忽略其余部分Fig.6.(a) Density matrix obtained by experiment; (b) the result of theoretical calculation; (c) the corresponding image of the matrix with the maximum singular value recovery.Since the auxiliary bit is | 1〉 and gives a meaningful result, only consider subspaces whose values are | 1〉 and the rest are ignored.

通过计算保真度F来定量判断实验结果的准确性:

其中 ρ 和 σ 分别是实验值和理论值, 得到保真度为99.84%.

实验结果与理论计算的误差主要来源于两方面.一方面是核磁共振谱仪有一定的涨落, 导致脉冲会有一定的误差, 这是由于实验设备造成的.除此之外, 在制备初始态、实现时间演化算符时, 均利用了梯度下降算法优化脉冲, 每个优化后的脉冲与目标存在很小的误差.综合以上两点, 实验结果在误差范围内验证了本文提出的量子矩阵低秩近似算法的正确性.

4 总 结

本文提出了一个基于奇异值分解的矩阵低秩近似量子算法.对于矩阵 A ∈Rp×q, 求解其低秩矩阵 A′∈Rp×q, 在经典计算机上计算耗时一般是O[poly(pq)], 而本文的量子算法, 在量子计算机上只需要耗时 O [log(pq)] 即可解决该问题, 对比经典计算机有指数加速.

本文第2节中提出的算法用到了相位估计, 在比特数较多的情况下, 量子主成分分析、量子推荐算法等量子算法都可以实现同样的目标, 在时间复杂度上基本相同.实验中利用了共振跃迁的算法,在这一类问题中等效替代了相位估计算法, 大大减少了对辅助比特的需求, 只通过1个辅助比特, 保留1个奇异值就较好地恢复了图像, 这是其他基于相位估计的算法目前不能实现的.

以数字1的灰度图像在Matlab中对应的矩阵为原始矩阵, 在核磁量子计算平台求解了该矩阵的低秩近似并展示了近似后的矩阵对应的图像.实验结果和理论值相比保真度为99.84%, 充分证明了本文算法的正确性.

在降维、图像处理等多个机器学习相关领域,矩阵的低秩近似问题具有重要意义.本文提出的量子算法只需要少量的辅助量子比特即可在O[log(pq)]时间内给出一个矩阵的低秩近似, 在解决经典算法耗时久的问题的同时大大减少了对量子比特资源的需求, 对量子计算机的实际应用有重要的意义.

感谢清华大学物理系龙桂鲁教授、南京师范大学计算机与电子信息学院段博佳老师提供的宝贵意见.