🧪 MLIR - Add QPE and Iterative Phase Estimation examples

[DRovara]

Background

We have manual MLIR sketches using scf/arith for classical control and memref<?xi1> for classical registers. See the first two comments below.

Problem

We need runnable, reviewed examples that exercise classical-quantum interplay and passes.

Proposal:

• Commit .mlir examples for QPE and IPE in MQTRef; include tests for parse/roundtrip and simple analysis checks; optional simulator hook to validate expected classical results.

Acceptance criteria:

• Examples and tests present; comments explain design choices.

This is one of the first steps in elaborating how to represent hybrid quantum-classical workflows.

[DRovara]

One potential thing to look at is how we would ideally envision e.g. OpenQASM 3 code in our MLIR dialects. For this, let us first look at the "Iterative Quantum Phase Estimation" algorithm implemented in OpenQASM 3 in its introductory paper:

/*
 * Iterative phase estimation
 */
OPENQASM 3;
include "stdgates.inc";

const n = 3; // number of iterations
const θ= 3 * π / 8; // phase angle on target qubit

qubit q; // phase estimation qubit
qubit r; // target qubit for the controlled-unitary gate
angle[n] c = 0; // phase estimation bits

// initialize
reset q;
reset r;

// prepare uniform superposition of eigenvectors of phase
h r;

// iterative phase estimation loop
for k in [1:n] { // implicitly cast val to int
  reset q;
  h q;
  ctrl @ pow(2**k) @ phase(θ) q, r;
  inv @ phase(c) q;
  h q;
  measure q -> c[0];
  // newest measurement outcome is associated to a π/2 phase shift
  // in the next iteration, so shift all bits of c left
  c <<= 1;
}

// Now c contains the n-bit estimate of ϕ in the
// eigenvalue e^{i*ϕ} and qreg r is projected to an
// approximate eigenstate of the phase gate.

One potential representation in MQTOpt might be:

