【问题标题】:Finding center point given distance matrix给定距离矩阵寻找中心点
【发布时间】:2018-03-31 20:20:05
【问题描述】:

我有一个矩阵(实际上是一个加载的图像),其中每个元素与某个未知中心点的距离为 L2。

这是一个简单的例子

A = [1.4142  1.0000  1.4142  2.2361]
    [1.0000  0.0000  1.0000  2.0000]
    [1.4142  1.0000  1.4142  2.2361]

在这种情况下,中心显然在坐标 (1,1) 处(索引为 A[1,1] 的索引为 0 的矩阵或二维数组)。

但是,在我的中心不受整数索引限制的情况下,它就不再那么明显了。比如给定这个矩阵B,我的中心坐标在哪里?

B = [3.0292  1.9612  2.8932  5.8252]
    [1.2292  0.1612  1.0932  4.0252]
    [1.4292  0.3612  1.2932  4.2252]

你怎么会发现在这种情况下的答案是在第 1.034 行和第 1.4 列?

我知道trilateration 解决方案(有provided MATLAB code to visualize that in 3D previously),但有没有更有效的方法(例如没有矩阵求逆的方法)?

这个问题与语言无关,因为我正在寻找更多算法帮助。如果您可以在解决方案中坚持使用 MATLAB、Python 或 C++,那就太好了;-)。

【问题讨论】:

  • A 或 B 中的值是精确值还是噪声?
  • 它们可以被视为相当精确
  • 我问是因为 - 正如 sascha 也提到的 - B 似乎是“错误的”。虽然A显然是正确的,但是,您给出的中心点的B应该是类似的东西[[1.74044707 1.10867308 1.90503431 1.90503438 0.60096256 1.60096256 1.60096256 1.60096256 1.040096201 1.04596201 1.04554121 1.04596201 1.040096291 1.040096201 1.040096201 1.040096201 1.86899861 1.70096256 1.70096256 1.60096256 1.60096201 1.
  • @JulianS.: B 绝对正确——我根据 (1.034, 1.4) 的真实值计算了这些距离。我认为混乱可能是你假设同一个中心?编辑:该死的。我有一个错字——我把它改成了正确的表述“我的中心不受限制......”
  • 我想我知道出了什么问题:B 中的值是平方的,所以你忘了取 sqrt()。在 A 中正确,但在 B 中不正确。

标签: matrix computer-vision distance mathematical-optimization numerical-methods


【解决方案1】:

虽然没有类似任务的经验,但我阅读了一些东西并尝试了一些东西。

当不熟悉这个主题时,似乎很难掌握,而且我找到的所有资源都有些混乱。

我仍然不清楚理论:

  • 是上述凸优化问题(局部最小值 = 全局最小值;意味着可以使用强大的求解器!)
    • 有更多关于更一般问题的资源(传感器网络 本地化),它们是非凸的,并且已经开发了极其复杂的方法
  • 您的三边测量方法是否能够利用 > 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

【讨论】:

    【解决方案2】:

    如果我理解正确,您有一个矩阵 A,其中 A[i,j] 保存从 (i,j) 到某个未知点 (y,x) 的距离。你可以像这样找到 (y,x):

    对 A 的每个元素求平方,得到一个矩阵 B。 然后我们想找到 (y,x) 所以

    (y-i)*(y-i) + (x-j)*(x-j) = B[i,j]
    

    从0,0方程中减去每个方程并重新排列:

    2*i*y + 2*j*x = B[0,0] + i*i + j*j - B[i,j]
    

    这可以通过线性最小二乘法解决。请注意,由于有 2 个未知数,所涉及的矩阵求逆(更好的因式分解)将在 2x2 矩阵上,因此不耗时。确实,只要给定 A 的维度,您就可以解析出所需的矩阵及其逆矩阵。

    【讨论】: