QSVM阅读笔记

P.Rebentrost等人在2014年提出的Quantum support vector machine(QSVM) 利用量子算法构造和训练SVM模型，并利用该模型对新的样本进行分类。QSVM的复杂度对比经典SVM有指数级别的缩小，因此适用于大规模数据的训练。QSVM的主体思想是将传统SVM这一二次规划问题转换为最小二乘问题，利用量子算法高效地求解最小二乘问题中的向量内积，基于HHL算法求解模型的参数，使用Swap-test算法对新的样本进行分类。

一、支持向量机(SVM)

SVM是定义在特征空间上的间隔最大的线性分类器，对于线性可分的数据集，SVM相当于求解凸二次规划问题；对于线性不可分的数据集，支持向量机引入了核技巧(kernel trick)，将数据集从低维线性不可分的输入空间映射到高维线性可分的特征空间，并在特征空间学习SVM线性分类器。因此，支持向量机既可以用来对线性可分的数据集进行分类，也可以对线性不可分的数据集进行分类，只是后者引入了核函数。

支持向量机的具体介绍可参阅李航的《统计学习方法》，这里只给出它的数学表达：
假设数据集${\left\{\left({\vec{x}}_j,y_j\right):{\vec{x}}_j\in {\mathbb{R}}^N,y_j=\pm 1\right\}}_{j=1\dots M}$中有$M$个带正负标签的样本，每个样本有$N$个特征，那么线性可分SVM的对偶最优化问题为，
$$\begin{array}{ll}\underset{\alpha }{\min }& \frac{1}{2}{\sum }_{i=1}^M{\sum }_{j=1}^M{\alpha }_i{\alpha }_j{y}_i{y}_j\left({\vec{x}}_i,{\vec{x}}_j\right)-{\sum }_{i=1}^M{\alpha }_i\\ \text{ s.t. }& {\sum }_{i=1}^M{\alpha }_i{y}_i=0\\ & {\alpha }_i⩾0,i=1,2,\cdots ,M\end{array}$$其中$({\vec{x}}_i,{\vec{x}}_j)$表示${\vec{x}}_i$和${\vec{x}}_j$的内积${\vec{x}}_j^T{\vec{x}}_i$。求解该最优化问题得到最优解${\alpha }^{\ast }={\left({\alpha }_1^{\ast },{\alpha }_2^{\ast },\dots ,{\alpha }_M^{\ast }\right)}^T$，进而得到SVM线性分类器
$$f(x)=\text{sign}\left(\sum _{i=1}^M{\alpha }_i^{\ast }{y}_i\left(\vec{x},{\vec{x}}_i\right)+b^{\ast }\right)$$其中$b^{\ast }=y_j-{\sum }_{i=1}^M{\alpha }_i^{\ast }{y}_i\left({\vec{x}}_i,{\vec{x}}_j\right)$且${\alpha }_i^{\ast }>0$。

论文中的数学表述与上式略有不同，在论文中作者令${\alpha }_i={\alpha }_i{y}_i$，并且将使目标函数最小化改为使目标函数最大化，因此线性可分SVM的对偶最优化问题变为，
$$\begin{array}{ll}\underset{\alpha }{\max }& {\sum }_{i=1}^M{\alpha }_i{y}_i-\frac{1}{2}{\sum }_{i=1}^M{\sum }_{j=1}^M{\alpha }_i{\alpha }_j\left({\vec{x}}_i,{\vec{x}}_j\right)\\ \text{ s.t. }& {\sum }_{i=1}^M{\alpha }_i=0\\ & {\alpha }_i{y}_i⩾0,i=1,2,\cdots ,M\end{array}$$SVM线性分类器变为
$$f(x)=\text{sign}\left(\sum _{i=1}^M{\alpha }_i^{\ast }\left(\vec{x},{\vec{x}}_i\right)+b^{\ast }\right)$$其中$b^{\ast }=y_j-{\sum }_{i=1}^M{\alpha }_i^{\ast }\left({\vec{x}}_i,{\vec{x}}_j\right)$且${\alpha }_i^{\ast }\ne 0$。

