基于子空间计算给定哈密顿量的期望值

基于子空间计算给定哈密顿量的期望值

在运行基于量子线路模型的量子算法时，系统的哈密顿量往往由多项构成，而每一子项可能只涉及单个或者几个量子比特，如果每次计算都构建完整的矩阵，会造成资源的浪费，尤其是当系统规模很大时，成本会很高。此时，我们可以分别在每个子项所涉及的那几个量子比特的子空间内直接对量子态进行求期望值，最后再对各子项求和，就可以快速完成对哈密顿量期望值计算。

我们首先考虑一个三比特的系统，并假设我们的哈密顿量表达式为
$$H=0.2\ast Y_0+Z_2-1.5\ast X_0{Y}_1$$
其计算结果为
$$E=0.2\ast ⟨\psi |Y_0|\psi ⟩+⟨\psi |Z_2|\psi ⟩-1.5\ast ⟨\psi |X_0|\psi ⟩\ast ⟨\psi |Y_1|\psi ⟩$$
我们针对其中一项，比如 $0.2\ast ⟨\psi |Y_0|\psi ⟩$，其系数为 0.2，而 $⟨\psi |Y_0|\psi ⟩$ 的计算可以拆分成两步：

1. 先计算 $|{\psi }^{\prime }⟩=Y_0|\psi ⟩$
2. 再计算 $⟨\psi |{\psi }^{\prime }⟩$

其中，步骤 1 的计算过程其实就是将一个 $Y$ 门作用在 $q_0$ 上，获得一个变换后的态，这一点可以通过我之前一篇博客中的方法 基于子空间完成对量子比特的操作 快速完成。而步骤 2 则更为简单，只需要求原来量子态 $|\psi ⟩$ 的转置共轭，然后再与 $|{\psi }^{\prime }⟩$ 做内积即可。考虑到我们这里用一维数组来表示 $|\psi ⟩$，其实只需要对 $|\psi ⟩$ 做一个共轭，然后再与 $|{\psi }^{\prime }⟩$ 做元素乘法再求和即可。将得到的结果再与系数相乘。而其它项则可以类似处理。最后总的哈密顿量的期望值即是所有项的和。

我们先来展示如何计算第一项 $0.2\ast ⟨\psi |Y_0|\psi ⟩$ 的期望值。

import copy
import numpy as np
from mindquantum import *

n_qubits = 3
acted_qubit = 1

random_states = random_state((int(2**n_qubits),), seed=0)
print('随机初始量子态为：\n')
print(random_states)

hams = ('Y0', 0.2) # 我们将每一项用一个数组表示
acted_qubit = int(hams[0][1])

ham = Hamiltonian(QubitOperator(hams[0], hams[1]))
sim = Simulator('mqvector', n_qubits)
sim.set_qs(random_states)
mq_value = sim.get_expectation(ham)
print('\nMindQuantum 计算的期望值：\n',mq_value)

applied_matrix = eval(hams[0][0]).matrix()

blank_state = np.zeros((int(2**n_qubits),), dtype=complex)

for i in range(int(2**n_qubits)):
    index_str = bin(i)[2:].zfill(n_qubits)
    reversed_index_str = index_str[::-1]
    if reversed_index_str[acted_qubit] == '0':
        blank_state[i] += applied_matrix[0,0]*random_states[i]
        new_index_str = reversed_index_str[:acted_qubit]+'1'+reversed_index_str[acted_qubit+1:]
        blank_state[int(new_index_str[::-1], base=2)] += applied_matrix[1,0]*random_states[i] 
        
    if reversed_index_str[acted_qubit] == '1':
        blank_state[i] += applied_matrix[1,1]*random_states[i]
        new_index_str = reversed_index_str[:acted_qubit]+'0'+reversed_index_str[acted_qubit+1:]
        blank_state[int(new_index_str[::-1], base=2)] += applied_matrix[0,1]*random_states[i]


our_value = np.sum(random_states.conj()*blank_state) * hams[1]
print('\n我们方法最终得到的结果为：')
print(our_value)

随机初始量子态为：

