量子退火与机器学习（1）：少量数据求解未知QUBO矩阵，以少见多

前言

主要是来自大関真之教授的直播课程: 【実践的量子ソリューション創出論・量子力学B・合同補講】第4回: 量子アニーリングによるブラックボックス最適化を実装する【東北大学全学教育・東北大学工学部】
这篇主要讲，怎么用少量数据去推定QUBO矩阵，然后迭代求解未知函数的方法。牵涉的知识如下：

• QUBO建模
• 压缩感知算法(Compressed Sensing)
• 稀疏建模(Sparse Modeling)
• ISTA算法(iterative shrinkage thresholding algorithm:软阈值迭代算法)

ー、复习QUBO：中药配伍的复杂性

提示：仅用公式进行问题描述，太难懂了，就举个例子，不用深究。

中药讲究配伍，即不同药材组合在一起能产生比单一药材更好的疗效，并且能减少副作用。但是，中药材之间的相互作用非常复杂，哪些药材组合在一起能更好地降血压、哪些药材组合会产生不良反应，这些都很难通过传统方式（例如人工经验）进行高效筛选。

1.QUBO 的介入：寻找最佳药材组合

QUBO 是一种数学优化技术，它特别适用于解决组合优化问题。我们可以将中药配伍问题转化为 QUBO 问题，然后利用量子退火或经典计算方法来寻找最佳的药材组合。

QUBO 如何应用于降血压中药配伍：

1. 定义二进制变量：

   • $$对于每一种可能用于降血压的中药材（比如，黄芪、决明子、菊花、钩藤、杜仲等），我们都定义一个二进制变量x_i。$$
   • $$如果x_i=1，则表示在最终的配伍中包含这种药材；如果x_i=0，则表示不包含这种药材。$$

2. 构建目标函数（成本函数）：

   • 目标函数需要反映出我们希望达成的疗效的综合打分，例如：
     • 疗效最大化： 包含能有效降低血压的药材组合。我们可以根据现有研究或实验数据，赋予每个药材一个 “降压能力” 的权重，然后尽可能选择权重高的药材组合。
     • 副作用最小化： 避免产生不良反应的药材组合。可以根据文献或实验数据，赋予每个药材一个 “副作用” 的权重，然后尽可能避免选择副作用权重高的药材组合。
     • 协同作用最大化： 鼓励选择有协同增效作用的药材组合。可以使用药材之间相互作用的实验数据来计算协同作用，并将其纳入目标函数。
   • 因此，目标函数会是这样的形式：
     $$E=\underset{对角元素}{\underbrace{\sum _i(Q_{ii}\ast x_i{x}_i)}}+\underset{上角元素}{\underbrace{\sum _{i,j(i<j)}(Q_{ij}\ast x_i{x}_j)}}$$
     • $$x_i是二进制变量，表示是否使用第i种药材。$$
     • $$Q_{ii}代表第i种药材的个体权重(例如，降压能力、副作用)。$$
     • $$Q_{ij}代表第i种和第j种药材的相互作用权重(例如，协同作用或不良反应)。$$
   • 目标是找到能使 Q 的值最小化的 xi 的组合。

3. 约束条件：

   • 有些情况下，我们可能需要加入一些约束条件，例如：
     • 配方中药材的总数不超过某个值（例如不超过5种）。
     • 必须包含某几种基础药材。
     • 必须避免某些药材同时出现。
   • 这些约束条件也会被转化为 QUBO 中的惩罚项（添加到目标函数中），以确保优化结果满足要求。

4. 优化：

   • 使用量子退火器，寻找使 QUBO 目标函数 E 最小化的二进制变量x 组合。
   • 计算出的 x_i 的值（0或1）就对应着最佳的配伍组合。

二、难题：QUBO矩阵未知的问题

1.为什么这么难？

1. 很多问题没有确定的QUBO矩阵
   比如，中药配伍的问题，你不能通过像TSP问题那样，已经知道地点位置，地点间距离，相应的约束条件。

2. 获得验证数据的周期太长或者难度太大。
   比如，中药配伍的话，你收集一个配方的实验数据，就需要很多人力物力，这样成本代价太高了，不能无限的验证下去。

已经有少量数据的情况下，怎么近似求解QUBO？

• 思路如下图:

数据足够多的话，是不是可以解方程。比如，中药配伍问题的情况，各变量的含义如下：

• x 变量就是用或者不用某位药,n维就代表有n种药。
• b变量就是每次不同中药组合的测量后的综合药效列表，假定有m个。
• a 就是每次不同的QUBO矩阵上三角里的元素n列表。a 是无数种可能的，但是里面肯定有一个是我们想要的接近现实的解。
  
  上面的式子有的难懂，给大家举个实例。