非线性支持向量机引入了核函数$K(\vec{x},\vec{z})={\varphi }^T(\vec{x})\varphi (\vec{z})$，$\varphi (\cdot )$是输入空间到特征空间的映射。将线性SVM的对偶问题以及分类器中的向量内积替换为核函数$K$便可得到非线性SVM的对偶最优化问题和分类器。所以，非线性SVM的对偶最优化问题的目标函数为
$$\begin{array}{ll}\underset{\alpha }{\max }& {\sum }_{i=1}^M{\alpha }_i{y}_i-\frac{1}{2}{\sum }_{i=1}^M{\sum }_{j=1}^M{\alpha }_i{\alpha }_{j}K({\vec{x}}_j,{\vec{x}}_i)\end{array}$$非线性SVM分类器为
$$f(x)=\text{sign}\left(\sum _{i=1}^M{\alpha }_i^{\ast }K({\vec{x}}_i,\vec{x})+b^{\ast }\right)$$注意，论文中作者虽然把由输入向量的内积构成的矩阵$K$称为kernel matrix，但并不是指的非线性SVM中的kernel。因此，作者提出的QSVM主要是将线性SVM用量子算法来实现，而对于非线性SVM的量子化，作者并没有做详细讨论。

二、最小二乘支持向量机(LSSVM)

Aram W. Harrow等人提出的HHL算法旨在利用量子计算的方法求解线性方程$A\vec{x}=\vec{b}$，而QSVM要利用HHL算法，首先需要把传统SVM的二次规划形式转换为$A\vec{x}=\vec{b}$的这种形式。实现这种转化的方法便是引入松弛变量$e_j$，将原始问题的不等式约束转换为等式约束
$$y_j\left(\vec{w}\cdot {\vec{x}}_j+b\right)\ge 1\to \left(\vec{w}\cdot {\vec{x}}_j+b\right)=y_j-y_j{e}_j$$，引入松弛变量后的原问题的对偶问题近似为
$$F\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)=\left(\begin{array}{cc}0& {\vec{1}}^T\\ \vec{1}& K+{\gamma }^{-1}I\end{array}\right)\left(\begin{array}{l}b\\ \vec{\alpha }\end{array}\right)=\left(\begin{array}{l}0\\ \vec{y}\end{array}\right)$$其中$\vec{1}=(1,\cdots ,1{)}^T$，$\vec{\alpha }={\left({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_M\right)}^T$，$K=\left[\begin{array}{cccc}({\vec{x}}_1,{\vec{x}}_1)& ({\vec{x}}_1,{\vec{x}}_2)& \cdots & ({\vec{x}}_1,{\vec{x}}_M)\\ ({\vec{x}}_2,{\vec{x}}_1)& ({\vec{x}}_2,{\vec{x}}_2)& \cdots & ({\vec{x}}_2,{\vec{x}}_M)\\ ⋮& ⋮& \ddots & ⋮\\ ({\vec{x}}_M,{\vec{x}}_1)& ({\vec{x}}_M,{\vec{x}}_2)& \cdots & ({\vec{x}}_M,{\vec{x}}_M)\end{array}\right]$，$\vec{y}=(y_1,y_2,\cdots ,y_M{)}^T$。
$b$和$\vec{\alpha }$即为待求参数，作者将利用HHL算法求解
$$\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)=F^{-1}\left(\begin{array}{l}0\\ \vec{y}\end{array}\right).$$

三、QSVM的训练

1. 非稀疏矩阵$\hat{F}$的指数化

回顾HHL算法，用量子算法求解矩阵$A$的逆，需要将$A$指数化为酉算子$e^{-iA\frac{t_0}{T}}$，并且要求
①$A^{†}=A$；
②$A$是稀疏的；
③$A$的特征值都大于0小于1，特征值越接近于0，指数化的效率越低。
$F$显然满足$F^{†}=F$，且$F$是正定的，令$\hat{F}=F/\mathrm{tr}F$，$\hat{F}$的特征值将处于0到1之间，但$F$一般不满足稀疏这一条件。对于非稀疏的埃尔米特矩阵的指数化，将采用一种近似的量子算法进行逼近。
先对$\hat{F}$进行一下拆分，$\hat{F}=\left(J+K+{\gamma }^{-1}I\right)/\mathrm{tr}F$，其中$J=\left(\begin{array}{cc}0& {\vec{1}}^T\\ \vec{1}& 0\end{array}\right)$，因此，
$$e^{-i\hat{F}\Delta t}=e^{-i\Delta t{\gamma }^{-1}I/\mathrm{t}\mathrm{r}\mathrm{F}}{e}^{-i\Delta tJ/\mathrm{t}\mathrm{r}\mathrm{F}}{e}^{-i\Delta tK/\mathrm{t}\mathrm{r}\mathrm{F}}+O\left(\Delta t^2\right)$$$\mathrm{tr}F$和$\mathrm{tr}K$是两个常数，不妨稍微改变一下HHL算法中quantum phase estimation的时间$\Delta t$，上式变为（后文中为了方便起见将$\Delta t^{\prime }$写作$\Delta t$）
$$e^{-i\hat{F}\Delta t^{\prime }}=e^{-i\Delta t^{\prime }{\gamma }^{-1}I/\mathrm{t}\mathrm{r}\mathrm{K}}{e}^{-i\Delta t^{\prime }J/\mathrm{t}\mathrm{r}\mathrm{K}}{e}^{-i\Delta t^{\prime }K/\mathrm{t}\mathrm{r}\mathrm{K}}+O\left(\Delta t^{\prime 2}\right)$$上式中无论${\gamma }^{-1}I/\mathrm{t}\mathrm{r}\mathrm{K}$还是$J/\mathrm{t}\mathrm{r}\mathrm{K}$都满足HHL算法的三条要求，因此都是可以高效地指数化的，只有$K/\mathrm{t}\mathrm{r}K$不满足第二条要求。所以如何高效地对$K/\mathrm{t}\mathrm{r}K$进行指数化，构造为酉算子$e^{-i\Delta tK/\mathrm{t}\mathrm{r}K}$是一个核心问题。

论文中作者采用了QPCA中的思想，近似地对非稀疏的埃尔米特矩阵进行指数化，$O\left(\Delta t^2\right)$就是近似误差。

令$\hat{K}=K/\mathrm{t}\mathrm{r}K$，显然$\hat{K}$的$trace$等于1且正定，因此可以把$\hat{K}$视density operator或density matrix，以下量子电路将产生状态$\hat{K}$，
Training data可以以量子态形式存储在quantum RAM中，
$$|{\vec{x}}_i⟩=1/\left|{\vec{x}}_i\right|\sum _{k=1}^N{\left({\vec{x}}_i\right)}_k|k⟩$$量子态的Training data和事先准备的量子态$\frac{1}{\sqrt{M}}{\sum }_{i=1}^M|i⟩$经过Oracle的作用后，形成的叠加态为
$$|\chi ⟩=1/\sqrt{N_{\chi }}\sum _{i=1}^M\left|{\vec{x}}_i\right||i⟩|{\vec{x}}_i⟩$$其中$N_{\chi }={\sum }_{i=1}^M{\left|{\vec{x}}_i\right|}^2$。利用density matrix描述量子系统状态，舍去Training data寄存器中的量子态，系统的状态为
ρ 1 = tr ⁡ 2 { ∣ χ ⟩ ⟨ χ ∣ } = tr ⁡ 2 { 1 N χ ∑ i , j = 1 M ∣ x ⃗ i ∣ ∣ x ⃗ j ∣ ∣ i ⟩ ∣ x ⃗ i ⟩ ⟨ j ∣ ⟨ x ⃗ j ∣ } = 1 N χ ∑ i , j = 1 M ⟨ x ⃗ j ∣ x ⃗ i ⟩ ∣ x ⃗ i ∣ ∣ x ⃗ j ∣ ∣ i ⟩ ⟨ j ∣ = K / t r K = K ^ \rho _{1}=\operatorname{tr}_{2}\{|\chi\rangle\langle\chi|\}=\operatorname{tr}_{2}\left\{\frac{1}{N_{\chi}} \sum_{i, j=1}^{M}|\vec{x}_{i}||\vec{x}_{j}|| i\rangle|\vec{x}_{i}\rangle\langle j|\langle\vec{x}_{j}|\right\}=\frac{1}{N_{\chi}} \sum_{i, j=1}^{M}\left\langle\vec{x}_{j}|\vec{x}_{i}\right\rangle|\vec{x}_{i}||\vec{x}_{j}|| i\rangle\langle j| = K/\mathrm{tr}K=\hat{K} ρ1​=tr2​{∣χ⟩⟨χ∣}=tr2​{Nχ​1​i,j=1∑M​∣x i​∣∣x j​∣∣i⟩∣x i​⟩⟨j∣⟨x j​∣}=Nχ​1​i,j=1∑M​⟨x j​∣x i​⟩∣x i​∣∣x j​∣∣i⟩⟨j∣=K/trK=K^这一步的计算复杂度为$O(\log MN)$，可以理解为$M$个$N$维向量${\vec{x}}_i$可以由$\log MN$个量子进行表示。

