HHL算法及其改进

Origin HHL

在 HHL 算法中，我们需要计算一个参数 $\delta$，该参数用于确定酉矩阵 $U$ 中旋转角度的大小，进而影响量子相位估计过程中测量结果的精度。

具体而言，在 HHL 算法中，我们需要构造一个酉矩阵 $U=e^{iAt}$，其中 $A$ 是待求解的厄米矩阵，$t$ 是待确定的时间参数。通过量子相位估计技术，我们可以将 $t$ 的取值范围缩小到 $[0,2\pi)$，并且使得相位估计结果精度达到 $\mathcal{O}\left(\frac{1}{\delta}\right)$。因此，我们希望选择一个合适的 $\delta$ 值，使得相位估计结果的精度可以满足我们的要求。

具体而言，HHL 算法中 $\delta$ 的计算公式为：

$$\delta=\frac{\text{polylog}(n)}{\sqrt{\lambda_{\min}\lambda_{\max}}}$$

其中 $n$ 是矩阵 $A$ 的维数，$\lambda_{\text{min}}$ 和 $\lambda_{\text{max}}$ 分别是 $A$ 的最小和最大特征值，$\text{polylog}(n)$ 表示多重对数函数。这个公式中，分母部分可以看做是矩阵的谱范数，也就是矩阵所有特征值绝对值的最大值。因此，$\delta$ 的值越小，矩阵的谱范数越大，量子相位估计的精度也就越高。

在代码中，计算 $\delta$ 的过程包含了以下几个步骤：

1. 将 $\lambda_{\text{min}}$ 缩放到 $[0,2^{n_l}-1]$ 的区间内，这里的 $n_l$ 是指用于表示二进制小数时使用的比特数。
2. 将 $\lambda_{\text{min}}$ 的二进制表示转换为十进制表示，并将小数点前的所有位数相加，得到 $\delta$ 的近似值。

需要注意的是，这里计算得到的是 $\delta$ 的近似值，而非精确值。因此，在实际使用中，可能需要根据具体的需求和精度要求，对 $\delta$ 进行进一步的调整和优化。

Qiskit Linear Solver源码解析

matrices

LinearSystemMatrix

线性问题矩阵的基类。

def __init__(
    self,
    num_state_qubits: int,
    tolerance: float,
    evolution_time: float,
    name: str = "ls_matrix",
) -> None:
    super().__init__(name=name)

    # 存储参数
    self._num_state_qubits = num_state_qubits
    self._reset_registers(num_state_qubits)
    self.tolerance = tolerance
    self.evolution_time = evolution_time

参数	含义
num_state_qubits	$|x\rangle$比特数
tolerance	精度
evolution_time	哈密顿模拟时间
name	名称

@abstractmethod
def eigs_bounds(self) -> Tuple[float, float]:
    """返回矩阵最小和最大特征值。"""
    raise NotImplementedError

@abstractmethod
def condition_bounds(self) -> Tuple[float, float]:
    """返回矩阵最小和最大条件数。"""
    raise NotImplementedError

@abstractmethod
def _reset_registers(self, num_state_qubits: int) -> None:
    """依据新的比特数重置量子寄存器。

    参数:
        num_state_qubits: 新的比特数。
    """
    raise NotImplementedError

@abstractmethod
def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit:
    """构建电路的“幂”。

    Args:
        power: 幂。
        matrix_power: 如果为True，电路会被转为矩阵并计算矩阵的幂。否则，若power为正，会重复相应次数的电路。           the implementation defaults to ``repeat``.

    Returns:
        酉幂的量子电路。
    """
    raise NotImplementedError

NumpyMatrix

接收一个NumPy数组进行初始化。

检查矩阵合法性。

def _check_configuration(self, raise_on_failure: bool = True) -> bool:
    valid = True

    if self.matrix.shape[0] != self.matrix.shape[1]:
        if raise_on_failure:
            raise AttributeError("Input matrix must be square!")
        return False
    if np.log2(self.matrix.shape[0]) % 1 != 0:
        if raise_on_failure:
            raise AttributeError("Input matrix dimension must be 2^n!")
        return False
    if not np.allclose(self.matrix, self.matrix.conj().T):
        if raise_on_failure:
            raise AttributeError("Input matrix must be hermitian!")
        return False

    return valid

重复电路，用 $U_A=e^{iAt}$ 生成酉变换。

def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit:
    """构建电路的“幂”。

    Args:
        power: 幂。
        matrix_power: 如果为True，电路会被转为矩阵并计算矩阵的幂。否则，若power为正，会重复相应次数的电路。           the implementation defaults to ``repeat``.

    Returns:
        酉幂的量子电路。
    """
    qc = QuantumCircuit(self.num_state_qubits)
    evolved = sp.linalg.expm(1j * self.matrix * self.evolution_time)
    # 构造酉电路
    qc.unitary(evolved, qc.qubits)
    return qc.power(power)

Solver

LinearSolver

仅包含一个solve方法。

