algo/qaoa.py

class QUBO

QUBO (Quadratic Unconstrained Binary Optimization) problems defined by a matrix Q to minimize the function:

$f(x)=\sum _{i}Q[i,i]\times x[i]+\sum _{i<j}Q[i,j]\times x[i]\times x[j]$

where x[i] are binary variables (0 or 1). The toHamiltonian method converts this QUBO representation into a Hamiltonian suitable for quantum algorithms like QAOA.

toHamiltonian

Convert a QUBO problem defined by matrix Q into a Hamiltonian representation.

Function iterates over Q constructs Hamiltonian in terms of Pauli Z ops. Linear terms are (where $i=j$) and quadratic terms are (where $i\ne j$). The resulting Hamiltonian is list[tuple], where

def toHamiltonian(Q: dict) -> Any [static]

Implementation

def toHamiltonian(Q: dict):
    ising = {}
    offset = 0.0
    for (i, j), val in Q.items():
        if i == j:
            ising[i] = ising.get(i, 0) - val / 2
            offset += val / 2
        else:
            key = tuple(sorted((i, j)))
            ising[key] = ising.get(key, 0) + val / 4
            ising[i] = ising.get(i, 0) - val / 4
            ising[j] = ising.get(j, 0) - val / 4
            offset += val / 4
    ham = []
    for keys, coeff in ising.items():
        if coeff == 0:
            continue
        if isinstance(keys, int):
            ham.append((coeff, 'Z', [keys]))
        else:
            ham.append((coeff, 'ZZ', list(keys)))
    return (ham, offset)

class QAOA

Quantum Approximate Optimization Algorithm (QAOA) implementation for qudits.

QAOA has a known structure, and therefore instead of being built on circuit, inherits directly from nn.Module. The forward method constructs the QAOA state based on the provided Hamiltonian, and the expectation method computes the expectation value of the Hamiltonian with respect to the current state.

name	type	default
d	int	None
device	str	None
hamiltonian	list	None
offset	float	0.0
qubo	dict	None
wires	int	None

init

No Definition provided

def __init__(self: Any, d: int, wires: int, qubo: dict, hamiltonian: list, offset: float, layers: Any, device: Any) -> Any

Implementation

def __init__(self, d: int, wires: int, qubo: dict=None, hamiltonian: list=None, offset: float=0.0, layers=1, device='cpu'):
    super().__init__()
    self.width = d ** wires
    self.wires = wires
    self.d = d
    self.device = device
    self.layers = layers
    self.qubo = qubo
    if hamiltonian is not None:
        self.hamiltonian = hamiltonian
        self.offset = offset
    elif qubo is not None:
        self.hamiltonian, self.offset = QUBO.toHamiltonian(qubo)
    else:
        raise ValueError("Either 'hamiltonian' or 'qubo' must be provided.")
    self._H_P_matrix = self.getMat()
    self.gammas = nn.Parameter(pt.rand(layers, device=device) * (2 * pt.pi))
    self.betas = nn.Parameter(pt.rand(layers, device=device) * pt.pi)
    self.gg = GG.Gategen(dim=d, device=device)
    self._H_factory = self.gg.U(self.gg.H, name='H')
    self._CX_factory = self.gg.U(self.gg.CX, name='CX')
    self.OpMap = {'Z': self._RZ, 'ZZ': self._RZZ}

gate

Helper function to apply a gate to the state x using the specified gate_class and index. Additional parameters for the gate can be passed via kwargs.

def gate(self: Any, x: Any, gate_class: Any, index: Any) -> Any

Implementation

def gate(self, x, gate_class, index, **kwargs):
    gate = gate_class(dim=self.d, wires=self.wires, index=index, device=self.device, **kwargs)
    return gate.forward(x)

forward

No Definition provided

def forward(self: Any) -> Any

Implementation

def forward(self):
    state = pt.zeros((self.width, 1), dtype=C64, device=self.device)
    state[0, 0] = 1.0
    for j in range(self.wires):
        state = self.gate(state, self._H_factory, index=[j])
    for i in range(self.layers):
        for coeff, gtype, indices in self.hamiltonian:
            angle = 2 * self.gammas[i] * coeff
            if gtype in self.OpMap:
                state = self.OpMap[gtype](state, angle, indices)
            else:
                raise NotImplementedError(f"Gate type '{gtype}' not supported.")
        for j in range(self.wires):
            state = self.gate(state, self.gg.RX, index=[j], angle=2 * self.betas[i])
    return state

getMat

Construct Hamiltonian matrix from self.hamiltonian.

def getMat(self: Any) -> Any

Implementation

def getMat(self):
    H_P = pt.zeros((self.width, self.width), dtype=C64, device=self.device)
    gg = GG.Gategen(dim=self.d, device=self.device)
    opmat = {'I': gg.I.tensor, 'Z': gg.Z.tensor, 'X': gg.X.tensor}
    for coeff, gtype, indices in self.hamiltonian:
        oplist = []
        if gtype == 'Z':
            oplist = [opmat['Z'] if i == indices[0] else opmat['I'] for i in range(self.wires)]
        elif gtype == 'ZZ':
            oplist = [opmat['Z'] if i in indices else opmat['I'] for i in range(self.wires)]
        else:
            raise NotImplementedError(f"Matrix for '{gtype}' not defined.")
        term = oplist[0]
        for k in range(1, len(oplist)):
            term = pt.kron(term, oplist[k])
        H_P += coeff * term
    return H_P

expectation

Compute the expectation value of the Hamiltonian with respect to the current state.