[0.21986706+0.38606503j 0.28652098+0.15361532j 0.2414806 +0.3171829j
 0.21829249+0.21188723j 0.16972566+0.22757146j 0.25875974+0.37081488j
 0.1753073 +0.02845865j 0.35726437+0.03490596j]

MindQuantum 计算的期望值：
 (-0.037964117246960305-5.204170427930421e-18j)

我们方法最终得到的结果为：
(-0.0379641172469603+0j)

接下来，我们展示如何计算
$$E=0.2\ast ⟨\psi |Y_0|\psi ⟩+⟨\psi |Z_2|\psi ⟩-1.5\ast ⟨\psi |X_0|\psi ⟩\ast ⟨\psi |Y_1|\psi ⟩$$
其中，为便于阅读，我们将第三项用一个数组表示为 (‘X0 Y1’, -1.5)。而在计算过程中，会自动合并掉字符串中间空格，以便于索引。

term = ('X0 Y1', -1.5)
term_str = term[0].replace(" ","")
print('\n原字符串为：\n', term_str)
print('\n修改后的字符串为：\n', term_str)
for i in range(0, len(term_str), 2):
    print('\n泡利算符为：')
    print(term_str[i])
    print('\n所作用的比特为：')
    print(term_str[i+1])

原字符串为：
 X0Y1

修改后的字符串为：
 X0Y1

泡利算符为：
X

所作用的比特为：
0

泡利算符为：
Y

所作用的比特为：
1

import copy
import numpy as np
from mindquantum import *

n_qubits = 3
acted_qubit = 1

random_states = random_state((int(2**n_qubits),), seed=0)
print('随机初始量子态为：\n', random_states)

hams = [('Y0', 0.2), ('Z2', 1), ('X0 Y1', -1.5)]

operators = QubitOperator()
for qubit_operator in hams:
    operators += QubitOperator(qubit_operator[0], qubit_operator[1])
    
ham = Hamiltonian(operators)
sim = Simulator('mqvector', n_qubits)
sim.set_qs(random_states)
mq_value = sim.get_expectation(ham)
print('\nMindQuantum 计算的期望值:\n',mq_value)

our_value = 0
for term in hams:
    term_value = 1
    term_str = term[0].replace(" ", "") # 消去空格，好方便索引
    for operator_index in range(0,len(term_str), 2):
        acted_qubit = int(term_str[operator_index+1])
        applied_matrix = eval(term_str[operator_index]).matrix()

        blank_state = np.zeros((int(2**n_qubits),), dtype=complex)

        for i in range(int(2**n_qubits)):
            index_str = bin(i)[2:].zfill(n_qubits)
            reversed_index_str = index_str[::-1]
            if reversed_index_str[acted_qubit] == '0':
                blank_state[i] += applied_matrix[0,0]*random_states[i]
                new_index_str = reversed_index_str[:acted_qubit]+'1'+reversed_index_str[acted_qubit+1:]
                blank_state[int(new_index_str[::-1], base=2)] += applied_matrix[1,0]*random_states[i] 

            if reversed_index_str[acted_qubit] == '1':
                blank_state[i] += applied_matrix[1,1]*random_states[i]
                new_index_str = reversed_index_str[:acted_qubit]+'0'+reversed_index_str[acted_qubit+1:]
                blank_state[int(new_index_str[::-1], base=2)] += applied_matrix[0,1]*random_states[i]

        operator_value = np.sum(random_states.conj()*blank_state)
        term_value *= operator_value

    our_value += term_value * term[1]
    
print('\n我们方法最终得到的结果为：\n')
print(our_value)

随机初始量子态为：
 [0.21986706+0.38606503j 0.28652098+0.15361532j 0.2414806 +0.3171829j
 0.21829249+0.21188723j 0.16972566+0.22757146j 0.25875974+0.37081488j
 0.1753073 +0.02845865j 0.35726437+0.03490596j]

MindQuantum 计算的期望值:
 (0.42190150611033195+1.3877787807814457e-17j)

我们方法最终得到的结果为：

(0.4744281907889403+0j)

可见，我们方法得到的结果和 MindQuantum 的结果是一样的。