当算法遇到物理思想——量子优化算法

作者：佚名浏览：次浏览日期：2024-05-06

很多问题都可以转化为最优化问题，传统最优化算法常常受到局部极值的约束而影响最优化效果。为了使算法跳出局部极值，许多科学家进行了改进和优化。引入量子力学的思想也是为了进一步改进已有的算法，提高其收敛速度和精度。这篇文章的内容最早是大三当助教的时候写的讨论课讲义的一篇，现在重新整理了一下顺便把代码和测试贴出来，算是把之前写的讲义全部给做完整了。转载请告知并得到允许，谢谢。

经典启发式优化算法
量子遗传算法
测试样例
代码

在最优化问题中，常常会提到几个著名的启发式算法——模拟退火算法、遗传算法。前者是由统计物理出发的启发式算法，基于正则分布构造接受拒绝算法来判断是否接受构型改变的提议。根据统计物理给出的直觉，当温度逐渐降低，系统的构型将逐渐“冻结”到最低能量状态，因此，将目标函数作为能量，进行退火即可在很大的搜索空间里找到能量最低点，也就对目标函数进行了最优化。但是同样的，我们也可以想像，当搜索空间的能量曲面十分复杂，有很多局部极小值时，如果当时温度使得系统的构型很难跳出这个小的谷，也就使得最优化到局部最优，而不是全局最优解。（由于退火温度逐渐下降、若当时很难跳出局部极小，之后也将更难跳出）如果学过计算物理的同学应该比较好理解，这和Ising模型的Monte Carlo算法非常类似，都是MCMC里的Metropolis-Hasting算法。

遗传算法则是由生物学启发的算法，将数字编码成不同进制的编码（编码进制会决定优化结果的精度），称为基因。而需要找的最优自变量往往是一个向量，因此，将各个分量分别编码之后得到的基因组合在一起，称为染色体，或者个体。在最优化过程中，我们构造一批自变量两相，将其编码后得到的个体的集合称为种群。遗传算法则是通过在种群中实施交叉、变异、选择、灾变等操作，对这一批种群进行不断的更新进化，并且进化方向与目标函数的最优化方向一致。所谓交叉，即两个染色体彼此交换一段编码片段。变异则是一个染色体的某个编码位点发生变化。选择是将种群中按照目标函数进行评测，将接近最优化目标的个体保留并复制，其他个体随机丢弃或减量。灾变则是随机生成新的个体加入种群并且随机去除种群中的已有个体，这一操作是为了避免被约束在局部最优解。所有操作都是依照事先设置好的概率规则进行。和模拟退火算法寻找能量最低点相反，遗传算法是找到目标函数的最大值。

1. 量子编码

量子遗传算法的着眼点在于将经典遗传算法的进制编码改为量子比特编码，使得在经典遗传算法中完全确定的染色体上的编码位点在量子遗传算法中成为态矢的线性组合。当我们采用二进制编码时，每个节点既可能为0，也可能为1。但是在用目标函数评测染色体的过程中，必须要采用确定的数值，那么量子化的染色体具体应该取值多少？我们引入一种新的操作，这是经典遗传算法中所没有的，也就是测量，我们对染色体上的每个编码位点都进行一次测量，之后染色体将塌缩到一个确定的值。再由这一测量结果代替染色体进行评测。

也就是现在，我们用张量积来表示具有N个编码位点的染色体 $|\\Psi\\rangle=|\\psi_1\\rangle \\otimes |\\psi_2\\rangle \\otimes \\cdots \\otimes|\\psi_N\\rangle$

测量操作就是通过均匀分布随机数来确定该位点的测量结果， $r \\sim U(0, 1), ~ \\mathrm{if}~ r < |\\alpha_k|^2, ~ \\mathrm{return}~ |0\\rangle, ~ \\mathrm{else ~ return}~ |1\\rangle$

可以想象，原本在经典遗传算法中代表确定数值的染色体，在量子遗传算法中由基矢的线性叠加，可以表示任意可能的数值。也就是说在某种意义下，量子遗传算法化“串行”为“并行”。

2. 量子门更新