基于QPCA中的思想，非稀疏的埃尔米特矩阵$\hat{K}$的指数化形成的酉算子$e^{-i\hat{K}\Delta t}$对某一状态$\rho$的作用，近似于swap matrix S（$S={\sum }_{m,n=1}^M|m⟩⟨n|\otimes |n⟩⟨m|$）的指数化形成的酉算子$e^{-iS\Delta t}$对$\hat{K}$和$\rho$形成的组合量子系统作用后，并舍弃$\hat{K}$变化后的状态而余下的状态，图解如上图所示。该部分的证明较为复杂，将在文章最后给出。由于swap matrix S的指数化是较为高效的，因此该方法在误差允许范围内实现了非稀疏埃尔米特矩阵的指数化。

2. 基于HHL算法的模型参数求解

至此，矩阵$\hat{F}$可以高效的指数化为酉算子$e^{-i\hat{F}\Delta t}$，因此可以利用HHL算法求解$\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)=F^{-1}\left(\begin{array}{l}0\\ \vec{y}\end{array}\right)$。对应的QSVM参数训练量子电路图如下图所示，

这里简要推导一下基于HHL算法怎样得到模型参数$\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)$。

假设矩阵$\hat{F}$的特征值为${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{M+1}$，对应的单位正交的特征向量为$|u_1⟩,|u_2⟩,\dots ,|u_{M+1}⟩$，这$M+1$个特征向量构成M+1维线性空间的一组标准正交基，将$|\tilde{y}⟩=\left(\begin{array}{c}0\\ \vec{y}\end{array}\right)$在这组基下分解得到$|\tilde{y}⟩={\sum }_{j=1}^{M+1}⟨u_j|\tilde{y}⟩|u_j⟩$，经过绿线框中的Improved quantum phase estimation，得到${\sum }_{j=1}^{M+1}⟨u_j|\tilde{y}⟩|{\lambda }_j⟩|u_j⟩$；引入ancilla qubit $|0⟩$，执行Controlled rotation，得到${\sum }_{j=1}^{M+1}⟨u_j|\tilde{y}⟩\left(\sqrt{1-\frac{1^2}{{\lambda }_j^2}}|0⟩+\frac{1}{{\lambda }_j}|1⟩\right)|{\lambda }_j⟩|u_j⟩$；对第一个寄存器中的量子态进行测量，当测量得到的量子态为$|1⟩$时，忽略第二个寄存器中的量子态，最终得到${\sum }_{j=1}^{M+1}\frac{1}{{\lambda }_j}⟨u_j|\tilde{y}⟩|u_j⟩$，即$\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)$。

注意，上述推导过程忽略了HHL算法中要求$|\tilde{y}⟩$为单位向量以及量子态的系数，导致某些量子态不是单位化的，并且HHL算法最终求解得出的结果是$\left(\begin{array}{c}b\\ \vec{\alpha }\end{array}\right)$的单位化。所以，最终得到的结果可以写为
$$|b,\vec{\alpha }⟩=\frac{1}{\sqrt{C}}\left(b|0⟩+\sum _{k=1}^M{\alpha }_k|k⟩\right)$$其中$C=b^2+{\sum }_{k=1}^M{\alpha }_k^2$。

四、QSVM分类器

如何利用以量子态形式存储的模型参数$|b,\vec{\alpha }⟩$对新的样本进行预测（分类）呢？作者利用了Swap-test算法。

