lloyd-max 最优标量量化算法分析

变限积分求导公式

假设有函数定义为:

\[K(x)=\int_{\phi(x)}^{\Psi(x)}f(t)dt \\ \frac{dK(x)}{dx} = f[\Psi(x)]\Psi(x)^{\prime} - f[\phi(x)]\phi(x)^{\prime} \]

量化失真与最优标量量化

对于N个量化区间的失真定义为：

\[ D=\sum_{i=1}^{N}(\int_{t_i}^{t_{i+1}}(x-\hat{x_i})^2f(x)dx) \]

当D能取得最小值时，即:

\[\frac{\partial(D)}{\partial{t_i}} = 0 \\ 令 M(x)= \int_{t_{i-1}}^{t_{i}}(x-\hat{x_{i-1}})^2f(x)dx + \int_{t_i}^{t_{i+1}}(x-\hat{x_i})^2f(x)dx \\ 则 \frac{\partial(D)}{\partial{t_i}} = \frac{\partial(M(x))}{\partial{t_i}} = 0 \\ (t_{i} - \hat{x_{i-1}})^2f(t_i)(t_i)^\prime - (t_{i-1} - \hat{x_{i-1}})^2f(t_{i-1})(t_{i-1})^{\prime} + (t_{i+1} - \hat{x_{i}})^2f(t_{i+1})(t_{i+1})^\prime - (t_{i} - \hat{x_{i}})^2f(t_{i})(t_{i})^\prime \\ = (t_{i} - \hat{x_{i-1}})^2f(t_i) - (t_{i} - \hat{x_{i}})^2f(t_{i}) = 0 \]

\(\because \qquad t_i < \hat{x_i} \\ \therefore\)

\[ t_i = \frac{\hat{x_{i-1}} + \hat{x_i}}{2} \qquad i = {1,2 ... N} \qquad (1) \]

令\(\frac{\partial{D}}{\partial{\hat{x_i}}} = 0\), 对该式展开求导，可得

\[令 \qquad M(x) = \int_{t_i}^{t_{i+1}}(x-\hat{x_i})^2f(x)dx \\ 则 \qquad M(X) =\int_{t_i}^{t_{i+1}}(x^2f(x))dx -\int_{t_i}^{t_{i+1}}(2x\hat{x_i})f(x)dx + \int_{t_i}^{t_{i+1}}(\hat{x_i}^2f(x))dx \\ \frac{\partial{M(x)}}{\partial{\hat{x_i}}} = \\ \int_{t_i}^{t_{i+1}}(2\hat{x_i})f(x)dx - \int_{t_i}^{t_{i+1}}(2xf(x))dx = \\ \int_{t_i}^{t_{i+1}}[2(\hat{x_i} - x)]f(x)dx = 0 \]

$\therefore $

\[\qquad\qquad \hat{x_i} = \frac{\int_{t_i}^{t_{i+1}}{xf(x)}dx}{\int_{t_i}^{t_{i+1}}{f(x)}dx} \qquad i = 1,2, ... N-1 \qquad （2） \]

式(1) 和 式（2）构成最优标量量化器的必要条件，即影响量化是真的所有变量的偏导数都为0，这变量为: 边界 量化估计值

其中式(2)和物理学中的质心公式形式相同，因此也可以理解成区间的质心

一般来说，直接求解这两个公式比较麻烦，现实中的算法步骤是：

代码实现

def interval_MSE(x,t1,t2):
    return integrate.quad(lambda t: ((t - x)**2) * f(t), t1, t2)[0]

该函数用于计算积分:

\[\int_{t1}^{t2}[(t-x)^2f(t)]dt \]

def MSE(t,x):
    s = interval_MSE(x[0], -float('Inf'), t[0]) + interval_MSE(x[-1], t[-1], float('Inf'))
    for i in range(1,len(x)-1):
        s = s + interval_MSE(x[i], t[i-1], t[i])
    return s

该函数用于计算积分:

\[MSE = \int_{-\infty}^{t_0}[(t-x_0)^2f(t)]dt + \int_{t_n}^{\infty}[(t-x_n)^2f(t)]dt \qquad 首尾特殊处理 \\ MSE += \sum_{i=1}^{N-1}[\int_{t_{i-1}^{t_i}}(t-x_i)^2f(t)dt] \]

其中n表示序列中的最后一个元素

def centroid(t1,t2):
    if integrate.quad(f, t1, t2)[0] == 0 or t1 == t2:
        return 0
    else:
        return integrate.quad(lambda t:t*f(t), t1, t2)[0] / integrate.quad(f, t1, t2)[0]

该函数用于计算区间质心：

\[\frac{\int_{t_1}^{t_2}(t*f(t))dt}{\int_{t_1}^{t_2}(f(t))dt} \]

# error_threshold 误差阈值
# 不过问题是，这个东西如何保证算法一定收敛
def maxlloyd(t,x,error_threshold):
    # 计算当前的MSE
    e = MSE(t,x)
    error = [e]
    c = 0
    # 控制300次以内，以及误差小于阈值
    while e > error_threshold and c < 300:
        c = c+1
        # 奇数
        if c%2 == 1:
            # adjust thresholds
            # 更新边界
            for i in range(len(t)):
                t[i] = 0.5 * ( x[i] + x[i+1] )
        else:
            # adjust levels
            # 更新量化估计值
            # 负∞和t_0 的质心
            x[0] = centroid(-float('Inf'), t[0])
            # t_n 和∞的质心
            x[-1] = centroid(t[-1], float('Inf'))
            # 其余每个区间的质心
            for i in range(1,len(x)-1):
                x[i] = centroid(t[i-1], t[i])
        # 计算误差
        e = MSE(t,x)
        # 增加误差记录
        error.append(e)
        print(e)
    return x,t,error