class LinearSolver(ABC):
    """线性求解器的抽象类。"""

    @abstractmethod
    def solve(
            self,
            matrix: Union[np.ndarray, QuantumCircuit],
            vector: Union[np.ndarray, QuantumCircuit],
            observable: Optional[
                Union[
                    LinearSystemObservable,
                    List[LinearSystemObservable],
                ]
            ] = None,
            observable_circuit: Optional[
                Union[QuantumCircuit, List[QuantumCircuit]]
            ] = None,
            post_processing: Optional[
                Callable[[Union[float, List[float]], int, float], float]
            ] = None,
    ) -> LinearSolverResult:

参数	含义
matrix	$A\vec{x}=\vec{b}$ 中的 $A$。
vector	$A\vec{x}=\vec{b}$ 中的 $b$。
observable	观测量，默认为HHL成功的概率。
observable_circuit	用于测量的电路，默认为None。
post_processing	对观测量的进一步处理，默认输出观测量原始值。

NumpyLinearSolver

这是一个经典求解器，能帮助我们理解HHL求解器的代码逻辑。

def solve(
    self,
    matrix: Union[np.ndarray, QuantumCircuit],
    vector: Union[np.ndarray, QuantumCircuit],
    observable: Optional[
        Union[
            LinearSystemObservable,
            List[LinearSystemObservable],
        ]
    ] = None,
    observable_circuit: Optional[
        Union[QuantumCircuit, List[QuantumCircuit]]
    ] = None,
    post_processing: Optional[
        Callable[[Union[float, List[float]], int, float], float]
    ] = None,
) -> LinearSolverResult:

    # 检查实例所属类（获取A和x）。
    if isinstance(vector, QuantumCircuit):
        vector = Statevector(vector).data
    if isinstance(matrix, QuantumCircuit):
        if hasattr(matrix, "matrix"):
            matrix = matrix.matrix
        else:
            matrix = Operator(matrix).data

    # 求解
    solution_vector = np.linalg.solve(matrix, vector)
    solution = LinearSolverResult()
    solution.state = solution_vector

    # 计算观测量
    if observable is not None:
        if isinstance(observable, list):
            solution.observable = []
            for obs in observable:
                solution.observable.append(
                    # 调用经典计算的测量
                    obs.evaluate_classically(solution_vector)
                )
        else:
            # 调用经典计算的测量
            solution.observable = observable.evaluate_classically(solution_vector)
    # 计算二范数
    solution.euclidean_norm = float(np.linalg.norm(solution_vector))
    return solution

HHL

利用HHL求解。

def __init__(
    self,
    epsilon: float = 1e-2,
    expectation: Optional[ExpectationBase] = None,
    quantum_instance: Optional[Union[Backend, QuantumInstance]] = None,
) -> None:

    super().__init__()

    self._epsilon = epsilon
    # Tolerance for the different parts of the algorithm as per [1]
    self._epsilon_r = epsilon / 3  # conditioned rotation
    self._epsilon_s = epsilon / 3  # state preparation
    self._epsilon_a = epsilon / 6  # hamiltonian simulation

    self._scaling = None  # scaling of the solution

    self._sampler = None
    self.quantum_instance = quantum_instance

    self._expectation = expectation

    # For now the default reciprocal implementation is exact
    self._exact_reciprocal = True
    # Set the default scaling to 1
    self.scaling = 1

• epsilon：要求精确解 $\vec{x}$ 与近似解 $\tilde{x}$ 之间的误差 $\|x-\tilde{x}\|\le \epsilon$。

• expectation：应用于期望值的期望转换器。如果为None，则使用PauliExpectation。

• QuantumInstance：量子计算环境的实例或后端，如果为None，就用Statevector。

可以看到HHL的solve函数逻辑和NumpyLinearSolver类似，都是给出解（NumpyLinearSolver是解向量，HHL是解电路），接着计算观测量。

构造电路

第一步，将向量 $\vec{b}$ 编码为 $|b\rangle$，注意这里 $\vec{b}$ 被标准化。我们直接调用iso或isometry函数

$$|0\rangle_{n_b}\to\sum_{i=0}^{2^{n_b}-1}b_i|i\rangle$$

第二步，构造矩阵 $A$ 的电路，已经有类NumpyMatrix；接着计算了受控旋转的参数 $C$，即代码中的 delta，还有估计特征值所需的量子比特数 $n_l$，并更新了哈密顿模拟的时间evolution_time和缩放比scaling。

第三步，构造受控旋转电路，已经有类ExactReciprocal，参数 $C$ 已经计算得到。

第四步，构造相位估计电路，已经有类PhaseEstimation。

第五步，拼接电路。

先构造一个空电路。

qc = QuantumCircuit(qb, ql, qf)

  q4: 

q5_0: 

q5_1: 

q5_2: 

  q6: 

拼上vector_circuit后，在 $q_4$ 上有一个模块。

qc.append(vector_circuit, qb[:])

      ┌──────────┐
  q4: ┤ Isometry ├
      └──────────┘
