八、TEBD 与 DMRG 算法

八、TEBD 与 DMRG 算法

1. TEBD 算法

---

       我们给定一个大小为 $W\times H$ 个张量组成的网络（设 H 为偶数），其图形表示如下：

该张量网络由三种不同的张量构成，在四个角上的为二阶张量，边上的为三阶张量，内部的为四阶张量，边上的数字表示第几个指标。很多物理问题可最终等价为类似的张量网络收缩计算问题，即将该张量网络收缩为一个数。如，经典模型热力学配分函数的计算，量子格点模型基态的计算等。严格收缩这种张量网络的计算，其代价随着 $W$ 和 $H$ 呈指数增大，属于 NP 难问题，因此必须采用近似方法收缩计算。

       TEBD 算法（time-evolving block decimation） 为一种基于矩阵乘积态的近似收缩张量网络的数值算法。其主要思路为：从处于边界的张量构成的 MPS 开始，一行一行地收缩张量网络。

       将张量网络从中间分成两部分，以下半部分为例，下边界的张量实际上组成了一个长度为 $W$ 的 MPS 态，记为 $|{\phi }^{D_0}⟩$ ，整个下半部分的收缩结果可记为另一个 MPS 态 $|{\phi }^D⟩$ ，中间每一层张量构成一个作用在 $|{\phi }^{D_0}⟩$ 上的算符，记为 $\hat{\rho }$ ，那么我们可以用符号表示下半部分的收缩为：
$$|{\phi }^D⟩={\hat{\rho }}^{\frac{H}{2}-1}|{\phi }^{D_0}⟩$$

我们可以用图形表示为：

类似的，张量的上半部分可以相应的定义为 $|{\phi }^{U_0}⟩$ 和 $|{\phi }^U⟩$ ，同样的，我们有：
$$|{\phi }^U⟩={\left({\hat{\rho }}^T\right)}^{\frac{H}{2}-1}|{\phi }^{U_0}⟩$$

那么整个张量网络的收缩结果为：
$$\mathrm{Z}=⟨{\phi }^U\mid {\phi }^D⟩=⟨{\phi }^{U_0}\left|{\hat{\rho }}^{H-2}\right|{\phi }^{D_0}⟩$$

所以，张量网络的 收缩变成了计算 $|{\phi }^U⟩$ 和 $|{\phi }^D⟩$ 。

具体的计算思路为：通过进行下图所示的张量收缩，计算 $(\frac{H}{2}-1)$ 次 $\hat{\rho }$ 对 MPS 的作用， $\hat{\rho }$ 被称为转移矩阵（transfer matrix），这里 $\hat{\rho }$ 具备的张量收缩的形式被称为矩阵乘积算符（matrix product operator , MPO）。

张量的每次收缩都会使 MPS 辅助指标的维数扩大 D 倍（D 为 $\hat{\rho }$ 水平方向指标的维数），也就是说，MPS 辅助指标的维数随收缩次数呈指数增大。所以我们引入 MPS 的最优裁剪，将维数限定为预先设定的截断维数 $\chi$ ，也就是对 MPS 收缩后依次每个张量进行中心正交化，然后对该张量进行奇异值分解，保留前 $\chi$ 个奇异值和奇异向量，保证每个辅助指标的大小都不超过 $\chi$。具体的算符步骤如下：

1. 用处于边界的对应张量初始化 $|{\phi }^{U_0}⟩$ 与 $|{\phi }^{D_0}⟩$
2. 计算 $|{\phi }^{U_{t+1}}⟩={\hat{\rho }}^T|{\phi }^{U_t}⟩,|{\phi }^{D_{t+1}}⟩=\hat{\rho }|{\phi }^{D_t}⟩$ ，得到新的 MPS 态
3. 如果 MPS 辅助指标的维数超过截断维数 $\chi$ ，那么进行最优化裁剪，使辅助指标截断到 $\chi$
4. 进行第二步和第三步 $\frac{H}{2}-1$ 次后，计算 $Z=⟨{\phi }^U|{\phi }^D⟩$

TEBD 算法计算一维格点模型的基态

       我们来考虑 N 个自旋构成的一维海森堡模型：
$$\hat{H}=\sum _{⟨i,j⟩}{\hat{H}}_{ij},{\hat{H}}_{ij}=\sum _{\alpha =x,y,z}{\hat{s}}_i^{\alpha }{\hat{s}}_j^{\alpha }$$

