单纯形投影算法

原创已于 2024-05-06 10:45:11 修改 · 814 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#人工智能

于 2024-04-25 17:30:49 首次发布

文章讨论了如何求解任意点到平移坐标轴面和标准单纯形的投影问题，涉及转换变量、导数分析、最小值求解以及Moreau分解等数学工具，提供了C++代码实现。

本文来自论文：Projection Onto A Simplex

一，任意点到平移坐标轴面的投影

1，求解目标

求：

几何意义就是，考虑平移坐标轴面：

求任意一点y到这个平移坐标轴面的投影z。

以n=2为例：

黑色线即坐标轴，红色线即平移坐标轴面，蓝色线展示了3个不同的投影例子。

2，转换变量

我们可以当t是常数，对于任意点y，求y的投影点z。

反过来，我们也可以把y当常数，即固定点，对于坐标轴面平移到不同的位置，即不同的t，有不同的投影点z。

3，求解结果

不防设y的n个分量是依次递增的，y1<=...<=yn，(下同)

那么上式的结果就是：

其中， $\hat{z}(t)_i=\left\{\begin{matrix} t,if\, y_i>t\\ y_i,if\, y_i<=t \end{matrix}\right.,1\leq i\leq n-1$ , $\hat{z}(t)_n=t$

4，f(t)的导数

f(t)在全体实数域上是可导的

5，f(t)的最小值

f(t)的最小值一定是在导数为0的点，且一定存在唯一的点。

当 $t=(\sum _{i=k}^ny_i-1)/(n-k+1),t>=y_{k-1}$ 时，f(t)取到最小值。

其中， $1\leq k\leq n,y_0=-\infty$

对于这个t有 $\hat{z}(t)_i=\left\{\begin{matrix} t,if\, y_i>t\\ y_i,if\, y_i<=t \end{matrix}\right.,1\leq i\leq n$

二，任意点到标准单纯形的投影

1，求解目标

标准单纯形：

对于任意点y，求投影：

2，公式变形

定义函数：

那么求解目标转化成：

3，Moreau 分解

参考原函数的近端映射函数和共轭函数的近端映射函数的对偶

即：

所以只需要求出

那么也就求出了

4，π函数的共轭函数

5，求解

其中，

所以，取 $t=(\sum _{i=k}^ny_i-1)/(n-k+1),t>=y_{k-1}$ ，其中， $1\leq k\leq n,y_0=-\infty$

再取 $\hat{z}(t)_i=\left\{\begin{matrix} t,if\, y_i>t\\ y_i,if\, y_i<=t \end{matrix}\right.,1\leq i\leq n$

得到的就是

6，总结

不防设y1<=...<=yn

则求

只需要2步：

（1）取唯一的 $t=(\sum _{i=k}^ny_i-1)/(n-k+1),t>=y_{k-1}$ ，其中， $1\leq k\leq n,y_0=-\infty$

（2） $x_i=max(y_i-t,0),1\leq i\leq n$

7，c++实现

class VectorOpt
{
public:
	//拓展数据域，加上id
	template<typename T>
	static inline vector<pair<T, int>>expandWithId(const vector<T>& v)
	{
		vector<pair<T, int>>ans;
		ans.resize(v.size());
		for (int i = 0; i < v.size(); i++)ans[i].first = v[i], ans[i].second = i;
		return ans;
	}
	//给vector拓展，加上id，但只按照原数据进行排序
	template<typename T>
	static inline vector<pair<T, int>> sortWithOrderMatch(const vector<T>& v)
	{
		vector<pair<T, int>>ans = expandWithId(v);
		sort(ans.begin(), ans.end(), cmpJustFirst<T, int>);
		return ans;
	}
private:
	template<typename T>
	static inline bool cmp(T a, T b)
	{
		return a < b;
	}
	template<typename T, typename T2>
	static inline bool cmpJustFirst(pair<T, T2> x, pair<T, T2> y)
	{
		return cmp(x.first, y.first);
	}
};

template<typename T>
vector<T> Simplexprojection(const vector<T> &x) //计算向量x到标准单纯形的投影
{
	vector<pair<T, int>> vp = VectorOpt::sortWithOrderMatch(x);
	int n = vp.size();
	T s = 0;
	for (int i = n - 1; i >= 0; i--) {
		s += vp[i].first;
		T t = (s - 1) / (n - i);
		if (i == 0 || t >= vp[i - 1].first) {
			vector<T> ans(vp.size());
			for (int i = n - 1; i >= 0; i--) {
				ans[vp[i].second] = max(vp[i].first - t, T{ 0 });
			}
			return ans;
		}
	}
	return x;
}