Qda

classical2quantum(A_c, b_c)

Convert a classical linear system Ax = b to quantum-compatible form

Returns:

Name	Type	Description
A_q	NDArray[float64]	Quantum-compatible matrix (Hermitian, power-of-2 dimension)
b_q	NDArray[float64]	Corresponding right-hand side vector
recover_x	Callable[[ndarray], ndarray]	Function to recover original solution from quantum solution

Source code in qalgo/qda/qram.py

def classical2quantum(
    A_c: np.ndarray | list, b_c: np.ndarray | list
) -> tuple[
    NDArray[np.float64], NDArray[np.float64], Callable[[np.ndarray], np.ndarray]
]:
    """
    Convert a classical linear system Ax = b to quantum-compatible form

    Returns:
        A_q: Quantum-compatible matrix (Hermitian, power-of-2 dimension)
        b_q: Corresponding right-hand side vector
        recover_x: Function to recover original solution from quantum solution
    """
    A_c = np.array(A_c, dtype=np.float64)
    b_c = np.array(b_c, dtype=np.float64)

    if A_c.shape[0] != A_c.shape[1]:
        raise ValueError("Input matrix A_c must be square.")
    if A_c.shape[0] != b_c.size:
        raise ValueError("Dimensions of A_c and b_c are incompatible.")

    original_dim = A_c.shape[0]
    hermitian_transform_done = False

    # Step 1: Hermitization (if necessary)
    # A simple check, though the problem implies it might not be
    # For robustness, we can always apply the embedding if not explicitly told it's Hermitian.
    # Or, if A_c IS Hermitian, we can skip. For now, let's assume we check.
    # A more robust check than `is_hermitian` for very small matrices or specific structures
    # might be needed, but `np.allclose` is generally good.
    if utils.is_hermitian(A_c):
        # print("Input A is already Hermitian.")
        A_herm = A_c.copy()
        b_herm = b_c.copy()
    else:
        print("Input A is not Hermitian. Applying transformation.")
        hermitian_transform_done = True
        n = A_c.shape[0]
        A_herm = np.zeros((2 * n, 2 * n), dtype=A_c.dtype)
        A_herm[:n, n:] = A_c
        A_herm[n:, :n] = A_c.conj().T

        b_herm = np.zeros(2 * n, dtype=A_c.dtype)
        b_herm[:n] = b_c

    # Step 2: Padding to power of 2
    herm_dim = A_herm.shape[0]
    padded_dim = utils.next_power_of_2(herm_dim)

    if padded_dim == herm_dim:
        # print("Dimension is already a power of 2.")
        A_q = A_herm
        b_q = b_herm
    else:
        print(f"Padding dimension from {herm_dim} to {padded_dim}")
        A_q = np.identity(padded_dim, dtype=A_herm.dtype)
        b_q = np.zeros(padded_dim, dtype=b_herm.dtype)

        A_q[:herm_dim, :herm_dim] = A_herm
        b_q[:herm_dim] = b_herm

    # Step 3: Normalize the matrix A_q and vector b_q
    b_q /= np.linalg.norm(b_q)
    A_q /= np.linalg.norm(A_q)

    # Step 4: Create the recovery function
    def recover_x(x_q: np.ndarray) -> np.ndarray:
        if x_q.size != padded_dim:
            print(x_q.size, padded_dim)
            raise RuntimeError(
                "Solution vector x_q has incorrect dimension for recovery."
            )

        x_herm = x_q[:herm_dim].copy()

        if hermitian_transform_done:
            # Original x was in the second half of the [0, x]^T solution vector
            # for the system [0 A; A_dag 0] [y; z] = [b; 0]
            # where solution is y=0, z=x. So x_herm = [0; x]
            # and herm_dime == 2 * original_dim in this case.
            if x_herm.size != 2 * original_dim:
                raise RuntimeError(
                    "Mismatch in dimensions during hermitian recovery logic."
                )
            return x_herm[original_dim:]
        else:
            # No hermitian transform was done, x_herm is directly the solution
            # (potentially padded, but x_herm[:original_dim] already handled that)
            # and herm_dim == original_dim in this case.
            if x_herm.size != original_dim:
                raise RuntimeError(
                    "Mismatch in dimensions during non-hermitian recovery logic."
                )
            return x_herm

    return np.array(A_q, dtype=np.float64), np.array(b_q, dtype=np.float64), recover_x

solve(A, b, kappa=None, p=1.3, step_rate=0.01)

Solves the system of linear equations Ax=b using the Quantum Discrete Adiabatic (QDA) algorithm.

This function serves as a high-level wrapper for the entire algorithmic workflow. It encapsulates the process from converting classical data into quantum states, constructing the necessary QRAM circuits, executing the quantum walk, to finally measuring the result and recovering the classical solution vector.

Parameters:

Name	Type	Description	Default
A	ndarray	The matrix A in the linear system Ax=b.	required
b	ndarray	The vector b in the linear system Ax=b.	required
kappa	Optional[float]	The condition number of matrix A. This value is critical for determining the number of steps in the adiabatic evolution. If set to None, the function will attempt to estimate it automatically via utils.condest(A). Defaults to None.	None
p	float	A key evolution parameter within the Quantum Walk sequence, which influences the construction of the Hamiltonian. Defaults to 1.3.	1.3
step_rate	float	A rate used to compute the total number of adiabatic evolution steps. A smaller value results in more steps, theoretically yielding higher accuracy at the cost of increased computation. Defaults to 0.01.	0.01