module {
  func.func @iterativeQuantumPhaseEstimation() -> f64 {
    %one = arith.constant 1 : i64 // a constant of value 1 used for some calculations
    %two_pi = arith.constant 6.283185307179586 : f64 // a constant of value 2 * pi used for angle calculations
    %two_f = arith.constant 2.0 : f64 // a constant of value 2.0 used for some calculations

    %n = arith.constant 3 : i64
    %n_f = arith.sitofp %n : i64 to f64 // convert the number of iterations to a floating-point value
    %theta = arith.constant 1.1780972450961724 : f64
    %angle_step = arith.divf %two_pi, %n_f : f64 // calculate the angle step for the quantum phase estimation algorithm

    %reg_0 = "mqtopt.allocQubitRegister"() <{size_attr = 2 : i64}> : () -> !mqtopt.QubitRegister
    
    %reg_1, %q_0 = "mqtopt.extractQubit"(%reg_0) <{index_attr = 0 : i64}> : (!mqtopt.QubitRegister) -> (!mqtopt.QubitRegister, !mqtopt.Qubit)
    %reg_2, %r_0 = "mqtopt.extractQubit"(%reg_1) <{index_attr = 1 : i64}> : (!mqtopt.QubitRegister) -> (!mqtopt.QubitRegister, !mqtopt.Qubit)
    %c_0 = arith.constant 0.0 : f64

    %r_1 = mqtopt.h() %r_0 : !mqtopt.Qubit

    %from = arith.constant 1 : index
    %to_i64 = arith.addi %n, %one : i64 // upper bound is exclusive, so we add 1
    %to = arith.index_cast %to_i64 : i64 to index // convert the number of iterations to an index type
    %step = arith.constant 1 : index // a constant step value of 1 for the loop
    %q_1, %r_2, %c_1 = scf.for %k = %from to %to step %step iter_args(%q_current = %q_0, %r_current = %r_1, %c_current = %c_0) -> (!mqtopt.Qubit, !mqtopt.Qubit, f64) {
        // >> reset q;
        %q_1 = "mqtopt.reset"(%q_current) : (!mqtopt.Qubit) -> (!mqtopt.Qubit)
        
        // >> h q;
        %q_2 = mqtopt.h() %q_1 : !mqtopt.Qubit

        // >> cp(2^k * theta) q, r
        %k_i = arith.index_cast %k : index to i64 // convert the loop index to an i64 type
        %angle_0 = arith.shli %one, %k_i : i64
        %angle_1 = arith.sitofp %angle_0 : i64 to f64
        %angle = arith.mulf %angle_1, %theta : f64
        %r_new, %q_3 = mqtopt.p(%angle) %r_current ctrl %q_2 : !mqtopt.Qubit ctrl !mqtopt.Qubit

        // >> inv @ phase(c) q;
        %c_neg = arith.negf %c_current : f64
        %q_4 = mqtopt.p(%c_neg) %q_3 : !mqtopt.Qubit

        // >> h q;
        %q_5 = mqtopt.h() %q_4 : !mqtopt.Qubit

        // >> measure q -> c[0]
        //      c[0] is the least-signinifcant bit of the current angle `c`.
        //      A single bit contributes 2pi/n = %angle_step to the angle.
        //      Therefore, if %m = 1, we just add %angle_step to the current angle.
        %q_6, %m = "mqtopt.measure"(%q_5) : (!mqtopt.Qubit) -> (!mqtopt.Qubit, i1)
        %c_1 = scf.if %m -> f64 {
          %c_if = arith.addf %c_current, %angle_step : f64
          scf.yield %c_if : f64
        } else {
          scf.yield %c_current : f64
        }

        // >> `c <<= 1`
        //      This effectively doubles the angle but gives the result modulo 2pi.
        %c_2 = arith.mulf %c_1, %two_f : f64
        %larger = arith.cmpf oge, %c_2, %two_pi : f64
        %c_3 = scf.if %larger -> f64 {
          // If the angle is larger than 2pi, we subtract 2pi to keep it in the range [0, 2pi).
          %c_if = arith.subf %c_2, %two_pi : f64
          scf.yield %c_if : f64
        } else {
          // If the angle is smaller than 2pi, we keep it as is.
          scf.yield %c_2 : f64
        }

        scf.yield %q_6, %r_new, %c_3 : !mqtopt.Qubit, !mqtopt.Qubit, f64
    }

    %reg_3 = "mqtopt.insertQubit"(%reg_2, %q_1) <{index_attr = 0 : i64}> : (!mqtopt.QubitRegister, !mqtopt.Qubit) -> !mqtopt.QubitRegister
    %reg_4 = "mqtopt.insertQubit"(%reg_3, %r_2) <{index_attr = 1 : i64}> : (!mqtopt.QubitRegister, !mqtopt.Qubit) -> !mqtopt.QubitRegister
    "mqtopt.deallocQubitRegister"(%reg_4) : (!mqtopt.QubitRegister) -> ()

    // Now c contains the n-bit estimate of ϕ in the
    // eigenvalue e^{i*ϕ} and qreg r is projected to an
    // approximate eigenstate of the phase gate.
    return %c_1 : f64
  }
}

or, using MQTRef:

