3

作为编写自己的弹性网络求解器的热身,我正在尝试使用坐标下降来获得足够快的普通最小二乘版本。

我相信我已经正确实现了坐标下降算法,但是当我使用“快速”版本(见下文)时,该算法非常不稳定,当特征数量为与样本数量相比,规模适中。

线性回归和 OLS

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

||b - 斧头||^2

如果 A[j] 是 A 的第 j 列并且 A[-j] 是没有列 j 的 A,并且对 A 的列进行归一化使得 ||A[j]||^2 = 1 对于所有 j,坐标然后是明智的更新

坐标下降:

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 的情况一致,我认为该算法在形式上是正确的,但在数值上不稳定。有没有人对这个猜测是否正确有任何想法,如果正确,如何在保持算法快速版本的(大部分)性能增益的同时纠正不稳定性?

4

1 回答 1

5

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
于 2016-12-16T22:26:43.640 回答