利用Oracle构造两个量子态$|\tilde{u}⟩$和$|\tilde{x}⟩$,
$$|\tilde{u}⟩=\frac{1}{\sqrt{N_{\tilde{u}}}}\left(b|0⟩|0⟩+\sum _{k=1}^M{\alpha }_k\left|{\vec{x}}_k\right||k⟩|{\vec{x}}_k⟩\right),N_{\tilde{u}}=b^2+\sum _{k=1}^M{\alpha }_k^2{\left|{\vec{x}}_k\right|}^2$$$$|\tilde{x}⟩=\frac{1}{\sqrt{N_{\tilde{x}}}}\left(|0⟩|0⟩+\sum _{k=1}^M|\vec{x}||k⟩|\vec{x}⟩\right),N_{\tilde{x}}=M|\vec{x}{|}^2+1$$对这两个量子执行Swap-test，测得ancilla qubit的量子态为$|1⟩$的概率为$P=\frac{1}{2}(1-|⟨\tilde{u}|\tilde{x}⟩{|}^2)$，其中
$$⟨\tilde{u}|\tilde{x}⟩=\frac{1}{\sqrt{N_{\tilde{x}}{N}_{\tilde{u}}}}\left(b+\sum _{k=1}^M{\alpha }_k\left|{\vec{x}}_k\right||\vec{x}|⟨{\vec{x}}_k|\vec{x}⟩\right)$$上式与线性SVM的分类决策函数$f(x)=\text{sign}\left({\sum }_{i=1}^M{\alpha }_i^{\ast }\left(\vec{x},{\vec{x}}_i\right)+b^{\ast }\right)$只是系数上的不同。因此，通过多次运行Swap-test算法，若测得$P<\frac{1}{2}$，说明新样本$\vec{x}$属于正类；若测得$P>\frac{1}{2}$，说明新样本$\vec{x}$属于负类。为了使最后的误差小于$ϵ$，需要运行Swap-test算法$O\left(P(1-P)/{ϵ}^2\right)$次。

参考文献
[1] Rebentrost P, Mohseni M, Lloyd S. Quantum support vector machine for big data classification[J]. Physical review letters, 2014, 113(13): 130503.
[2] Lloyd S, Mohseni M, Rebentrost P. Quantum principal component analysis[J]. Nature Physics, 2014, 10(9): 631-633.
[3] Buhrman H, Cleve R, Watrous J, et al. Quantum fingerprinting[J]. Physical Review Letters, 2001, 87(16): 167902.
[4] Harrow A W, Hassidim A, Lloyd S. Quantum algorithm for linear systems of equations[J]. Physical review letters, 2009, 103(15): 150502.
[5] Suykens J A K, Vandewalle J. Least squares support vector machine classifiers[J]. Neural processing letters, 1999, 9(3): 293-300.
[6] 量子支持向量机QSVM（一）
[7] 量子支持向量机QSVM（二）

公式证明

1. 证明$S={\sum }_{m,n=1}^M|m⟩⟨n|\otimes |n⟩⟨m|$能够交换两个量子的状态：
   设有量子$|\alpha ⟩={\sum }_{i=1}^M{a}_i|i⟩$，$|\beta ⟩={\sum }_{j=1}^M{b}_j|j⟩$，
   $$\begin{array}{rl}S\left(|\alpha ⟩\otimes |\beta ⟩\right)& =\sum _{m,n=1}^M|m⟩⟨n|\otimes |n⟩⟨m|\left(\sum _{i,j=1}^M{a}_i|i⟩\otimes b_j|j⟩\right)\\ & =\sum _{m,n=1}^M\sum _{i,j=1}^M{a}_i{b}_j\left(|m⟩⟨n|i⟩\right)\otimes \left(|n⟩⟨m|j⟩\right)\\ & =\sum _{m,n=1}^M{b}_m|m⟩\otimes a_n|n⟩\\ & =|\beta ⟩\otimes |\alpha ⟩.\end{array}$$
