一个利用Tensorflow求解几何问题的例子

原创
2018/03/02 09:44
阅读数 7.5K

知乎上有一个问题,内容是已知空间三个点的坐标,求三个点所构成的圆的圆心坐标(编程实现)?

根据圆的定义,这道题的核心就是找到一个点,到已知的三个点的距离相等,利用数学知识可以求解如下:

例如 :给定a(x1,y1) b(x2,y2) c(x3,y3)求外接圆心坐标O(x,y)
1. 首先,外接圆的圆心是三角形三条边的垂直平分线的交点,我们根据圆心到顶点的距离相等,可以列出以下方程:       
(x1-x)*(x1-x)+(y1-y)*(y1-y)=(x2-x)*(x2-x)+(y2-y)*(y2-y);       
(x2-x)*(x2-x)+(y2-y)*(y2-y)=(x3-x)*(x3-x)+(y3-y)*(y3-y);
2.化简得到:        
2*(x2-x1)*x+2*(y2-y1)y=x2^2+y2^2-x1^2-y1^2;        
2*(x3-x2)*x+2*(y3-y2)y=x3^2+y3^2-x2^2-y2^2;        
令:A1=2*(x2-x1);            
B1=2*(y2-y1);            
C1=x2^2+y2^2-x1^2-y1^2;            
A2=2*(x3-x2);            
B2=2*(y3-y2);            
C2=x3^2+y3^2-x2^2-y2^2;        
即:A1*x+B1y=C1;            
A2*x+B2y=C2;
3.最后根据克拉默法则:          
x=((C1*B2)-(C2*B1))/((A1*B2)-(A2*B1));          
y=((A1*C2)-(A2*C1))/((A1*B2)-(A2*B1));

当然,我们今天不是来学习数学公式和数学推导的。Tensorflow是google开源的一款深度学习的工具,其实我们可以利用Tensoflow提供了强大的数学计算能力来求解类似的数学问题。

这道题,我们可以利用梯度下降算法,因为圆心是一个最优解,任何其它点都不满条件。(前提是这三个点不在一条直线上,否则是没有解的)

好了,我们先看代码先,然后在解释。

import tensorflow as tf
import numpy

# Parameters
learning_rate = 0.1
training_epochs = 3000
display_step = 50

# Training Data, 3 points that form a triangel
train_X = numpy.asarray([3.0,6.0,9.0])
train_Y = numpy.asarray([7.0,9.0,7.0])

# tf Graph Input
X = tf.placeholder("float")
Y = tf.placeholder("float")

# Set vaibale for center
cx = tf.Variable(3, name="cx",dtype=tf.float32)
cy = tf.Variable(3, name="cy",dtype=tf.float32)

# Caculate the distance to the center and make them as equal as possible
distance = tf.pow(tf.add(tf.pow((X-cx),2),tf.pow((Y-cy),2)),0.5)
mean = tf.reduce_mean(distance)
cost = tf.reduce_sum(tf.pow((distance-mean),2)/3)
# Gradient descent
optimizer = tf.train.GradientDescentOptimizer(learning_rate).minimize(cost)

# Initialize the variables (i.e. assign their default value)
init = tf.global_variables_initializer()

# Start training
with tf.Session() as sess:
    sess.run(init)

    # Fit all training data
    for epoch in range(training_epochs):
        sess.run(optimizer, feed_dict={X: train_X, Y: train_Y})
        c = sess.run(cost, feed_dict={X: train_X, Y:train_Y})
        if (c - 0) < 0.0000000001:
            break
        #Display logs per epoch step
        if (epoch+1) % display_step == 0:
            c = sess.run(cost, feed_dict={X: train_X, Y:train_Y})
            m = sess.run(mean, feed_dict={X: train_X, Y:train_Y})
            print "Epoch:", '%04d' % (epoch+1), "cost=", "{:.9f}".format(c), \
                "CX=", sess.run(cx), "CY=", sess.run(cy), "Mean=", "{:.9f}".format(m)

    print "Optimization Finished!"
    training_cost = sess.run(cost, feed_dict={X: train_X, Y: train_Y})
    print "Training cost=", training_cost, "CX=", round(sess.run(cx),2), "CY=", round(sess.run(cy),2), "R=", round(m,2), '\n'

