numpy奇怪的结果是解一个不等式

Question

我用两种格式写了一个不等式，第一种形式是转化为多项式

多项式和矩阵乘法是完全相同的

V0j=yj.T @ P @ yj=12137.5 * (y1**2) + 11438.7 * (y2**2) +(26.381 * 2)*y1*y2

还有下一个常数

lambdaprimezero*nrmedoe0+minelement=8.920678

这意味着多项式不等式和下一个不等式是相同的，并且必须有相同的答案

yj.T @ P @ yj-lambdaprimezero*nrmedoe0<minelement

在朋友的帮助下，我可以提取多项式不等式有效的点（第一个程序）

当我在矩阵一中使用firs程序（多项式不等式）的答案时，问题就出现了

由于它们是相同的，send程序必须在这些点返回足够小的值，但它不是

#first part  
import numpy as np

# 12137.5x^2 + 11438.7y^2 + (26.381*2)xy = 0.000731

Y1 = np.linspace(-0.003, 0.003, num=100)
Y2 = np.linspace(-0.003, 0.003, num=100)

pointsInsideEllipse = []

for y1 in Y1:
    for y2 in Y2:
        if 12137.5 * (y1**2) + 11438.7 * (y2**2) +(26.381 * 2)*y1*y2 < 8.920678:
            pointsInsideEllipse.append([y1, y2])
#print(pointsInsideEllipse)
y=[y1, y2]
P=np.array([[12137.5, 26.381], [26.381,11438.7]])
yj=np.array(y)
pointset=np.array(pointsInsideEllipse)
def msquarefunc(yj):
    VALUE=yj.T @ P @ yj
    return VALUE
point=pointsInsideEllipse
for point in pointset:
     if msquarefunc(point)<8.920678:
      
############################################
###########################
       
# second program(matrix form)
       V0j=yj.T @ P @ yj
       testfeaturej=V0j-lambdaprimezero*nrmedoe0
#print(V0j)
#print(lambdaprimezero*nrmedoe0)
       cj=abs(testfeaturej)
print(cj)
if cj<=minelement:
   print('small enough')

有什么问题

Answer 1

您使用了y=[y1, y2]，而它应该是y=[Y1, Y2]。此外，请注意point=pointsInsideEllipse是无用的，因为point是由以下循环设置的。最大的问题在于第一个假设实际上是错误的。事实上，yj.T @ P @ yj并不等于12137.5 * (y1**2) + 11438.7 * (y2**2) +(26.381 * 2)*y1*y2。我不明白为什么这是真的。在一种情况下，乘以的y1值与第二种情况下的值似乎不一样。 请注意，代码的其余部分不清楚、不完整，因此很难测试/运行

您可以使用矢量化调用，使用Numpy使代码更快、更清晰。此外，您还可以使用以下代码轻松测试等式假设是否错误：

Y1 = np.linspace(-0.003, 0.003, num=100)
Y2 = np.linspace(-0.003, 0.003, num=100)

yj = np.vstack([Y1, Y2])

V0j_1 = yj.T @ P @ yj
V0j_2 = 12137.5 * (Y1.reshape(-1,1)**2) + 11438.7 * (Y2**2) + (26.381*2)*np.outer(Y1,Y2)

# Another (less efficient) way of computing V0j_2:
#V0j_2_bis = np.array([[12137.5 * (y1**2) + 11438.7 * (y2**2) +(26.381 * 2)*y1*y2 for y2 in Y2] for y1 in Y1])
#print(np.allclose(V0j_2, V0j_2_bis)) # True

distance = (V0j_1 - V0j_2)**2
print(np.allclose(V0j_1, V0j_2)) # False
print(distance)

# This shows the two matrices have different symmetries (and so the assumption is wrong)
print(np.allclose(V0j_1, V0j_1.T)) # True
print(np.allclose(V0j_2, V0j_2.T)) # False

distance矩阵显示出巨大的差异（与浮点错误无关）。实际上V0j_1和V0j_2的一些值甚至是相同符号的注释