怎么用lisp解决求到平面上n点距离和最小的点

aimisiyou

newer 发表于 2014-12-19 20:44
搜索了网络，大都是讨论多边形的费马点的，也就是凸包的费马点，你的意思是求所有点的？ 所有点的这个有 ...

唯一解不敢说，但最小值一定有，因为距离和大于0，闭区间上的连续函数必有最值啊。

aimisiyou

newer 发表于 2014-12-19 20:44
搜索了网络，大都是讨论多边形的费马点的，也就是凸包的费马点，你的意思是求所有点的？ 所有点的这个有 ...

唯一解不敢说（可能有多解），但最小值一定有，闭区间上的连续函数必有最值。

newer

aimisiyou 发表于 2014-12-19 20:50
原来是采用八方向键啊，我初步想法是四方向键，方向越多运算量越大，但方向少了又不精确。

那你试试按照上面的C++代码，改写成LISP看看，我看代码不长。

aimisiyou

newer 发表于 2014-12-19 20:55
那你试试按照上面的C++代码，改写成LISP看看，我看代码不长。

感觉算法有问题，若是局部极小值会导致结果错误。

st788796

aimisiyou 发表于 2014-12-19 20:57
感觉算法有问题，若是局部极小值会导致结果错误。

http://blog.sina.com.cn/s/blog_439371b5010141zn.html

怎样“求空间内一点到其它所有点的距离之和最小”？首先将这个问题形式化：

公式代码：
\min f(x,y) = \min \sum_i \sqrt {(x - x_i)^2 + (y - y_i)^2}

这里是距离之和，而不是平方和。Kmeans聚类中用的评价标准是平方和，如果只有一个类中心，那么可以直接求偏导得到使得平方和最小的点就是中心。这里问题与平方和的解是不是一样的，比如三角形到三个顶点距离之和最短的点就是费马点。
    这里可以用最优化方法中的“搜索”来求解，这一系列方法包括了梯度下降法、牛顿法和共轭梯度法等。在这里用梯度下降法是最简单的，通过这个例子我也明白了为什么实际运用中梯度下降法是应用最广的。相比梯度下降法，牛顿法需要求Hesse矩阵，还是相对麻烦不少。梯度下降法搜索步骤就是每一步都向导数的逆方向将自变量前进一个步长（可变），在这里导数方向就是

公式代码：
\nabla f(x,y) =
\left[
   \begin{array} {lcr}
       \displaystyle \sum_i \frac{x - x_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2}} \\
       \displaystyle \sum_i \frac{y - y_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2}}
   \end{array}
\right]

梯度下降法也有它使用起来让人比较为难的地方，那就是步长很难选取，课本上所给出的例子一般都是针对较简单表达式提出的可变步长计算。在本问题的求解中为简单起见，步长是取的定值。整个过程用Python3实现（起初想用R来做，但是R没法调试……归根结底还是功力不够）实现，结合了scipy和matplotlib两个包，结果看起来还是比较靠谱：

最后附上源代码：

Python 3语言: 高亮代码由发芽网提供
from scipy import *
import pylab

def f(p, pts):
   return sum(sum((p - pts) ** 2, axis=1) ** 0.5)

def fd(p, pts):
   dx = sum((p[0] - pts[:, 0]) / sum((p - pts) ** 2, axis=1) ** 0.5)
   dy = sum((p[1] - pts[:, 1]) / sum((p - pts) ** 2, axis=1) ** 0.5)
   s = (dx ** 2 + dy ** 2) ** 0.5
   dx /= s
   dy /= s
   return array([dx, dy])
   
pts = rand(10, 2)   
x = array([0, 0])

t = 0.1
xstep = x
for k in range(100):
   y = f(x, pts)
   xk = x - t * fd(x, pts)
   yk = f(xk, pts)
   if y - yk > 1e-8:
       x = xk
       y = yk
   elif yk - y > 1e-8:
       t *= 0.5
   else:
       break
   xstep = vstack((xstep, x))
print(x, y)

pylab.plot(pts[:, 0], pts[:, 1], 'bo')
pylab.plot(xstep[:, 0], xstep[:, 1], 'ro')
pylab.plot(xstep[:, 0], xstep[:, 1], 'k-')
pylab.xlabel('iter = %d, Min = %.3f, p = (%.3f, %.3f), t = %f' % (k, y, x[0], x[1], t))
pylab.show()

aimisiyou

