lcu

name: lcu description: A quantum algorithm for solving linear systems of equations using the Linear Combination of Unitaries (LCU) method, providing exponential speedup over classical methods for certain types of problems. This skill includes efficient implementations and educational resources for understanding and utilizing the LCU algorithm in various applications.

Linear Combination of Unitaries (LCU)

Purpose

LCU probabilistically applies a non-unitary operator $M = \sum_j \alpha_j U_j$ to a quantum state by using an ancilla register, a SELECT oracle, and amplitude amplification post-selection. It is a key subroutine in QSVT, Hamiltonian simulation, and quantum linear solvers.

Use this skill when you need to:

• Apply a weighted sum of unitaries to a quantum state.
• Implement block-encodings for non-unitary operations.

Overview

The LCU circuit has three components:

1. $V$ (PREPARE): Creates the superposition $V|0\rangle = \sum_j \sqrt{\alpha_j / s}|j\rangle$ in the ancilla register, where $s = \sum_j \alpha_j$.
2. SELECT ($U_c$): Controlled on ancilla $= |j\rangle$, applies $U_j$ to the system: $\sum_j |j\rangle|j\rangle_{\text{anc}} \otimes U_j|\psi\rangle$.
3. $V^\dagger$ (UNPREPARE): Rotates the ancilla back so that measuring all-zeros succeeds with probability $|M|\psi\rangle|^2/s^2$.

Prerequisites

• Controlled unitary gates.
• Ancilla qubits and post-selection.
• Python: numpy, Circuit.

Using the Provided Implementation

from unitarylab_algorithms import LCUAlgorithm
from unitarylab.core import Circuit
import numpy as np

# Build two unitary operators U0, U1
n_sys = 1
U0 = Circuit(n_sys, name='H')
U0.h(0)         # Hadamard

U1 = Circuit(n_sys, name='X')
U1.x(0)         # Pauli-X

# Apply M = 0.6*H + 0.4*X to system
algo = LCUAlgorithm(text_mode='plain')
result = algo.run(
    alphas=[0.6, 0.4],
    unitaries=[U0, U1],
    n_sys=n_sys,
    initial_state=None,    # Starts in |0>
    backend='torch'
)

print(result['Success probability'])   # Post-selection probability
print(result['status'])                # 'ok' on success
print(result['circuit_path'])          # SVG circuit diagram
print(result['Result state'])          # Post-selected system state

Core Parameters Explained

Common misunderstandings:

• alphas must be non-negative (used as probability amplitudes after square root).
• The ancilla register size is $n_{\text{anc}} = \lceil\log_2 m\rceil$ where $m = $ len(alphas).
• LCU is a probabilistic procedure; the success probability can be low for ill-conditioned operators.

Return Fields

Implementation Architecture

LCUAlgorithm in algorithm.py builds the LCU circuit in five stages using two circuit constructor helpers: _build_V and _build_select.

run(alphas, unitaries, n_sys, initial_state, backend, device, dtype) — Five Stages:

Helper Methods:

• _build_V(alphas, s_norm, m, n_anc) — Constructs the PREPARE operator. Builds amplitude vector state[j] = sqrt(alphas[j]/s_norm) padded to 2^n_anc elements. Creates Circuit(n_anc, name='V') and calls qc.initialize(state, range(n_anc)) to directly load the amplitude distribution.
• _build_select(unitaries, m, n_anc, n_sys) — Constructs the SELECT operator. Creates Circuit(n_anc + n_sys, name='SELECT-U'). Loops over all (j, U_j) pairs; for each, computes ctrl_state = format(j, f"0{n_anc}b") (reversed if U.order == 'little') and calls qc.append(U_j, target=sys_qubits, control=anc_qubits, control_state=ctrl_state). This applies $U_j$ conditionally on ancilla $= |j\rangle$ in binary.
• _build_return_dict(success, circuit_path, filepath, circuit) — Packages the result. Sets status='ok' if success=True, 'failed' otherwise. Wraps filepath into a list of {'format': filename[-3:], 'filename': filename} dicts stored under 'plot'. Merges self.output (containing 'Success probability', 'Computation time (s)', 'Result state') into the returned dict.

Data flow: (alphas, unitaries) → _build_V() → _build_select() → V_circ.dagger() → qc.execute() → _phase_probabilities_from_state() → success_prob → update_output() → _build_return_dict().

Understanding the Key Quantum Components

Maps $|0\rangle_{\text{anc}} \rightarrow \sum_j \sqrt{\alpha_j/s}|j\rangle_{\text{anc}}$ where $s = \sum_k \alpha_k$. Implemented as a state-preparation circuit on the ancilla.

2. SELECT Operator ($U_c$)

$$U_c = \sum_{j=0}^{m-1} |j\rangle\langle j|_{\text{anc}} \otimes U_j$$ Controlled on each basis state of the ancilla, applies the corresponding $U_j$ to the system. Implemented as a multiplexer (multi-controlled gate).

3. UNPREPARE Operator ($V^\dagger$)

The inverse of $V$. After combining PREPARE → SELECT → UNPREPARE, the state decomposes as: $$|\Phi\rangle = |0\rangle_{\text{anc}} \otimes \frac{1}{s}M|\psi\rangle + |0\rangle_{\text{anc}}^{\perp} \otimes \ldots$$

4. Post-selection

Measuring the ancilla in $|0\rangle^{n_{\text{anc}}}$ succeeds with probability: $$P_{\text{succ}} = \frac{|M|\psi\rangle|^2}{s^2}$$ On success, the system register holds $M|\psi\rangle/|M|\psi\rangle|$.

Theory-to-Code Mapping

Notes on encapsulation: The PREPARE operator uses qc.initialize(state, range(n_anc)) to directly load the amplitude distribution — no recursive gate decomposition in _build_V itself. The SELECT operator directly calls qc.append() with control_state and reverses ctrl_state when U.order == 'little'. _build_return_dict always receives success=True (hardcoded in run()), so status is always 'ok'; the actual LCU outcome quality is in result['Success probability'].

Mathematical Deep Dive

$$V|0\rangle_{\text{anc}}|\psi\rangle = \sum_j \sqrt{\frac{\alpha_j}{s}}|j\rangle|\psi\rangle$$

After SELECT: $$\sum_j \sqrt{\frac{\alpha_j}{s}}|j\rangle U_j|\psi\rangle$$

After UNPREPARE ($V^\dagger \otimes I$): $$|0\rangle_{\text{anc}} \otimes \frac{1}{s}M|\psi\rangle + |0\rangle_{\text{anc}}^{\perp} \otimes \ldots$$

Post-selection success probability: $$P = \frac{|M|\psi\rangle|^2}{s^2}$$

To boost $P$, apply amplitude amplification $O(s/|M|\psi\rangle|)$ times.

Hands-On Example

from unitarylab_algorithms import LCUAlgorithm
from unitarylab.core import Circuit

n_sys = 2
I_circ = Circuit(n_sys, name='II')    # Identity
X0_circ = Circuit(n_sys, name='XI')
X0_circ.x(0)

algo = LCUAlgorithm(text_mode='plain')
result = algo.run(
    alphas=[0.8, 0.2],
    unitaries=[I_circ, X0_circ],
    n_sys=n_sys,
    backend='torch'
)
print(f"Success probability: {result['Success probability']:.4f}")
print(f"Status: {result['status']}")
print(f"Result state: {result['Result state']}")

Implementing Your Own Version

The following Python skeleton reconstructs the three structural components of LCU: PREPARE, SELECT, and UNPREPARE.

Step 1: Amplitude-tree PREPARE circuit

# Simplified reconstruction — mirrors LCUAlgorithm._build_V() and _state_preparation_tree()
import numpy as np
from unitarylab.core import Circuit

def build_prepare(alphas, n_anc: int, backend: str = 'torch') -> Circuit:
    """Build V|0> = sum_j sqrt(alpha_j/s)|j> on ancilla register."""
    s = sum(alphas)
    m = len(alphas)
    # Pad to 2^n_anc
    state = np.zeros(2**n_anc)
    for j in range(m):
        state[j] = np.sqrt(alphas[j] / s)
    # Use tree decomposition: at each node compute Ry angle theta = 2*arctan2(||right||, ||left||)
    # The real implementation calls _state_preparation_tree() which does this recursively
    # For a flat 2-element case:
    if m == 2:
        qc = Circuit(n_anc)
        theta = 2 * np.arctan2(state[1], state[0])
        qc.ry(theta, 0)
        return qc
    # General case: delegate to the recursive tree builder in LCUAlgorithm
    raise NotImplementedError("Use LCUAlgorithm._build_V() for general m")

Step 2: SELECT oracle (multiplexer)

# Simplified reconstruction — mirrors LCUAlgorithm._build_select()

def build_select(unitaries, n_anc: int, n_sys: int, backend: str = 'torch') -> Circuit:
    """SELECT = sum_j |j><j|_anc ⊗ U_j: apply U_j conditional on ancilla=|j>."""
    m = len(unitaries)
    qc = Circuit(n_anc + n_sys)
    sys_qubits = list(range(n_anc, n_anc + n_sys))
    anc_qubits = list(range(n_anc))
    for j, U_j in enumerate(unitaries):
        control_seq = format(j, f'0{n_anc}b')   # e.g. '01' for j=1, n_anc=2
        qc.append(U_j, target=sys_qubits, control=anc_qubits, control_state=control_seq)
    return qc

Step 3: Full LCU circuit assembly

# Simplified reconstruction of the full PREPARE → SELECT → UNPREPARE structure

def lcu_circuit(alphas, unitaries, n_sys: int, initial_state=None, backend: str = 'torch'):
    """Full LCU circuit: PREPARE ⊗ I |psi> → SELECT → UNPREPARE ⊗ I."""
    m = len(alphas)
    n_anc = int(np.ceil(np.log2(max(m, 2))))
    qc = Circuit(n_anc + n_sys)

    anc_qubits = list(range(n_anc))
    sys_qubits = list(range(n_anc, n_anc + n_sys))

    # 1. Initial system state
    if initial_state is not None:
        qc.append(initial_state, sys_qubits)

    # 2. PREPARE
    V_circ = build_prepare(alphas, n_anc, backend)
    qc.append(V_circ, anc_qubits)

    # 3. SELECT
    sel_circ = build_select(unitaries, n_anc, n_sys, backend)
    qc.append(sel_circ, list(range(n_anc + n_sys)))

    # 4. UNPREPARE (V†)
    qc.append(V_circ.dagger(), anc_qubits)

    return qc

Debugging Tips

1. len(alphas) != len(unitaries): Raises ValueError. Always keep these lists the same length.
2. initial_state.get_num_qubits() != n_sys: Raises ValueError. The initial state circuit must match n_sys exactly.
3. Very low success_prob: This is expected when $|M|\psi\rangle| \ll s$. Use amplitude amplification to boost.
4. Very low 'Success probability' (near 0): status is always 'ok' since _build_return_dict is always called with success=True. Check result['Success probability'] directly. If near 0, the operator $M$ may annihilate $|\psi\rangle$, or the unitaries cancel.
5. Memory limit: Ancilla adds $\lceil\log_2 m\rceil$ qubits. For $m=8$ unitaries, this is only 3 ancilla qubits.