[DEMO] QML for Phase Detection

[SaverioMonaco]

Title:
Supervised and Unsipervised Quantum Machine Learning models for the phase detection of the ANNNI spin model

Summary:
Showcase of two foundational models in QML for Phase Detection for Spin Systems. One follows a supervised approach, the other follows an unsupervised approach.

Relevant references:

Possible Drawbacks:
QCNN models is considered now outdated.

Related GitHub Issues:

---

If you are writing a demonstration, please answer these questions to facilitate the marketing process.

• GOALS — Why are we working on this now?

  Eg. Promote a new PL feature or show a PL implementation of a recent paper.
  Show a working example of the founding models in QML for Phase Detection. Show how MPS and QML can be used together.

• AUDIENCE — Who is this for?

  Eg. Chemistry researchers, PL educators, beginners in quantum computing.
  Beginners in Quantum Computing, Chemistry researchers

• KEYWORDS — What words should be included in the marketing post?
  ANNNI, QCNN, QAD, Phase Detection

• Which of the following types of documentation is most similar to your file?
  (more details here)

• [ x] Tutorial
• Demo
• How-to

[sc-80710]

[SaverioMonaco]

• The image annni_sketch.png is a rough sketch for the thumbnail. It needs to be updated by replacing $J2$ with $\kappa$ and $J$ with $J_1$ to ensure consistency with the Hamiltonian described in Eq. (1).

• The other two images, annni_pd.png and annni_pd_analytical.png, were originally generated using matplotlib.pyplot. I imported them as images to reduce code volume. I will add the code to reproduce these images shortly.

[SaverioMonaco]

• annni_pd.png

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm

# Kosterlitz-Thouless transition line
def kt_transition(k):
    return 1.05 * np.sqrt((k - 0.5) * (k - 0.1))

# Ising transition line
def ising_transition(k):
    return np.where(k == 0, 1, (1 - k) * (1 - np.sqrt((1 - 3 * k + 4 * k**2) / (1 - k))) / np.maximum(k, 1e-9))

# Floating Phase transition line
def bkt_transition(k):
    return 1.05 * (k - 0.5)
    
# Get the phase from the DMRG transition lines
def get_phase(k, h):
    if k < .5 and h < ising_transition(k):
        return 0
    elif k > .5 and h < kt_transition(k):
        return 1
    return 2

# Generate the phase diagram data from DMRG transition lines
ks100, hs100 = np.linspace(0, 1, 100), np.linspace(0, 2, 100)
K100, H100 = np.meshgrid(ks100, hs100)
img_dmrg = np.vectorize(get_phase)(K100, H100)

colors = ['#80bfff', '#fff2a8',  '#80f090', '#da8080',]
phase_labels = ["Ferromagnetic", "Antiphase", "Paramagnetic", "Trash Class",]
cmap = ListedColormap(colors)

bounds = [-0.5, 0.5, 1.5, 2.5, 3.5]
norm = BoundaryNorm(bounds, cmap.N)

# Plot the phase diagram
plt.figure(figsize=(4, 4))
plt.imshow(img_dmrg, cmap=cmap, aspect="auto", origin="lower", extent=[0, 1, 0, 2], norm=norm)

# Plot the transition lines.
k_vals1 = np.linspace(0.0, 0.5, 50)
k_vals2 = np.linspace(0.5, 1.0, 50)
plt.plot(k_vals1, ising_transition(k_vals1), 'k', lw=2)
plt.plot(k_vals2, kt_transition(k_vals2), 'k', lw=2)
plt.plot(k_vals2, bkt_transition(k_vals2), 'k', ls = '--', lw=2)

# Create legend entries
for color, phase in zip(colors, phase_labels[:-1]):
    plt.scatter([], [], color=color, label=phase, edgecolors='black')
