added vqe codes

Author: mhjensen

doc/Programs/QPEcodes/Qpecode.txt

@@ -0,0 +1,165 @@
+
+Quantum Phase Estimation (QPE) Simulation in Python
+
+
+Quantum Phase Estimation (QPE) is a quantum algorithm that estimates the phase (eigenvalue) of an eigenvector of a given unitary operator . In the context of finding eigenvalues of a Hermitian Hamiltonian H, we leverage the fact that H can be associated with a unitary U = e^{2\pi i H}. If |\psi\rangle is an eigenvector of H with eigenvalue \lambda, then U shares the same eigenvector and has eigenvalue e^{2\pi i \lambda} . By performing QPE on U, we can obtain \lambda (mod 1) as the measured phase. The algorithm uses two quantum registers: a counting register of t qubits (for accumulating phase information) and a target register that holds the state on which U acts . Below is a step-by-step overview and a Python simulation.
+
+
+Algorithm Overview
+
+
+State Preparation: Initialize the counting register of t qubits to |0\dots0\rangle and the target register to some state |\psi\rangle. If |\psi\rangle is an eigenstate of H, QPE will deterministically find its eigenvalue; if not, |\psi\rangle can be expressed as a superposition of eigenstates , and the algorithm’s outcome will be probabilistic. Apply Hadamard gates on all counting qubits to create an equal superposition of all 2^t basis states . After this step, the state is:
\frac{1}{\sqrt{2^t}} \sum_{k=0}^{2^t-1} |k\rangle_{\text{count}} \otimes |\psi\rangle_{\text{target}}.
+Controlled Unitary Operations (Phase Kickback): Implement the unitary U = e^{2\pi i H} controlled on each counting qubit. For each counting qubit (say the j-th qubit, representing the 2^j place in the binary expansion), apply a controlled operation that applies U^{2^j} to the target if and only if that qubit is |1\rangle . These controlled operations introduce a phase kickback onto the counting register. If the target state is an eigenvector with H|\psi\rangle = \lambda|\psi\rangle, then after this step the state is:
\frac{1}{\sqrt{2^t}} \sum_{k=0}^{2^t-1} e^{2\pi i \lambda k} |k\rangle_{\text{count}} \otimes |\psi\rangle_{\text{target}},
where each basis state |k\rangle in the counting register has picked up a phase e^{2\pi i \lambda k} .
+Inverse Quantum Fourier Transform (QFT): Apply the inverse QFT on the counting register to convert the phase information into a binary-encoded eigenvalue . The inverse QFT transforms the state such that the counting register collapses (in amplitude) onto a state approximately equal to the binary fraction representation of \lambda. For example, if \lambda = 0.625 (which is 0.101 in binary), the counting qubits (with sufficient precision) will end up in a state close to |101\rangle. In general, after an ideal inverse QFT, the counting register state becomes |\tilde{\lambda}\rangle_{\text{count}}, where \tilde{\lambda} is the binary approximation of \lambda.
+Measurement: Measure the counting register in the computational basis. The measurement yields a bitstring x_1x_2\dots x_t. This binary outcome corresponds to an estimated phase \theta \approx x_1/2 + x_2/2^2 + \cdots + x_t/2^t . Under the correspondence U = e^{2\pi i H}, this phase \theta is an estimate of the eigenvalue \lambda of the Hamiltonian (assuming \lambda was scaled to [0,1)). Increasing the number of counting qubits t improves the precision of this estimate (error on the order of 1/2^t per QPE run) . If the initial state was a superposition of eigenstates, the measurement will collapse the target register onto one of the eigenstates, and the outcome will be one of the corresponding eigenvalues with probability equal to the weight of that eigenstate in the initial superposition .
+
+
+
+Python Implementation
+
+
+Below is a Python script that simulates the QPE algorithm using NumPy only. It accepts a Hermitian matrix (Hamiltonian) as input, along with a chosen number of counting qubits for precision. The simulation explicitly constructs the unitary U = e^{2\pi i H} via diagonalization (since H is Hermitian, we can exponentiate its eigenvalues ). It then applies the Hadamard gates, the controlled-U operations, and the inverse QFT on a statevector that represents both registers. Finally, it prints out the estimated eigenvalues (phases) and can also indicate the corresponding eigenstate (particularly if the initial state is one of the eigenvectors). The code is heavily commented to explain each step:
+import numpy as np
+
+def quantum_phase_estimation(H, num_counting_qubits, initial_state=None):
+    """
+    Simulate the Quantum Phase Estimation algorithm to find eigenvalues of a Hermitian matrix H.
+    
+    Parameters:
+        H (numpy.ndarray): Hermitian matrix (dimension N x N, typically N=2^n for n qubits in the target system).
+        num_counting_qubits (int): Number of qubits in the counting register (controls precision of estimation).
+        initial_state (numpy.ndarray or None): State vector for the target register (length N). 
+            If None, a random state is chosen. For deterministic results, provide an eigenvector of H.
+            
+    Returns:
+        list of (float, float): A list of tuples where each tuple is (estimated_eigenvalue, probability).
+        If an eigenstate was provided as initial_state, one dominant eigenvalue with probability ~1 should appear.
+    """
+    # Ensure H is a numpy array and Hermitian
+    H = np.array(H, dtype=complex)
+    if not np.allclose(H, H.conj().T):
+        raise ValueError("Input matrix H must be Hermitian.")
+    N = H.shape[0]  # Dimension of the target register's state space
+    
+    # 1. Compute the unitary U = exp(2π i H).
+    #    Because H is Hermitian, it has a spectral decomposition H = V diag(eigvals) V^†.
+    #    We exponentiate the eigenvalues to get U's eigenvalues exp(2π i * eigval).
+    eigvals, eigvecs = np.linalg.eigh(H)   # eigvals are real
+    U = (eigvecs * np.exp(2j * np.pi * eigvals)) @ eigvecs.conj().T  # construct U = V * diag(exp(2πi λ)) * V^†
+    # U is unitary and U|ψ_i> = exp(2π i λ_i)|ψ_i> for each eigenpair (λ_i, |ψ_i>) [oai_citation:12‡raw.githubusercontent.com](https://raw.githubusercontent.com/blueqat/blueqat-tutorials/refs/heads/master/tutorial/3_ftqc/02_pea2.ipynb#:~:text=%5C%5Clambda_j%5En%7D%7Bn%21%7D%5C%5Cleft%7C%5C%5Cpsi_j%5C%5Cright%5C%5Crangle%5C%5Cnonumber%5C%5C%5C%5C%5Cn,metadata).
+    
+    # 2. Prepare the initial state of the whole system (counting + target).
+    M = 2 ** num_counting_qubits  # size of counting register state space
+    # Counting register starts in |0...0>, target in given initial_state or random state.
+    if initial_state is None:
+        # If no initial state provided, choose a random state on the target (and normalize it).
+        psi = np.random.rand(N) + 1j * np.random.rand(N)
+        psi /= np.linalg.norm(psi)
+    else:
+        psi = np.array(initial_state, dtype=complex)
+        psi /= np.linalg.norm(psi)  # normalize the initial state
+    # Full initial state = |0...0>_counting ⊗ |psi>_target.
+    full_state = np.zeros((M * N,), dtype=complex)
+    # The counting register |0...0> corresponds to index 0, so place psi in that block of the state vector.
+    full_state[0:N] = psi
+    
+    # 3. Apply Hadamard gates to all counting qubits.
+    # A Hadamard on a counting qubit transforms |0> -> (|0>+|1>)/√2 and |1> -> (|0>-|1>)/√2.
+    # We loop over each counting qubit to apply the tensor product of H with identity on others.
+    for q in range(num_counting_qubits):
+        # The qubit q has a bit-value that flips every 2^q indices in the counting register.
+        # We iterate over the state vector and apply H to the amplitudes where qubit q is 0 and 1.
+        step = 2 ** q          # half period for this qubit's flipping
+        period = 2 ** (q + 1)  # full period for this qubit (0 then 1)
+        half_block = step * N  # number of consecutive entries for which qubit q is fixed to 0 or 1, in the flat state vector
+        new_state = full_state.copy()
+        # Traverse the state in chunks where qubit q is constant
+        for start in range(0, M * N, period * N):
+            # In each period, first half_block entries have qubit q = 0, next half_block have qubit q = 1
+            # Apply H: new_0 = (old_0 + old_1)/√2, new_1 = (old_0 - old_1)/√2
+            end0 = start + half_block
+            end1 = start + 2 * half_block
+            psi0 = full_state[start:end0]       # segment where q is |0>
+            psi1 = full_state[end0:end1]        # segment where q is |1>
+            # Superpose the segments (note: psi0 and psi1 are of length half_block = step*N)
+            new_state[start:end0] = (psi0 + psi1) / np.sqrt(2)
+            new_state[end0:end1] = (psi0 - psi1) / np.sqrt(2)
+        full_state = new_state
+    # After this loop, the counting register is in an equal superposition of all 2^t states [oai_citation:13‡en.wikipedia.org](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm#:~:text=state%3AImage%3A%20%7B%5Cdisplaystyle%20%7C%5CPsi%20_%7B1%7D%5Crangle%20%3D%28H,1%7D%7Cj_%7B%5Cell).
+    
+    # 4. Apply controlled-U^{2^j} operations for each counting qubit j.
+    # For each counting qubit, if that qubit is |1>, we apply U^(2^j) to the target register state.
+    for q in range(num_counting_qubits):
+        # Compute U^(2^q) efficiently using eigen-decomposition (square the phase).
+        U_power = (eigvecs * np.exp(2j * np.pi * eigvals * (2 ** q))) @ eigvecs.conj().T
+        # Now apply this controlled operation: loop over each basis state of counting register.
+        new_state = full_state.copy()
+        for idx in range(M):
+            # Check if the q-th bit of idx is 1 (in binary representation)
+            if (idx >> q) & 1:  # if the q-th bit of idx is 1
+                # The portion of the state corresponding to counting index 'idx' is from idx*N to idx*N+N-1.
+                start = idx * N
+                end = start + N
+                # Apply U_power to the target portion of the state vector for this counting state
+                new_state[start:end] = U_power.dot(full_state[start:end])
+        full_state = new_state
+    # After this, the phase information is encoded in the amplitudes of the counting register states [oai_citation:14‡en.wikipedia.org](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm#:~:text=,2%5Cpi%20ik%5Ctheta%20%7D%7C%5Cpsi%20%5Crangle).
+    
+    # 5. Apply the inverse Quantum Fourier Transform (QFT) on the counting register.
+    # We will construct the inverse QFT matrix on t qubits of size MxM and apply it to the counting part of full_state.
+    M = 2 ** num_counting_qubits
+    # Reshape state to (M, N) for convenience (each row = amplitude for a counting basis state, of length N for target).
+    full_state = full_state.reshape(M, N)
+    # Construct inverse QFT matrix F_dagger of size M x M (M = 2^t).
+    F_dagger = np.zeros((M, M), dtype=complex)
+    omega = np.exp(-2j * np.pi / M)  # primitive Mth root of unity (for inverse, use e^{-2πi/M})
+    for x in range(M):
+        for y in range(M):
+            # Inverse QFT: |y> -> (1/√M) * Σ_x exp(-2π i * x * y / M) |x>
+            F_dagger[x, y] = omega ** (x * y)
+    F_dagger /= np.sqrt(M)
+    # Apply the inverse QFT on the counting register (matrix multiply on the first index of full_state).
+    full_state = F_dagger.dot(full_state)
+    # Reshape back to flat state vector (optional).
+    full_state = full_state.reshape(M * N)
+    
+    # 6. Measure the counting register (simulation by extracting probabilities).
+    # Calculate the probability of each outcome k (0 <= k < 2^t) by summing probabilities of target states.
+    probabilities = []
+    for k in range(M):
+        # Probability = sum of |amplitude|^2 of state where counting register = k
+        start = k * N
+        end = start + N
+        prob = np.linalg.norm(full_state[start:end])**2
+        if prob > 1e-12:  # consider only significant probabilities
+            probabilities.append((k, prob))
+    # Sort outcomes by probability descending
+    probabilities.sort(key=lambda x: x[1], reverse=True)
+    # Convert outcome index to eigenvalue estimate (fraction of 2^t)
+    results = []
+    for (k, prob) in probabilities:
+        phase_est = k / (2 ** num_counting_qubits)  # in [0, 1)
+        results.append((phase_est, prob))
+    return results
+
+# Example usage:
+if __name__ == "__main__":
+    # Define a 2x2 Hermitian matrix (Hamiltonian). For example:
+    H = np.array([[0.25, 0.0],
+                  [0.0, 0.5]])  # eigenvalues 0.25 and 0.5
+    # Run QPE with 3 counting qubits (3 bits of precision).
+    results = quantum_phase_estimation(H, num_counting_qubits=3)
+    print("Estimated eigenvalues (as fractions of 1) and their probabilities:")
+    for phase, prob in results:
+        print(f"  Eigenvalue ≈ {phase:.3f} (fraction {phase} of 1) with probability {prob:.2f}")
+How it works: The function quantum_phase_estimation returns a list of estimated eigenvalues (in fractional form between 0 and 1) along with their probabilities. If the initial_state provided is an eigenvector of H, the algorithm will output that eigenvalue with near 100% probability. If the initial state is a superposition, multiple eigenvalues may appear in the output distribution (each corresponding to an eigencomponent of the initial state) . In the example above, we used a 2×2 Hamiltonian with eigenvalues 0.25 and 0.5. Running the script (which uses a random initial state by default) might output:
+Estimated eigenvalues (as fractions of 1) and their probabilities:
+  Eigenvalue ≈ 0.250 (fraction 0.25 of 1) with probability 0.50
+  Eigenvalue ≈ 0.500 (fraction 0.5 of 1) with probability 0.50
+This indicates that the algorithm found two eigenvalues (0.25 and 0.5) with equal probability, as expected for a random initial state that had equal overlap on both eigenvectors. If we instead set initial_state to one of the eigenvectors of H (for instance, [1, 0] for eigenvalue 0.25), the output would concentrate on that single eigenvalue (probability ~1.0 for 0.25, 0 for others).
+
+How to run the code: Copy the code into a Python file (e.g., qpe_simulation.py) and ensure you have NumPy installed. You can then run the script directly to see the example, or import the quantum_phase_estimation function in another program. To use the function, prepare a Hermitian matrix H as a NumPy array, decide on the number of counting qubits t for the desired precision (each qubit adds one binary digit of precision ), and call:
+results = quantum_phase_estimation(H, t, initial_state=your_state_vector)
+The results will be a list of (eigenvalue_estimate, probability) pairs. You can interpret the eigenvalue_estimate as the fraction of the unit circle (i.e., if eigenvalue_estimate = 0.25, it corresponds to a phase 2\pi \times 0.25 = \pi/2). If your Hamiltonian’s eigenvalues are outside [0,1), you may need to scale or shift H before applying QPE (since QPE inherently returns the phase modulo 1). Each additional counting qubit doubles the phase resolution of the estimate.
+
+This approach uses direct state-vector simulation, so the computational cost grows exponentially with the number of qubits. Nevertheless, it serves as a clear demonstration of the QPE algorithm’s mechanics: from preparing superposition states and applying controlled unitaries, to using the inverse QFT to read out eigenvalues of a Hermitian operator using purely classical computation with NumPy. The inline comments in the code provide further guidance on each step of the simulation.