simplex法(单纯形) 并在python实现简单的应用

我最近在学习动物遗传育种GBLUP的一个算法-dfreml

一般做gwas,这些传统软件如asreml dmu gmat(我老师写的,有需要可以与我联系) 等等,计算方差组分都是运用求导 用EM和AI迭代算法求出

有没有一种不用求导的算法计算方差组分呢,那就是非求导的df-reml

Df-reml基本包括两步

第一步即为构建似然函数

第二步为搜索未知值使似然函数极大化-simplex法

simplex法的主要步骤

1. Simplex法求似然函数的最大值

假设有p个需要估计的参数:  

                                                  

设置 p+1 组参数的初值:    

                   

通过第一步求取的F即可算出p+1组值即:

                                      

假设未知参数的个数为2,即设置3组初值,通过带入F 得 F0,F1, F2

                              设置三组初值t0,t1,t2

 

 以上即为simplex的一个迭代的过程

说的通俗点就是有两个未知的变量,设置三组两个未知变量的初值,将每个初值往设定好的函数里带入并且计算,每一次迭代舍弃一组数,并用保留的几组数产生新的一组数,经过多次迭代即可得到最为希望想要的结果

由于专业原因,本文跳过第一步,直接介绍simplex法在df-reml的应用

                          参考了张勤教授编写的<<动物遗传育种中的计算方法>>动物遗传育种中的计算方法本书全面系统的介绍了在动物遗传育种中常用的主要计算方法，内容包括混合模型方程组得相关计算基数、Monte Carlo 方法、MCMC算法三部分。https://book.sciencereading.cn/shop/book/Booksimple/show.do?id=BF8D8CC9939C441D39698611F5B7F7F7A000

为了验证simplex法,还是需要构建一个函数

假设我构建的函数为两个自变量即

def lu(x, y):
    f = pow(x, 2) + pow(y, 2)
    return f

这个函数有一个最小值即为(0,0)

为了实现simplex法,我认为首先需要对给予的初值进行一个编号排序,以下函数可以实现它

def newt(t0, t1, t2):
    f0 = lu(t0[0], t0[1])
    f1 = lu(t1[0], t1[1])
    f2 = lu(t2[0], t2[1])
    fmax = max(f0, f1, f2)
    fmin = min(f0, f1, f2)
    if f0 == fmax:
        ret2 = t0
        if f1 == fmin:
            ret0 = t1
            ret1 = t2
            fmid = f2
        elif f2 == fmin:
            ret0 = t2
            ret1 = t1
            fmid = f1
    elif f1 == fmax:
        ret2 = t1
        if f0 == fmin:
            ret0 = t0
            ret1 = t2
            fmid = f2
        elif f2 == fmin:
            ret0 = t2
            ret1 = t0
            fmid = f0
    elif f2 == fmax:
        ret2 = t2
        if f0 == fmin:
            ret0 = t0
            ret1 = t1
            fmid = f1
        elif f1 == fmin:
            ret0 = t1
            ret1 = t0
            fmid = f0
    var = np.var(np.array([f1, f2, f0]))
    return ret0, ret1, ret2, fmin, fmid, fmax, var

这个函数返回的元组从左到右即为排序的初值,初值带入函数计算的值,以及计算结果的方差

即  fmin = f(ret0)  fmid = f(ret1)

完整的python代码:

import numpy as np
import datetime


star = datetime.datetime.now()


def lu(x, y):
    f = pow(x, 2) + pow(y, 2)
    return f


np.random.seed(1)
wg = np.random.random((3, 2))


def newt(t0, t1, t2):
    f0 = lu(t0[0], t0[1])
    f1 = lu(t1[0], t1[1])
    f2 = lu(t2[0], t2[1])
    fmax = max(f0, f1, f2)
    fmin = min(f0, f1, f2)
    if f0 == fmax:
        ret2 = t0
        if f1 == fmin:
            ret0 = t1
            ret1 = t2
            fmid = f2
        elif f2 == fmin:
            ret0 = t2
            ret1 = t1
            fmid = f1
    elif f1 == fmax:
        ret2 = t1
        if f0 == fmin:
            ret0 = t0
            ret1 = t2
            fmid = f2
        elif f2 == fmin:
            ret0 = t2
            ret1 = t0
            fmid = f0
    elif f2 == fmax:
        ret2 = t2
        if f0 == fmin:
            ret0 = t0
            ret1 = t1
            fmid = f1
        elif f1 == fmin:
            ret0 = t1
            ret1 = t0
            fmid = f0
    var = np.var(np.array([f1, f2, f0]))
    return ret0, ret1, ret2, fmin, fmid, fmax, var


for item in range(500):
    t_0 = wg[0, :]
    t_1 = wg[1, :]
    t_2 = wg[2, :]
    xuz = newt(t_0, t_1, t_2)
    if xuz[-1] <= 1.0e-8:
        print("grid search done")
        break
    tm = (xuz[0] + xuz[1]) / 2
    tr = tm + 0.5 * (tm - xuz[2])
    fr = lu(tr[0], tr[1])
    if fr >= xuz[5]:
        print("overextended")
        tc = tm + 0.5 * (xuz[2] - tm)
        fc = lu(tc[0], tc[1])
        if fc < min(fr, xuz[5]):
            print("The extension was overdone, but the contraction was successful")
            tn = tc
        else:
            print("#########################################################err#######################################")
            wg = np.array((xuz[0], (xuz[0] + xuz[1]) / 2, (xuz[0] + xuz[2]) / 2))
            continue
    else:
        if xuz[4] <= fr < xuz[5]:
            print("extension success")
            tn = tm + 0.5 * (tr - tm)
        elif xuz[3] <= fr < xuz[4]:
            tn = tr
        else:
            print("Extension and its success")
            te = tm + 2 * (tr - tm)
            fe = lu(te[0], te[1])
            if fe < fr:
                tn = te
            else:
                tn = tr
    wg = np.array((tn, xuz[0], xuz[1]))
    print(wg)

整体的代码思路是给定函数规则,随机出现3×2的初值,迭代搜索(每一次的迭代需要重新排序)

代码迭代终止的条件是迭代的过程中计算三组数 F0 F1 F2 的方差小于10的负8次方

经过运算的结果为

大致能够计算出取极点的解

将函数更换为

                               

 即

def lu(x, y):
    f = pow(pow(x - 1, 2)-1, 2) + pow(pow(y + 1, 2)-1, 2)
    return f

 得出的结果为

                                    是不完整的解

由此可以反映出一些问题,函数的构建时如果函数是一个无限递增的函数,再使用simplex求取最大值出现时的未知数的值,这可能会出现问题,可能本人水平有限未能找出解决的办法,再一个如果函数的取极值的点有多个,最后的结果可能也会出现问题

本文只是一个simplex法的最基本的使用方法,仅仅分享一下学习的过程,我感觉这是一个非求导的较为良好的方法,提供了一个很好的想法和思路