运行以上的python代码,结果如下:

Epoch: 0050 cost= 0.290830940 CX= 5.5859795 CY= 2.6425467 Mean= 5.657848835
Epoch: 0100 cost= 0.217094064 CX= 5.963002 CY= 3.0613017 Mean= 5.280393124
Epoch: 0150 cost= 0.173767462 CX= 5.997781 CY= 3.5245996 Mean= 4.885882378
Epoch: 0200 cost= 0.126330480 CX= 5.9999194 CY= 4.011508 Mean= 4.485837936
Epoch: 0250 cost= 0.078660280 CX= 5.9999976 CY= 4.4997787 Mean= 4.103584766
Epoch: 0300 cost= 0.038911112 CX= 5.9999976 CY= 4.945466 Mean= 3.775567770
Epoch: 0350 cost= 0.014412695 CX= 5.999998 CY= 5.2943544 Mean= 3.535865068
Epoch: 0400 cost= 0.004034557 CX= 5.999998 CY= 5.5200934 Mean= 3.390078306
Epoch: 0450 cost= 0.000921754 CX= 5.999998 CY= 5.6429324 Mean= 3.314131498
Epoch: 0500 cost= 0.000187423 CX= 5.999998 CY= 5.7023263 Mean= 3.278312683
Epoch: 0550 cost= 0.000035973 CX= 5.999998 CY= 5.7292333 Mean= 3.262284517
Epoch: 0600 cost= 0.000006724 CX= 5.999998 CY= 5.7410445 Mean= 3.255288363
Epoch: 0650 cost= 0.000001243 CX= 5.999998 CY= 5.746154 Mean= 3.252269506
Epoch: 0700 cost= 0.000000229 CX= 5.999998 CY= 5.7483506 Mean= 3.250972748
Epoch: 0750 cost= 0.000000042 CX= 5.999998 CY= 5.749294 Mean= 3.250416517
Epoch: 0800 cost= 0.000000008 CX= 5.999998 CY= 5.749697 Mean= 3.250178576
Epoch: 0850 cost= 0.000000001 CX= 5.999998 CY= 5.749871 Mean= 3.250076294
Epoch: 0900 cost= 0.000000000 CX= 5.999998 CY= 5.7499437 Mean= 3.250033140
Optimization Finished!
Training cost= 9.8869656e-11 CX= 6.0 CY= 5.75 R= 3.25 

经过900多次的迭代,圆心位置是(6.0,5.75),半径是3.25。

# Parameters
learning_rate = 0.1
training_epochs = 3000
display_step = 50
  • learning_rate 是梯度下降的速率,这个值越大,收敛的越快,但也有可能会错过最优解
  • training_epochs是学习迭代的次数
  • display_step是每多少次迭代显示当前的计算结果
# Training Data, 3 points that form a triangel
train_X = numpy.asarray([3.0,6.0,9.0])
train_Y = numpy.asarray([7.0,9.0,7.0])

# tf Graph Input
X = tf.placeholder("float")
Y = tf.placeholder("float")

# Set vaibale for center
cx = tf.Variable(3, name="cx",dtype=tf.float32)
cy = tf.Variable(3, name="cy",dtype=tf.float32)
  • train_X,train_Y是三个点的x,y坐标,这里我们选了(3,7)(6,9)(9,7)三个点
  • X,Y是计算的输入,在计算过程中我们会使用训练数据输入X,Y
  • cx,cy是我们想要找的圆心点,初始值设置为(3,3),一般的学习算法会使用随机的初始值,这里我选了三角形中的一个点,这样做一般会减少迭代的次数。
# Caculate the distance to the center and make them as equal as possible
distance = tf.pow(tf.add(tf.pow((X-cx),2),tf.pow((Y-cy),2)),0.5)
mean = tf.reduce_mean(distance)
cost = tf.reduce_sum(tf.pow((distance-mean),2)/3)
# Gradient descent
optimizer = tf.train.GradientDescentOptimizer(learning_rate).minimize(cost)

这几行代码是算法的核心。

  • distance是利用两个点的距离公式算出三个点到圆心的距离
  • mean是三个距离的平均值
  • cost是三个距离的方差,我们的目标是让三个点到圆心的距离一样,也就是方差最小(cx/cy为圆心的时候,这个方差为零)
  • optimizer是梯度下降的训练函数,目标是使得cost(方差)最小

