From e3f0d63765096ada1161bb5185fcf8a005151bb8 Mon Sep 17 00:00:00 2001 From: Yutaro Iiyama Date: Wed, 26 Apr 2023 01:26:39 +0900 Subject: [PATCH] reordered operators + using SparsePauliOp --- source/ja/dynamics_simulation.md | 32 ++++----- source/ja/more_dynamics.md | 38 +++++------ source/ja/spectrum_estimation.md | 26 ++++---- source/qc_workbook/dynamics.py | 66 +++++++++++++++++-- source/qc_workbook/hamiltonian.py | 105 ------------------------------ 5 files changed, 111 insertions(+), 156 deletions(-) delete mode 100644 source/qc_workbook/hamiltonian.py diff --git a/source/ja/dynamics_simulation.md b/source/ja/dynamics_simulation.md index 492dc167..3260cd76 100644 --- a/source/ja/dynamics_simulation.md +++ b/source/ja/dynamics_simulation.md @@ -315,7 +315,7 @@ $$ ```{math} :label: heisenberg -H = -J \sum_{j=0}^{n-2} (\sigma^X_j\sigma^X_{j+1} + \sigma^Y_j\sigma^Y_{j+1} + \sigma^Z_j \sigma^Z_{j+1}) +H = -J \sum_{j=0}^{n-2} (\sigma^X_{j+1}\sigma^X_{j} + \sigma^Y_{j+1}\sigma^Y_{j} + \sigma^Z_{j+1} \sigma^Z_{j}) ``` です。ここで、$\sigma^{[X,Y,Z]}_j$は第$j$スピンに作用するパウリ演算子です。 @@ -323,7 +323,7 @@ H = -J \sum_{j=0}^{n-2} (\sigma^X_j\sigma^X_{j+1} + \sigma^Y_j\sigma^Y_{j+1} + \ ただし、式{eq}`heisenberg`の和の記法には実は若干の省略があります。例えば第$j$項をより正確に書くと、 $$ -I_0 \otimes \dots \otimes I_{j-1} \otimes \sigma^X_j \otimes \sigma^X_{j+1} \otimes I_{j+2} \otimes \dots I_{n-1} +I_{n-1} \otimes \dots \otimes I_{j+2} \otimes \sigma^X_{j+1} \otimes \sigma^X_{j} \otimes I_{j-1} \otimes \dots I_{0} $$ です。ここで$\otimes$は線形演算子間の「テンソル積」を表しますが、聞き慣れない方は掛け算だと思っておいてください。重要なのは、式{eq}`heisenberg`の各項が、上で触れたように$n$個の基底演算子の積になっているということです。さらに、この系では隣接スピン間の相互作用しか存在しないため、ハミルトニアンが$n-1$個の項に分解できています。 @@ -335,27 +335,27 @@ $$ 時間発展演算子は $$ -U_H(t) = \exp \left[ \frac{i\omega t}{2} \sum_{j=0}^{n-2} (\sigma^X_j\sigma^X_{j+1} + \sigma^Y_j\sigma^Y_{j+1} + \sigma^Z_j \sigma^Z_{j+1}) \right] +U_H(t) = \exp \left[ \frac{i\omega t}{2} \sum_{j=0}^{n-2} (\sigma^X_{j+1}\sigma^X_{j} + \sigma^Y_{j+1}\sigma^Y_{j} + \sigma^Z_{j+1} \sigma^Z_{j}) \right] $$ ですが、ハミルトニアンの各項が互いに可換でないので、シミュレーションでは鈴木・トロッター分解を用いて近似します。各時間ステップ$\Delta t$での近似時間発展は $$ -\tilde{U}_{H;\Delta t} = \prod_{j=0}^{n-2} \exp\left( \frac{i \omega \Delta t}{2} \sigma^X_j\sigma^X_{j+1} \right) \exp\left( \frac{i \omega \Delta t}{2} \sigma^Y_j\sigma^Y_{j+1} \right) \exp\left( \frac{i \omega \Delta t}{2} \sigma^Z_j\sigma^Z_{j+1} \right) +\tilde{U}_{H;\Delta t} = \prod_{j=0}^{n-2} \exp\left( \frac{i \omega \Delta t}{2} \sigma^X_{j+1}\sigma^X_{j} \right) \exp\left( \frac{i \omega \Delta t}{2} \sigma^Y_{j+1}\sigma^Y_{j} \right) \exp\left( \frac{i \omega \Delta t}{2} \sigma^Z_{j+1}\sigma^Z_{j} \right) $$ です。 ### 量子ゲートでの表現 -これを回転ゲートと制御ゲートで表します。まず$\exp(\frac{i \omega \Delta t}{2} \sigma^Z_j\sigma^Z_{j+1})$について考えてみましょう。この演算子の$j$-$(j+1)$スピン系の4つの基底状態への作用は +これを回転ゲートと制御ゲートで表します。まず$\exp(\frac{i \omega \Delta t}{2} \sigma^Z_{j+1}\sigma^Z_{j})$について考えてみましょう。この演算子の$j$-$(j+1)$スピン系の4つの基底状態への作用は $$ \begin{align} -\upket_j \upket_{j+1} \rightarrow e^{i \omega \Delta t / 2} \upket_j \upket_{j+1} \\ -\upket_j \downket_{j+1} \rightarrow e^{-i \omega \Delta t / 2} \upket_j \downket_{j+1} \\ -\downket_j \upket_{j+1} \rightarrow e^{-i \omega \Delta t / 2} \downket_j \upket_{j+1} \\ -\downket_j \downket_{j+1} \rightarrow e^{i \omega \Delta t / 2} \downket_j \downket_{j+1} +\upket_{j+1} \upket_{j} \rightarrow e^{i \omega \Delta t / 2} \upket_{j+1} \upket_{j} \\ +\upket_{j+1} \downket_{j} \rightarrow e^{-i \omega \Delta t / 2} \upket_{j+1} \downket_{j} \\ +\downket_{j+1} \upket_{j} \rightarrow e^{-i \omega \Delta t / 2} \downket_{j+1} \upket_{j} \\ +\downket_{j+1} \downket_{j} \rightarrow e^{i \omega \Delta t / 2} \downket_{j+1} \downket_{j} \end{align} $$ @@ -386,16 +386,16 @@ $$ \end{align} $$ -と変換するので（確認してください）、まさに$\exp(\frac{i \omega \Delta t}{2} \sigma^Z_j\sigma^Z_{j+1})$の表現になっています。 +と変換するので（確認してください）、まさに$\exp(\frac{i \omega \Delta t}{2} \sigma^Z_{j+1}\sigma^Z_{j})$の表現になっています。 -残りの2つの演算子も同様にパリティに対する回転で表せますが、CNOTで表現できるのは$Z$方向のパリティだけなので、先にスピンを回転させる必要があります。$\exp(\frac{i \omega \Delta t}{2} \sigma^X_j\sigma^X_{j+1})$による変換は +残りの2つの演算子も同様にパリティに対する回転で表せますが、CNOTで表現できるのは$Z$方向のパリティだけなので、先にスピンを回転させる必要があります。$\exp(\frac{i \omega \Delta t}{2} \sigma^X_{j+1}\sigma^X_{j})$による変換は $$ \begin{align} -\rightket_j \rightket_{j+1} \rightarrow e^{i \omega \Delta t / 2} \rightket_j \rightket_{j+1} \\ -\rightket_j \leftket_{j+1} \rightarrow e^{-i \omega \Delta t / 2} \rightket_j \leftket_{j+1} \\ -\leftket_j \rightket_{j+1} \rightarrow e^{-i \omega \Delta t / 2} \leftket_j \rightket_{j+1} \\ -\leftket_j \leftket_{j+1} \rightarrow e^{i \omega \Delta t / 2} \leftket_j \leftket_{j+1} +\rightket_{j+1} \rightket_{j} \rightarrow e^{i \omega \Delta t / 2} \rightket_{j+1} \rightket_{j} \\ +\rightket_{j+1} \leftket_{j} \rightarrow e^{-i \omega \Delta t / 2} \rightket_{j+1} \leftket_{j} \\ +\leftket_{j+1} \rightket_{j} \rightarrow e^{-i \omega \Delta t / 2} \leftket_{j+1} \rightket_{j} \\ +\leftket_{j+1} \leftket_{j} \rightarrow e^{i \omega \Delta t / 2} \leftket_{j+1} \leftket_{j} \end{align} $$ @@ -415,7 +415,7 @@ circuit.h(1) circuit.draw('mpl') ``` -最後に、$\exp(\frac{i \omega \Delta t}{2} \sigma^Y_j\sigma^Y_{j+1})$に対応する回路は +最後に、$\exp(\frac{i \omega \Delta t}{2} \sigma^Y_{j+1}\sigma^Y_{j})$に対応する回路は ```{code-cell} ipython3 :tags: [remove-input] diff --git a/source/ja/more_dynamics.md b/source/ja/more_dynamics.md index 7e15ed13..f9908e3b 100644 --- a/source/ja/more_dynamics.md +++ b/source/ja/more_dynamics.md @@ -56,10 +56,10 @@ $\newcommand{\minusket}{\ket{-}}$ import numpy as np import matplotlib.pyplot as plt from qiskit import QuantumCircuit, transpile +from qiskit.quantum_info import SparsePauliOp from qiskit_aer import AerSimulator # このワークブック独自のモジュール -from qc_workbook.dynamics import plot_heisenberg_spins, bit_expectations_sv, bit_expectations_counts -from qc_workbook.hamiltonian import make_hamiltonian, diagonalized_evolution +from qc_workbook.dynamics import plot_heisenberg_spins, bit_expectations_sv, bit_expectations_counts, diagonalized_evolution ``` ```{code-cell} ipython3 @@ -266,7 +266,7 @@ $$ と表現できることがわかります。一方、$\Phi^{\dagger}_j\Phi_{j+1}$に関しては、やや込み入った議論{cite}`PhysRevD.13.1043`の末、 $$ -\Phi^{\dagger}_j\Phi_{j+1} \rightarrow i \sigma^+_j \sigma^-_{j+1} +\Phi^{\dagger}_j\Phi_{j+1} \rightarrow i \sigma^-_{j+1} \sigma^+_{j} $$ が正しい表現であることがわかっています。ここで、 @@ -275,16 +275,16 @@ $$ \sigma^{\pm} = \frac{1}{2}(\sigma^X \pm i \sigma^Y) $$ -です。ハミルトニアンには$\Phi_j\Phi^{\dagger}_{j+1} \rightarrow i \sigma^-_j \sigma^+_{j+1}$も登場するので、二つの項を合わせると +です。また、このワークブックでの約束に従って、右辺の$j$と$j+1$の順番をひっくり返してあります。ハミルトニアンには$\Phi_j\Phi^{\dagger}_{j+1} \rightarrow i \sigma^+_{j+1}\sigma^-_{j}$も登場するので、二つの項を合わせると $$ -\Phi^{\dagger}_{j} \Phi_{j+1} + \Phi_j \Phi^{\dagger}_{j+1} \rightarrow \frac{i}{2} (\sigma^X_j \sigma^X_{j+1} + \sigma^Y_j \sigma^Y_{j+1}) +\Phi^{\dagger}_{j} \Phi_{j+1} + \Phi_j \Phi^{\dagger}_{j+1} \rightarrow \frac{i}{2} (\sigma^X_{j+1}\sigma^X_j + \sigma^Y_{j+1}\sigma^Y_j) $$ となります。まとめると、 $$ -H \rightarrow \frac{1}{4a} \left\{ \sum_{j=0}^{n-2} (\sigma^X_j \sigma^X_{j+1} + \sigma^Y_j \sigma^Y_{j+1}) + J \sum_{j=1}^{n-2} (n - j - 1) \sum_{k=0}^{j-1} \sigma^Z_k \sigma^Z_j + \sum_{j=0}^{n-1} \left[ (-1)^{j+1} \mu - J \flhalf{n-j} \right] \sigma^Z_j \right\} +H \rightarrow \frac{1}{4a} \left\{ \sum_{j=0}^{n-2} (\sigma^X_{j+1}\sigma^X_j + \sigma^Y_{j+1}\sigma^Y_j) + J \sum_{j=1}^{n-2} (n - j - 1) \sum_{k=0}^{j-1} \sigma^Z_j \sigma^Z_k + \sum_{j=0}^{n-1} \left[ (-1)^{j+1} \mu - J \flhalf{n-j} \right] \sigma^Z_j \right\} $$ です。ただし、計算過程で現れる定数項（恒等演算子に比例する項）は時間発展において系の状態に全体位相をかける作用しか持たないため、無視しました。 @@ -309,18 +309,18 @@ $$ **ヒント**: -上のハミルトニアンのパラメターの値は参考文献{cite}`Martinez_2016`と同一です。したがって、$n=4$, $\omega \Delta t = \pi/8$とすれば、論文中の図3aを再現できるはずです。答え合わせに使ってください。 +上のハミルトニアンのパラメターの値は参考文献{cite}`Martinez_2016`と同一です（ただし$\sigma$積の順番は逆です）。したがって、$n=4$, $\omega \Delta t = \pi/8$とすれば、論文中の図3aを再現できるはずです。答え合わせに使ってください。 -また、問題を解くためのヒントではありませんが、ハイゼンベルグモデルと同様にこのモデルでも対角化による厳密解を比較的簡単にプロットできるように道具立てがしてあります。下のコードのテンプレートでは、シミュレーション回路と厳密解を計算するためのハミルトニアンのパウリ行列分解だけ指定すれば、`plot_heisenberg_spins`と同様のプロットが作成されるようになっています。パウリ行列分解を指定するには、`paulis`と`coeffs`という二つのリストを作ります。これらのリストの長さはハミルトニアンの項数で、`paulis`の各要素は対応する項のパウリ行列のリスト、`coeffs`の各要素はその項にかかる係数にします。例えば +また、問題を解くためのヒントではありませんが、ハイゼンベルグモデルと同様にこのモデルでも対角化による厳密解を比較的簡単にプロットできるように道具立てがしてあります。下のコードのテンプレートでは、シミュレーション回路と厳密解を計算するためのハミルトニアンのパウリ行列分解だけ指定すれば、`plot_heisenberg_spins`と同様のプロットが作成されるようになっています。パウリ行列分解を指定するには、`paulis`と`coeffs`という二つのリストを作り、Qiskitの`SparsePauliOp`というクラスに渡します。これらのリストの長さはハミルトニアンの項数で、`paulis`の各要素は対応する項のパウリ行列のリスト、`coeffs`の各要素はその項にかかる係数にします。例えば $$ -H = 0.5 \sigma^X_0 \sigma^Y_1 I_2 + I_0 \sigma^Z_1 \sigma^X_2 +H = 0.5 \sigma^X_2 \sigma^Y_1 I_0 + I_2 \sigma^Z_1 \sigma^X_0 $$ というハミルトニアンに対しては、 ```{code-block} python -paulis = [['x', 'y', 'i'], ['i', 'z', 'x']] +paulis = ['XYI', 'IZX'] coeffs = [0.5, 1.] ``` @@ -378,25 +378,25 @@ simulator = AerSimulator() circuits = transpile(circuits, backend=simulator) sim_job = simulator.run(circuits, shots=shots) sim_counts_list = sim_job.result().get_counts() +``` -## Numerical solution through diagonalization +```{code-cell} ipython3 +:tags: [raises-exception, remove-output] -# Construct the Hamiltonian -paulis = [] -coeffs = [] +## Numerical solution through diagonalization ################## ### EDIT BELOW ### ################## -paulis = [['i'] * n] -coeffs = None +paulis = ['I' * n] +coeffs = [1.] ################## ### EDIT ABOVE ### ################## -hamiltonian = make_hamiltonian(paulis, coeffs) +hamiltonian = SparsePauliOp(paulis, coeffs).to_matrix() # Initial state as a statevector initial_state = np.zeros(2 ** n, dtype=np.complex128) @@ -417,11 +417,11 @@ plt.plot(time_points, number_density(bit_exp)) initial_probs = np.square(np.abs(initial_state)) fmt = f'{{:0{n}b}}' initial_counts = dict((fmt.format(idx), prob) for idx, prob in enumerate(initial_probs) if prob != 0.) -sim_counts_list = [initial_counts] + sim_counts_list +sim_counts_list_with_init = [initial_counts] + sim_counts_list # Plot the simulation results time_points = np.linspace(0., omegadt * M, M + 1, endpoint=True) -_, bit_exp = bit_expectations_counts(time_points, sim_counts_list, n) +_, bit_exp = bit_expectations_counts(time_points, sim_counts_list_with_init, n) plt.plot(time_points, number_density(bit_exp), 'o') ``` diff --git a/source/ja/spectrum_estimation.md b/source/ja/spectrum_estimation.md index 531d1323..9fde24e6 100644 --- a/source/ja/spectrum_estimation.md +++ b/source/ja/spectrum_estimation.md @@ -49,7 +49,7 @@ $\newcommand{\ket}[1]{|#1\rangle}$ 前回登場したハイゼンベルグモデルのハミルトニアンは $$ -H = -J \sum_{j=0}^{n-2} (\sigma^X_j\sigma^X_{j+1} + \sigma^Y_j\sigma^Y_{j+1} + \sigma^Z_j \sigma^Z_{j+1}) \quad (J > 0) +H = -J \sum_{j=0}^{n-2} (\sigma^X_{j+1}\sigma^X_{j} + \sigma^Y_{j+1}\sigma^Y_{j} + \sigma^Z_{j+1} \sigma^Z_{j}) \quad (J > 0) $$ というものでした。このハミルトニアンが表しているのは、空間中で一列に並んだスピンを持つ粒子が、隣接粒子間で相互作用を及ぼしているような系でした。ここで相互作用は、スピンの向きが揃っているときにエネルギーが低くなるようなものでした。したがって、全てののスピンが同じ方向を向いているときにエネルギーが最も低くなることが予想されました。 @@ -57,7 +57,7 @@ $$ 今回は、このハミルトニアンに外部からの磁場の影響を入れます。外部磁場がある時は、スピンが磁場の方向を向いているときにエネルギーが低くなります。したがって、外部磁場を$+Z$方向にかけるとすれば、ハミルトニアンは $$ -H = -J \sum_{j=0}^{n-1} (\sigma^X_j\sigma^X_{j+1} + \sigma^Y_j\sigma^Y_{j+1} + \sigma^Z_j \sigma^Z_{j+1} + g \sigma^Z_j) +H = -J \sum_{j=0}^{n-1} (\sigma^X_{j+1}\sigma^X_{j} + \sigma^Y_{j+1}\sigma^Y_{j} + \sigma^Z_{j+1} \sigma^Z_{j} + g \sigma^Z_j) $$ となります。このハミルトニアンにはもう一点前回と異なる部分があります。前回はスピンに関する和を$j=0$から$n-2$まで取ることで、両端のスピンは「内側」のスピンとしか相互作用をしないような境界条件を採用していました。今回は和を$n-1$まで取っています。$\sigma^{X,Y,Z}_n$を$\sigma^{X,Y,Z}_0$と同一視することで、これは「周期境界条件」（一列ではなく環状に並んだスピン）を表します。 @@ -71,10 +71,10 @@ $$ import numpy as np import matplotlib.pyplot as plt from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, transpile +from qiskit.quantum_info import SparsePauliOp from qiskit.visualization import plot_histogram from qiskit_aer import AerSimulator # ワークブック独自のモジュール -from qc_workbook.hamiltonian import make_hamiltonian from qc_workbook.show_state import show_state print('notebook ready') @@ -91,18 +91,22 @@ g = 0. # Construct the Hamiltonian matrix paulis = list() coeffs = list() + +xx_template = 'I' * (n_s - 2) + 'XX' +yy_template = 'I' * (n_s - 2) + 'YY' +zz_template = 'I' * (n_s - 2) + 'ZZ' + for j in range(n_s): - paulis.append(list('x' if k in (j, (j + 1) % n_s) else 'i' for k in range(n_s))) - coeffs.append(-J) - paulis.append(list('y' if k in (j, (j + 1) % n_s) else 'i' for k in range(n_s))) - coeffs.append(-J) - paulis.append(list('z' if k in (j, (j + 1) % n_s) else 'i' for k in range(n_s))) - coeffs.append(-J) + paulis.append(xx_template[j:] + xx_template[:j]) + paulis.append(yy_template[j:] + yy_template[:j]) + paulis.append(zz_template[j:] + zz_template[:j]) + coeffs += [-J] * 3 + if g != 0.: - paulis.append(list('z' if k == j else 'i' for k in range(n_s))) + paulis.append('I' * (n_s - j - 1) + 'Z' + 'I' * j) coeffs.append(-J * g) -hamiltonian = make_hamiltonian(paulis, coeffs) +hamiltonian = SparsePauliOp(paulis, coeffs).to_matrix() # Diagonalize and obtain the eigenvalues and vectors eigvals, eigvectors = np.linalg.eigh(hamiltonian) diff --git a/source/qc_workbook/dynamics.py b/source/qc_workbook/dynamics.py index 27737da4..cf4ba243 100644 --- a/source/qc_workbook/dynamics.py +++ b/source/qc_workbook/dynamics.py @@ -2,8 +2,8 @@ import matplotlib as mpl import matplotlib.pyplot as plt from qiskit import QuantumCircuit +from qiskit.quantum_info import SparsePauliOp -from .hamiltonian import tensor_product, make_hamiltonian, diagonalized_evolution def make_heisenberg_circuits(n_spins, M, omegadt): circuits = [] @@ -71,6 +71,45 @@ def bit_expectations_sv(time_points, statevectors): return x, y +def diagonalized_evolution(hamiltonian, initial_state, time, num_steps=100): + """Diagonalize the given reduced Hamiltonian and evolve the initial state by exp(-i time*hamiltonian). + + Args: + hamiltonian (np.ndarray(shape=(D, D), dtype=np.complex128)): Hamiltonian matrix divided by hbar. + initial_state (np.ndarray(shape=(D,), dtype=np.complex128)): Initial state vector. + time (float): Evolution time. + num_steps (int): Number of steps (T) to divide time into. + + Returns: + np.ndarray(shape=(T,), dtype=float): Time points. + np.ndarray(shape=(D, T), dtype=np.complex128): State vector as a function of time. + """ + + num_dim = hamiltonian.shape[0] + num_qubits = np.round(np.log2(num_dim)).astype(int) + + # Create the array of time points + time_points = np.linspace(0., time, num_steps, endpoint=True) + + ## Diagonalize the Hamiltonian + eigvals, eigvectors = np.linalg.eigh(hamiltonian) + + ## Decompose the initial state vector into a linear combination of eigenvectors + # Matrix eigvectors has the form [v_0 v_1 v_2 ..], where v_i . v_j = delta_ij + # -> eigvectors^dagger @ initial_state = coefficients for the eigenvector decomposition of the initial state vector + initial_coeff = eigvectors.T.conjugate() @ initial_state + # Initial state as a matrix [c_0 v_0, c_1 v_1, ...] (shape (D, D)) + initial_state_matrix = eigvectors * initial_coeff + + ## Time-evolve the initial state to each time point + # Phase at each time point (shape (D, T)) + phase = np.outer(-1.j * eigvals, time_points) + phase_factor = np.exp(phase) + statevectors = initial_state_matrix @ phase_factor # shape (D, T) + + return time_points, statevectors + + def bit_expectations_counts(time_points, counts_list, num_bits): """Compute the bit expectation values from experiment results. @@ -134,11 +173,11 @@ def plot_heisenberg_spins(counts_list, num_spins, initial_state, omegadt, add_th # Construct the numerical Hamiltonian matrix from a list of Pauli operators paulis = list() for j in range(num_spins - 1): - paulis.append(list('x' if k in (j, j + 1) else 'i' for k in range(num_spins))) - paulis.append(list('y' if k in (j, j + 1) else 'i' for k in range(num_spins))) - paulis.append(list('z' if k in (j, j + 1) else 'i' for k in range(num_spins))) + paulis.append('I' * (num_spins - j - 2) + 'XX' + 'I' * j) + paulis.append('I' * (num_spins - j - 2) + 'YY' + 'I' * j) + paulis.append('I' * (num_spins - j - 2) + 'ZZ' + 'I' * j) - hamiltonian = make_hamiltonian(paulis) + hamiltonian = SparsePauliOp(paulis).to_matrix() # Compute the statevector as a function of time from Hamiltonian diagonalization time_points, statevectors = diagonalized_evolution(-0.5 * hamiltonian, initial_state, omegadt * num_steps) @@ -198,3 +237,20 @@ def plot_heisenberg_spins(counts_list, num_spins, initial_state, omegadt, add_th ax.set_xlabel(r'$\omega t$') ax.set_ylabel(r'$\langle S_z \rangle$') + + +def tensor_product(ops): + """Recursively apply np.kron to construct a tensor product of the operators. + + Args: + ops (List): List of (2, 2) arrays. + + Returns: + np.ndarray(shape=(2 ** nops, 2 ** nops), dtype=np.complex128): Tensor product of ops. + """ + + prod = 1. + for op in ops: + prod = np.kron(op, prod) + + return prod diff --git a/source/qc_workbook/hamiltonian.py b/source/qc_workbook/hamiltonian.py deleted file mode 100644 index 57ce34e2..00000000 --- a/source/qc_workbook/hamiltonian.py +++ /dev/null @@ -1,105 +0,0 @@ -import sys -import numpy as np - -def tensor_product(ops): - """Recursively apply np.kron to construct a tensor product of the operators. - - Args: - ops (List): List of (2, 2) arrays. - - Returns: - np.ndarray(shape=(2 ** nops, 2 ** nops), dtype=np.complex128): Tensor product of ops. - """ - - prod = 1. - for op in ops: - prod = np.kron(op, prod) - - return prod - - -def make_hamiltonian(paulis, coeffs=None): - """Compute the Hamiltonian matrix from its Pauli decomposition. - - Args: - paulis (List(List(str))): Terms in Pauli decomposition. Form [['i', 'x', 'z', ..], ['x', 'y', 'i', ..], ..] - All inner lists must be of the same length (number of qubits n). - coeffs (None or List(float)): If not None, the list must be of the same length as the outer list of paulis - and should specify the coefficient of each term. - - Returns: - np.ndarray(shape=(2 ** n, 2 ** n), dtype=np.complex128): The numerical Hamiltonian matrix. The first qubit - corresponds to the least significant digit. - """ - - if len(paulis) == 0: - return np.array([[0.]], dtype=np.complex128) - - qubit_nums = set(len(term) for term in paulis) - assert len(qubit_nums) == 1, 'List of paulis must all have the same length.' - - if coeffs is None: - coeffs = [1.] * len(paulis) - - num_qubits = qubit_nums.pop() - - # Basis matrices - basis_matrices = { - 'i': np.array([[1., 0.], [0., 1.]], dtype=np.complex128), - 'x': np.array([[0., 1.], [1., 0.]], dtype=np.complex128), - 'y': np.array([[0., -1.j], [1.j, 0.]], dtype=np.complex128), - 'z': np.array([[1., 0.], [0., -1.]], dtype=np.complex128) - } - - # Start with an empty matrix - hamiltonian = np.zeros((2 ** num_qubits, 2 ** num_qubits), dtype=np.complex128) - - for iterm, term in enumerate(paulis): - try: - ops = list(basis_matrices[op.lower()] for op in term) - except KeyError as err: - sys.stderr.write('Invalid operator {} in term {}\n'.format(err.args[0], iterm)) - raise - - hamiltonian += coeffs[iterm] * tensor_product(ops) - - return hamiltonian - - -def diagonalized_evolution(hamiltonian, initial_state, time, num_steps=100): - """Diagonalize the given reduced Hamiltonian and evolve the initial state by exp(-i time*hamiltonian). - - Args: - hamiltonian (np.ndarray(shape=(D, D), dtype=np.complex128)): Hamiltonian matrix divided by hbar. - initial_state (np.ndarray(shape=(D,), dtype=np.complex128)): Initial state vector. - time (float): Evolution time. - num_steps (int): Number of steps (T) to divide time into. - - Returns: - np.ndarray(shape=(T,), dtype=float): Time points. - np.ndarray(shape=(D, T), dtype=np.complex128): State vector as a function of time. - """ - - num_dim = hamiltonian.shape[0] - num_qubits = np.round(np.log2(num_dim)).astype(int) - - # Create the array of time points - time_points = np.linspace(0., time, num_steps, endpoint=True) - - ## Diagonalize the Hamiltonian - eigvals, eigvectors = np.linalg.eigh(hamiltonian) - - ## Decompose the initial state vector into a linear combination of eigenvectors - # Matrix eigvectors has the form [v_0 v_1 v_2 ..], where v_i . v_j = delta_ij - # -> eigvectors^dagger @ initial_state = coefficients for the eigenvector decomposition of the initial state vector - initial_coeff = eigvectors.T.conjugate() @ initial_state - # Initial state as a matrix [c_0 v_0, c_1 v_1, ...] (shape (D, D)) - initial_state_matrix = eigvectors * initial_coeff - - ## Time-evolve the initial state to each time point - # Phase at each time point (shape (D, T)) - phase = np.outer(-1.j * eigvals, time_points) - phase_factor = np.exp(phase) - statevectors = initial_state_matrix @ phase_factor # shape (D, T) - - return time_points, statevectors