如何在 Python 函数最小化中加快速度以找到椭圆方程的解答案

【问题标题】：How to gain speed in Python function minimization to find solutions to Ellipsoid equation如何在 Python 函数最小化中加快速度以找到椭圆方程的解
【发布时间】：2018-07-07 09:50:09
【问题描述】：

简介

使用 Python，我想检索一组满足以下描述椭圆体的方程的解：

其中 H 是一个正定矩阵。

为了检索满足任意正定矩阵 30 维矩阵 H 和任意 y 条件的向量 x，我将最小化以下函数：

from scipy.optimize import minimize
import numpy as np

N = 30
H = np.identity(N)
y = 5

def fct(x):
    return np.abs(np.dot(x, np.dot(H, x))-y)

x0 = np.ones(N)/ N
solution = minimize(fct, x0, method='Nelder-Mead')['x']

根据我的研究目标，我想检索满足大量 (=500) 不同正定矩阵 H 条件的许多不同解决方案。因此，我将使用不同的起始值多次最小化函数每个矩阵 H 用 x0 表示。我想调整解决问题的方法，以便从问题的特殊特征中获得一些性能（在计算时间方面）。

问题分析

我对数值优化知之甚少（计划有时间时深入研究），但在阅读了一篇有趣的帖子here 和一些scipy lectures 之后，我得出以下结论：

鉴于 H 是正定矩阵，我的问题是凸优化之一。
我相信我的问题是无约束优化之一。
我不确定我的问题是好还是坏。解决方案确实会随着起始值的变化而变化，但我认为我的问题是有条件的。

问题

如果 H 是正定矩阵，我的问题是无约束凸优化问题是否正确？

Nelder-Mead 算法确实会返回“有意义”的解决方案。我可以通过利用问题的特征来加快速度吗？或者可能通过使用不同的优化算法？我知道有某些方法可以提供雅可比和粗麻布显着提高性能，但我不确定这是否适用于我的问题。

非常感谢您的帮助。

【问题讨论】：

为什么选择 Nelder-mead？只需使用 (L-)BFGS 并提供 jacobian/hessian 信息，或者仅使用 QP 求解器，您不需要任何初始点。查看实际上与公式不同的代码，您可能会考虑您的 l1 与 l2 惩罚（您正在做 l1），在后一种情况下，它也可以通过最小二乘求解器来解决。
如果H是半定的，你对H的逆定义是什么？此外，方程使用了 H 的逆，但代码创建并使用了一个称为 H 的矩阵。在方程中使用 H 的逆真的是你想要的吗？
@sascha，我最初没有使用 (L-)BFGS 的原因是它似乎只返回等坐标解，即所有值在每个维度上重合的解向量 x。我正在寻找许多解决方案，这些解决方案将作为进一步分析的输入，包括非等式解决方案。但是谢谢你的回答。我猜 Nelder-Mead 比（L-BFGS）慢一点。另外：我的代码与公式有何不同？我纠正了正定矩阵的逆。
@WarrenWeckesser 我实际上正在处理正定矩阵。我在问题中对此进行了更正！我确实需要 H 的逆，因为单位矩阵的逆是单位矩阵。

标签： python performance scipy minimization convex-optimization

【解决方案1】：

您要解决的问题是一个可以表示为 F(x, params) = 0 的方程，因此使用 root-finding function 而不是最小化函数更自然。在这种情况下，您有 30 个变量和一个方程，因此问题是高度欠定的，如果您使用求根函数，您将不得不处理这个问题。

该方程定义了正定二次形式的level set，因此您知道（正如您所指出的）解集是一个（超）椭球体。您可以利用这种结构轻松地在椭球上生成点，而无需借助数值求解器。

不难生成的三组点分别是：椭球与坐标轴相交的点；椭球相对于坐标轴的极值（即边界框与椭球相切的点）；以及椭球主轴与椭球相交的点。您还可以通过首先在单位超球面上生成随机点，然后将这些点转换为椭球来在椭球上生成随机点。这也体现在脚本中。最后，如果x 是一个解决方案，那么-x 也是。（脚本中没有使用这个事实。）

这是一个创建上述点集的脚本。我使用 10 个变量而不是 30 个变量来限制打印输出的数量。