相比于经典遗传算法中的交叉变异，量子遗传算法由于采用了量子编码，所以需要用量子门对染色体进行更新。和量子力学测量之后的塌缩不同，算法中的测量操作只是使染色体塌缩到确定态进行数字运算，并不改变其线性叠加的形式。一般来说，可以通过旋转变化来对染色体的某一位点进行更新，定义

$U(\ heta)=\\left(\\begin{array}{cc}\\cos \ heta & - \\sin\ heta\\\\ \\sin \ heta & \\cos \ heta \\end{array}\\right)$

旋转更新后的叠加系数为 $\\left( \\begin{array}{c}\\alpha_k' \\\\ \\beta_k'\\end{array}\\right)=U(\ heta) \\left( \\begin{array}{c}\\alpha_k \\\\ \\beta_k\\end{array}\\right)$

下面的问题就是如何确定旋转角，为了保证当前种群中的最优染色体进化，我们根据如下表格进行更新，使得叠加系数尽可能接近最优染色体。这里 $x_i$ 表示当前染色体第i个位点测量后的值， $best_i$ 表示最优染色体测量后的第i位点， $f$ 是目标函数。这里的 $0.01\\pi$ 是进化的步长，可以调节。

3. 其他量子操作

一般来说，由旋转更新就够了，因为量子编码已经对最优化引入了随机性，可以一定程度上避免系统被约束在局部极大。

（1）量子交叉：交换两个染色体的进化方向，即旋转变换。

（2）量子变异：交换某一位点的叠加系数。

（3）量子灾变：算法收敛之后，对种群实施扰动，加入新的个体，防止系统约束在局部极大。

4. 模拟退火算法的量子版本

我记得量子退火算法应该是加入了隧穿，也就缓解了上面图示中的被局部极小约束的情况。

1. 只有两个局部极大

这里用了函数 $f(x_1, x_2)=30 e^{- ((x_1-2)^2 + (x_2-2)^2)/10}+ 20 e^{- ((x_1+2)^2 + (x_2+2)^2)/5}, ~ -5 \\leq x_1, x_2 \\leq 5$ 来做测试，有两个局部极大，分别是 $(-2, -2, 20.5495), (2, 2, 30.0182)$ ，为了让轨迹进入第一个局部极大的basin，在初始化的时候叠加系数没有初始化成相等的值。可以看到，轨迹刚开始进入第一个basin，之后跳出进入第二个basin并搜索到最大值。

2. 局部极大非常多

我选了一个局部极大非常多的函数来测试， $f(x_1, x_2)=x_1 \\sin(4\\pi x_1) + x_2 \\sin(20 \\pi x_2)$ ，并且限制在 $-3 \\leq x_1 \\leq 12, ~ 4.1 \\leq x_2 \\leq 5.8$ ，下面画了函数的轮廓图，从轮廓图可以看出这个函数的局部极大多的有多么丧心病狂。蓝色点是搜索轨迹，可以看到，基本上十步左右就收敛了，后面基本上都是在附近搜索，收敛在 $(x_1, x_2, f(x_1, x_2))=(11.62559, 5.325093, 16.950273)$

具体这个函数的最大值是多少，用画轮廓图的时候的网格数组（步长0.01）取最大值， $(x_1, x_2, f(x_1, x_2))=(11.6700, 5.7300, 17.0566)$

效果还算可以吧，可能再调一下旋转步长会搜索的更准确一点。

### classqga.py
### definition of Quantum Genetic Algorithm class
### Created by Yulong Dong, current PhD student at UC Berkeley, dongyl at berkeley.edu

import numpy as np
from evaluate import EvalFitness

# ---------- Parameters
N_Var = 2 		# number of variables
LenVar = 20		# length of each variable after compiling
BoundVar = [[-3,12.1],[4.1,5.8]]	# boundaries of domain of each variable

class chrom_quantum :
    def __init__(self):
        self.prob = np.ones([N_Var*LenVar,2], dtype=float)
        self.prob /= np.sqrt(2)
        self.value = [0 for i in range(N_Var)]
        self.fitness = 0

    def collapse(self):
        self.bin = (np.random.rand(N_Var*LenVar) > self.prob[:,0]**2)*int(1)

    def bin2deg(self):
        self.collapse()
        deg_list = [int("".join(map(str, self.bin.tolist()[i*LenVar:(i+1)*LenVar])), 2) for i in range(N_Var)]
        self.value = [BoundVar[i][0] + float(deg_list[i])/(2**LenVar-1) * (BoundVar[i][1]-BoundVar[i][0]) for i in range(N_Var)]

    def QuantumGate(self, bestbin, bestfitness):
        for i in range(N_Var*LenVar):
            alpha = self.prob[i,0]
            beta = self.prob[i,1]
            binary = self.bin[i]
            best_b = bestbin[i]
            if ((binary==0)and(best_b==0)) or ((binary==1)and(best_b==1)):
                (delta, sign) = (0, 0)
            elif (binary==0)and(best_b==1)and(self.fitness<bestfitness):
                delta = 0.01*np.pi
                if alpha*beta>0: sign = 1
                elif alpha*beta<0: sign = -1
                elif alpha==0: sign = 0
                elif beta==0: sign = np.sign(np.random.randn())
            elif (binary==0)and(best_b==1)and(self.fitness>=bestfitness):
                delta = 0.01*np.pi
                if alpha*beta>0: sign = -1
                elif alpha*beta<0: sign = 1
                elif alpha==0: sign = np.sign(np.random.randn())
                elif beta==0: sign = 0
            elif (binary==1)and(best_b==0)and(self.fitness<bestfitness):
                delta = 0.01*np.pi
                if alpha*beta>0: sign = -1
                elif alpha*beta<0: sign = 1
                elif alpha==0: sign = np.sign(np.random.randn())
                elif beta==0: sign = 0
            elif (binary==1)and(best_b==0)and(self.fitness>=bestfitness):
                delta = 0.01*np.pi
                if alpha*beta>0: sign = 1
                elif alpha*beta<0: sign = -1
                elif alpha==0: sign = 0
                elif beta==0: sign = np.sign(np.random.randn())
            rotdeg = sign*delta
            gate = np.array([[np.cos(rotdeg), -np.sin(rotdeg)],[np.sin(rotdeg), np.cos(rotdeg)]])
            self.prob[count,:] = gate.dot(self.prob[count,:])

def GetFitness(obj):
    obj.bin2deg()
    obj.fitness = EvalFitness(obj.value)
    return obj.fitness

def copybin(bin):
    return (bin==1)*int(1)

### main.py
### Created by Yulong Dong, current PhD student at UC Berkeley, dongyl at berkeley.edu

from classqga import *
# ---------- Parameters
MAX_GEN = 200		# max generation
SIZE_POP = 20		# size of population
trace = open("trace.txt",'w')
#trace.write("gen\	fitness\	\	value\
")
# ---------- Initialize Population
chrom_list = [chrom_quantum() for i in range(SIZE_POP)]
fitness_list = map(GetFitness, chrom_list)
best = []		# 0: bestfitness, 1: best index, 2: best value, 3: best binary
best.append(max(fitness_list))
best.append(fitness_list.index(best[0]))
best.append(chrom_list[best[1]].value)
best.append(copybin(chrom_list[best[1]].bin))
trace.write("%d\	%f\	%f\	%f\
"%(1, best[0], best[2][0], best[2][1]))
print "1 complete of %d\
"%MAX_GEN
for i in range(1, MAX_GEN):
    # ---------- Upadate Data
    for j in range(SIZE_POP): chrom_list[j].QuantumGate(best[3], best[0])
    fitness_list = map(GetFitness, chrom_list)
    # ---------- Record new best data
    newbest = max(fitness_list)
    if newbest > best[0]:
        best[0] = newbest
        best[1] = fitness_list.index(best[0])
        best[2] = chrom_list[best[1]].value
        best[3] = copybin(chrom_list[best[1]].bin)
    trace.write("%d\	%f\	%f\	%f\
"%(i+1, best[0], best[2][0], best[2][1]))
    print "%d complete of %d\
"%(i+1, MAX_GEN)
trace.flush()
trace.close()

### evaluate.py
### Description: target function to find maximum
### Created by Yulong Dong, current PhD student at UC Berkeley, dongyl at berkeley.edu

from math import sin, pi
def EvalFitness(value):
    return sin(4*pi*value[0])*value[0]+sin(20*pi*value[1])*value[1]

平台注册入口