下面就是训练的过程了:

# Start training
with tf.Session() as sess:
    sess.run(init)

    # Fit all training data
    for epoch in range(training_epochs):
        sess.run(optimizer, feed_dict={X: train_X, Y: train_Y})
        c = sess.run(cost, feed_dict={X: train_X, Y:train_Y})
        if (c - 0) < 0.0000000001:
            break
        #Display logs per epoch step
        if (epoch+1) % display_step == 0:
            c = sess.run(cost, feed_dict={X: train_X, Y:train_Y})
            m = sess.run(mean, feed_dict={X: train_X, Y:train_Y})
            print "Epoch:", '%04d' % (epoch+1), "cost=", "{:.9f}".format(c), \
                "CX=", sess.run(cx), "CY=", sess.run(cy), "Mean=", "{:.9f}".format(m)

    print "Optimization Finished!"
    training_cost = sess.run(cost, feed_dict={X: train_X, Y: train_Y})
    print "Training cost=", training_cost, "CX=", round(sess.run(cx),2), "CY=", round(sess.run(cy),2), "R=", round(m,2), '\n'
  • 初始化tf的session
  • 开始迭代
  • 计算cost值,当cost小于一定的值的时候,推出迭代,说明我们已经找到了圆心
  • 最后打印出训练的结果

原题目是空间上的点,我的例子是平面上的点,其实没有本质差别。可以加一个Z轴的数据。这个题,三维其实是多余的,完全可以把空间上的三个点投影到平面上来解决。

我用JS实现过另一个算法,但是不是总收敛。大家可以参考着看一看:

其中,绿色三个点条件点,红色的圆是最终的学习结果,黄色的中心点学习的轨迹。

利用这个例子,我想说的是:

  • Tensorflow不仅仅是一个深度学习的工具,它提供了强大的数据计算能力,可以用于解决很多的数学问题
  • 机器学习的本质是通过一组数据来找到答案,这也是数学的作用,所以很多的数学问题都可以用机器学习的思路来解决。

 

展开阅读全文
打赏
7
1 收藏
分享
加载中
楼主,你好,可以用thensorflow求解目标规划问题吗
2018/03/07 11:45
回复
举报
if __name__ =="__main__":
points = [[3, 7],[6, 9],[9, 7]]
point = [3, 3]

lam = 0.1
count = 1
while DifSum(point, points) > math.exp(-2):
ori = point
#梯度下降
point0 = ori0 - lam * DeltaDifSum_x(ori, points)
point1 = ori1 - lam * DeltaDifSum_y(ori, points)
count += 1

print("count=", count)
print("circle:", point)
2018/03/07 11:02
回复
举报
#目标函数X偏导
def DeltaDifSum_x(point, points):
d1 = Dif(point, points0)
d2 = Dif(point, points1)
d3 = Dif(point, points2)
delta_d1 = Delta_x(point, points0)
delta_d2 = Delta_x(point, points1)
delta_d3 = Delta_x(point, points2)
delta_f1 = (d1 - d3)/abs(d1 - d3) * (delta_d1 - delta_d3)
delta_f2 = (d2 - d3)/abs(d2 - d3) * (delta_d2 - delta_d3)
return delta_f1 + delta_f2

#目标函数Y偏导
def DeltaDifSum_y(point, points):
d1 = Dif(point, points[0])
d2 = Dif(point, points[1])
d3 = Dif(point, points[2])
delta_d1 = Delta_y(point, points[0])
delta_d2 = Delta_y(point, points[1])
delta_d3 = Delta_y(point, points[2])
delta_f1 = (d1 - d3)/abs(d1 - d3) * (delta_d1 - delta_d3)
delta_f2 = (d2 - d3)/abs(d2 - d3) * (delta_d2 - delta_d3)
return delta_f1 + delta_f2
2018/03/07 11:02
回复
举报
该评论暂时无法显示,详情咨询 QQ 群:912889742
不错
2018/03/07 10:56
回复
举报
不错不错,好例子,收藏。
2018/03/02 13:53
回复
举报
更多评论
打赏
6 评论
1 收藏
7
分享
返回顶部
顶部