说明

L1算法中的核心公式简明易懂，但是从核心的优化式向复杂的迭代式的推导却不是那么简单。本文将针对相关公式进行一次推导，最后根据验证的迭代式进行代码实现。

公式推导

L1核心公式如下：
$\underset{X}{\operatorname{argmin}} \sum_{i \in I} \sum_{j \in J}\left\|x_{i}-q_{j}\right\| \theta\left(\left\|x_{i}-q_{j}\right\|\right)+R(X) \tag{1} \\ \theta(r)=e^{-r^{2} /(h / 2)^{2}} , R(X)=\sum_{i \in I} \gamma_{i} \sum_{i^{\prime} \in I \backslash\{i\}} \frac{\theta\left(\left\|x_{i}-x_{i^{\prime} }\right\|\right)}{\sigma_{i}\left\|x_{i} -x_{i^{\prime}}\right\|}$

(1)中 $X, x_i$ 代表sample point集合和单个点， $Q,q_j$ 代表原本的点云点集和单个点， $\gamma_i$ 代表 $x_i$ 对应的平衡常数(balancing constants)， $\sigma_i$ 代表 $x_i$ 的有向度(directionality degree)。
将 $x_i$ 看做函数中的变量，则有以下式子

$\underset{X}{arg\,min}\:f(x_i)= \underset {j \in J}{\sum} f_1(||x_i-q_j||) + \gamma _i \underset {i^, \in I/ \{ i \}}{\sum} f_2(||x_i-x_{i^,}||) \tag {1.1}$

$f_1$ 偏微分：
$\begin{align} \frac {\partial f_1}{\partial x_i}&= \frac {\partial f_1}{\partial ||x_i-q_i||} \times \, \frac {\partial ||x_i-q_i||}{\partial x_i}\\ &= (1-\frac {2{(x_i-q_j)}^2}{{(h/2)}^2}) \theta (||x_i-q_j||) \times \frac {x_i-q_i}{||x_i-q_i||}\\ &=(1-\frac {2{(x_i-q_j)}^2}{{(h/2)}^2}) (x_i-q_j) \alpha _{ij} \tag{1.2} \end{align}$

$f_2$ 偏微分：
$\begin{align} \frac {\partial f_2}{\partial x_i}&= \frac {\partial f_2}{\partial ||x_i-x_{i^,}||} \times \, \frac {\partial ||x_i-x_{i^,}||}{\partial x_i}\\ &=-(1+\frac {2{(x_i-x_{i^,})}^2}{(h/2)^2})\frac {\theta (||x_i-x_{i^,}||)}{\sigma _i(x_i-x_{i^,})^2}\times \frac {x_i-x_{i^,}}{||x_i-x_{i^,}||}\\&= -(1+\frac {2{(x_i-x_{i^,})}^2}{(h/2)^2}) \frac {x_i-x_{i^,}}{\sigma _i} \beta _{ii^,} \tag {1.3} \end{align}$

由于式(1.2)和(1.3)中的系数 $(1-\frac {2{(x_i-q_j)}^2}{{(h/2)}^2})$ 和 $-(1+\frac {2{(x_i-x_{i^,})}^2}{(h/2)^2})$ 在实际运算中近似为1和-1，所以系数可以直接去掉，直接保留符号即可。

最终去掉系数的化简式如下：
$\sum_{j \in J}\left(x_{i}-q_{j}\right) \alpha_{i j}-\gamma_{i} \sum_{i^{\prime} \in I \backslash\{i\}} \frac{x_{i}-x_{i^{\prime}}}{\sigma_{i}} \beta_{i i^{\prime}}=0, i \in I \tag{2} \\ \alpha_{i j}=\frac{\theta\left(\left\|x_{i}-q_{j}\right\|\right)} {\left\|x_{i}-q_{j}\right\|},\: \beta _{ii^,}=\frac {\theta (||x_i-x_{i^,}||)} {||x_i-x_{i^,}||^3}$
P.S. 原论文 $\beta _{ii^,}=\frac {\theta (||x_i-x_{i^,}||)} {||x_i-x_{i^,}||^2}$ ，但是按照(1.3)计算 $\beta _{ii^,}=\frac {\theta (||x_i-x_{i^,}||)} {||x_i-x_{i^,}||^3}$ ，此处采用的是3次方的这个定义。

