最前沿・量子退火建模方法(1) : subQUBO讲解和python实现

前言

量子退火机在小规模问题上的效果得到了有效验证,但是由于物理量子比特的大规模制备以及噪声的影响,还没有办法再大规模的场景下应用。
这时候就需要我们思考,如何通过软件的方法怎么样把大的问题分解成小的问题,以便通过现在小规模的量子退火机解决。主要思路就是,同样的QUBO建模,怎么使用更少的量子比特。

下面的文章中,量子退火机伊辛机会混用。


一、subQUBO的创新点

先行的研究中,使用启发式方法将大型问题划分为较小的问题,并使用伊辛机进行求解,但划分后的问题的答案与原始大型问题的答案并不相同。达成协议的理论条件仍不清楚。早稻田大学的研究者开发出了subQUBO算法在保证分解后的小问题也能保证在原始大问题上的理论上做出了突破。

Atobe, Yuta, Masashi Tawada, and Nozomu Togawa. "Hybrid annealing method based on subQUBO model extraction with multiple solution instances." IEEE Transactions on Computers 71.10 (2021): 2606-2619.

subQUBO的创新点

  1. 首先研究将大规模组合优化问题划分为较小问题而不失去最优性的条件。该条件成立的话就证明,如果用伊辛机来解决一个满足条件的小问题,它就会和原来的大规模问题的答案相匹配。
  2. 还提出了一种新算法,成功地从大规模问题中提取出此类条件,将原始大规模问题缩小到伊辛机可以解决的问题规模,并迭代求解。所提出的算法通过基于理论支持将大规模问题分解为更小的问题来解决它,使得以比传统技术更高的精度解决原始大规模问题成为可能。
    在这里插入图片描述

二、subQUBO的详细思路

1. 怎么把大规模问题分解成小问题

1.1 逻辑前提:挑出错误后,回炉重造

  • 大规模组合优化问题的QUBO建模中,最终的答案由多个量子比特集合组成。
  • 如果你创建一个小规模问题,其中包括最终解的量子比特集合中的,所有不正确的量子比特集合
  • 并使用伊辛机解决该问题,则所有最终解的不正确的量子比特集合都将被纠正为正确的量子比特集合作为解。

1.2 具体实现:

实现方法: 可以创建一个大致包含所有不正确的量子比特集合的小问题,并使用伊辛机重复解决它。

  • 不正确的量子比特集合创建:
    – 我们使用传统的经典计算器来准备问题的多个候选答案。这些候选答案不一定是正确的,但在比较经典计算器求解得到的多个答案的量子比特集合的最终值。
    – 多个候选中匹配一致的就是正确的量子比特集合
    – 答案不匹配且不同的就是不正确的量子比特集合

  • 通过仅提取不正确的量子比特集合,并使用真实的伊辛机进行求解,最终可以获得整体的正确答案。

1.3 业界影响:

传统上,伊辛机很难解决大规模问题,因为可用位数受到硬件限制,但通过使用这种方法,可以使用伊辛机进行计算。因此,人们认为可以使用伊辛机(包括量子退火机)扩展现实世界组合优化问题的用例。此外,本研究尝试将经典计算机与伊辛机相结合来解决问题,这将大大扩展伊辛机的使用范围。

最新成果,参考以下新闻:
Quanmatic Co., Ltd.利用量子计算技术解决方案规模突破1亿比特

https://prtimes.jp/main/html/rd/p/000000015.000117406.html

三、subQUBO的python实现

  1. 导入库
import random
import itertools
import numpy as np
from dataclasses import dataclass
  1. 设置subQUBO所需参数
N_I = 20 # instance池
N_E = 10 # subQUBO的抽取次数
N_S = 10 # N_I个instance池中抽取的解的个数
sub_qubo_size = 5 # subQUBO的量子比特数
  1. QUBO建模

# 为了简单,使用TSP作为例子
NUM_CITY = 4
ALPHA = 1
np.random.seed(0)
num_spin = NUM_CITY ** 2

distance_mtx = np.random.randint(1, 100, (NUM_CITY, NUM_CITY))
distance_mtx = np.tril(distance_mtx) + np.tril(distance_mtx).T - 2 * np.diag(distance_mtx.diagonal())