2. 显然有$S^{†}=S$，所以$⟨\alpha |⟨\beta |S={\left(S|\alpha ⟩|\beta ⟩\right)}^{†}={\left(|\beta ⟩|\alpha ⟩\right)}^{†}=⟨\beta |⟨\alpha |$.
3. 因为$S^2\left(|\alpha ⟩\otimes |\beta ⟩\right)=S\left(|\beta ⟩\otimes |\alpha ⟩\right)=|\alpha ⟩\otimes |\beta ⟩$，所以$S^2=I$.
4. $e^{-iS\Delta t}=\cos (\Delta tS)-i\sin (\Delta tS)=I\cos \Delta t-iS\sin \Delta t$.
5. 证明$e^{-i\hat{K}\Delta t}\rho e^{i\hat{K}\Delta t}\approx {\mathrm{tr}}_1\left\{e^{-iS\Delta t}\hat{K}\otimes \rho e^{iS\Delta t}\right\}=\rho -i\Delta t[\hat{K},\rho ]+O\left(\Delta t^2\right)$.
   设$\hat{K}={\sum }_{i=1}^M{\lambda }_i|k_i⟩⟨k_i|$，$\rho ={\sum }_{j=1}^M{p}_j|{\psi }_j⟩⟨{\psi }_j|$，其中${\sum }_{i=1}^M{\lambda }_i=1$，${\sum }_{j=1}^M{p}_j=1$。注意，density matrix的展开并不一定是谱分解，因此${\lambda }_i$和$p_j$并不一定是特征值，$|k_i⟩$和$|{\psi }_j⟩$并不一定是标准正交的特征向量。
   $$\begin{array}{rl}{e}^{-i\hat{K}\Delta t}\rho e^{i\hat{K}\Delta t}& \approx \left(I-i\hat{K}\Delta t\right)\rho \left(I+i\hat{K}\Delta t\right)\\ & =\sum _{j=1}^M\left[\left(I-i\Delta t\sum _{i=1}^M{\lambda }_i|k_i⟩⟨k_i|\right)p_j|{\psi }_j⟩⟨{\psi }_j|\left(I+i\Delta t\sum _{i=1}^M{\lambda }_i|k_i⟩⟨k_i|\right)\right]\\ & =\sum _{j=1}^M\left(p_j|{\psi }_j⟩⟨{\psi }_j|-i\Delta t\sum _{i=1}^M{\lambda }_i{p}_j⟨k_i|{\psi }_j⟩|k_i⟩⟨{\psi }_j|+i\Delta t\sum _{i=1}^M{\lambda }_i{p}_j⟨{\psi }_j|k_i⟩|{\psi }_j⟩⟨k_i|+\Delta t^2\sum _{i,\gamma =1}^M{p}_j{\lambda }_i{\lambda }_{\gamma }{\left|⟨k_i|{\psi }_j⟩\right|}^2|k_i⟩⟨k_{\gamma }|\right)\\ & =\rho -i\Delta t\left(\sum _{i,j=1}^M{\lambda }_i{p}_j⟨k_i|{\psi }_j⟩|k_i⟩⟨{\psi }_j|-\sum _{i,j=1}^M{\lambda }_i{p}_j⟨{\psi }_j|k_i⟩|{\psi }_j⟩⟨k_i|\right)+\Delta t^2\sum _{i,j,\gamma =1}^M{p}_j{\lambda }_i{\lambda }_{\gamma }{\left|⟨k_i|{\psi }_j⟩\right|}^2|k_i⟩⟨k_{\gamma }|\\ & =\rho -i\Delta t\left[\hat{K},\rho \right]+O\left(\Delta t^2\right)\end{array}$$ tr ⁡ 1 { e − i S Δ t K ^ ⊗ ρ e i S Δ t } = tr ⁡ 1 { ( I cos ⁡ Δ t − i S sin ⁡ Δ t ) ( ∑ i , j = 1 M λ i p j ∣ k i ⟩ ⟨ k i ∣ ⊗ ∣ ψ j ⟩ ⟨ ψ j ∣ ) ( I cos ⁡ Δ t + i S sin ⁡ Δ t ) } = ∑ i , j = 1 M λ i p j tr ⁡ 1 { ( I cos ⁡ Δ t − i S sin ⁡ Δ t ) ( ∣ k i ⟩ ∣ ψ j ⟩ ⊗ ⟨ k i ∣ ⟨ ψ j ∣ ) ( I cos ⁡ Δ t + i S sin ⁡ Δ t ) } = ∑ i , j = 1 M λ i p j tr ⁡ 1 { cos ⁡ 2 Δ t ∣ k i ⟩ ∣ ψ j ⟩ ⟨ k i ∣ ⟨ ψ j ∣ − i sin ⁡ Δ t cos ⁡ Δ t ∣ ψ j ⟩ ∣ k i ⟩ ⟨ k i ∣ ⟨ ψ j ∣ + i sin ⁡ Δ t cos ⁡ Δ t ∣ k i ⟩ ∣ ψ j ⟩ ⟨ ψ j ∣ ⟨ k i ∣ + sin ⁡ 2 Δ t ∣ ψ j ⟩ ∣ k i ⟩ ⟨ ψ j ∣ ⟨ k i ∣ } = ∑ i , j = 1 M λ i p j tr ⁡ 1 { cos ⁡ 2 Δ t ∣ k i ⟩ ⟨ k i ∣ ⊗ ∣ ψ j ⟩ ⟨ ψ j ∣ − i sin ⁡ Δ t cos ⁡ Δ t ∣ ψ j ⟩ ⟨ k i ∣ ⊗ ∣ k i ⟩ ⟨ ψ j ∣ + i sin ⁡ Δ t cos ⁡ Δ t ∣ k i ⟩ ⟨ ψ j ∣ ⊗ ∣ ψ j ⟩ ⟨ k i ∣ + sin ⁡ 2 Δ t ∣ ψ j ⟩ ⟨ ψ j ∣ ⊗ ∣ k i ⟩ ⟨ k i ∣ } = ∑ i , j = 1 M λ i p j { cos ⁡ 2 Δ t ∣ ψ j ⟩ ⟨ ψ j ∣ − i sin ⁡ Δ t cos ⁡ Δ t ( ⟨ k i ∣ ψ j ⟩ ∣ k i ⟩ ⟨ ψ j ∣ − ⟨ ψ j ∣ k i ⟩ ∣ ψ j ⟩ ⟨ k i ∣ ) + sin ⁡ 2 Δ t ∣ k i ⟩ ⟨ k i ∣ } = ρ cos ⁡ 2 Δ t − i sin ⁡ Δ t cos ⁡ Δ t ∑ i , j = 1 M λ i p j ( ⟨ k i ∣ ψ j ⟩ ∣ k i ⟩ ⟨ ψ j ∣ − ⟨ ψ j ∣ k i ⟩ ∣ ψ j ⟩ ⟨ k i ∣ ) + K ^ sin ⁡ 2 Δ t = ρ − i sin ⁡ Δ t cos ⁡ Δ t [ K ^ , ρ ] + O ( Δ t 2 ) \begin{aligned} \operatorname{tr}_{1}\left\{e^{-i S \Delta t} \hat{K} \otimes \rho e^{i S \Delta t}\right\} &= \operatorname{tr}_{1}\left\{\left(I\cos\Delta t-iS\sin\Delta t\right)\left(\sum_{i,j=1}^{M}\lambda_{i}p_{j}|k_{i}\rangle\langle k_{i}|\otimes|\psi_{j}\rangle\langle \psi_{j}|\right)\left(I\cos\Delta t+iS\sin\Delta t\right)\right\} \\ &= \sum_{i,j=1}^{M}\lambda_{i}p_{j}\operatorname{tr}_{1}\left\{\left(I\cos\Delta t-iS\sin\Delta t\right)\left(|k_{i}\rangle|\psi_{j}\rangle\otimes\langle k_{i}|\langle \psi_{j}|\right)\left(I\cos\Delta t+iS\sin\Delta t\right)\right\} \\ &= \sum_{i,j=1}^{M}\lambda_{i}p_{j}\operatorname{tr}_{1}\left\{\cos^{2}\Delta t|k_{i}\rangle|\psi_{j}\rangle\langle k_{i}|\langle \psi_{j}|-i\sin\Delta t\cos\Delta t|\psi_{j}\rangle|k_{i}\rangle\langle k_{i}|\langle \psi_{j}|+i\sin\Delta t\cos\Delta t|k_{i}\rangle|\psi_{j}\rangle\langle \psi_{j}|\langle k_{i}|+\sin^{2}\Delta t|\psi_{j}\rangle|k_{i}\rangle\langle \psi_{j}|\langle k_{i}|\right\} \\ &= \sum_{i,j=1}^{M}\lambda_{i}p_{j}\operatorname{tr}_{1}\left\{\cos^{2}\Delta t|k_{i}\rangle\langle k_{i}|\otimes|\psi_{j}\rangle\langle \psi_{j}|-i\sin\Delta t\cos\Delta t|\psi_{j}\rangle\langle k_{i}|\otimes|k_{i}\rangle\langle \psi_{j}|+i\sin\Delta t\cos\Delta t|k_{i}\rangle\langle \psi_{j}|\otimes|\psi_{j}\rangle\langle k_{i}|+\sin^{2}\Delta t|\psi_{j}\rangle\langle \psi_{j}|\otimes|k_{i}\rangle\langle k_{i}|\right\} \\ &= \sum_{i,j=1}^{M}\lambda_{i}p_{j}\left\{\cos^{2}\Delta t|\psi_{j}\rangle\langle \psi_{j}|-i\sin\Delta t\cos\Delta t\left(\langle k_{i}|\psi_{j}\rangle|k_{i}\rangle\langle\psi_{j}|-\langle \psi_{j}|k_{i}\rangle|\psi_{j}\rangle\langle k_{i}|\right)+\sin^{2}\Delta t|k_{i}\rangle\langle k_{i}|\right\} \\ &= \rho\cos^{2}\Delta t-i\sin\Delta t\cos\Delta t \sum_{i,j=1}^{M}\lambda_{i}p_{j}\left(\langle k_{i}|\psi_{j}\rangle|k_{i}\rangle\langle\psi_{j}|-\langle \psi_{j}|k_{i}\rangle|\psi_{j}\rangle\langle k_{i}|\right)+\hat{K}\sin^{2}\Delta t \\ &= \rho-i\sin\Delta t\cos\Delta t\left[\hat{K},\rho\right]+O\left(\Delta t^{2}\right) \end{aligned} tr1​{e−iSΔtK^⊗ρeiSΔt}​=tr1​{(IcosΔt−iSsinΔt)(i,j=1∑M​λi​pj​∣ki​⟩⟨ki​∣⊗∣ψj​⟩⟨ψj​∣)(IcosΔt+iSsinΔt)}=i,j=1∑M​λi​pj​tr1​{(IcosΔt−iSsinΔt)(∣ki​⟩∣ψj​⟩⊗⟨ki​∣⟨ψj​∣)(IcosΔt+iSsinΔt)}=i,j=1∑M​λi​pj​tr1​{cos2Δt∣ki​⟩∣ψj​⟩⟨ki​∣⟨ψj​∣−isinΔtcosΔt∣ψj​⟩∣ki​⟩⟨ki​∣⟨ψj​∣+isinΔtcosΔt∣ki​⟩∣ψj​⟩⟨ψj​∣⟨ki​∣+sin2Δt∣ψj​⟩∣ki​⟩⟨ψj​∣⟨ki​∣}=i,j=1∑M​λi​pj​tr1​{cos2Δt∣ki​⟩⟨ki​∣⊗∣ψj​⟩⟨ψj​∣−isinΔtcosΔt∣ψj​⟩⟨ki​∣⊗∣ki​⟩⟨ψj​∣+isinΔtcosΔt∣ki​⟩⟨ψj​∣⊗∣ψj​⟩⟨ki​∣+sin2Δt∣ψj​⟩⟨ψj​∣⊗∣ki​⟩⟨ki​∣}=i,j=1∑M​λi​pj​{cos2Δt∣ψj​⟩⟨ψj​∣−isinΔtcosΔt(⟨ki​∣ψj​⟩∣ki​⟩⟨ψj​∣−⟨ψj​∣ki​⟩∣ψj​⟩⟨ki​∣)+sin2Δt∣ki​⟩⟨ki​∣}=ρcos2Δt−isinΔtcosΔti,j=1∑M​λi​pj​(⟨ki​∣ψj​⟩∣ki​⟩⟨ψj​∣−⟨ψj​∣ki​⟩∣ψj​⟩⟨ki​∣)+K^sin2Δt=ρ−isinΔtcosΔt[K^,ρ]+O(Δt2)​