module {
  func.func @iterativeQuantumPhaseEstimation() -> f64 {
    %one = arith.constant 1 : i64 // a constant of value 1 used for some calculations
    %two_pi = arith.constant 6.283185307179586 : f64 // a constant of value 2 * pi used for angle calculations
    %two_f = arith.constant 2.0 : f64 // a constant of value 2.0 used for some calculations

    %n = arith.constant 3 : i64
    %n_f = arith.sitofp %n : i64 to f64 // convert the number of iterations to a floating-point value
    %theta = arith.constant 1.1780972450961724 : f64
    %angle_step = arith.divf %two_pi, %n_f : f64 // calculate the angle step for the quantum phase estimation algorithm

    %reg = "mqtref.allocQubitRegister"() <{size_attr = 2 : i64}> : () -> !mqtref.QubitRegister
    
    %q = "mqtref.extractQubit"(%reg) <{index_attr = 0 : i64}> : (!mqtref.QubitRegister) -> !mqtref.Qubit
    %r = "mqtref.extractQubit"(%reg) <{index_attr = 1 : i64}> : (!mqtref.QubitRegister) -> !mqtref.Qubit
    %c_0 = arith.constant 0.0 : f64

    mqtref.h() %r

    %from = arith.constant 1 : index
    %to_i64 = arith.addi %n, %one : i64 // upper bound is exclusive, so we add 1
    %to = arith.index_cast %to_i64 : i64 to index // convert the number of iterations to an index type
    %step = arith.constant 1 : index // a constant step value of 1 for the loop
    %c_1 = scf.for %k = %from to %to step %step iter_args(%c_current = %c_0) -> f64 {
        // >> reset q;
        "mqtref.reset"(%q) : (!mqtref.Qubit) -> ()
        
        // >> h q;
        mqtref.h() %q

        // >> cp(2^k * theta) q, r
        %k_i = arith.index_cast %k : index to i64 // convert the loop index to an i64 type
        %angle_0 = arith.shli %one, %k_i : i64
        %angle_1 = arith.sitofp %angle_0 : i64 to f64
        %angle = arith.mulf %angle_1, %theta : f64
        mqtref.p(%angle) %r ctrl %q

        // >> inv @ phase(c) q;
        %c_neg = arith.negf %c_current : f64
        mqtref.p(%c_neg) %q

        // >> h q;
        mqtref.h() %q

        // >> measure q -> c[0]
        //      c[0] is the least-signinifcant bit of the current angle `c`.
        //      A single bit contributes 2pi/n = %angle_step to the angle.
        //      Therefore, if %m = 1, we just add %angle_step to the current angle.
        %m = "mqtref.measure"(%q) : (!mqtref.Qubit) -> (i1)
        %c_1 = scf.if %m -> f64 {
          %c_if = arith.addf %c_current, %angle_step : f64
          scf.yield %c_if : f64
        } else {
          scf.yield %c_current : f64
        }

        // >> `c <<= 1`
        //      This effectively doubles the angle but gives the result modulo 2pi.
        %c_2 = arith.mulf %c_1, %two_f : f64
        %larger = arith.cmpf oge, %c_2, %two_pi : f64
        %c_3 = scf.if %larger -> f64 {
          // If the angle is larger than 2pi, we subtract 2pi to keep it in the range [0, 2pi).
          %c_if = arith.subf %c_2, %two_pi : f64
          scf.yield %c_if : f64
        } else {
          // If the angle is smaller than 2pi, we keep it as is.
          scf.yield %c_2 : f64
        }

        scf.yield %c_3 : f64
    }

    "mqtref.deallocQubitRegister"(%reg) : (!mqtref.QubitRegister) -> ()

    // Now c contains the n-bit estimate of ϕ in the
    // eigenvalue e^{i*ϕ} and qreg r is projected to an
    // approximate eigenstate of the phase gate.
    return %c_1 : f64
  }
}

[burgholzer]

Thanks for creating this issue and providing a first potential solution!
I tried to not look at your solution before proposing my own version based on some discussions I had with Claude Sonnet 4. I purposely tried to make this rather generic so that it also works if n is actually a runtime value to ensure that this is scalable.
Those discussions ended up in