# <<< Objective term >>>
qubo_obj = np.zeros((NUM_CITY**2, NUM_CITY**2), dtype=np.int32)
for t_u_v in itertools.product(range(NUM_CITY), repeat=3):
    t, u, v = t_u_v[0], t_u_v[1], t_u_v[2]
    idx_i = NUM_CITY * t + u
    if t < NUM_CITY - 1:
        idx_j = NUM_CITY * (t + 1) + v
    elif t == NUM_CITY - 1:
        idx_j = v
    qubo_obj[idx_i, idx_j] += distance_mtx[u, v]
qubo_obj = np.triu(qubo_obj) + np.tril(qubo_obj).T - np.diag(np.diag(qubo_obj))

# <<< Constraint term >>>
qubo_constraint = np.zeros((NUM_CITY**2, NUM_CITY**2), dtype=np.int32)
# Calculate constraint term1 : 1-hot of horizontal line
for t in range(NUM_CITY):
    for u in range(NUM_CITY - 1):
        for v in range(u + 1, NUM_CITY):
            qubo_constraint[NUM_CITY*t+u, NUM_CITY*t+v] += ALPHA * 2
# Linear term
for t_u in itertools.product(range(NUM_CITY), repeat=2):
    qubo_constraint[NUM_CITY*t_u[0]+t_u[1], NUM_CITY*t_u[0]+t_u[1]] += ALPHA * (-1)
const_constraint = ALPHA * NUM_CITY

# Calculate constraint term2 : 1-hot of vertical line
# Quadratic term
for u in range(NUM_CITY):
    for t1 in range(NUM_CITY - 1):
        for t2 in range(t1+1, NUM_CITY):
            qubo_constraint[NUM_CITY*t1+u, NUM_CITY*t2+u] += ALPHA * 2
# Linear term
for u_t in itertools.product(range(NUM_CITY), repeat=2):
    qubo_constraint[NUM_CITY*u_t[1]+u_t[0], NUM_CITY*u_t[1]+u_t[0]] += ALPHA * (-1)
const_constraint += ALPHA * NUM_CITY
  1. 创建instance池

@dataclass
class Solution():
    """
    Solution information.

    Attributes:
        x (np.ndarray): n-sized solution composed of binary variables
        energy_all (float): energy value obtained from QUBO-matrix of all term
        energy_obj (float): energy value obtained from QUBO-matrix of objective term
        energy_constraint (float): energy value obtained from QUBO-matrix of constraint term
        constraint (bool): flag whether the solution satisfies the given constraint
    """
    x: np.ndarray
    energy_all: float = 0
    energy_obj: float = 0
    energy_constraint: float = 0
    constraint: bool = True

    @classmethod
    def energy(cls, qubo:np.ndarray, x: np.ndarray, const=0) -> float:
        """
        Calculate the enrgy from the QUBO-matrix & solution x

        Args:
            qubo (np.ndarray): n-by-n QUBO-matrix
            x (np.ndarray): n-sized solution composed of binary variables
            const (int, optional): _description_. Defaults to 0.

        Returns:
            float: Energy value.
        """
        return float(np.dot(np.dot(x, qubo), x) + const)
    @classmethod
    def check_constraint(cls, qubo: np.ndarray, x: np.ndarray, const=0) -> bool:
        """
        Check whether the solution satisfies the constraints.

        Args:
            qubo (np.ndarray): QUBO-model of the constraint term.
            x (np.ndarray): solution that you want to check.
            const (int, optional): constant of the constraint term. Defaults to 0.

        Returns:
            bool: Return True if the solution satisfy.
                  Return False otherwise.
        """
        return True if cls.energy(qubo, x, const) == 0 else False

  1. subQUBO Hybrid Annealing Algorithm
# https://ieeexplore.ieee.org/document/9664360

# <<< Line 2-4 >>>
# Initialize the Instance Pool
pool = []
for i in range(N_I):
    # ====================
    # 实验时改动此参数
    x = np.random.randint(0, 2, num_spin) # 生成随机解
    # ====================
    energy_obj = Solution.energy(qubo_obj, x)
    energy_constraint = Solution.energy(qubo=qubo_constraint, x=x, const=const_constraint)
    pool.append(
        Solution(
            x = x,
            energy_all = energy_obj + energy_constraint,
            energy_obj = energy_obj,
            energy_constraint = energy_constraint,
            constraint = Solution.check_constraint(qubo=qubo_constraint, x=x, const=const_constraint)
        )
    )
