如何优化包含“k-最大”运算符的 objective 函数？
How to optimize objective function that contains the “k-biggest” operator?
我想解决以下优化问题：
min_{x,y} f(xA+yB)
s.t: {x,y} >= 0 , (x+y)=1
其中，{A,B}是方阵，{x,y}是标量。
函数f(A)=sum(s_i)，而s_i表示最大的和[=31= A 的第 i 行中的 ]k 个条目。
例如如果
A= [1 2 3
5 1 7
3 4 2]
然后 f(A) 对于 k=2 将是
(2+3)+(5+7)+(4+3)=24
如果我能计算 f 关于 (x,y) 的导数，那么问题就很容易解决了。
我还知道在 objective 中处理“max(a)”的一种方法是添加堆栈变量 t 到 objective 并添加 a<=t 到约束。但是我找不到上述问题那么简单。
为了比较解的最优性，我使用Matlab通用求解器解决了上述问题，但我需要知道如何自己优化它。
（这个问题有点不清楚。我不知道你是想学习如何将其表述为凸优化问题（作业？）还是以某种方式解决它。你没有提供你的工作方法和很多缺少更多。请在下一个问题中记住这一点)
使用强大的建模语言进行凸优化
cvxpy (python)（可能还有 cvxopt (matlab) 和 convex.jl (julia)）支持这种约束。所以整个问题看起来像：
代码：
from cvxpy import *
import numpy as np
np.random.seed(1)
# Constants
A = np.random.normal(size=(5,5))
B = np.random.normal(size=(5,5))
k = 2
# Vars
x = Variable(2)
# Constraints
constraints = [sum_entries(x) == 1,
x >= 0.0]
# Problem
func_inner = x[0] * A + x[1] * B
func = sum([sum_largest(func_inner[i, :], 2) for i in range(A.shape[1])])
objective = Minimize(func)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
这些工具还通过构造来证明凸性，这对于我们以后可能使用的其他方法（例如投影梯度下降）非常重要。这是一个凸优化问题！
手动凸优化
让我们关注 sum_largest(x, k) 的公式，x 是变量向量，k 是正常数。
这可以通过引入一些辅助变量来制定。让我们修复 x 并重新制定：
sum_largest([1, 0, 3], 2):
x = [1, 0, 3] # input
q = Variable(1) # aux-var: single-dim var
t = Variable(vec_len) # aux-vars: dim = 3
Min sum(t) + k * q
st. x[0] <= t[0] + q,
x[1] <= t[1] + q,
x[2] <= t[2] + q,
t[0] >= 0,
t[1] >= 0,
t[2] >= 0
此重新表述取自 cvxpy's sources!
测试（就 cvxpy 的功能而言不是最好的代码）：
from cvxpy import *
import numpy as np
# Constants
x = np.array([-1, -2, -0.5, 1, 0, 3])
k = 4
# Vars
q = Variable()
t = Variable(x.shape[0])
# Constraints
constraints = [x[0] <= t[0] + q,
x[1] <= t[1] + q,
x[2] <= t[2] + q,
x[3] <= t[3] + q,
x[4] <= t[4] + q,
x[5] <= t[5] + q,
t >= 0.0]
objective = Minimize(sum(t) + k * q)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
print(problem.value)
输出：
3.4999999942200737 (IPM-based solver; therefore only approximating optimum)
投影梯度下降
可能会想到使用一些基于梯度下降的方法。在这种情况下，需要选择使用哪种梯度计算：
- 手动梯度计算
- 自动微分
- 数值微分
我不会详细介绍（请参阅下一段了解原因），但由于缺少对多维排序的梯度支持（我在函数中使用了它），我使用 autograd 的尝试失败了-def).
但是： 投影梯度下降效果很好（并且比所有那些针对大规模问题的凸优化求解器快得多），即使使用数值微分（调用功能很多！）。我测试了一些非单调的 PGD 方法，效果很好。
人们可以使用对概率单纯形的有效线性时间投影（Duchi 等人，“用于高维学习的 l1-Ball 的有效投影”）！请记住：这是一个受约束的优化问题并且纯GD 将不起作用（不重新制定）。
省略代码。
就解决了？要聪明！使用标量最小化！
你刚得到 2 variables x and y，都在 [0,1] 中，约束条件 x + y == 1。
这两个变量合并为一个，因此我们可以使用 1d 优化/标量优化，这更简单！
而不是优化x和y，我们只是优化x和y=1.0 -x是隐含的！
当然可以使用有界方法！
代码：
import numpy as np
from scipy.optimize import minimize_scalar
from time import perf_counter
np.random.seed(1)
# Constants
# ---------
k = 2
N = 500
A = np.random.random(size=(N,N))
A_zeros = np.random.choice([0, 1], size=(N, N)) # more interesting data
A *= A_zeros # """
B = np.random.random(size=(N,N)) # """
B_zeros = np.random.choice([0, 1], size=(N, N)) # """
B *= B_zeros # """
# Function
# --------
# remark: partition / incomplete sort might be better for large-scale
def func(x):
return np.sum(np.sort(x*A + (1.-x)*B, axis=1)[:, -k:]) # y=1.0-x !!!
start = perf_counter()
res = minimize_scalar(func, method='bounded', bounds=(0., 1.))
end = perf_counter()
print(res)
print('used: ', end-start)
输出：
fun: 930.68327527871179
message: 'Solution found.'
nfev: 20
status: 0
success: True
x: 0.50511029941339591
used: 0.20979814249560497
请记住，所有呈现的内容都使用了一个事实，即它是一个凸问题，其中每个局部最小值都是全局最小值，cvxpy 证明了这一点！
我想解决以下优化问题：
min_{x,y} f(xA+yB)
s.t: {x,y} >= 0 , (x+y)=1
其中，{A,B}是方阵，{x,y}是标量。
函数f(A)=sum(s_i)，而s_i表示最大的和[=31= A 的第 i 行中的 ]k 个条目。
例如如果
A= [1 2 3
5 1 7
3 4 2]
然后 f(A) 对于 k=2 将是
(2+3)+(5+7)+(4+3)=24
如果我能计算 f 关于 (x,y) 的导数，那么问题就很容易解决了。
我还知道在 objective 中处理“max(a)”的一种方法是添加堆栈变量 t 到 objective 并添加 a<=t 到约束。但是我找不到上述问题那么简单。
为了比较解的最优性，我使用Matlab通用求解器解决了上述问题，但我需要知道如何自己优化它。
（这个问题有点不清楚。我不知道你是想学习如何将其表述为凸优化问题（作业？）还是以某种方式解决它。你没有提供你的工作方法和很多缺少更多。请在下一个问题中记住这一点)
使用强大的建模语言进行凸优化
cvxpy (python)（可能还有 cvxopt (matlab) 和 convex.jl (julia)）支持这种约束。所以整个问题看起来像：
代码：
from cvxpy import *
import numpy as np
np.random.seed(1)
# Constants
A = np.random.normal(size=(5,5))
B = np.random.normal(size=(5,5))
k = 2
# Vars
x = Variable(2)
# Constraints
constraints = [sum_entries(x) == 1,
x >= 0.0]
# Problem
func_inner = x[0] * A + x[1] * B
func = sum([sum_largest(func_inner[i, :], 2) for i in range(A.shape[1])])
objective = Minimize(func)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
这些工具还通过构造来证明凸性，这对于我们以后可能使用的其他方法（例如投影梯度下降）非常重要。这是一个凸优化问题！
手动凸优化
让我们关注 sum_largest(x, k) 的公式，x 是变量向量，k 是正常数。
这可以通过引入一些辅助变量来制定。让我们修复 x 并重新制定：
sum_largest([1, 0, 3], 2):
x = [1, 0, 3] # input
q = Variable(1) # aux-var: single-dim var
t = Variable(vec_len) # aux-vars: dim = 3
Min sum(t) + k * q
st. x[0] <= t[0] + q,
x[1] <= t[1] + q,
x[2] <= t[2] + q,
t[0] >= 0,
t[1] >= 0,
t[2] >= 0
此重新表述取自 cvxpy's sources!
测试（就 cvxpy 的功能而言不是最好的代码）：
from cvxpy import *
import numpy as np
# Constants
x = np.array([-1, -2, -0.5, 1, 0, 3])
k = 4
# Vars
q = Variable()
t = Variable(x.shape[0])
# Constraints
constraints = [x[0] <= t[0] + q,
x[1] <= t[1] + q,
x[2] <= t[2] + q,
x[3] <= t[3] + q,
x[4] <= t[4] + q,
x[5] <= t[5] + q,
t >= 0.0]
objective = Minimize(sum(t) + k * q)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
print(problem.value)
输出：
3.4999999942200737 (IPM-based solver; therefore only approximating optimum)
投影梯度下降
可能会想到使用一些基于梯度下降的方法。在这种情况下，需要选择使用哪种梯度计算：
- 手动梯度计算
- 自动微分
- 数值微分
我不会详细介绍（请参阅下一段了解原因），但由于缺少对多维排序的梯度支持（我在函数中使用了它），我使用 autograd 的尝试失败了-def).
但是： 投影梯度下降效果很好（并且比所有那些针对大规模问题的凸优化求解器快得多），即使使用数值微分（调用功能很多！）。我测试了一些非单调的 PGD 方法，效果很好。
人们可以使用对概率单纯形的有效线性时间投影（Duchi 等人，“用于高维学习的 l1-Ball 的有效投影”）！请记住：这是一个受约束的优化问题并且纯GD 将不起作用（不重新制定）。
省略代码。
就解决了？要聪明！使用标量最小化！
你刚得到 2 variables x and y，都在 [0,1] 中，约束条件 x + y == 1。
这两个变量合并为一个，因此我们可以使用 1d 优化/标量优化，这更简单！
而不是优化x和y，我们只是优化x和y=1.0 -x是隐含的！
当然可以使用有界方法！
代码：
import numpy as np
from scipy.optimize import minimize_scalar
from time import perf_counter
np.random.seed(1)
# Constants
# ---------
k = 2
N = 500
A = np.random.random(size=(N,N))
A_zeros = np.random.choice([0, 1], size=(N, N)) # more interesting data
A *= A_zeros # """
B = np.random.random(size=(N,N)) # """
B_zeros = np.random.choice([0, 1], size=(N, N)) # """
B *= B_zeros # """
# Function
# --------
# remark: partition / incomplete sort might be better for large-scale
def func(x):
return np.sum(np.sort(x*A + (1.-x)*B, axis=1)[:, -k:]) # y=1.0-x !!!
start = perf_counter()
res = minimize_scalar(func, method='bounded', bounds=(0., 1.))
end = perf_counter()
print(res)
print('used: ', end-start)
输出：
fun: 930.68327527871179
message: 'Solution found.'
nfev: 20
status: 0
success: True
x: 0.50511029941339591
used: 0.20979814249560497
请记住，所有呈现的内容都使用了一个事实，即它是一个凸问题，其中每个局部最小值都是全局最小值，cvxpy 证明了这一点！