plt.plot([], [], 'k', label='Transition Lines')
plt.text(0.15, 0.5, s=r'$h_I$', fontsize=12)
plt.text(0.70, 0.5, s=r'$h_C$', fontsize=12)
plt.text(0.8, 0.2, s=r'$h_{BKT}$', fontsize=12)

plt.legend(), plt.xlabel("k"), plt.ylabel("h"), plt.title("Phase diagram"), plt.show();

• annni_pd_analytical.png

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm

# Kosterlitz-Thouless transition line
def kt_transition(k):
    return 1.05 * np.sqrt((k - 0.5) * (k - 0.1))

# Ising transition line
def ising_transition(k):
    return np.where(k == 0, 1, (1 - k) * (1 - np.sqrt((1 - 3 * k + 4 * k**2) / (1 - k))) / np.maximum(k, 1e-9))

# Floating Phase transition line
def bkt_transition(k):
    return 1.05 * (k - 0.5)
    
# Get the phase from the DMRG transition lines
def get_phase(k, h):
    if k < .5 and h < ising_transition(k):
        return 0
    elif k > .5 and h < kt_transition(k):
        return 1
    return 2

# Generate the phase diagram data from DMRG transition lines
ks100, hs100 = np.linspace(0, 1, 100), np.linspace(0, 2, 100)
K100, H100 = np.meshgrid(ks100, hs100)
img_dmrg = np.vectorize(get_phase)(K100, H100)

colors = ['#80bfff', '#fff2a8',  '#80f090', '#da8080',]
phase_labels = ["Ferromagnetic", "Antiphase", "Paramagnetic", "Trash Class",]
cmap = ListedColormap(colors)

bounds = [-0.5, 0.5, 1.5, 2.5, 3.5]
norm = BoundaryNorm(bounds, cmap.N)

# Plot the phase diagram
plt.figure(figsize=(4, 4))
plt.imshow(img_dmrg, cmap=cmap, aspect="auto", origin="lower", extent=[0, 1, 0, 2], norm=norm)

# Plot the transition lines.
k_vals1 = np.linspace(0.0, 0.5, 50)
k_vals2 = np.linspace(0.5, 1.0, 50)
plt.plot(k_vals1, ising_transition(k_vals1), 'k', lw=2)
plt.plot(k_vals2, kt_transition(k_vals2), 'k', lw=2)
plt.plot(k_vals2, bkt_transition(k_vals2), 'k', ls = '--', lw=2)

# Create legend entries
for color, phase in zip(colors, phase_labels[:-1]):
    plt.scatter([], [], color=color, label=phase, edgecolors='black')
plt.plot([], [], 'k', label='Transition Lines')

# Highlight analytical points (k = 0 and h = 0)
plt.scatter(np.zeros_like(hs100), hs100, color='red', s=20, label="Analytical Points")  # y-axis points
plt.scatter(ks100, np.zeros_like(ks100), color='red', s=20)  # x-axis point

plt.text(0.15, 0.5, s=r'$h_I$', fontsize=12)
plt.text(0.70, 0.5, s=r'$h_C$', fontsize=12)
plt.text(0.8, 0.2, s=r'$h_{BKT}$', fontsize=12)

plt.legend(), plt.xlabel("k"), plt.ylabel("h"), plt.title("Phase diagram"), plt.show();

[github-actions]

Thank you for opening this pull request.

You can find the built site at this link.

Deployment Info:

• Pull Request ID: 1335
• Deployment SHA: 557ab675881bec491fe4a31814d7d0ca67da83e1
  (The Deployment SHA refers to the latest commit hash the docs were built from)

Note: It may take several minutes for updates to this pull request to be reflected on the deployed site.

[daniela-angulo]

This was a joy to read. I added a few comments/suggestions.

[daniela-angulo]

I am just curious about the minus sign in the front. The Hamiltonian from the paper has a positive sign.

[SaverioMonaco]