• x

• a

• b

上面的式子变换一下：

下面解释一下变换后的式子中各个变量的维度:

1. 向量 b 是 m 维向量: b ∈ ℝᵐ

2. 矩阵 A = [a⁽¹⁾, …, a⁽ⁿ⁾] 的维度是:

   • 每个 a⁽ⁱ⁾ 是 m 维向量
   • 一共有 n 个这样的向量
   • 所以 A 的维度是 m × n

3. x是 n(n+1)/2 维向量（QUBO矩阵的上三角里所有元素）: $$x\in R^{n(n+1)/2}$$

4. 通过矩阵乘法 Ax:

   • A(m×n) × x(n×1) = b(m×1)
   • 结果 b 是 m 维向量，与原始定义一致

三、稀疏建模(Sparse Modeling)

线性方程组大家都知道，学完线性代数，也都知道可以换成矩阵形式。我就直接贴上wiki截图了。
https://zh.wikipedia.org/zh-cn/%E7%BA%BF%E6%80%A7%E6%96%B9%E7%A8%8B%E7%BB%84

• 一般情况下，1个方程解1个未知数，2个方程解2个未知数，这是我们平时接触较多的求解线性系统的情况，称之为适定系统。
  那如果，一个方程有两个未知数呢？这种情况就是欠定系统了。
  在压缩感知理论中，一般用下列式子来表示一个欠定系统：
  $$\mathbf{b}=\mathbf{Ax}$$
  $$其中{\mathbb{R}}^{M\times N},X\in {\mathbb{R}}^M,b\in {\mathbb{R}}^N.且当M<N时，系统维欠定系统.$$

方程组的数量不足意味着决定解的条件不足。由于条件不足，如果再增加一些条件就可以确定解。
例如，如果预先知道解，通过将其代入，就可以有效地减少N。现在假设已知解，且x的各分量几乎为0。
在这种情况下，可以从方程组中删除值为0的分量。如果将非零项的数量记为K，那么从M个方程实际上就是求解K个非零分量，即使M < N，只要M > K，就可以求解。

这种大部分分量为零或预期为零的性质称为"稀疏性"，具有这种性质的解称为"稀疏解"。

1. 欠定系统中的稀疏解

下面的所有截图都在这个日文资料里：https://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/Presentation/lecturenote20160727.pdf
对于N维的未知向量x，M维的实数值向量b和M × N的观测矩阵A，假设满足以下关系：

这里即使M < N，当x的分量中大部分为零（具有稀疏性）时，如果非零分量的数量K满足M > K，就可以求得解。
然而，这K个非零分量究竟在哪里？这是未知的。那么如何求解呢？
虽然遗憾，但没有决定性的方法，只能从N个分量中选取K个分量，寻找满足y = Ax的解。从N个中选取K个的组合数，随着N的增大会呈指数级增长。对于高维问题，进行这样的计算在现实中是不可行的。而且虽然说是K个非零项，但K这个数字真的已知吗？这也不一定知道。

因此，当这些非零分量的数量也未知时，应该采取什么样的策略来寻找满足b = Ax的解呢？
其实就是用各种正则化L0,L1,L2正则。

• L0正则：
  
  $$\Vert \mathbf{x}{\Vert }_0，代表非0解的个数。越小越稀疏。$$

• L1正则：
  
  $$\Vert \mathbf{x}{\Vert }_1=|x_1|+|x_2|+\cdots +|x_N|,代表x的绝对值总和。0越多，\Vert \mathbf{x}{\Vert }_1越小越稀疏。$$

• L2正则：
  
  $$\Vert \mathbf{x}{\Vert }_2=\sqrt{x_1^2+x_2^2+\cdots +x_N^2}$$

2. L1和L2的选择：

下面的图是一个L1和L2求解的结果，明显L1成功获得了真实解，L2失败。

三、压缩感知算法(Compressed Sensing)

根据正则的性质，我们已经知道可以获得这样的解选择技术。

• 这时需要思考的问题是：我们真的需要稀疏解吗？真正的解是稀疏解吗？
  • 前者关注的是变量选择的问题。当我们对方程的真实解不感兴趣，而只是在寻找能满足方程的最少变量组合时，这是一个重要的问题。
  • 至于后者，当我们不是要选择变量而是要寻找真实解时，就需要考虑稀疏解是否合适。对于本质上具有稀疏解的方程问题，有选择性地找出稀疏解会产生巨大的效果。

