最小绝对偏差（LAD）

最新推荐文章于 2024-12-09 20:34:38 发布

颹蕭蕭

最新推荐文章于 2024-12-09 20:34:38 发布

阅读量8.1k

点赞数 2

分类专栏： # 机器学习 # 优化问题

我们不生产知识，我们只是互联网的搬运工

本文链接：https://blog.csdn.net/itnerd/article/details/105935302

版权

机器学习同时被 2 个专栏收录

136 篇文章

订阅专栏

优化问题

32 篇文章

订阅专栏

最小绝对偏差 (Least Absolute Deviations, LAD) 与最小二乘法（假设误差服从高斯分布）类似：当假设线性回归的误差服从拉普拉斯分布时，最小绝对偏差回归是对参数的最大似然估计。

问题描述

$\min_{x} \quad || Wx-y||_1$
等价于
$\begin{array}{ll} \min & \nu \\ s.t. & \nu = || Wx-y||_1 \end{array}$
等价于
$\begin{array}{ll} \min & \nu \\ s.t. & || Wx-y||_1 \leq \nu \end{array}$
等价于
$\begin{array}{ll} \min & \nu \\ s.t. & Wx-y \leq \nu, \\ & -Wx+y \leq \nu, \end{array}$
即
$\begin{array}{ll} \min & \nu \\ s.t. & Wx-\nu \leq y, \\ & -Wx -\nu \leq -y, \end{array}$
即
$\begin{array}{ll} \min & \left[\begin{array}{ll} 0 &I\end{array}\right]^\top \left[\begin{array}{ll} x \\ \nu \end{array}\right] \\\\ s.t. & \left[\begin{array}{cc} W & -I \\ -W & -I\end{array}\right] \left[\begin{array}{c} x \\ \nu\end{array}\right] \leq \left[\begin{array}{c} y \\ -y\end{array}\right] \end{array}$

可见，该问题可以转化成标准的线性规划问题！！！

cvxopt 求解器

#Sources: http://cvxopt.org/examples/mlbook/l1.html?highlight=l1

from cvxopt import blas, lapack, solvers
from cvxopt import matrix, spdiag, mul, div, sparse 
from cvxopt import spmatrix, sqrt, base

def l1(P, q):
    P,q = matrix(P), matrix(q)
    m, n = P.size
    c = matrix(n*[0.0] + m*[1.0])
    h = matrix([q, -q])

    def Fi(x, y, alpha = 1.0, beta = 0.0, trans = 'N'):    
        if trans == 'N':
            u = P*x[:n]
            y[:m] = alpha * ( u - x[n:]) + beta*y[:m]
            y[m:] = alpha * (-u - x[n:]) + beta*y[m:]
        else:
            y[:n] =  alpha * P.T * (x[:m] - x[m:]) + beta*y[:n]
            y[n:] = -alpha * (x[:m] + x[m:]) + beta*y[n:]

    def Fkkt(W): 
        d1, d2 = W['d'][:m], W['d'][m:]
        D = 4*(d1**2 + d2**2)**-1
        A = P.T * spdiag(D) * P
        lapack.potrf(A)

        def f(x, y, z):
            x[:n] += P.T * ( mul( div(d2**2 - d1**2, d1**2 + d2**2), x[n:]) 
                + mul( .5*D, z[:m]-z[m:] ) )
            lapack.potrs(A, x)

            u = P*x[:n]
            x[n:] =  div( x[n:] - div(z[:m], d1**2) - div(z[m:], d2**2) + 
                mul(d1**-2 - d2**-2, u), d1**-2 + d2**-2 )

            z[:m] = div(u-x[n:]-z[:m], d1)
            z[m:] = div(-u-x[n:]-z[m:], d2)
        return f

    uls =  +q
    lapack.gels(+P, uls)
    rls = P*uls[:n] - q 

    x0 = matrix( [uls[:n],  1.1*abs(rls)] ) 
    s0 = +h
    
    Fi(x0, s0, alpha=-1, beta=1) 

    if max(abs(rls)) > 1e-10:  
        w = .9/max(abs(rls)) * rls
    else: 
        w = matrix(0.0, (m,1))
    z0 = matrix([.5*(1+w), .5*(1-w)])

    dims = {'l': 2*m, 'q': [], 's': []}
    
    s0 = np.array(s0)
    s0[s0<=0]=1e-8
    s0 =  matrix(s0)
    
    sol = solvers.conelp(c, Fi, h, dims, kktsolver = Fkkt,  
        primalstart={'x': x0, 's': s0}, dualstart={'z': z0})
    return sol['x'][:n]

测试

A = np.random.random((10,10))
x = np.random.randn((10))
y = A @ x
n = np.random.laplace(0, 0.01, 10)
y += n

res = l1(A,y)
res = np.squeeze(res)

plt.subplot(1,2,1)
plt.bar(np.arange(len(res)),res, alpha=0.5)
plt.title('solution')
plt.subplot(1,2,2)
plt.bar(np.arange(len(x)),x, alpha=0.5)
plt.title('ground truth')