关于12380网站建设文件,南京企业网站设计建设,桂林北站怎么去阳朔,精准引流怎么推广文章目录 前言ー、复习QUBO#xff1a;中药配伍的复杂性1.QUBO 的介入#xff1a;寻找最佳药材组合 二、难题#xff1a;QUBO矩阵未知的问题1.为什么这么难#xff1f; 三、稀疏建模(Sparse Modeling)1. 欠定系统中的稀疏解2. L1和L2的选择#xff1a; 三、压缩感知算法(C… 文章目录 前言ー、复习QUBO中药配伍的复杂性1.QUBO 的介入寻找最佳药材组合 二、难题QUBO矩阵未知的问题1.为什么这么难 三、稀疏建模(Sparse Modeling)1. 欠定系统中的稀疏解2. L1和L2的选择 三、压缩感知算法(Compressed Sensing)1. 压缩感知的性质2. ISTA算法 四、Python实现1. 代码和结果解释 四、总结 前言
主要是来自大関真之教授的直播课程: 【実践的量子ソリューション創出論・量子力学B・合同補講】第4回: 量子アニーリングによるブラックボックス最適化を実装する【東北大学全学教育・東北大学工学部】 这篇主要讲怎么用少量数据去推定QUBO矩阵然后迭代求解未知函数的方法。牵涉的知识如下
QUBO建模压缩感知算法(Compressed Sensing)稀疏建模(Sparse Modeling)ISTA算法(iterative shrinkage thresholding algorithm:软阈值迭代算法) ー、复习QUBO中药配伍的复杂性
提示仅用公式进行问题描述太难懂了就举个例子不用深究。
中药讲究配伍即不同药材组合在一起能产生比单一药材更好的疗效并且能减少副作用。但是中药材之间的相互作用非常复杂哪些药材组合在一起能更好地降血压、哪些药材组合会产生不良反应这些都很难通过传统方式例如人工经验进行高效筛选。
1.QUBO 的介入寻找最佳药材组合
QUBO 是一种数学优化技术它特别适用于解决组合优化问题。我们可以将中药配伍问题转化为 QUBO 问题然后利用量子退火或经典计算方法来寻找最佳的药材组合。
QUBO 如何应用于降血压中药配伍 定义二进制变量 对于每一种可能用于降血压的中药材比如黄芪、决明子、菊花、钩藤、杜仲等我们都定义一个二进制变量 x i 。 对于每一种可能用于降血压的中药材比如黄芪、决明子、菊花、钩藤、杜仲等我们都定义一个二进制变量 x_i。 对于每一种可能用于降血压的中药材比如黄芪、决明子、菊花、钩藤、杜仲等我们都定义一个二进制变量xi。 如果 x i 1 则表示在最终的配伍中包含这种药材如果 x i 0 则表示不包含这种药材。 如果 x_i 1则表示在最终的配伍中包含这种药材如果 x_i 0则表示不包含这种药材。 如果xi1则表示在最终的配伍中包含这种药材如果xi0则表示不包含这种药材。 构建目标函数成本函数 目标函数需要反映出我们希望达成的疗效的综合打分例如 疗效最大化 包含能有效降低血压的药材组合。我们可以根据现有研究或实验数据赋予每个药材一个 “降压能力” 的权重然后尽可能选择权重高的药材组合。副作用最小化 避免产生不良反应的药材组合。可以根据文献或实验数据赋予每个药材一个 “副作用” 的权重然后尽可能避免选择副作用权重高的药材组合。协同作用最大化 鼓励选择有协同增效作用的药材组合。可以使用药材之间相互作用的实验数据来计算协同作用并将其纳入目标函数。 因此目标函数会是这样的形式 E ∑ i ( Q i i ∗ x i x i ) ⏟ 对角元素 ∑ i , j ( i j ) ( Q i j ∗ x i x j ) ⏟ 上角元素 E \underbrace{\sum_i(Q_{ii} * x_ix_i)}_{对角元素} \underbrace{\sum_{i,j(ij)}(Q_{ij} * x_ix_j)}_{上角元素} E对角元素 i∑(Qii∗xixi)上角元素 i,j(ij)∑(Qij∗xixj) x i 是二进制变量表示是否使用第 i 种药材。 x_i是二进制变量表示是否使用第i种药材。 xi是二进制变量表示是否使用第i种药材。 Q i i 代表第 i 种药材的个体权重 ( 例如降压能力、副作用 ) 。 Q_{ii}代表第i种药材的个体权重 (例如降压能力、副作用)。 Qii代表第i种药材的个体权重(例如降压能力、副作用)。 Q i j 代表第 i 种和第 j 种药材的相互作用权重 ( 例如协同作用或不良反应 ) 。 Q_{ij}代表第i种和第j种药材的相互作用权重 (例如协同作用或不良反应)。 Qij代表第i种和第j种药材的相互作用权重(例如协同作用或不良反应)。 目标是找到能使 Q 的值最小化的 xi 的组合。 约束条件 有些情况下我们可能需要加入一些约束条件例如 配方中药材的总数不超过某个值例如不超过5种。必须包含某几种基础药材。必须避免某些药材同时出现。 这些约束条件也会被转化为 QUBO 中的惩罚项添加到目标函数中以确保优化结果满足要求。 优化 使用量子退火器寻找使 QUBO 目标函数 E 最小化的二进制变量x 组合。计算出的 x_i 的值0或1就对应着最佳的配伍组合。
二、难题QUBO矩阵未知的问题
1.为什么这么难 很多问题没有确定的QUBO矩阵 比如中药配伍的问题你不能通过像TSP问题那样已经知道地点位置地点间距离相应的约束条件。 获得验证数据的周期太长或者难度太大。 比如中药配伍的话你收集一个配方的实验数据就需要很多人力物力这样成本代价太高了不能无限的验证下去。
已经有少量数据的情况下怎么近似求解QUBO
思路如下图: 数据足够多的话是不是可以解方程。比如中药配伍问题的情况各变量的含义如下
x 变量就是用或者不用某位药,n维就代表有n种药。b变量就是每次不同中药组合的测量后的综合药效列表假定有m个。a 就是每次不同的QUBO矩阵上三角里的元素n列表。a 是无数种可能的但是里面肯定有一个是我们想要的接近现实的解。 上面的式子有的难懂给大家举个实例。x
x_typevaluex₁1x₂0x₃1x₁x₂0x₁x₃1x₂x₃0
a
a⁽¹⁾a⁽²⁾a⁽³⁾a⁽⁴⁾a⁽⁵⁾a⁽⁶⁾a₁0.5-0.30.8-0.40.6-0.7a₂-0.60.7-0.20.5-0.80.2a₃0.4-0.50.9-0.60.3-0.4
b
bvalueb₁1.9b₂-1.6b₃1.6
上面的式子变换一下 下面解释一下变换后的式子中各个变量的维度: 向量 b 是 m 维向量: b ∈ ℝᵐ 矩阵 A [a⁽¹⁾, …, a⁽ⁿ⁾] 的维度是: 每个 a⁽ⁱ⁾ 是 m 维向量一共有 n 个这样的向量所以 A 的维度是 m × n x是 n(n1)/2 维向量QUBO矩阵的上三角里所有元素: x ∈ R n ( n 1 ) / 2 x ∈ ℝ^{n(n1)/2} x∈Rn(n1)/2 通过矩阵乘法 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个未知数这是我们平时接触较多的求解线性系统的情况称之为适定系统。 那如果一个方程有两个未知数呢这种情况就是欠定系统了。 在压缩感知理论中一般用下列式子来表示一个欠定系统 b A x \mathbf{b} \mathbf{A} \mathbf{x} bAx 其中 R M × N , X ∈ R M , b ∈ R N . 且当 M N 时系统维欠定系统 . 其中\mathbb{R}^{M \times N}, X \in \mathbb{R}^M, b \in \mathbb{R}^N. 且当M N时系统维欠定系统. 其中RM×N,X∈RM,b∈RN.且当MN时系统维欠定系统.
方程组的数量不足意味着决定解的条件不足。由于条件不足如果再增加一些条件就可以确定解。 例如如果预先知道解通过将其代入就可以有效地减少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维的未知向量xM维的实数值向量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正则 ∥ x ∥ 0 代表非 0 解的个数。越小越稀疏。 \|\mathbf{x}\|_0代表非0解的个数。越小越稀疏。 ∥x∥0代表非0解的个数。越小越稀疏。 L1正则 ∥ x ∥ 1 ∣ x 1 ∣ ∣ x 2 ∣ ⋯ ∣ x N ∣ , 代表 x 的绝对值总和。 0 越多 ∥ x ∥ 1 越小越稀疏。 \|\mathbf{x}\|_1 |x_1| |x_2| \cdots |x_N|, 代表x的绝对值总和。0越多\|\mathbf{x}\|_1越小越稀疏。 ∥x∥1∣x1∣∣x2∣⋯∣xN∣,代表x的绝对值总和。0越多∥x∥1越小越稀疏。 L2正则 ∥ x ∥ 2 x 1 2 x 2 2 ⋯ x N 2 \|\mathbf{x}\|_2 \sqrt{x_1^2 x_2^2 \cdots x_N^2} ∥x∥2x12x22⋯xN2
2. L1和L2的选择 下面的图是一个L1和L2求解的结果明显L1成功获得了真实解L2失败。
三、压缩感知算法(Compressed Sensing)
根据正则的性质我们已经知道可以获得这样的解选择技术。
这时需要思考的问题是我们真的需要稀疏解吗真正的解是稀疏解吗 前者关注的是变量选择的问题。当我们对方程的真实解不感兴趣而只是在寻找能满足方程的最少变量组合时这是一个重要的问题。至于后者当我们不是要选择变量而是要寻找真实解时就需要考虑稀疏解是否合适。对于本质上具有稀疏解的方程问题有选择性地找出稀疏解会产生巨大的效果。
压缩感知这个框架是利用正则的特性从欠定方程组中获得稀疏解从而更准确地确定我们想要了解的内容。它就像信息科学中的名侦探。 特别是通过L1范数最小化来估计原始信息的方法被称为基追踪Basis Pursuit。
1. 压缩感知的性质
当观测矩阵A的各分量从均值为0、方差为1的高斯分布生成时以下列曲线为边界在α较大且ρ较小的区域内通过L1正则最小化可以以极高的概率成功恢复原始信号。其中α M/Nρ K/NPQ(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)\} α112π te2t2{1−2Q(t)} ρ 1 − ρ 2 ( e − t 2 2 t 2 π − Q ( t ) ) \frac{\rho}{1-\rho} 2\left(\frac{e^{-\frac{t^2}{2}}}{t\sqrt{2\pi}}-Q(t)\right) 1−ρρ2(t2π e−2t2−Q(t)) Q ( t ) ∫ t ∞ e − x 2 2 2 π d x Q(t) \int_t^{\infty}\frac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi}}dx Q(t)∫t∞2π e−2x2dx α M N , ρ K N \alpha \frac{M}{N}, \quad \rho \frac{K}{N} αNM,ρNK 上图展示了压缩感知中L1正则最小化重构的可行性边界。让我详细解释一下
坐标轴含义
横轴 ρ K/N表示稀疏度信号中非零元素的比例纵轴 α M/N表示测量数与信号维度的比值压缩比
图中的区域
蓝色区域这是L1范数最小化能够成功重构原始信号的区域 当(ρ,α)点落在这个区域内时我们可以以很高的概率通过L1最小化重构出原始信号特别是在α较大即测量数较多且ρ较小即信号较稀疏的情况下重构成功率最高
分界线
实线曲线表示L1重构的理论边界虚线 α ρ这条对角线表示测量数等于非零元素个数的情况
实际意义
这个图帮助我们理解在给定信号稀疏度ρ的情况下需要多少测量值由α决定才能成功重构在蓝色区域内压缩感知是有效的即可以用少量测量重构出原始信号区域外则表示测量数不足无法保证信号重构的成功
这个图对于实际应用压缩感知非常有用它可以帮助我们确定所需的最小测量数以保证可以成功重构具有特定稀疏度的信号。下面这句话很重要我说三遍。
压缩感知重要的不仅仅是选择稀疏解关键在于不能仅仅是选择差不多解还需要其中包含正确答案。压缩感知重要的不仅仅是选择稀疏解关键在于不能仅仅是选择差不多解还需要其中包含正确答案。压缩感知重要的不仅仅是选择稀疏解关键在于不能仅仅是选择差不多解还需要其中包含正确答案。
2. ISTA算法
ISTA是一个通过L1正则化迭代求解欠定系统的算法流程如下证明自己网上可查 令t 0初始化x[0]。例如可以设置 x [ 0 ] A T y x[0] A^T y x[0]ATy 通过平方完成法求解g(x)的二次函数近似的顶点 v [ t ] x [ t ] ( 1 / L λ ) A T ( y − A x [ t ] ) v[t] x[t] (1/Lλ)A^T(y - Ax[t]) v[t]x[t](1/Lλ)AT(y−Ax[t]) 应用软阈值函数 x [ t 1 ] S 1 / L ( v [ t ] ) x[t1] S_{1/L}(v[t]) x[t1]S1/L(v[t]) 重复步骤2-4直到满足终止条件。
四、Python实现
import numpy as np
import matplotlib.pyplot as plt
from openjij import SASampler
from IPython.display import clear_outputdef grad_comp(y, A, x):计算梯度Args:y: 观测值向量A: 测量矩阵x: 当前解向量Returns:grad: 梯度向量grad -np.dot(A.T, (y - A.dot(x)))return graddef 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] thrreturn zdef opt_qvec(x, x0, y, A, Tall10, p10.0, flagTrue):使用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):# 更新xx Atemp.dot(y) (np.eye(Nvec) - Atemp.dot(A)).dot(z u)# 更新z软阈值步骤z SoftThr(x - u, 1/p)# 更新对偶变量uu u (z - x)# 如果需要绘制优化过程if flag:clear_output(True)plt.plot(x)plt.plot(x0)plt.show()return xdef Xmat_make(x):构造QUBO问题的特征向量Args:x: 输入向量Returns:Xvec: 包含一阶项和二阶项的特征向量Ns len(x)# 向量长度为一阶项数量加上二阶项数量Xvec np.zeros(Ns Ns*(Ns-1)//2)# 填充一阶项t 0for 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 1return Xvecdef ycomp(Xvec, Qvec):计算QUBO问题的能量Args:Xvec: 特征向量Qvec: QUBO系数向量Returns:Ene: 能量值Ene np.dot(Xvec, Qvec)return Enedef 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 0for i in range(Ns):for j in range(Ns):if i j:QUBO[i,j] Qoff[t]t t 1return 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, Tall100)# 构造学习到的QUBO矩阵
QUBO QUBO_create(Qinf, Ns)# 使用量子退火采样器求解QUBO问题
sampler SASampler()
sampleset sampler.sample_qubo(QUBO, num_reads1)# 迭代优化过程
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, Tall10, flagFalse)QUBO QUBO_create(Qinf, Ns)# 使用量子退火采样器获得新解sampleset sampler.sample_qubo(QUBO, num_reads1)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)breakxlist.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 总体思路回顾
目标 使用模拟退火算法SA或者量子退火算法QA来找到一个QUBO问题的最优解但QUBO矩阵本身是未知的“黑盒”。难点 QUBO矩阵是未知的我们无法直接使用标准的退火方法。解决方法 使用压缩感知算法逐步猜测和逼近真实的QUBO矩阵并在这个过程中利用退火算法进行优化。关键 从客户黑盒那里获得数据然后用这些数据来推断Q矩阵。
1.3 压缩感知算法的应用
使用压缩感知算法的核心体现在opt_qvec函数内部和整个迭代过程中它的思想是
稀疏性假设 假设QUBO矩阵是稀疏的即有很多元素为零。数据采集 通过不断询问比如做问卷问专家黑盒获取数据可以理解为通过不断迭代模拟退火算法来寻找更好的01向量。逐步逼近 使用采集到的数据反推拟合出一个稀疏的QUBO矩阵。更新和迭代 然后使用这个推导出的Q矩阵进行退火并继续这个采样拟合的过程直到找到一个比较好的Q矩阵来推断。
1.4 最后的输出结果解读
最终图像部分 x轴表示退火优化的迭代步骤y轴表示能量值。蓝色曲线表示模拟退火算法在尝试优化寻找更低的能量过程中每个采样点所对应的能量值橘色曲线真实情况的能量值用来对比模拟退火算法找到的解和真实解之间的差距。 解读 数据与优化协同作用 这种蓝色线和黄色线的同步下降生动地展示了压缩感知算法的核心——通过模拟退火或量子退火算法的优化搜索不断引导QUBO矩阵的逼近同时利用新的01向量的数据使推导的矩阵越来越精确最终在黑盒优化问题中找到好的解。蓝色尖峰出现 蓝色线的尖峰通常表示模拟退火算法在搜索过程中随机尝试到了一个能量比较高的状态。这是退火算法的探索性的一部分它会尝试从当前的局部最优解“跳出”看看是否有更低的能量值。这种尖峰通常表示对目前解的否定。
四、总结
这个教程真的是很直观地讲解了最先进的QUBO建模技术以少见多。有什么问题欢迎指正。 下一篇更深入的讲解ISTA算法的升级版ADMM算法