压缩感知这个框架是利用正则的特性从欠定方程组中获得稀疏解，从而更准确地确定我们想要了解的内容。它就像信息科学中的名侦探。
特别是，通过L1范数最小化来估计原始信息的方法被称为基追踪（Basis Pursuit）。

1. 压缩感知的性质

当观测矩阵A的各分量从均值为0、方差为1的高斯分布生成时，以下列曲线为边界，在α较大且ρ较小的区域内，通过L1正则最小化可以以极高的概率成功恢复原始信号。其中α = M/N，ρ = K/NP，Q(t)是标准正态分布的尾部概率积分。
1 α = 1 + π 2 t e t 2 2 { 1 − 2 Q ( t ) } \frac{1}{\alpha} = 1 + \sqrt{\frac{\pi}{2}}te^{\frac{t^2}{2}}\{1-2Q(t)\} α1​=1+2π​ ​te2t2​{1−2Q(t)}
$$\frac{\rho }{1-\rho }=2\left(\frac{e^{-\frac{t^2}{2}}}{t\sqrt{2\pi }}-Q(t)\right)$$
$$Q(t)={\int }_t^{\infty }\frac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi }}dx$$
$$\alpha =\frac{M}{N},\rho =\frac{K}{N}$$

上图展示了压缩感知中L1正则最小化重构的可行性边界。让我详细解释一下：

1. 坐标轴含义：

• 横轴 ρ = K/N：表示稀疏度（信号中非零元素的比例）
• 纵轴 α = M/N：表示测量数与信号维度的比值（压缩比）

2. 图中的区域：

• 蓝色区域：这是L1范数最小化能够成功重构原始信号的区域
  • 当(ρ,α)点落在这个区域内时，我们可以以很高的概率通过L1最小化重构出原始信号
  • 特别是在α较大（即测量数较多）且ρ较小（即信号较稀疏）的情况下，重构成功率最高

3. 分界线：

• 实线曲线：表示L1重构的理论边界
• 虚线 α = ρ：这条对角线表示测量数等于非零元素个数的情况

4. 实际意义：

• 这个图帮助我们理解在给定信号稀疏度ρ的情况下，需要多少测量值（由α决定）才能成功重构
• 在蓝色区域内，压缩感知是有效的，即可以用少量测量重构出原始信号
• 区域外则表示测量数不足，无法保证信号重构的成功

这个图对于实际应用压缩感知非常有用，它可以帮助我们确定所需的最小测量数，以保证可以成功重构具有特定稀疏度的信号。下面这句话很重要，我说三遍。

• 压缩感知，重要的不仅仅是选择稀疏解，关键在于不能仅仅是选择"差不多"解，还需要其中包含正确答案。
• 压缩感知，重要的不仅仅是选择稀疏解，关键在于不能仅仅是选择"差不多"解，还需要其中包含正确答案。
• 压缩感知，重要的不仅仅是选择稀疏解，关键在于不能仅仅是选择"差不多"解，还需要其中包含正确答案。

2. ISTA算法

ISTA是一个通过L1正则化，迭代求解欠定系统的算法，流程如下（证明自己网上可查）：

1. 令t = 0，初始化x[0]。例如可以设置$$x[0]=A^{T}y$$

2. 通过平方完成法求解g(x)的二次函数近似的顶点：
   $$v[t]=x[t]+(1/L\lambda )A^T(y-Ax[t])$$

3. 应用软阈值函数：
   $$x[t+1]=S_{1/L}(v[t])$$

4. 重复步骤2-4直到满足终止条件。

四、Python实现

import numpy as np
import matplotlib.pyplot as plt
from openjij import SASampler
from IPython.display import clear_output

def grad_comp(y, A, x):
    """
    计算梯度
    Args:
        y: 观测值向量
        A: 测量矩阵
        x: 当前解向量
    Returns:
        grad: 梯度向量
    """
    grad = -np.dot(A.T, (y - A.dot(x)))
    return grad

def SoftThr(v, thr):
    """
    软阈值函数实现
    Args:
        v: 输入向量
        thr: 阈值
    Returns:
        z: 经过软阈值处理的向量
    """
    z = np.zeros(len(v))
    # 处理大于阈值的元素
    itemp = np.where(v > thr)
    z[itemp] = v[itemp] - thr
    # 处理小于-阈值的元素
    itemp = np.where(v <= -thr)
    z[itemp] = v[itemp] + thr
    return z