Raises:

Type	Description
ValueError	Raised if the dimension of matrix A, after being processed by the internal classical2quantum function, is not a power of 2. This is a common requirement for quantum algorithms operating on qubit-based registers.

Returns:

Type	Description
ndarray	np.ndarray: The calculated solution vector x for the linear system.

Source code in qalgo/qda/qram.py

def solve(
    A: np.ndarray,
    b: np.ndarray,
    kappa: Optional[float] = None,
    p: float = 1.3,
    step_rate: float = 0.01,
) -> np.ndarray:
    """Solves the system of linear equations Ax=b using the Quantum Discrete Adiabatic (QDA) algorithm.

    This function serves as a high-level wrapper for the entire algorithmic workflow. It encapsulates the process from converting classical data into quantum states, constructing the necessary QRAM circuits, executing the quantum walk, to finally measuring the result and recovering the classical solution vector.

    Args:
        A (np.ndarray): The matrix A in the linear system Ax=b.
        b (np.ndarray): The vector b in the linear system Ax=b.
        kappa (Optional[float], optional): The condition number of matrix A. This value is critical for determining the number of steps in the adiabatic evolution. If set to None, the function will attempt to estimate it automatically via `utils.condest(A)`. Defaults to None.
        p (float, optional): A key evolution parameter within the Quantum Walk sequence, which influences the construction of the Hamiltonian. Defaults to 1.3.
        step_rate (float, optional): A rate used to compute the total number of adiabatic evolution steps. A smaller value results in more steps, theoretically yielding higher accuracy at the cost of increased computation. Defaults to 0.01.

    Raises:
        ValueError: Raised if the dimension of matrix A, after being processed by the internal `classical2quantum` function, is not a power of 2. This is a common requirement for quantum algorithms operating on qubit-based registers.

    Returns:
        np.ndarray: The calculated solution vector x for the linear system.
    """
    A = np.array(A, dtype=np.float64)
    b = np.array(b, dtype=np.float64)

    if kappa is None:
        kappa = utils.condest(A)
    print(f"{kappa = }")

    steps = compute_step_rate(step_rate, kappa)
    print(f"{steps = }")

    A, b, recover_x = classical2quantum(A, b)

    data_size = 50
    rational_size = 51

    exponent = 15
    log_column_size = int(np.ceil(np.log2(A.shape[0])))
    if A.shape[0] != 2**log_column_size:
        raise ValueError(
            "Matrix dimension is not a power of 2. Call 'classical2quantum' first."
        )

    conv_A = scale_and_convert_vector(A.flatten(order="F"), exponent, data_size)
    conv_b = scale_and_convert_vector(b, exponent, data_size)
    data_tree_A = make_vector_tree(conv_A, data_size)
    data_tree_b = make_vector_tree(conv_b, data_size)
    addr_size = log_column_size * 2 + 1

    qram_A = sq.QRAMCircuit_qutrit(addr_size, data_size, data_tree_A)
    qram_b = sq.QRAMCircuit_qutrit(log_column_size + 1, data_size, data_tree_b)
    print("QRAMCircuit ok")

    state = sq.SparseState()

    main_reg = sq.AddRegister(
        "main_reg", sq.StateStorageType.UnsignedInteger, log_column_size
    )(state)
    anc_UA = sq.AddRegister(
        "anc_UA", sq.StateStorageType.UnsignedInteger, log_column_size
    )(state)
    anc_4 = sq.AddRegister("anc_4", sq.StateStorageType.Boolean, 1)(state)
    anc_3 = sq.AddRegister("anc_3", sq.StateStorageType.Boolean, 1)(state)
    anc_2 = sq.AddRegister("anc_2", sq.StateStorageType.Boolean, 1)(state)
    anc_1 = sq.AddRegister("anc_1", sq.StateStorageType.Boolean, 1)(state)

    sq.State_Prep_via_QRAM(qram_b, "main_reg", data_size, rational_size)(state)
    WalkSequence_via_QRAM_Debug(
        qram_A,
        qram_b,
        A,
        b,
        "main_reg",
        "anc_UA",
        "anc_1",
        "anc_2",
        "anc_3",
        "anc_4",
        steps,
        kappa,
        p,
        data_size,
        rational_size,
    )(state)

    # Calculate the total probability of the subspace where anc_UA, anc_2, anc_3 are 0
    prob_inv0 = sq.PartialTraceSelect({anc_UA: 0, anc_2: 0, anc_3: 0})(state)
    prob0 = (1.0 / prob_inv0) ** 2

    print("Success probability after walk sequence:", prob0)

    sol, _ = sq.PartialTraceSelect(
        {anc_UA: 0, anc_1: 1, anc_2: 0, anc_3: 0, anc_4: 0}
    ).get_projected_full(state)
    sol = np.array(sol, np.complex128)
    sol = recover_x(sol.real)

    return sol