func.func @phase_estimation(%n: i64, %theta: f64) -> memref<?xi1> {
  %c0 = arith.constant 0 : index
  %c1 = arith.constant 1 : index
  %c2_i32 = arith.constant 2 : i32

  // Convert n to index for loop bounds
  %n_index = arith.index_cast %n : i64 to index
  
  // For inclusive range [1:n], we need to go from 1 to n+1 (exclusive upper bound)
  %n_plus_1 = arith.addi %n_index, %c1 : index

  // Static qubits
  %q = mqtref.qubit 0
  %r = mqtref.qubit 1
  
  // Use memref for classical angle storage with dynamic size
  %c_angle = memref.alloc(%n_index) : memref<?xi1>

  // Initialize qubits
  mqtref.reset %q
  mqtref.reset %r

  // Prepare uniform superposition
  mqtref.h %r

  // Iterative phase estimation loop for k in [1:n] (inclusive)
  scf.for %k = %c1 to %n_plus_1 step %c1 {
    mqtref.reset %q
    mqtref.h %q

    // Controlled phase gate with power 2^k using floating point theta
    %k_i32 = arith.index_cast %k : index to i32
    %power = math.ipowi %c2_i32, %k_i32 : i32
    pow(%power) @ mqtref.p(%theta) %r ctrl %q

    // Compute accumulated phase from all bits
    %pi = arith.constant 3.141592653589793 : f64
    %two_f64 = arith.constant 2.0 : f64
    %zero_f64 = arith.constant 0.0 : f64

    %accumulated_phase = scf.for %i = %c0 to %n_index step %c1 iter_args(%acc = %zero_f64) -> f64 {
      %bit = memref.load %c_angle[%i] : memref<?xi1>
      %bit_f64 = arith.uitofp %bit : i1 to f64

      // Calculate π/2^(i+1) for the phase contribution
      %exp = arith.addi %i, %c1 : index
      %exp_i32 = arith.index_cast %exp : index to i32
      %power_of_2 = math.ipowi %c2_i32, %exp_i32 : i32
      %power_f64 = arith.sitofp %power_of_2 : i32 to f64

      %phase_unit = arith.divf %pi, %power_f64 : f64
      %phase_contrib = arith.mulf %bit_f64, %phase_unit : f64
      %new_acc = arith.addf %acc, %phase_contrib : f64

      scf.yield %new_acc : f64
    }

    inv @ mqtref.p(%accumulated_phase) %q

    mqtref.h %q

    // Measure into c[0], then shift left
    %measurement = mqtref.measure %q

    // Store measurement at position 0
    memref.store %measurement, %c_angle[%c0] : memref<?xi1>

    // Perform left shift: move all bits one position right
    // Start from the end to avoid overwriting values we still need
    %n_minus_1 = arith.subi %n_index, %c1 : index
    scf.for %i = %n_minus_1 to %c1 step %-1 {
      %prev_i = arith.subi %i, %c1 : index
      %bit = memref.load %c_angle[%prev_i] : memref<?xi1>
      memref.store %bit, %c_angle[%i] : memref<?xi1>
    }
  }

  // Return the resulting angle register
  return %c_angle : memref<?xi1>
}

// Example usage for the original case (n=3, θ=3π/8)
func.func @phase_estimation_example() -> memref<?xi1> {
  %n = arith.constant 3 : i64
  
  // Compute θ = 3π/8 explicitly
  %pi = arith.constant 3.141592653589793 : f64
  %three = arith.constant 3.0 : f64
  %eight = arith.constant 8.0 : f64
  %theta_num = arith.mulf %three, %pi : f64
  %theta = arith.divf %theta_num, %eight : f64

  %result = func.call @phase_estimation(%n, %theta) : (i64, f64) -> memref<?xi1>
  return %result : memref<?xi1>
}

I generally agree that this gets complicated rather quickly. Especially due to the angle computations as well as the bit shifting of the registers. Maybe the iterative QFT is a better example because it is simpler (but probably not).

The key things I see as potentially interesting:

• the use of the memref MLIR builtin for referring to classical values. Seems to be a common practice and also what QIRO is doing.
• the definition of inv and pow modifiers in MLIR to more directly mirror the OpenQASM capabilities. We already have the inv functionality on the QuantumComputation side and had the power functionality planned for unitary hack. However, we generally evaluate these modifiers eagerly instead of lazily carrying them around until they need to be resolved (like how we do this for controls).
• the rest, to me, seems like code changes that are more related to the way that angle computations and shift expressions are realised than to the actual representation of hybrid programs. Naturally, it would be interesting to cover the respective translations at some point, but I suppose this is not a short term goal. One could imagine extending the existing OpenQASM3 parser to directly translate to MLIR (similar for Qiskit).
• I believe MLIR itself should also be able to perform quite some optimisations on the classical code; performing things like replacing the floating point pow in my code snippet with an integer bit shift as in your example.

Anything that comes to your mind when looking at the code snippet above?

[burgholzer]

Just adding this here to summarize some of the offline discussion:
Several aspects involved in this have been split out into separate issues

• 🔁 MLIR - Gate modifiers framework #1129 🔁 MLIR - inv (adjoint) modifier (op- and region-level) #1130 🔁 MLIR - pow(r) modifier (integer r first; floats deferred) #1131
• 🧮 MLIR - Decide and document classical result semantics (SSA vs memory via memref) #1123
• 🔁 MLIR - Import classic-controlled operations from QuantumComputation #1109
• 📐 MLIR - Optionally represent angles as fixed-point (quantized) type #1128

Once these have been figured out, we should add a dedicated example representing iterative quantum phase estimation to the documentation and the tests.
I correspondingly reworded the issue description and title.