其主要思路为：利用基于 TRotter - Suzuki 分解的退火算法，将随机初始化的 MPS 态演化到基态，我们将一次一次的演化过程转换为张量网络的收缩问题，并使用 TEBD 算法求解该收缩。其具体步骤为

1. 获得演化算符 $e^{-\tau {\hat{H}}_{ij}}$ 的系数张量（$\tau$ 为小的正实数），当哈密顿算符确定后，我们就能用下式计算出每个演化算符
   $$G_{s_1^{\prime }{s}_2^{\prime }{s}_1{s}_2}=⟨s_1^{\prime }{s}_2^{\prime }\left|e^{-\tau {\hat{H}}_{ij}}\right|s_1{s}_2⟩$$
   用图形可以表示为：
   通过计算，我们可以将基态演化算符
   $$\hat{P}=\underset{K\to \infty }{\lim }{\left(e^{-\tau {\sum }_{⟨i,j⟩}{\hat{H}}_{ij}}\right)}^K$$

   表示为由四阶张量 G 构成的张量网络，如下图
   

2. 随机初始化一个 MPS 态，放置于张量底部，整个张量网络变为下图所示
   

3. 将靠近MPS的半层张量与MPS中相关的张量进行收缩，部分近邻的两个张量被收缩成为一个四阶张量，如下图
   

4. 利用规范变换，将 MPS 变换为中心正交形式，正交中心为第一个四阶张量。

5. 对中心张量进行奇异值分解，如果奇异谱的维数大于预设的截断维数 $\chi$，则仅保留 $\chi$ 个最大的奇异值与相应的奇异向量
   

6. 移动正则中心到下一个四阶张量，重复第五步，直到所有四阶张量被分解为三阶张量

7. 进行完上述步骤后，MPS 被变换称为演化前的形式，即都是三阶张量的形式，重复第三步到第六步，收缩下一个半层张量直到MPS收敛。

TEBD 算法的具体实施方案并不唯一，我们可以每次只演化一个局域门，将一个局域门演化为三阶张量后再去演化下一个，如下图所示：

2. DMRG 密度矩阵重整化群算法

---

       TEBD 算法计算基态的思路是退火，密度矩阵重整化群算法（DMRG）采用基于最大本征态求解对应的最优化问题，即求解下面的极小化问题：
$$E_g=\underset{⟨\phi \mid \phi ⟩=1}{\min }⟨\phi |\hat{H}|\phi ⟩$$

当哈密顿量为 $\hat{H}={\sum }_{⟨i,j⟩}{\hat{H}}_{ij}$ ，有：
$$E_g=\underset{⟨\phi \mid \phi ⟩=1}{\min }\sum _{⟨i,j⟩}⟨\phi \left|{\hat{H}}_{ij}\right|\phi ⟩$$

他们每一项对应的张量网络图如下所示，以 $⟨\phi \left|{\hat{H}}_{23}\right|\phi ⟩$ 为例可以用图表示为：
DMRG 的策略是更新各个张量，使能量达到最小，具体的更新策略也不唯一。

       在 单点 DMRG 中，每次更新 MPS 中的一个张量，其余张量看成是给定的参数，一起来更新一个张量，该方法称为交替最小二乘算法（ALS）。

下面我们将单个张量的优化问题等价为局域矩阵最大本征值的问题，以更新第二个张量为例，定义有效哈密顿量 $\widehat{h_2}$ 为下面各项的求和：

其中每一项收缩之后的结果为六阶张量，$\widehat{h_2}$ 为多个六阶张量的求和，所以也是六阶张量，我们表示为

上面的过程就是我们先计算跟第二个张量无关的收缩，得到张量，下面我们将第二个张量和得到的六阶张量收缩，得到系统的能量变为：
$$E=⟨A^{(2)}\left|{\hat{h}}_2\right|A^{(2)}⟩$$

可以用图形表示为：

最优化问题变为了：
$$\underset{⟨\phi \mid \phi ⟩=1}{\min }⟨A^{(2)}\left|{\hat{h}}_2\right|A^{(2)}⟩$$

此时我们可以通过规范变换，将正则中心移动至 $A^{(2)}$ ，有：
$$⟨\phi \mid \phi ⟩=1\iff \left|A^{(2)}\right|=1$$

所以最优化问题最终变成了：
$$\underset{|A^{(2)}|=1}{\min }⟨A^{(2)}\left|{\hat{h}}_2\right|A^{(2)}⟩$$