q5_0: ────────────

q5_1: ────────────

q5_2: ────────────

  q6: ────────────

接着拼上QPE电路。

qc.append(phase_estimation, ql[:] + qb[:])

      ┌──────────┐┌──────┐
  q4: ┤ Isometry ├┤3     ├
      └──────────┘│      │
q5_0: ────────────┤0     ├
                  │  QPE │
q5_1: ────────────┤1     ├
                  │      │
q5_2: ────────────┤2     ├
                  └──────┘
  q6: ────────────────────

然后是受控旋转电路。

$$|x\rangle |0\rangle \mapsto \cos\frac{s}{x}|x\rangle|0\rangle + \sin\frac{s}{x}|x\rangle |1\rangle$$

qc.append(reciprocal_circuit, ql[::-1] + [qf[0]])

      ┌──────────┐┌──────┐        
  q4: ┤ Isometry ├┤3     ├────────
      └──────────┘│      │┌──────┐
q5_0: ────────────┤0     ├┤2     ├
                  │  QPE ││      │
q5_1: ────────────┤1     ├┤1     ├
                  │      ││  1/x │
q5_2: ────────────┤2     ├┤0     ├
                  └──────┘│      │
  q6: ────────────────────┤3     ├
                          └──────┘

最后是 $\text{QPE}^{\dagger}$。

qc.append(phase_estimation.inverse(), ql[:] + qb[:])

      ┌──────────┐┌──────┐        ┌─────────┐
  q4: ┤ Isometry ├┤3     ├────────┤3        ├
      └──────────┘│      │┌──────┐│         │
q5_0: ────────────┤0     ├┤2     ├┤0        ├
                  │  QPE ││      ││  QPE_dg │
q5_1: ────────────┤1     ├┤1     ├┤1        ├
                  │      ││  1/x ││         │
q5_2: ────────────┤2     ├┤0     ├┤2        ├
                  └──────┘│      │└─────────┘
  q6: ────────────────────┤3     ├───────────
                          └──────┘           

observables

observables是关于测量的组件。

LinearSystemObservable

这是线性观测的抽象类。

class LinearSystemObservable(ABC):

    @abstractmethod
    def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]:
        """The observable operator.

        Args:
            num_qubits: The number of qubits on which the observable will be applied.

        Returns:
            The observable as a sum of Pauli strings.
        """
        raise NotImplementedError

    @abstractmethod
    def observable_circuit(
        self, num_qubits: int
    ) -> Union[QuantumCircuit, List[QuantumCircuit]]:
        """The circuit implementing the observable.

        Args:
            num_qubits: The number of qubits on which the observable will be applied.

        Returns:
            The observable as a QuantumCircuit.
        """
        raise NotImplementedError

    @abstractmethod
    def post_processing(
        self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1
    ) -> float:
        """Evaluates the given observable on the solution to the linear system.

        Args:
            solution: The probability calculated from the circuit and the observable.
            num_qubits: The number of qubits where the observable was applied.
            scaling: Scaling of the solution.

        Returns:
            The value of the observable.
        """
        raise NotImplementedError

    @abstractmethod
    def evaluate_classically(
        self, solution: Union[np.ndarray, QuantumCircuit]
    ) -> float:
        """Calculates the analytical value of the given observable from the solution vector to the
         linear system.

        Args:
            solution: The solution to the system as a numpy array or the circuit that prepares it.

        Returns:
            The value of the observable.
        """
        raise NotImplementedError

Hybrid HHL

定义

特征值编码

已知实数 $\lambda\in (0,1)$ 为哈密顿矩阵 $A$ 的特征值，其二进制编码为

$$\lambda=0.b_1b_2\cdots$$

其中 $b_k\in\{0,1\}$ 是二进制编码的第 $k$ 位。$\forall n\in\mathbb{N}$，$n$ 比特估计为

$$\lambda(n)=b_1b_2\cdots b_n\approx 2^n \lambda$$

第 $k$ 特征均值

令 $\left\{\lambda_j\right\}_{j=1}^l$ 为哈密顿矩阵 $A$ 的非零特征值，可对 $k\in\mathbb{N}$ 定义第 $k$ 特征均值

$$\bar{m}_k=\frac{1}{l}\sum_{j=1}^{l} b_k^j$$

其中 $b_k^j$ 是 $\lambda_j$ 的第 $k$ 位比特。如果 $\bar{m}_k=0,1$，则称其为固定的；在这种情况下， $\forall 1\leqslant j_1,j_2\leqslant l$ 显然有 $b_k^{j_1}=b_k^{j_2}$。

$n$ 完美估计

$\lambda$ 是哈密顿矩阵 $A$ 的特征值，且 $n\in\mathbb{N}$。

1. 若 $2^n\lambda=\lambda(n)$，则称 $\lambda$ 是 $n$ 完美估计。

2. 若所有 $A$ 的特征值都是 $n$ 完美估计，则称 $A$ 是 $n$ 完美估计。

Improved HHL