2014年西安网络赛 Ellipsoid 题解 模拟退火

网赛那天我正在珠海旅游，由于旅游安排的景点比较少，到中午的时候就把所有的景点看完了，于是下午就回来了。

大概四点多到公司的吧。

然后看看榜，Ellipsoid 这道题过得人还是比较多的，于是看了一下题意。

题意

说的是给你一个椭球，求椭球上的点到原点最近的点。

第一感觉是这不是一道高等数学上的题吗？

但是想想ACM不会出这个没意思的题，所以肯定还有其他方法的。

看看告诉你的椭球方程：

Ax^2 + By^2 + Cz^2 + Dyz + Exz + Fxy = 1.

这个椭球是个特殊的椭球，因为原点总是在椭球内。

这里需要先介绍一个小知识点：

对于球形方程，将某点带入方程，根据值得正负可以判断这个点是在球内还是球外。

• 小于0 球内
• 等于0 球上
• 大于0 球外

既然在球内，自己在脑子里想想它的三维图形，发现球上的点到原点的距离满足中间小，两边大的，也就是凸性的。

这样就很容易想到使用三分来做。

不过三分一般用来解决二维的，对于三维的，使用模拟退火更方便。

模拟退火

什么是模拟退火呢？

假设当前凸性满足凸性，我们可以随意找一个点。

如果当前点不是最优点，则当前点到最优点的那个方向上的点应该都比当前点更优，也就是找到当前点周围的一些更有的点。

之后我们把找到的更优的点来替换当前点。

然后继续上面的，在周围找是否有更有的点。

对于上面的方法，我们需要解决几个问题。

1. 周围的点怎么选。
2. 什么时候停止上面的迭代。

对于第一个问题，一般会指定一个半径R.

然后得到到原点距离大概为R的在球上的点。

至于具体取那些店，这个就需要自己定义一个取点的公式了。

公式定义的好，代码实现简单，定义的不好，很容易就写残了。

看了第一个问题的说明，对于第二个问题，实际上也清楚了。

由于有个半径，我们没找到比当前点更有的点时，要做的是半径减半，然后继续找。

由于这些最优值都是浮点类型，所以半径小于要求的误差即可停止了。

到这里，我们又有一个问题：半径R的初始值该怎么取？

这个需要根据实际情况来取，最终目的是可以快速找到那个最优值，而且还要避免陷入死循环。

实际上只要注意一点就不会陷入死循环：只在找到更优值时才更新当前最优值.也就是找到一个和当前一样优的值时，一定不要更新，且半径一定要减半。

至于更新一个最优值后，半径要不要减半这个问题，我一直感觉不减半比较好，不过当时小舟学长当时给我的建议是减半。

假设每次都减半，则可以搜索的区间是:R(1 + 1/2 + 1/2^2 + 1/2^3 + …) = 2R.

所以如果你确定最优值到当前点的距离，那就可以取一个合适的R了。

对于这个模拟退火，可以使用球求抛物线的最低点来模拟一下就理解了。

假设抛物线方程是 y = Ax^2 + Bx + C ，求最小值y,误差小于esp=10e-5即可。

我们随意取抛物线上一个(x0,y0), 去一个半径R=1000000.

则可以得到新的两个点 (x0+R,y1),(x0-R, y2).

如果两个值有一个比(x0,y0)更优(即比y0小)，则更新(x0,y0),否则R=R/2;

当R<esp 时，我们就找到答案了。

解决思路

好了，不多说了，大概就是这个意思。

现在对于这道题，我们同样面临着上面说的几个问题。

• 随机选一个点，选哪个点呢？

我选的是这三个点的最优点:(0,0,z),(0,y,0),(x,0,0).

• 选多大的半径。

由于这是个椭球，我们通过三维转换可以大概写成下面的标准椭球的形式。

A(x +a)^2 + B(y + b)^2 + C(z + c)^2 = r

而这里我假设这个椭球就是标准椭球。

则三个半径差不多为 1/A, 1/B, 1/C.

然后取一个最大的作为半径即可。

其实这里算是随便取得一个值吧，什么三个半径也只是估计的公式。

半径取得小了，跑的时间长点，只要不是太小，都不影响的。

• 怎么取周围的点。

需要保证点在椭球上。

这个说简单也简单，说难也难。

主要看想的点的公式了。

由于这里是三维空间，我们不知道更优点在那个方位，所以四面八方都需要取一些点才可靠。

于是我想到了立方体，我们目前的点算是中心那个点，半径为立方体边长的一般，这样可以快速定位其余的26个点。

但是我们的目标不是那26个立方体上的点，而是椭球上的点。

所以这里选择先确定两个维度的值，代入方程解出第三维度的值。

这样就可以保证点一定在椭球上了。

这样可以得到几个点呢？

两维确定一条线，一条可能与椭球有两个焦点，所以最多有54个点。

当然实际情况并没有这么多，但是这些都是常数级的，都不是问题。

看下面的代码就知道了。

for(int i=-1; i<=1; i++) {
    for(int j=-1; j<=1; j++) {
        getZ(x+r*i,y+r*j);
        getX(y+r*i,z+r*j);
        getY(x+r*i,z+r*j);
    }
}

详见代码

这样模拟退火的一些小障碍都得到了解决，接下来就看看AC的代码吧。

详见我的github ( tiankonguse ):https://github.com/tiankonguse/ACM/blob/master/hdu/5040.cpp

致小舟学长