虽然没有类似任务的经验,但我阅读了一些东西并尝试了一些东西。
当不熟悉这个主题时,似乎很难掌握,而且我找到的所有资源都有些混乱。
我仍然不清楚理论:
- 是上述凸优化问题(局部最小值 = 全局最小值;意味着可以使用强大的求解器!)
- 有更多关于更一般问题的资源(传感器网络
本地化),它们是非凸的,并且已经开发了极其复杂的方法
- 您的三边测量方法是否能够利用 > 3 个点(三边测量与多边测量;至少这段代码看起来不像它可以,这意味着:性能不佳,有噪音!)
这里有一些具有两种方法的示例代码:
-
A: 凸优化:SOCP-Relaxation
-
B: 非线性规划优化
- 使用 scipy.optimize 实现
- 在我的合成实验中非常完美;即使在嘈杂的情况下也能取得良好的效果;尽管我们使用的是数值微分(这里很难使用自动差异)
一些补充说明:
- 在我看来,您的示例 B 肯定有一些(非常糟糕的)噪音或其他问题,因为我的方法完全不可行;而方法 B 尤其适合我的合成数据(至少这是我的印象)
代码:
import numpy as np
import cvxpy as cvx
from scipy.spatial.distance import cdist
from scipy.optimize import minimize
np.random.seed(1)
""" Create noise-free (not anymore!) fake-problem """
real_x = np.random.random(size=2) * 3
M, N = 5, 10
NOISE_DISTS = 0.1
pos = np.array([(i,j) for i in range(M) for j in range(N)]) # ugly -> tile/repeat/stack
real_x_stacked = np.vstack([real_x for i in range(pos.shape[0])])
Y = cdist(pos, real_x[np.newaxis])
Y += np.random.normal(size=Y.shape)*NOISE_DISTS # Let's add some noise!
print('-----')
print('PROBLEM')
print('-------')
print('real x: ', real_x)
print('dist mat: ', np.round(Y,3).T)
""" Helper """
def cost(x, Y, pos):
res = np.linalg.norm(pos - x, ord=2, axis=1) - Y.ravel()
return np.linalg.norm(res, 2)
print('cost with real_x (check vs. noisy): ', cost(real_x, Y, pos))
""" SOLVER SOCP """
def solve_socp_relax(pos, Y):
x = cvx.Variable(2)
y = cvx.Variable(pos.shape[0])
fake_stack = [x for i in range(pos.shape[0])] # hacky
objective = cvx.sum_entries(cvx.norm(y - Y))
x_stacked = cvx.reshape(cvx.vstack(*fake_stack), pos.shape[0], 2) # hacky
constraints = [cvx.norm(pos - x_stacked, 2, axis=1) <= y]
problem = cvx.Problem(cvx.Minimize(objective), constraints)
problem.solve(solver=cvx.ECOS, verbose=False)
return x.value.T
""" SOLVER NLP """
def solve_nlp(pos, Y):
sol = minimize(cost, np.zeros(pos.shape[1]), args=(Y, pos), method='BFGS')
# print(sol)
return sol.x
""" TEST """
print('-----')
print('SOLVE')
print('-----')
socp_relax_sol = solve_socp_relax(pos, Y)
print('SOCP RELAX SOL: ', socp_relax_sol)
nlp_sol = solve_nlp(pos, Y)
print('NLP SOL: ', nlp_sol)
输出:
-----
PROBLEM
-------
real x: [ 1.25106601 2.16097348]
dist mat: [[ 2.444 1.599 1.348 1.276 2.399 3.026 4.07 4.973 6.118 6.746
2.143 1.149 0.412 0.766 1.839 2.762 3.851 4.904 5.734 6.958
2.377 1.432 0.856 1.056 1.973 2.843 3.885 4.95 5.818 6.84
2.711 2.015 1.689 1.939 2.426 3.358 4.385 5.22 6.076 6.97
3.422 3.153 2.759 2.81 3.326 4.162 4.734 5.627 6.484 7.336]]
cost with real_x (check vs. noisy): 0.665125233772
-----
SOLVE
-----
SOCP RELAX SOL: [[ 1.95749275 2.00607253]]
NLP SOL: [ 1.23560791 2.16756168]
编辑: 使用非线性最小二乘法而不是更通用的 NLP 方法可以实现进一步的加速(尤其是大规模)!我的结果仍然相同(如果问题是凸的,则如预期的那样)。 NLP/NLS 之间的时间可能看起来像 9 秒对 0.5 秒!
这是我推荐的方法!
def solve_nls(pos, Y):
def res(x, Y, pos):
return np.linalg.norm(pos - x, ord=2, axis=1) - Y.ravel()
sol = least_squares(res, np.zeros(pos.shape[1]), args=(Y, pos), method='lm')
# print(sol)
return sol.x
特别是第二种方法 (NLP) 也将运行在更大的实例上(cvxpy 的开销很痛苦;这不是 SOCP 求解器的缺点,它应该可以更好地扩展!)。
这里是M, N = 500, 1000 的一些输出,还有一些噪音:
-----
PROBLEM
-------
real x: [ 12.51066014 21.6097348 ]
dist mat: [[ 24.706 23.573 23.693 ..., 1090.29 1091.216
1090.817]]
cost with real_x (check vs. noisy): 353.354267797
-----
SOLVE
-----
NLP SOL: [ 12.51082419 21.60911561]
used: 5.9552763315495625 # SECONDS
所以在我的实验中它是有效的,但我不会给出任何全局收敛保证或重构保证(仍然缺少一些理论)。
起初我虽然关于在 NLP 求解器中使用松弛 SOCP 问题的全局最优作为初始点,但我没有找到任何需要这样做的示例!
一些只是为了好玩的视觉效果:
M, N = 20, 30
NOISE_DISTS = 0.2
...
import matplotlib.pyplot as plt
plt.imshow(Y.reshape(M, N), cmap='viridis', interpolation='none')
plt.colorbar()
plt.scatter(nlp_sol[1], nlp_sol[0], color='red', s=20)
plt.xlim((0, N))
plt.ylim((0, M))
plt.show()
还有一些超级嘈杂的情况(表现不错!):
M, N = 50, 100
NOISE_DISTS = 5
-----
PROBLEM
-------
real x: [ 12.51066014 21.6097348 ]
dist mat: [[ 22.329 18.745 27.588 ..., 94.967 80.034 91.206]]
cost with real_x (check vs. noisy): 354.527196716
-----
SOLVE
-----
NLP SOL: [ 12.44158986 21.50164637]
used: 0.01050068340320306