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 回答