目录
Background
Analysis
1. python 符号运算
2. python嵌套优化问题
Reference
Background
信息论project， 求LDPC的variable nodes和check nodes的degrees，也即λ和ρ，使得code rate最大， 抛开背景，就是个优化问题，形式如下：
是给定常数，假设是自己设定好的（这里设置为, 即），只需要求.
Analysis
这个优化问题的第一条constraint非常特别，首先以一个多项式带入另一个多项式的x，自己算是没法算的，阶数非常高，因此要先用符号运算求出constraint的表达式，并且要想在一个区间内都满足，必须转换成在这个区间内的最小值大于等于0的一个求最小值的优化问题；
其次它在x的连续的区间内都要满足，引入了新的变量x，容易想到嵌套优化问题的思路，但直接scipy.optimize.minimize这个库会因为额外变量x而报错（具体原因我就不写了，试过就知道）， google了很久都没搜到，最后想出了一个解决方案，用优化问题定义另一个优化问题以作为前者的constraint。
1. python 符号运算
sympy库 用来符号运算，定义好符号后，用这些符号作为多项式的系数，整个多项式就可以做运算了。
这里用scipy.poly1d 生成多项式，系数为p1,p2,p3,p4（scipy.poly1d就是numpy.poly1d），但用poly1d库生成的多项式虽然可以带入另一个多项式的x中，但不能print出表达式，看不到运算后的结果，因此这里选择用符号加法运算自己定义多项式。
from sympy import symbols
from scipy import poly1d
# 定义符号
x, p1, p2, p3 = symbols('x p1 p2 p3')
# 一次定义多个符号,返回tuple
vrs = symbols('p:5') # (p0, p1, p2, p3, p4)
vrs = symbols('p1:5') # (p1, p2, p3, p4)
# 用poly1d生成多项式，可拥有多项式任何属性，如获取系数、获取x为某个值时的结果等，
# 虽然可以带入另一个多项式的x，但不能print(G)
ρ = ploy1d(list(vrs))
print(ρ) # p0*x**4 + p1*x**3 + p2*x**2 + p3*x + p4
print(ρ.c) # 获取系数
print(ρ(1)) # x = 1时多项式的系数
# 自己定义多项式， sum运用了符号的加法运算来生成多项式
def polyn(x,vrs):
n = len(vrs)
result = 0 + sum([v * x**(n - i - 1) for i, v in enumerate(vrs)])
return result
e = 0.4 # constant
l = [0, 0, 0, 1, 0, 0] # 假设λ多项式的系数为l3 = 1, 其他系数为0
λ = poly1d(l) # x**2
ρ = polyn(x, vrs)
G = 1 - e * λ(1 - ρ)
print(G) # 1 - 0.4*(-p0*x**4 - p1*x**3 - p2*x**2 - p3*x - p4 + 1)**2
2. python嵌套优化问题
scipy.optimize库中，求解多元优化问题，可用fmin；如果带定义域（局部最小值），可用fminbound，只需设定x的取值区间；如果还要有constraints，就要用minimize。
这三个方法的优化目标函数都只能输入一个list作为唯一输入，这个list包含多个变量，即目标函数f(x), x[0], x[1], x[2]等都是变量，因此可以是多元的。
但是，如果一个优化问题的约束条件也是一个优化问题（如求某个函数的最小值，这个函数由原来优化的变量来确定），一定不要用fmin/fminbound这些封装好的库，因为没法重写源代码。
可以考虑用一个优化问题定义另一个优化问题，注意内层优化问题的返回值要取优化目标函数的值res.fun。
minimize的优化方法的选择有很多，可能迭代的速度有所不同，这里选了'SLSQP', 其他方法也都可以试一下。'SLSQP'的官方文档介绍：
Method SLSQP uses Sequential Least SQuares Programming to minimize a function of several variables with any combination of bounds, equality and inequality constraints.The method wraps the SLSQP Optimization subroutine originally implemented by Dieter Kraft [12]. Note that the wrapper handles infinite values in bounds by converting them into large floating values.
import numpy as np
from scipy.optimize import fminbound, minimize
# objective function: f
def f(l):
n = len(l)
l = [l[i] / (n - i + 1) for i in range(n)]
return sum(l)
# inner_cons
def cons(l):
# 用outer_minimize的变量l作为系数定义inner_minimize
# 这里的多项式就是对应上面符号运算的结果，只把p换成l
# G=1 - 0.4*(-p0*x**4 - p1*x**3 - p2*x**2 - p3*x - p4 + 1)**2
funbnd = lambda x: 1 - 0.4*(-l[0]*x**3 - l[1]*x**2 - l[2]*x - l[3] + 1)**2
bnds2 = ((1 - e, 1),)
# optimize f under cons
res = minimize(funbnd, x0=(1 - e,), method='SLSQP', bounds=bnds2)
# print("cons2:", res.fun)
# print(res)
return res.fun
# outer_cons
#* ineq means to be non-negative, eq means to be zero
Cons = (
{'type': 'ineq', 'fun': cons}, # inner_cons
{'type': 'eq', 'fun': lambda l: sum(l) - 1} # l1 + l2 + ... + ln = 1
)
# initial value of variables：array, tuple or list
initial_x = np.array([0 for _ in range(4)])
# initial_x = (0, 0, 0, 0)
# initial_x = [0, 0, 0, 0]
# bounds： array, tuple or List
bnds = ((0, 1), (0, 1),(0, 1), (0,1))
# outer_ominimize f under Cons
res = minimize(f, x0=initial_x, method='SLSQP', bounds=bnds, constraints=Cons)
print(res)
result:
x :array([1,0,0,0])是优化变量的最终取值，分别对应p1=1,p2=p3=p4=0,对应的多项式是
nfev:35是迭代次数；fun:0.2是优化目标函数的最终值。
Reference
scipy.optimize.fminbound — SciPy v1.7.1 Manual
scipy.optimize.minimize — SciPy v1.7.1 Manual
用Python学数学之Sympy代数符号运算 - 知乎
Python-Scipy 求函数局部极值/最大值/最小值 — 浮云的博客
Numpy中的多项式表示及拟合_蜗牛、Gray-程序员资料 - 程序员资料
python 3.x - Scipy minimize: How to pass args to both the objective and the constraint - Stack Overflow
Sympy常见多个变量【一行代码创建】_肥宅Sean-CSDN博客