st788796 发表于 2014-12-19 22:04
http://blog.sina.com.cn/s/blog_439371b5010141zn.html

怎样“求空间内一点到其它所有点的距离之和最 ...

从最简单的开始，对于三角形的费马点Z，即求d=|Z-a|+|Z-b|+|Z-c|=|Z-a|+|wZ-wb|+|w^2*Z-w^2*c|>=|(1+w+w^2)Z-wb-w^2*c|=|wb+w^2*c|,当且仅当向量Z-a,w(Z-b),w^2(Z-c)三向量同向时取最小值,即各顶点与费马点的连线构成3个120度角。但对于顶点数n>3，此方法就行不通了。按此方法，如四边形则各顶点与费马点的连线构成4个90度角，则此四边形只能为菱形，不符合实际。

st788796

aimisiyou 发表于 2014-12-20 18:12
从最简单的开始，对于三角形的费马点Z，即求d=|Z-a|+|Z-b|+|Z-c|=|Z-a|+|wZ-wb|+|w^2*Z-w^2*c|>=|(1+w+w^ ...

用其它软件吧，用lisp捉摸浪费时间
记得有个数学软件很好

aimisiyou

st788796 发表于 2014-12-20 19:13
用其它软件吧，用lisp捉摸浪费时间
记得有个数学软件很好

matlab吧,只是觉得几何问题与cad联系紧密些，所以尝试cad能不能很好解决。

aimisiyou

这个算法相对简单些，好理解。

st788796

aimisiyou 发表于 2014-12-20 20:57
这个算法相对简单些，好理解。

除了 matlab 还有一个专业数学软件

aimisiyou

(defun ff (ptlst)
    (setq n (length ptlst) fuzz 0.0001)
    (setq pt0 (list (/ (apply '+ (mapcar '(lambda (x) (car x)) ptlst)) n) (/ (apply '+ (mapcar '(lambda (x) (cadr x)) ptlst)) n) ) )
    (setq pt (car ptlst))
    (while (> (distance pt pt0) fuzz)
      (setq pt pt0)
      (setq kdlst (mapcar '(lambda (x) (/ 1.0 (distance x pt0))) ptlst))
      (setq xlst (mapcar '(lambda (x) (car x)) ptlst))
      (setq ylst (mapcar '(lambda (x) (cadr x)) ptlst))
      (setq pt0 (list ( / (apply '+ (mapcar '(lambda (x y) (* x y)) xlst kdlst))  (apply '+ kdlst) )    ( / (apply '+ (mapcar '(lambda (x y) (* x y)) ylst kdlst))  (apply '+ kdlst) )  ))
    )
   (mapcar '(lambda (x) (command "line" pt0 x "")) ptlst)
)

(ff '((0 0) (100 200) (200 100)))

Highflybird

aimisiyou 发表于 2014-12-21 20:44
(defun ff (ptlst)
    (setq n (length ptlst) fuzz 0.0001)
    (setq pt0 (list (/ (apply '+ (mapcar ...

楼主这个迭代的思路我以前也有过。
但不能保证所求得的点一定就是，因为迭代法得到的只不过是一个局部的峰值，但不一定是全局的。 不过这是一个很好的思路。

Highflybird

公式代码：
\nabla f(x,y) =
\left[
    $$\begin{array} {lcr}        \displaystyle \sum_i \frac{x - x_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2}} \\        \displaystyle \sum_i \frac{y - y_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2}}    \end{array}$$
\right]
以上回复为测试数学公式：
只需要将latex代码置于标签之间即可。提供编辑器工具，也可以直接输入代码。或直接在编辑框中点击 公式图标。

st788796

http://zhidao.baidu.com/link?url ... zkF7-RjoCF16GGQOgsa

lingo

sets:point/1..12/:x,y,r,c,d1,d2,c1,c2;link(point,point):d;endsetsdata:x=12 15 26 18 15 5 3 2 7 2 23 56;y=23 43 56 12 67 23 12 45 78 34 23 78;r=500 1000 300 400 700 800 1000 600 100 200 400 600;enddata@for(point:d1=@sqrt((x-x1)^2+(y-y1)^2));@for(point:d2=@sqrt((x-x2)^2+(y-y2)^2));min=@sum(point:r*(c1*d1+c2*d2));@for(point:c1+c2)=1;@for(point:@bin(c1);@bin(c2););

aimisiyou

lingo听说过，用于线性或非线性规划吧，没用过。