ascending_order_idx = np.argsort(np.array(list(map(lambda sol: sol.energy_all, pool))))
pool = [pool[i] for i in ascending_order_idx]

# <<< Line 5 >>>
# Find the best solution
ascending_order_idx = np.argsort(np.array(list(map(lambda sol: sol.energy_all, pool))))
x_best = pool[ascending_order_idx[0]]

for _ in range(1): # <<< Line 6 >>>
    # <<< Line 7-8 >>>
    # Obtain a quasi-ground-state solution for every N_I solution instance by a classical QUBO solver
    for solution_i in pool:
        # ====================
        # 实验时改动此参数
        x = np.random.randint(0, 2, num_spin) # 生成随机解
        # ====================

        # Update the solution info
        solution_i.x = x
        energy_obj = solution_i.energy(qubo_obj, x)
        energy_constraint = solution_i.energy(qubo_constraint, x, const_constraint)
        solution_i.energy_all = energy_obj + energy_constraint
        solution_i.energy_obj = energy_obj
        solution_i.energy_constraint = energy_constraint
        solution_i.constraint = solution_i.check_constraint(qubo=qubo_constraint, x=x, const=const_constraint)

    for i in range(N_E): # <<< Line 9 >>>
        # <<< Line 10 >>>
        # Select N_S solution instance randomly from the pool
        n_s_pool = random.sample(pool, N_S)

        # <<< Line 11-14 >>>
        # Calculate variance of each spin x_i in N_S instance poolSolution.check_constraint(qubo_constraint, x, const_constraint)
        vars_of_x = np.array([sum(n_s_pool[k].x[j] for k in range(N_S)) - N_S/2 for j in range(num_spin)])

        # <<< Line 15 >>>
        # Select a solution randomly from N_S solution instance pool as a tentative solution
        solution_tmp = random.choice(n_s_pool)

        # Extract a subQUBO
        extracted_spin_idx = np.argsort(vars_of_x)[:sub_qubo_size]
        non_extracted_spin_idx = np.argsort(vars_of_x)[sub_qubo_size:]
        subqubo_obj = np.array([[qubo_obj[j, k] for k in extracted_spin_idx] for j in extracted_spin_idx])
        subqubo_constraint = np.array([[qubo_constraint[j, k] for k in extracted_spin_idx] for j in extracted_spin_idx])
        for idx_i in range(sub_qubo_size):
            subqubo_obj[idx_i, idx_i] += sum(qubo_obj[idx_i, idx_j] * solution_tmp.x[idx_j] for idx_j in non_extracted_spin_idx)
            subqubo_constraint[idx_i, idx_i] += sum(qubo_constraint[idx_i, idx_j] * solution_tmp.x[idx_j] for idx_j in non_extracted_spin_idx)

        # <<< Line 16 >>>
        # Optimize the subQUBO using an Ising machine
        # ====================
        # 实验时改动此参数
        x_sub = np.random.randint(0, 2, sub_qubo_size) # 生成随机解
        # ====================
        # Combine the quasi-ground-state solution from the subQUBO with the tentative solution X_t(solution_tmp)
        for idx, val in enumerate(extracted_spin_idx):
            solution_tmp.x[idx] = x_sub[idx]

        # <<< Line 17 >>>
        # Add the solution into the pool
        pool.append(solution_tmp)

    # <<< Line 18 >>>
    # Find the best soliution
    ascending_order_idx = np.argsort(np.array(list(map(lambda sol: sol.energy_all, pool))))
    x_best = pool[ascending_order_idx[0]]

    # <<< Line 19 >>>
    # Arrange the N_I instance pool
    sorted_pool = [pool[i] for i in ascending_order_idx]
    pool = sorted_pool[:N_I]

pool, x_best

总结

subQUBO思路很简单,希望大家可以看着代码,理解如果实现。这个算法已经被早稻田大学申请专利了。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/540025.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

哪些因素影响阻抗控制?网格铜的妙用

原文来自微信公众号&#xff1a;工程师看海&#xff0c;与我联系&#xff1a;chunhou0820 看海原创视频教程&#xff1a;《运放秘籍》 大家好&#xff0c;我是工程师看海&#xff0c;原创文章欢迎点赞分享&#xff01; 前文介绍了传输线、特性阻抗以及信号的反射概念&#xff…

【漏洞复现】用友 NC PaWfm SQL注入漏洞

0x01 产品简介 用友NC是用友网络科技股份有限公司开发的一款大型企业数字化平台。它主要用于企业的财务核算、成本管理、资金管理、固定资产管理、应收应付管理等方面的工作&#xff0c;致力于帮助企业建立科学的财务管理体系&#xff0c;提高财务核算的准确性和效率。 0x02 …

在线批量生成URL HTML单页网页程序

输入前缀、开始数字、结束数字、后缀 即可快速生成 几万、十万、百万 条链接。 支持 一键复制、 一键导出本地 txt 文件。 源码免费下载地址抄笔记 (chaobiji.cn)

Spring MVC应用分层(三层架构)

该片文章主要是对 Spring MVC应用分层&#xff08;三层架构&#xff09;进行简单的介绍和学习。 一、介绍 1、什么是应用分层 应用分层 是一种 软件开发设计思想 , 它将应用程序分成N个层次, 这N个层次分别负责各自的职责, 多个 层次之间协同提供完整的功能. 根据项目的复杂…

正则表达式:量词(三)

正则表达式中的量词有以下几种:1. *: 匹配前面的字符0次或多次。2. : 匹配前面的字符1次或多次。3.?: 匹配前面的字符0次或1次。4. {n}: 匹配前面的字符恰好n次。5. {n,}: 匹配前面的字符至少n次。6. {n,m}:匹配前面的字符至少n次&#xff0c;但不超过m次。 以下是使用Python的…

0.1 + 0.2 不等于 0.3 ?这是为什么?一篇讲清楚!!!

0.1 0.2 不等于 0.3 &#xff1f;这是为什么&#xff1f;一篇讲清楚&#xff01;&#xff01;&#xff01; 分类 编程技术 在很多编程语言中&#xff0c;我们都会发现一个奇怪的现象&#xff0c;就是计算 0.1 0.2&#xff0c;它得到的结果并不是 0.3&#xff0c;比如 C、C、…

C语言--结构体大小

