Julia中的坐标下降算法用于最小二乘法不收敛

Question

作为编写我自己的弹性网络解算器的热身,我正在尝试使用坐标下降来实现足够快的普通最小二乘法.

我相信我已经正确地实现了坐标下降算法,但是当我使用"快速"版本(见下文)时,该算法非常不稳定,输出的回归系数通常在特征数量为6时浮动64位浮点数.与样品数量相比,尺寸适中.

线性回归和OLS

如果b = A*x,其中A是矩阵,未知回归系数的xa向量,y是输出,我想找到最小化的x

|| b - Ax || ^ 2

如果A [j]是A的第j列,而A [-j]是没有列j的A,并且A的列被归一化,所以对于所有j,|| A [j] || ^ 2 = 1,坐标然后是逐步更新

坐标下降:

x[j]  <--  A[j]^T * (b - A[-j] * x[-j])

我正在关注这些注释(第9-10页),但推导是简单的微积分.

它指出,不是一直重新计算A [j] ^ T(b - A [-j]*x [-j]),更快的方法是

快速坐标下降:

x[j]  <--  A[j]^T*r + x[j]

其中总残差r = b - Ax是在环路坐标之外计算的.这些更新规则的等效性来自注意到Ax = A [j]*x [j] + A [-j]*x [-j]并重新排列术语.

我的问题是,虽然第二种方法确实更快,但只要特征数量与样本数量相比不小,它在数值上就会非常不稳定.我想知道是否有人可能会对为什么会出现这种情况有所了解.我应该注意到,第一种更稳定的方法仍然开始不同意更多标准方法,因为特征数量接近样本数量.

朱莉娅代码

以下是两个更新规则的一些Julia代码:

function OLS_builtin(A,b)
    x = A\b
    return(x)
end

function OLS_coord_descent(A,b)    
    N,P = size(A)
    x = zeros(P)
    for cycle in 1:1000
        for j = 1:P 
            x[j] = dot(A[:,j], b - A[:,1:P .!= j]*x[1:P .!= j])
        end    
    end
    return(x)
end

function OLS_coord_descent_fast(A,b) 
    N,P = size(A)
    x = zeros(P)
    for cycle in 1:1000
        r = b - A*x
        for j = 1:P
            x[j] += dot(A[:,j],r)
        end    
    end
    return(x)
end

问题的例子

我使用以下内容生成数据:

n = 100
p = 50
? = 0.1
?_nz = float([i*(-1)^i for i in 1:10])

? = append!(?_nz,zeros(Float64,p-length(?_nz)))
X = randn(n,p); X .-= mean(X,1); X ./= sqrt(sum(abs2(X),1))
y = X*? + ?*randn(n); y .-= mean(y);

在这里我使用p = 50,并且我在OLS_coord_descent(X,y)和之间得到了很好的一致OLS_builtin(X,y),而OLS_coord_descent_fast(X,y)回归系数则返回指数大的值.

当p小于约20时,OLS_coord_descent_fast(X,y)与其他两个一致.

推测

由于事情同意p << n的制度,我认为该算法在形式上是正确的,但在数值上是不稳定的.有没有人对这个猜测是否正确有任何想法,如果有的话如何纠正不稳定性同时保留算法快速版本的(大部分)性能提升？

Answer 1

快速回答:您r在每次x[j]更新后忘记更新.以下是固定功能,其行为如下OLS_coord_descent:

function OLS_coord_descent_fast(A,b) 
    N,P = size(A)
    x = zeros(P)
    for cycle in 1:1000
        r = b - A*x
        for j = 1:P
            x[j] += dot(A[:,j],r)
            r -= A[:,j]*dot(A[:,j],r)   # Add this line
        end    
    end
    return(x)
end