[Bug fix] Update fixed depth hamiltonian simulation demo (#1344)

The problem is that the computed horizontal CSA $\mathfrak{h}$ is not
necessarily composed of just Pauli words, but soms of Paulis (Pauli
sentence). The old version was assuming everything was Pauli, did the
decomposition in terms of Pauli words, and accordingly counted elements.

Updating this here to decompose in the actual basis of the horizontal
CSA $\mathfrak{h}$ and show explicitly that the decomposition is stricly
in it.

edit: this works all well locally with correct jax, optax and pl
version, but for some reason fails in CI 🅰️ 🅰️ 🅰️
edit2: increasing the tolerance in `horizontal_cartan_subalgebra` fixes
this. A bit problematic that this cannot be reproduced locally 😢

---------

Co-authored-by: Yushao Chen (Jerry) <chenys13@outlook.com>
Co-authored-by: David Wierichs <davidwierichs@gmail.com>

Author: 3 people

demonstrations/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition.metadata.json

@@ -6,7 +6,7 @@
         }
     ],
     "dateOfPublication": "2024-12-19T00:00:00+00:00",
-    "dateOfLastModification": "2025-03-26T00:00:00+00:00",
+    "dateOfLastModification": "2025-04-11T00:00:00+00:00",
     "categories": [
         "Quantum Computing",
         "Algorithms"

demonstrations/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition.py

@@ -169,7 +169,7 @@
 # :func:`~.pennylane.liealg.horizontal_cartan_subalgebra` returns some additional information, which we will
 # not use here.
 
-g, k, mtilde, h, _ = horizontal_cartan_subalgebra(k, m)
+g, k, mtilde, h, _ = horizontal_cartan_subalgebra(k, m, tol=1e-8)
 len(g), len(k), len(mtilde), len(h)
 
 ##############################################################################
@@ -220,10 +220,10 @@
 #
 
 def run_opt(
-    value_and_grad,
+    loss,
     theta,
-    n_epochs=500,
-    lr=0.1,
+    n_epochs=1000,
+    lr=0.05,
 ):
     """Boilerplate JAX optimization"""
     value_and_grad = jax.jit(jax.value_and_grad(loss))
@@ -274,7 +274,7 @@ def loss(theta):
     A = K_m @ v_m @ K_m.conj().T
     return jnp.trace(A.conj().T @ H_m).real
 
-theta0 = jnp.ones(len(k), dtype=float)
+theta0 = jax.random.normal(jax.random.PRNGKey(0), shape=(len(k),), dtype=float)
 
 thetas, energy, _ = run_opt(loss, theta0, n_epochs=1000, lr=0.05)
 plt.plot(energy - np.min(energy))
@@ -297,13 +297,18 @@ def loss(theta):
 # .. math:: h_0 = K_c^\dagger H K_c.
 
 h_0_m = Kc_m.conj().T @ H_m @ Kc_m
-h_0 = qml.pauli_decompose(h_0_m)
 
-print(len(h_0))
+# decompose h_0_m in terms of the basis of h
+basis = [qml.matrix(op, wire_order=range(n_wires)) for op in h]
+coeffs = qml.pauli.trace_inner_product(h_0_m, basis)
+
+# ensure that decomposition is correct, i.e. h_0_m is truely an element of just h
+h_0_m_recomposed = np.sum([c * op for c, op in zip(coeffs, basis)], axis=0)
+print("Decomposition of h_0 is faithful: ", np.allclose(h_0_m_recomposed, h_0_m, atol=1e-10))
+
+# sanity check that the horizontal CSA is Abelian, i.e. all its elements commute
+print("All elements in h commute with each other: ", qml.liealg.check_abelian(h))
 
-# assure that h_0 is in \mathfrak{h}
-h_vspace = qml.pauli.PauliVSpace(h)
-not h_vspace.is_independent(h_0.pauli_rep)
 
 ##############################################################################
 #
@@ -331,6 +336,7 @@ def loss(theta):
 t = 1.
 U_exact = qml.exp(-1j * t * H)
 U_exact_m = qml.matrix(U_exact, wire_order=range(n_wires))
+h_0  = qml.dot(coeffs, h)
 
 def U_kak(theta_opt, t):
     qml.adjoint(K)(theta_opt, k)