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

我相信我已经正确地实现了坐标下降算法,但是当我使用"快速"版本(见下文)时,该算法非常不稳定,输出的回归系数通常在特征数量为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) …

我正在尝试编写一个快速坐标下降算法来解决普通最小二乘回归问题.以下Julia代码有效,但我不明白它为什么要分配这么多内存

function OLS_cd{T<:Float64}(A::Array{T,2}, b::Array{T,1}, tolerance::T=1e-12) 
    N,P = size(A)
    x = zeros(P)
    r = copy(b)
    d = ones(P)
    while sum(d.*d) > tolerance
        @inbounds for j = 1:P
            d[j] = sum(A[:,j].*r)
            x[j] += d[j]
            r -= d[j]*A[:,j]
        end
    end
    return(x)
end

关于我生成的数据

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

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

用@benchmark OLS_cd(X, y) …