The signs in the ANNNI Hamiltonian have always been a source of confusion. The important thing is getting the relative sign between the nearest-neighbor and next-nearest-neighbor interactions right, otherwise, there is no competition between them, which leads to no actual phase transition (I know this firsthand, it cost me a week of debugging at the time).

The original (1988) paper defines the hamiltonian as (Eq. 2.1):
$$H = - J_1\Big(\sum_i \sigma_x^i \sigma_x^{i+1} - \big( \underbrace{-\frac{J_2 }{J_1}}_{\kappa} \big) \sigma_x^i \sigma_x^{i+1} ... \Big) = - J_1\Big(\sum_i \sigma_x^i \sigma_x^{i+1} - \kappa \sigma_x^i \sigma_x^{i+1} ... \Big)$$
For $J_1 > 0$ and $J_2 \leq 0$ hence $\kappa \geq 0$

In our paper we defined it as (Eq. 1)
$$H = J_1 \sum_i \sigma_x^i\sigma_x^{i+1} - \kappa \sigma_x^{i}\sigma_x^{i+2} + h \sigma_z^i$$

I assume in our paper we forgot to add the minus sign before the nearest-neighbor term. Thank you for pointing that out!

I would personally leave as it is the equation in the notebook, as it is the same formulation of the original 1988 paper.

[daniela-angulo]

Since these two methods, VQE and MPS, weren't used in this demo, it would be best not to go into details about them here. We can mention them and add references if needed, but it's best to not emphasize them.

[daniela-angulo]

it would be helpful to add a sentence to define \tilde{p}_j as well

[SaverioMonaco]

I added in the commit a better explanation of the loss function and all its terms

[daniela-angulo]

The notation here looks tricky. One of the states has the \otimes K, which is clear, but for \phi is not so clear what the superscripts N-K mean.

[SaverioMonaco]

I rephrased the whole paragraph in the commit, I hope it is now clearer!

[daniela-angulo]

Suggested change

	# After training the circuit to optimally compress the (0,0)(0,0) state, we evaluate the compression score for all other input states using the learned parameters.
	# After training the circuit to optimally compress the (0,0) state, we evaluate the compression score for all other input states using the learned parameters.

[daniela-angulo]

Suggested change

	# In this section, we prepare the ground states of the system, which will serve as inputs for both QML models. Several methods can be used:
	#
	# * **Variational Quantum Eigensolver (VQE)**
	#
	# VQE is a technique introduced in [#Peruzzo]_, which leverages the Rayleigh-Ritz variational principle to approximate the lowest-energy states of a given Hamiltonian, as demonstrated in the :doc:`demo on VQE </demos/tutorial_vqe>`.
	#
	# * **Matrix Product States (MPS)**
	#
	# MPS can be efficiently computed on classical hardware at a low cost and provide accurate approximations for local quantum systems. Several techniques are being developed to optimally embed MPS into quantum circuits as seen for example in the :doc:`demo on Constant-depth preparation of MPS with dynamic circuits </demos/tutorial_constant_depth_mps_prep>`.
	#
	#
	# For simplicity, in this demo, we compute the ground state directly by finding the *eigenvector* corresponding to the lowest eigenvalue of the Hamiltonian. The resulting states are then loaded into the quantum circuits using PennyLane’s :class:`~pennylane.StatePrep`.
	# In this section, we prepare the ground states of the system, which will serve as inputs for both QML models. Several methods can be used, including **Variational Quantum Eigensolver (VQE)**, introduced in [#Peruzzo]_ and demonstrated in the :doc:`demo on VQE </demos/tutorial_vqe>`, and **Matrix Product States (MPS)**, illustrated in the :doc:`demo on Constant-depth preparation of MPS with dynamic circuits </demos/tutorial_constant_depth_mps_prep>`.
	# For simplicity, in this demo, we compute the ground state directly by finding the *eigenvector* corresponding to the lowest eigenvalue of the Hamiltonian. The resulting states are then loaded into the quantum circuits using PennyLane’s :class:`~pennylane.StatePrep`.