如何在 Python 函数最小化中获得速度以找到椭球方程的解
How to gain speed in Python function minimization to find solutions to Ellipsoid equation
简介
使用 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 表示。我想调整我解决问题的方法，以便从我的问题的特殊特征中获得一些性能（在计算时间方面）。
问题分析
我对数值优化知之甚少（计划在有空的时候深入研究）但是在阅读了有趣的 post and some scipy lectures 之后，我得出以下结论：
- 鉴于 H 是正定矩阵，我的问题是凸优化之一。
- 我相信我的问题是无约束优化之一。
- 我不确定我的问题是好是坏。解决方案确实会随着起始值的变化而变化，但我认为我的问题是有条件的。
问题
如果 H 是正定矩阵，我的问题是无约束凸优化问题，这样说是否正确？
Nelder-Mead 算法执行 return 解决 'make sense'。我能否通过利用问题的特征来加快速度？或者可能通过使用不同的优化算法？我知道有一些方法可以提供 jacobian 和 hessian 显着提高性能，但我不确定这是否适用于我的问题。
非常感谢您的帮助。
如果H是正定的，你提供的问题肯定是无约束凸的。这是因为正定矩阵的逆也是正定的，所以当你取雅可比行列式时，你得到 H 逆加上 H 逆转置，这本身就是半正定的（见 here and here). Because of this convexity, I'd recommend you use a convex solver like CVXPY. CVXPY has a bunch of different solver options too, so you can try out different QP options from this list 并看看哪个得到最佳运行时（尽管这可能是一个依赖于输入的问题）。
您要求解的问题是一个可以表示为 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.]
简介
使用 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 表示。我想调整我解决问题的方法，以便从我的问题的特殊特征中获得一些性能（在计算时间方面）。
问题分析
我对数值优化知之甚少（计划在有空的时候深入研究）但是在阅读了有趣的 post
- 鉴于 H 是正定矩阵，我的问题是凸优化之一。
- 我相信我的问题是无约束优化之一。
- 我不确定我的问题是好是坏。解决方案确实会随着起始值的变化而变化，但我认为我的问题是有条件的。
问题
如果 H 是正定矩阵，我的问题是无约束凸优化问题，这样说是否正确？
Nelder-Mead 算法执行 return 解决 'make sense'。我能否通过利用问题的特征来加快速度？或者可能通过使用不同的优化算法？我知道有一些方法可以提供 jacobian 和 hessian 显着提高性能，但我不确定这是否适用于我的问题。
非常感谢您的帮助。
如果H是正定的，你提供的问题肯定是无约束凸的。这是因为正定矩阵的逆也是正定的，所以当你取雅可比行列式时，你得到 H 逆加上 H 逆转置，这本身就是半正定的（见 here and here). Because of this convexity, I'd recommend you use a convex solver like CVXPY. CVXPY has a bunch of different solver options too, so you can try out different QP options from this list 并看看哪个得到最佳运行时（尽管这可能是一个依赖于输入的问题）。
您要求解的问题是一个可以表示为 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.]