基本数据类型占用的字节数分别为:char(1),short(2),int(4),long(4),long long(8),float(4),double(8)。 分析一下下面结构体占用的字节数。 struct A { int a; }; struct B { char a; int b; }; int main() { printf("sizeof(struct A)%d\n", sizeof(struct A));//测…

云计算:OVS 集群 使用 Geneve 流表

目录 一、实验 1.环境 2.OVS 集群 使用 Geneve 流表 二、问题 1.VXLAN与Geneve区别 一、实验 1.环境 (1) 主机 表1 宿主机 主机架构软件IP网卡备注ovs_controller控制端 karaf 0.7.3 192.168.204.63 1个NAT网卡 &#xff08;204网段&#xff09; 已部署ovs_server01服务…

推荐一款轻量级的hosts文件编辑器(免安装版)

在管理和编辑hosts文件时&#xff0c;一款简单而有效的工具是非常重要的。下面推荐一款免安装版的轻量级hosts文件编辑器&#xff0c;让你轻松管理你的hosts文件。 windows系统默认hosts文件位置 下载地址&#xff1a;https://www.alipan.com/s/8kSns9eAi9f

Day23_学点儿Java_多态复习

1 做错的选择题 Java中的多态性是通过以下哪个机制实现的&#xff1f;&#xff08;&#xff09; A. 方法重载 B. 方法覆盖 C. 抽象类 D. 接口2 多态复习 2.1 学点儿Java_Day7_继承、重载、重写、多态、抽象类 2.2 面向对象四大基本特征 封装、抽象、继承、多态 封装 面向…

视频号小店究竟有什么秘密,值得商家疯狂入驻,商家必看!

大家好&#xff0c;我是电商花花。 我们都知道视频号和抖音本身都是一个短视频平台&#xff0c;但是随着直播电商的发展&#xff0c;背后的流量推动逐步显露出强大的红利市场和变现机会。 视频号小店流量大和赚钱之外&#xff0c;还非常适合普通人创业。 这也使得越来越多的…

久菜盒子|可实现需求|Advanced Economic Evaluation(作业,资源免费获取)

数据资源下载地址&#xff1a;https://download.csdn.net/download/weixin_68126662/89101333 R代码参考&#xff1a; ############################################################################# ### ### GLBH0046: Advanced Economic Evaluation ### ### Practical 1…

BUUCTF刷题十一道(12)附-SSTI专题二

文章目录 学习文章[CISCN2019 华东南赛区]Web11【存疑】[RootersCTF2019]I_<3_Flask 学习文章 SSTI-服务端模板注入漏洞 flask之ssti模板注入从零到入门 CTFSHOW SSTI篇-yu22x SSTI模板注入绕过&#xff08;进阶篇&#xff09;-yu22x SSTI模板注入学习-竹言笙熙 全部总结看最…

六:ReentrantLock —— 可重入锁

目录 1、ReentrantLock 入门2、ReentrantLock 结构2.1、构造方法&#xff1a;默认为非公平锁2.2、三大内部类2.2、lock()&#xff1a;加锁【不可中断锁】2.2.1、acquire() 方法 —— AQS【模板方法】2.2.2.1 tryAcquire() 方法 —— AQS&#xff0c;由子类去实现2.2.2.2. addWa…

引人共鸣的情感视频素材在哪找?今天看这五个网站

朋友们好啊&#xff0c;最近是不是不少人都在发愁啊&#xff0c;优秀创作者做视频用的视频素材哪来的啊&#xff1f;今天我为朋友们准备了几个优秀的视频素材网站&#xff0c;让你们做视频不再缺少素材&#xff0c;然后还有几个辅助创作的工具&#xff0c;都是你们需要的&#…

了解虚拟路由器冗余协议(VRRP)

虚拟路由器冗余协议&#xff08;VRRP&#xff09;是一种被广泛使用的网络协议&#xff0c;旨在增强网络的可靠性和可用性。对于网络管理员和工程师来说&#xff0c;了解VRRP是确保网络能够实现无缝故障转移和保持不间断连接的关键。本文将深入探讨VRRP的基础知识&#xff0c;包…

贪心算法|968.监控二叉树

力扣题目链接 class Solution { private:int result;int traversal(TreeNode* cur) {// 空节点&#xff0c;该节点有覆盖if (cur NULL) return 2;int left traversal(cur->left); // 左int right traversal(cur->right); // 右// 情况1// 左右节点都有覆盖if (le…

GPT与Python结合应用于遥感降水数据处理、ERA5大气再分析数据的统计分析、干旱监测及风能和太阳能资源评估等大气科学关键场景

如何结合最新AI模型与Python技术处理和分析气候数据。介绍包括GPT-4等先进AI工具&#xff0c;旨在帮助大家掌握这些工具的功能及应用范围。 内容覆盖使用GPT处理数据、生成论文摘要、文献综述、技术方法分析等实战案例&#xff0c;使大家能够将AI技术广泛应用于科研工作。特别关…

摩天大楼为什么建不成

小小学校让搞什么生活中的数学&#xff0c;推荐主题各种高大上&#xff0c;而我独爱简单&#xff0c;昨天讲了大概&#xff0c;仅从电梯开销说明摩天大楼为什么不能无限高&#xff0c;今天作文记下。不过最终&#xff0c;我还是没有选择这个题目&#xff0c;而是帮助小小讲区块…

【Git教程】(十)版本库之间的依赖 —— 项目与子模块之间的依赖、与子树之间的依赖 ~

Git教程 版本库之间的依赖 1️⃣ 与子模块之间的依赖2️⃣ 与子树之间的依赖&#x1f33e; 总结 在 Git 中&#xff0c;版本库是发行单位&#xff0c;代表的是一个版本&#xff0c;而分支或标签则只能被创建在版本库这个整体中。如果一个项目中包含了若干个子项目&#xff0c;…