ExpVal $⟨\psi |H|\psi ⟩$, where $|\psi ⟩$ is the state obtained from the forward method, and H is the Hamiltonian matrix constructed in getMat. The offset is added to the computed expectation value to account for any constant terms in the Hamiltonian.

def expectation(self: Any) -> Any

Implementation

def expectation(self):
    final_state = self.forward()
    exp_val = pt.vdot(final_state.squeeze(), (self._H_P_matrix @ final_state).squeeze()).real
    return exp_val + self.offset

solve

No Definition provided

def solve(self: Any, func: Energy, optimizer: Any, steps: Any, lr: Any, hook: Any) -> Any

Implementation

def solve(self, func: Energy=None, optimizer=None, steps=100, lr=0.1, hook=devnull):
    if optimizer is None:
        optimizer = pt.optim.Adam(self.parameters(), lr=lr)
    for step in range(steps):
        optimizer.zero_grad()
        loss = self.expectation()
        loss.backward()
        optimizer.step()
        hook(loss, step)
    with pt.no_grad():
        final_state = self.forward()
        Pi = (pt.abs(final_state) ** 2).squeeze()
        maxP = pt.argmax(Pi).item()
        solution = format(maxP, f'0{self.wires}b')
        solution = [int(bit) for bit in solution]
    out = {'solution': solution, 'probabilities': Pi.cpu().flatten()}
    if func is not None and self.qubo is not None:
        soltensr = [pt.tensor(bit) for bit in solution]
        out['value'] = func(self.qubo, soltensr)
    return out

class ClockSolver

Variational qudit optimizer using direct matrix exponentiation (clock/Potts model).

Works for any local dimension $d\ge 2$. Each layer applies a phase separator $U_P(\gamma )=\exp (-i\gamma H_P)$ and a mixer $U_B(\beta )=\exp (-i\beta H_B)$ where $H_B=-\sum _i(X_d^{(i)}+X_d^{(i)†})$.

Initial state: $H_d|0{⟩}^{\otimes n}$. Parameters $\gamma ,\beta$ are optimized with Adam.

name	type	default
d	int	None
device	str	None
hamiltonian	list	None
layers	int	None
offset	float	None
qubo	dict	None
width	int	None
wires	int	None

init

No Definition provided

def __init__(self: Any, d: int, wires: int, qubo: dict, hamiltonian: list, offset: float, layers: int, device: str) -> Any

Implementation

def __init__(self, d: int, wires: int, qubo: dict=None, hamiltonian: list=None, offset: float=0.0, layers: int=1, device: str='cpu'):
    super().__init__()
    self.d = d
    self.wires = wires
    self.width = d ** wires
    self.layers = layers
    self.device = device
    self.qubo = qubo
    if hamiltonian is not None:
        self.hamiltonian = hamiltonian
        self.offset = offset
    elif qubo is not None:
        self.hamiltonian, self.offset = QUBO.toHamiltonian(qubo)
    else:
        raise ValueError("Either 'hamiltonian' or 'qubo' must be provided.")
    self.gammas = nn.Parameter(pt.rand(layers, device=device) * (2 * pt.pi))
    self.betas = nn.Parameter(pt.rand(layers, device=device) * pt.pi)
    with pt.no_grad():
        self._H_P = self._buildHP()
        self._H_B = self._buildHB()

forward

Return the QAOA state $|\psi (\gamma ,\beta )⟩$ as a complex vector of length $d^n$.

def forward(self: Any) -> pt.Tensor

Implementation

def forward(self) -> pt.Tensor:
    gg = GG.Gategen(dim=self.d, device=self.device)
    ket0 = pt.zeros(self.d, dtype=C64, device=self.device)
    ket0[0] = 1.0
    plus = gg.H.tensor @ ket0
    state = plus
    for _ in range(self.wires - 1):
        state = pt.kron(state, plus)
    for i in range(self.layers):
        U_P = pt.linalg.matrix_exp(-1j * self.gammas[i].to(C64) * self._H_P)
        state = U_P @ state
        U_B = pt.linalg.matrix_exp(-1j * self.betas[i].to(C64) * self._H_B)
        state = U_B @ state
    return state

expectation

Expectation value $⟨\psi |H_P|\psi ⟩+\mathrm{offset}$.

def expectation(self: Any) -> pt.Tensor

Implementation

def expectation(self) -> pt.Tensor:
    state = self.forward()
    exp_val = pt.vdot(state, self._H_P @ state).real
    return exp_val + self.offset

solve

Optimise gammas/betas using PyTorch Adam (or a provided optimizer).

def solve(self: Any, func: 'Energy | None', optimizer: Any, steps: int, lr: float, hook: 'Callable') -> dict

Implementation

def solve(self, func: 'Energy | None'=None, optimizer=None, steps: int=200, lr: float=0.1, hook: 'Callable'=devnull) -> dict:
    if optimizer is None:
        optimizer = pt.optim.Adam(self.parameters(), lr=lr)
    for step in range(steps):
        optimizer.zero_grad()
        loss = self.expectation()
        loss.backward()
        optimizer.step()
        hook(loss, step)
    with pt.no_grad():
        final_state = self.forward()
        Pi = pt.abs(final_state) ** 2
        maxP = int(pt.argmax(Pi).item())
        solution = self._decode(maxP)
    out: dict = {'solution': solution, 'probabilities': Pi.cpu()}
    if func is not None and self.qubo is not None:
        out['value'] = func(self.qubo, [pt.tensor(x) for x in solution])
    return out