def opt_qvec(x, x0, y, A, Tall=10, p=10.0, flag=True):
    """
    使用ADMM算法优化QUBO向量
    Args:
        x: 初始解向量
        x0: 目标解向量
        y: 观测值向量
        A: 测量矩阵
        Tall: 最大迭代次数
        p: ADMM惩罚参数
        flag: 是否显示优化过程图像
    Returns:
        x: 优化后的解向量
    """
    N = A.shape[0]
    # 计算A的伪逆相关矩阵
    Atemp = A.dot(A.T)
    Ainv = np.linalg.inv(Atemp)
    Atemp = A.T.dot(Ainv)
    Nvec = len(x)

    # ADMM算法的辅助变量
    z = np.zeros(Nvec)
    u = np.zeros(Nvec)
    
    # ADMM迭代
    for t in range(Tall):
        # 更新x
        x = Atemp.dot(y) + (np.eye(Nvec) - Atemp.dot(A)).dot(z + u)
        # 更新z（软阈值步骤）
        z = SoftThr(x - u, 1/p)
        # 更新对偶变量u
        u = u + (z - x)

        # 如果需要，绘制优化过程
        if flag:
            clear_output(True)
            plt.plot(x)
            plt.plot(x0)
            plt.show()

    return x

def Xmat_make(x):
    """
    构造QUBO问题的特征向量
    Args:
        x: 输入向量
    Returns:
        Xvec: 包含一阶项和二阶项的特征向量
    """
    Ns = len(x)
    # 向量长度为一阶项数量加上二阶项数量
    Xvec = np.zeros(Ns + Ns*(Ns-1)//2)

    # 填充一阶项
    t = 0
    for i in range(Ns):
        Xvec[t] = x[i]
        t = t + 1

    # 填充二阶项（交互项）
    for i in range(Ns):
        for j in range(Ns):
            if i < j:
                Xvec[t] = x[i]*x[j]
                t = t + 1

    return Xvec

def ycomp(Xvec, Qvec):
    """
    计算QUBO问题的能量
    Args:
        Xvec: 特征向量
        Qvec: QUBO系数向量
    Returns:
        Ene: 能量值
    """
    Ene = np.dot(Xvec, Qvec)
    return Ene

def QUBO_create(Qvec, Ns):
    """
    从向量形式构造QUBO矩阵
    Args:
        Qvec: QUBO系数向量
        Ns: 系统大小
    Returns:
        QUBO: QUBO矩阵
    """
    # 计算二阶项数量
    Noff = (Ns*(Ns-1))//2
    # 提取对角项（一阶项）和非对角项（二阶项）
    Qdiag = Qvec[:Ns]
    Qoff = Qvec[Ns:]

    # 构造QUBO矩阵
    QUBO = np.diag(Qdiag)

    # 填充非对角元素
    t = 0
    for i in range(Ns):
        for j in range(Ns):
            if i < j:
                QUBO[i,j] = Qoff[t]
                t = t + 1

    return QUBO

# 主程序开始

# 设置系统大小
Ns = 20

# 生成随机的QUBO问题
# 生成对角项
Qdiag = np.random.randn(Ns)
QUBO = np.diag(Qdiag)

# 生成稀疏的非对角项
Noff = (Ns*(Ns-1))//2
Qoff = np.random.randn(Noff)
rho = 0.2  # 稀疏度参数
mask = (np.random.rand(Noff) < rho)
Qoff = mask*Qoff

# 合并对角项和非对角项
Qvec = np.concatenate((Qdiag,Qoff))

# 生成训练数据
M = 100  # 训练样本数
Adata = []  # 特征矩阵
ydata = []  # 能量值

# 随机生成训练样本
for d in range(M):
    # 生成随机二值向量
    x = (np.random.rand(Ns) > 0.5)
    x = x.astype(np.int16)
    # 计算特征向量和对应能量
    Xvec = Xmat_make(x)
    Ene = ycomp(Xvec,Qvec)
    Adata.append(Xvec)
    ydata.append(Ene)

# 将数据转换为numpy数组
y = np.array(ydata)
A = np.array(Adata)

# 使用ADMM算法学习QUBO参数
Nvec = Noff + Ns
Qinf = np.zeros(Nvec)
Qinf = opt_qvec(Qinf, Qvec, y, A, Tall=100)

# 构造学习到的QUBO矩阵
QUBO = QUBO_create(Qinf, Ns)

# 使用量子退火采样器求解QUBO问题
sampler = SASampler()
sampleset = sampler.sample_qubo(QUBO, num_reads=1)

# 迭代优化过程
Ns = 20
Ndata = 5  # 初始数据点数
Nall = 195  # 总迭代次数

# 初始化数据集
Adata = []
ydata = []

for d in range(Ndata):
    x = (np.random.rand(Ns) > 0.5)
    x = x.astype(np.int16)
    Xvec = Xmat_make(x)
    Ene = ycomp(Xvec,Qvec)
    Adata.append(Xvec)
    ydata.append(Ene)

# 记录优化过程中的能量
Enelist = []
Eneminlist = []
xlist = []
Qinf = np.dot(A.T,y)

# 主优化循环
for d in range(Nall):
    # 更新QUBO参数
    y = np.array(ydata)
    A = np.array(Adata)
    Qinf = opt_qvec(Qinf, Qvec, y, A, Tall=10, flag=False)
    QUBO = QUBO_create(Qinf, Ns)
    
    # 使用量子退火采样器获得新解
    sampleset = sampler.sample_qubo(QUBO, num_reads=1)
    x = sampleset.record[0][0]

    # 检查是否重复解
    for xtemp in xlist:
        if np.array_equal(x,xtemp):
            x = (np.random.rand(Ns) > 0.5)
            x = x.astype(np.int16)
            break

    xlist.append(x)

    # 计算新解的能量
    Xvec = Xmat_make(x)
    Ene = ycomp(Xvec,Qvec)
    Enelist.append(Ene)
    Enemin = np.min(Enelist)
    Eneminlist.append(Enemin)

    # 更新数据集
    ydata.append(Ene)
    Adata.append(Xvec)

    # 绘制优化过程
    clear_output(True)
    plt.plot(Enelist)
    plt.plot(Eneminlist)
    plt.show()

1. 代码和结果解释

1.1 代码细节

代码其实挺简化，但我们也可以从中看到一些细节：

• Noff = int((Ns * (Ns - 1)) / 2) 计算了QUBO矩阵中非对角线的个数。
• mask 的作用是只考虑稀疏的那些Qij。
• np.random.rand(Ns) 在模拟实验中用于产生随机的01向量。
• opt_qvec 是关键的函数，里面通过数据拟合Q矩阵，并用此Q矩阵进行退火优化。

1.2 总体思路回顾：

1. 目标： 使用模拟退火算法（SA）或者量子退火算法（QA）来找到一个QUBO问题的最优解，但QUBO矩阵本身是未知的（“黑盒”）。
2. 难点： QUBO矩阵是未知的，我们无法直接使用标准的退火方法。
3. 解决方法： 使用压缩感知算法，逐步猜测和逼近真实的QUBO矩阵，并在这个过程中利用退火算法进行优化。
4. 关键： 从客户（黑盒）那里获得数据，然后用这些数据来推断Q矩阵。

1.3 压缩感知算法的应用：

使用压缩感知算法的核心体现在opt_qvec函数内部和整个迭代过程中，它的思想是：

1. 稀疏性假设： 假设QUBO矩阵是稀疏的（即有很多元素为零）。
2. 数据采集： 通过不断询问（比如做问卷，问专家）黑盒获取数据，可以理解为通过不断迭代模拟退火算法来寻找更好的01向量。
3. 逐步逼近： 使用采集到的数据，反推（拟合）出一个稀疏的QUBO矩阵。
4. 更新和迭代： 然后使用这个推导出的Q矩阵进行退火，并继续这个采样拟合的过程，直到找到一个比较好的Q矩阵来推断。

1.4 最后的输出结果解读：

• 最终图像部分：
  • x轴表示退火优化的迭代步骤，y轴表示能量值。
  • 蓝色曲线：表示模拟退火算法在尝试优化（寻找更低的能量）过程中，每个采样点所对应的能量值
  • 橘色曲线：真实情况的能量值，用来对比模拟退火算法找到的解和真实解之间的差距。
• 解读：
  • 数据与优化协同作用： 这种蓝色线和黄色线的同步下降，生动地展示了压缩感知算法的核心——通过模拟退火（或量子退火）算法的优化搜索，不断引导QUBO矩阵的逼近，同时利用新的01向量的数据，使推导的矩阵越来越精确，最终在黑盒优化问题中找到好的解。
  • 蓝色尖峰出现： 蓝色线的尖峰，通常表示模拟退火算法在搜索过程中，随机尝试到了一个能量比较高的状态。这是退火算法的探索性的一部分，它会尝试从当前的局部最优解“跳出”，看看是否有更低的能量值。这种尖峰通常表示对目前解的否定。

四、总结

这个教程，真的是很直观地讲解了最先进的QUBO建模技术，以少见多。有什么问题，欢迎指正。
下一篇，更深入的讲解ISTA算法的升级版ADMM算法