import numpy as np
from scipy.linalg import sqrtm
from scipy.stats import ortho_group


#--------------------------------------------
# Generate example input
#--------------------------------------------

dim = 10

# Create a positive definite matrix H with shape (dim, dim).
# evals is the vector of eigenvalues of H.  This can be replaced
# with any vector of length `dim` containing positive values.
evals = (np.arange(1, dim + 1)/2)**2
# Use a random orthogonal matrix to generate H.
R = ortho_group.rvs(dim)
H = R.T.dot(np.diag(evals).dot(R))

# y determines the level set to be computed.
y = 3.0

#--------------------------------------------
# Compute various points on the ellipsoid
#--------------------------------------------

Hinv = np.linalg.inv(H)

# Ellipsoid intercepts of the coordinate axes
xintercepts = np.diag(np.sqrt(y / np.diag(Hinv)))

# Ellipsoid coordinate extrema
xextrema = H.dot(np.diag(np.sqrt(y / np.diag(H))))

# Ellipsoid intersections with principal axes
evals, evecs = np.linalg.eigh(Hinv)
xprincipal = np.sqrt(y/evals)*evecs

# Generate random points on the ellipsoid.
# The distribution of the random points is not uniform on the ellipsoid.
nsample = 5
r = np.random.randn(H.shape[0], nsample)
# u contains random points on the hypersphere with radius 1.
u = r / np.linalg.norm(r, axis=0)
xrandom = sqrtm(H).dot(np.sqrt(y)*u)

#--------------------------------------------
# Print results
#--------------------------------------------

np.set_printoptions(precision=4, linewidth=120)

print("Columns are the ellipsoid coordinate axis intercepts:")
print(xintercepts)
print("Check:", np.array([x.dot(Hinv.dot(x)) for x in xintercepts.T]))

print()

print("Columns are points where the ellipsoid is tangent to its bounding box:")
print(xextrema)
print("Check:", np.array([x.dot(Hinv.dot(x)) for x in xextrema.T]))

print()

print("Columns are points where the principal axes of the ellipsoid intersect the ellipsoid:")
print(xprincipal)
print("Check:", np.array([x.dot(Hinv.dot(x)) for x in xprincipal.T]))

print()

print("Columns are random points on the ellipsoid:")
print(xrandom)
print("Check:", np.array([x.dot(Hinv.dot(x)) for x in xrandom.T]))

这是脚本运行时的输出。 H 是随机的，因此每次运行脚本时输出都会不同。每个矩阵后以Check: 开头的行显示了为每一列计算x.dot(Hinv.dot(x)) 的结果。 Check: 之后显示的值都应该是 y（有一个小的公差以允许正常的数值不精确）。

Columns are the ellipsoid coordinate axis intercepts:
[[ 2.9341  0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      2.752   0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      2.0081  0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      4.3061  0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      3.4531  0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      1.7946  0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      3.5877  0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      1.8925  0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      1.2435  0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.      2.9075]]
Check: [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]

Columns are points where the ellipsoid is tangent to its bounding box:
[[ 4.2205 -1.6619  0.9486 -1.0971  0.0672  0.0988  2.0569 -0.2758  0.4345 -1.1246]
 [-2.392   6.0745 -0.1126  0.262   0.1952  0.7275 -1.829  -3.3036  0.3353 -0.4161]
 [ 1.2413 -0.1023  5.5224 -1.9892  1.3277  0.6662  0.6003 -2.6505 -1.0728  0.65  ]
 [-1.3493  0.2239 -1.8697  5.1908  0.0367 -0.9623  0.1469  1.5594 -0.3419  0.4431]
 [ 0.0917  0.1851  1.385   0.0407  5.7611 -2.3986  0.884   0.1803 -2.3425 -1.4966]
 [ 0.1385  0.7082  0.7134 -1.0963 -2.4621  5.9136 -2.6088 -0.8365  4.7242 -0.2406]
 [ 2.5176 -1.5554  0.5615  0.1462  0.7926 -2.2789  5.1657  0.0131 -1.8862 -0.0804]
 [-0.3829 -3.1874 -2.8128  1.7606  0.1834 -0.829   0.0149  5.8607 -1.2844  2.1692]
 [ 0.4579  0.2456 -0.8642 -0.293  -1.8089  3.5539 -1.6244 -0.9749  4.4486 -1.8706]
 [-1.4002 -0.3599  0.6185  0.4485 -1.3651 -0.2138 -0.0818  1.945  -2.2096  5.2549]]
Check: [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]

Columns are points where the principal axes of the ellipsoid intersect the ellipsoid:
[[-0.4701  1.1006 -3.008   1.1252 -0.8771  0.3136 -1.0045  2.0293  0.2615 -0.0509]
 [ 3.2968  1.5837  4.3875 -0.4139 -0.9073 -0.2182 -1.6669  0.5303 -0.3532  0.2126]
 [ 0.9165  3.6613 -1.9066 -2.764   1.0091  1.727   0.7535  0.1335 -0.5318  0.3281]
 [-1.5272 -1.9226  2.4674  1.0623 -0.8445  3.5199  0.6734  0.3716  0.018  -0.0609]
 [-2.2736  3.3975  1.2729  1.3181  3.3424  0.6659 -1.0468 -0.2749  0.5697 -0.0946]
 [ 4.6168 -2.0728 -1.9233 -0.4141  1.2233  1.3196 -1.3411 -0.447  -0.3202 -0.3887]
 [-2.7789  1.8796 -1.6604  0.6303 -2.7613  0.7791 -1.534  -1.2874 -0.1734  0.0563]
 [-3.7189 -3.8936 -0.3987 -0.0928  1.7785 -0.2127 -1.1247  0.2619 -0.7321  0.3339]
 [ 3.2782 -1.6147 -1.2554  1.8212  0.203   0.5629 -0.1255 -0.4354  0.8151  0.5622]
 [-1.3981 -1.6701  0.4087 -4.573  -0.423   0.279  -0.8141  0.0621  0.9307 -0.0114]]
Check: [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]

Columns are random points on the ellipsoid:
[[ 1.5979 -1.8849  0.4878 -1.8069 -0.8968]
 [-0.6583 -0.304   0.5729 -0.6411  1.1796]
 [-0.857  -1.0776  1.4855 -2.2664 -1.6846]
 [ 0.1828  0.3555  2.4142  2.7053  3.0825]
 [-0.4522  0.6217  3.8429 -2.4424  0.0587]
 [-0.6646 -2.4867 -0.3544  2.0967 -2.9529]
 [-1.0759 -1.2378  0.7088 -0.561   0.2613]
 [ 2.6662  2.862   0.3385  0.9958  0.821 ]
 [-1.0947 -3.0513 -0.9832  2.4314 -1.6179]
 [ 0.9955  2.3737  0.14   -0.1343 -1.0845]]
Check: [ 3.  3.  3.  3.  3.]

【讨论】：

非常感谢@WarrenWeckesser 的详尽而有见地的回答。您的解释清楚地表明了我如何利用问题结构轻松/有效地在椭圆体上生成点，而无需求助于数值求解器。非常感谢！
一个问题@WarrenWeckesser。 xrandom = sqrtm(H).dot(np.sqrt(y)*u) 是否不应该用 H 的倒数定义：xrandom = sqrtm(Hinv).dot(np.sqrt(y)*u) ？
为了从单位超球面映射到超椭球体，我们乘以二次型矩阵的平方根倒数。在你的问题中，二次形式的矩阵是 Hinv，它的逆矩阵是 H。（你也可以试试，你会发现它不起作用。）

【解决方案2】：

如果 H 是正定的，那么您提供的问题肯定是无约束凸的。这是因为正定矩阵的逆矩阵也是正定矩阵，所以当你取雅可比矩阵时，你会得到 H 逆加 H 逆转置，这本身就是半正定矩阵（参见 here 和 here）。由于这种凸性，我建议您使用像 CVXPY 这样的凸求解器。 CVXPY 也有许多不同的求解器选项，因此您可以从 this list 尝试不同的 QP 选项，看看哪个获得了最佳运行时间（尽管这可能是一个与输入相关的问题）。

【讨论】：

感谢@cmitch 回答了所提供的问题是否为无约束凸问题。