化简式又可以提取为：
$(\underset {j \in J}{\sum} \alpha _{ij}-\frac {\gamma _i \sum _{i^, \in I/\{ i\}} \beta _{ii^,}} {\sigma _i})x_i=\underset {j \in J}{\sum} \alpha _{ij}q_{j}-\frac {\gamma _i \sum _{i^, \in I/\{ i\}} \beta _{ii^,}x_{i^,}} {\sigma _i} \: \tag {2.1}$

根据新的斥力系数定义 $\mu=\frac{\gamma_{i} \sum_{i^{\prime} \in I \backslash\{i\}} \beta_{i i^{\prime}}}{\sigma_{i}^2 \sum_{j \in J} \alpha_{i j}}, \forall i \in I$ (刊误：原论文里使用公式为 $\mu=\frac{\gamma_{i} \sum_{i^{\prime} \in I \backslash\{i\}} \beta_{i i^{\prime}}}{\sigma_{i} \sum_{j \in J} \alpha_{i j}}, \forall i \in I$ ，经分析 $\sigma_i$ 应为笔误)，可将(2.2)化简为：

$\left(1-\mu \sigma_{i}\right) x_{i}=\frac{\sum_{j \in J} q_{j} \alpha_{i j}}{\sum_{j \in J} \alpha_{i j}} - \mu \sigma_{i} \frac{\sum_{i^{\prime} \in I \backslash\{i\}} x_{i^{\prime}} \beta_{i i^{\prime}}}{\sum_{i^{\prime} \in I \backslash\{i\} } \beta_{i i^{\prime}}} \tag {3}$

简单移项后可得:
$x_{i} = \frac{\sum_{j \in J} q_{j} \alpha_{i j}}{\sum_{j \in J} \alpha_{i j}} + \mu \sigma_{i} \frac{\sum_{i^{\prime} \in I \backslash\{i\}} (x_i - x_{i^{\prime} }) \beta_{i i^{\prime}}}{\sum_{i^{\prime} \in I \backslash\{i\}} \beta_{i i^{\prime}}} \tag {3.1}$

使用不动点迭代^[1]推导出迭代公式：
$x_{i}^{k+1}=\frac{\sum_{j \in J} q_{j} \alpha_{i j}^{k}}{\sum_{j \in J} \alpha_{i j}^{k}}+\mu \sigma_{i}^{k} \frac{\sum_{i^{\prime} \in I \backslash\{i\}}\left(x_{i}^{k}-x_{i^{\prime}}^{k}\right) \beta_{i i^{\prime}}^{k}} {\sum_{i^{\prime} \in I \backslash\{i\}} \beta_{i i^{\prime}}^{k}} \\ \alpha_{i j}=\frac{\theta\left(\left\|x_{i}-q_{j}\right\|\right)} {\left\|x_{i}-q_{j}\right\|},\: \beta _{ii^,}=\frac {\theta (||x_i-x_{i^,}||)} {||x_i-x_{i^,}||^3} \tag {4}$
可以看到，(3.1)和迭代公式(4)已经基本一毛一样了。

代码实现

整个迭代的流程包括:

计算邻居(self neighbor 和origin neighbor)
计算 $f_1$ 项(computeAlphasTerms)和 $f_2$ 项(computeBetasTerms)
根据 $f_1$ 和 $f_2$ 计算新的迭代坐标（公式4）
重复上述步骤直到满足停止条件

迭代函数

double L1median::iterateReturnError()
{
    /*位置迭代更新*/
    paraPtr->add("neighborhood_size", h);
    for (int i = 0; i < sampleInfo.size(); ++i) {
        sampleInfo[i].kind = (sampleInfo[i].kind == pi::Candidate) ? pi::Sample : sampleInfo[i].kind;
    }

    computeAlphasTerms();
    computeBetasTerms();
    double error = updateSamplePos();
    emit infoSignal("RMS:" + QString::number(error, 'f', 6)+'\n');

    /*位置同步*/
    if (synSampleWithInfo()) {
        emit iterateSignal();
    } else {
        emit errorSignal(QString("fail to synchronize sample info"));
        return 0.0;
    }
    /*位置更新后的预备信息计算*/
    computeNeighbors();
    computeSigmas();
    return error;
}

其它具体实现部分过于复杂，此处按下不表，有时间可能另开一章进行说明。

https://blog.csdn.net/jbb0523/article/details/52459797 ↩

L1算法分析和实现(2) 从优化公式到迭代公式

L1算法分析和实现(2) 从优化公式到迭代公式

说明

公式推导

代码实现