作为编写自己的弹性网络求解器的热身,我正在尝试使用坐标下降来获得足够快的普通最小二乘版本。
我相信我已经正确实现了坐标下降算法,但是当我使用“快速”版本(见下文)时,该算法非常不稳定,当特征数量为与样本数量相比,规模适中。
线性回归和 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 的情况一致,我认为该算法在形式上是正确的,但在数值上不稳定。有没有人对这个猜测是否正确有任何想法,如果正确,如何在保持算法快速版本的(大部分)性能增益的同时纠正不稳定性?