所以 $A^{(2)}$ 就是算符 $\widehat{h_2}$ 的最低本征态，我们应该将 $A^{(2)}$ 更新为 $\widehat{h_2}$ 的最低本征态，求解 $\widehat{h_2}$ 的最低本征态的计算复杂度是可控的，也就是求解维数为 $d{\chi }^2\times d{\chi }^2$ 的厄密矩阵的本征态。我们直接将 $\widehat{h_2}$ 转换为矩阵，进行本征值分解，得到最小的本征态替换为 $A^{(2)}$ 。

我们可以用上面的方式，不断地重复更新 $A^{(1)}、A^{(2)}、A^{(3)}、A^{(4)}、A^{(5)}、A^{(6)}、A^{(7)}$ ，直到 MPS 态收敛。

得到基态MPS后，我们可利用中心正交形式简化对算符观测量的计算，当正交中心与观测量所在格点重合时，我们可以将下面左图的计算简化为右图：

当正交中心不与观测量所在格点重合时，在不移动正交中心的情况下，计算可做例如下图的简化：

综上所述，DMRG 算法，对第 n 个张量 $A^{(n)}$ 的更新步骤为：

1. 将MPS正交中心移动到第 n 个张量
2. 计算第 n 个张量对应的有效哈密顿量 $\widehat{h_n}$
3. 将 $A^{(n)}$ 更新为 $\widehat{h_n}$ 的最低本征态

整个 DMRG 算法的步骤为：

1. 随机初始化 MPS 态中的各个张量 $\left\{A^{(n)}\right\}$
2. 按上面的步骤一次更新第 1 到 n 个张量，再依次更新第 n 到第 1 个张量
3. 如果 MPS 收敛，计算完成，否则重复进行第二步

DMRG算法的关键在于将整个MPS的优化问题化成多个有效哈密顿量的最小本征问题，循环优化直至MPS收敛。

3. 基于自动微分的基态变分算法

---

       自动微分技术用来在机器学习中计算模型变分参数关于损失函数的梯度，我们称之为反向传播算法。

       这里我们可以用自动微分算法求解实对称矩阵的最大本征向量，即求解如下的极大值问题
$$\underset{|v|=1}{\max }\left|v^{\mathrm{T}}Mv\right|$$

       我们可以在上面的式子前添加一个符号将求解最大值问题转换为求解极小值问题，我们定义损失函数
$$f=-\frac{v^{\mathrm{T}}Mv}{|v{|}^2}$$

其中的 $\frac{v}{|v|}$ 是对向量的归一化，因此该式满足上面的 $|v|=1$ 的条件，该损失函数的最小本征向量就是 $\frac{v}{|v|}$ 。

       为了计算 v ，我们可以计算 $v$ 关于 $f$ 的梯度 $\frac{df}{dv}$ ，要使得 $f$ 减小，需要将 $v$ 沿负梯度方向进行更新，即每次将 $v$ 做如下改变：
$$v\Rightarrow v-\eta \frac{df}{dv}$$

其中，$\eta$ 为认为给定的常数，被称为更新步长或学习率。

       进行多次迭代更新之后，得到收敛的 $v$ ，那么最大本征向量 $\tilde{v}$ 与最大本征值 $\lambda$ 满足：
$$\tilde{v}=\frac{v}{|v|},\lambda ={\tilde{v}}^{\mathrm{T}}M\tilde{v}$$

我们可以用程序实现上述算法，上述算法中的梯度 $\frac{df}{dv}$ 可以直接使用自动微分技术计算，可以使用优化器让程序自适应地控制学习率，以达到较好的稳定性和收敛速率。

       根据最大本征值问题与基态问题地等价性，显然，自动微分方法可用于求解多体系统地基态，对应的优化问题可以写为：
E = min ⁡ { A ( n ) } ⟨ φ ∣ H ^ ∣ φ ⟩ ⟨ φ ∣ φ ⟩ E=\min _{\{A^{(n)}\}} \frac{\langle\varphi|\widehat{H}| \varphi\rangle}{\langle\varphi \mid \varphi\rangle} E={A(n)}min​⟨φ∣φ⟩⟨φ∣H ∣φ⟩​

其中，变分参数为 MPS 中地各个张量 { A ( n ) } \{A^{(n)}\} {A(n)} ，梯度更新地公式为：
$$A^{(n)}\Rightarrow A^{(n)}-\eta \frac{\partial E}{\partial A^{(n)}}$$

这里我们可以更加直接地将 MPS 看作是对基态地一种特殊的参数化形式，能量就是变分的损失函数。之后我们就可以使用 BP 算法和各种优化器进行梯度的更新了。