Алгоритм сглаживания с координатами в Юлии для наименьших квадратов, не сходящихся

Question

Как разогрев для написания моего собственного решеточного сетевого решателя, я пытаюсь получить достаточно быструю версию обычных наименьших квадратов, реализованных с использованием конического спуска.

Я считаю, что я реализовал алгоритм координат спуска правильно, но когда я использую «быстрый» вариант (см ниже), алгоритма душевнобольно нестабильный, выводя коэффициенты регрессии, которые обычно переполнение 64-битный поплавок когда количество функций имеет умеренный размер по сравнению с количеством выборок.

Линейная регрессия и МНК

Если б = А * х, где А представляет собой матрицу, ха вектор неизвестных коэффициентов регрессии, и у есть выход, я хочу найти х, что сводит к минимуму

|| b - Ax ||^2

Если A [j] - это j-й столбец A и A [-j], то A без столбца j, а столбцы A нормированы так, что || A [j ] ||^2 = 1 для всех j, тогда будет координатное обновление

Координаты Спуск:

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

Я следую вместе с these notes (page 9-10), но вывод является простым исчислением.

Он отметил, что вместо того, чтобы повторно вычислив A [J]^T (B - A [-j] * х [-j]) все время, более быстрый способ сделать это с

Быстрая Координировать Спуск:

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

где общий остаток r = b - Ax вычисляется вне цикла по координатам. Эквивалентность этих правил обновления следует из того, что Ax = A [j] * x [j] + A [-j] * x [-j] и перестраивает термины.

Моя проблема заключается в том, что, хотя второй метод действительно быстрее, он дико численно неустойчив для меня, когда количество функций невелико по сравнению с количеством выборок. Мне было интересно, может ли кто-нибудь понять, почему это так. Следует отметить, что первый более стабильный метод по-прежнему не согласуется с более стандартными методами, так как число признаков приближается к числу выборок.

Юлия Код

Ниже некоторые Джулии код для двух правил обновления:

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) соглашается с двумя другими.

ГИПОТЕЗА

Поскольку все согласуется для режима р < < п, я думаю, что алгоритм формально правильно, но численно нестабильна. У кого-нибудь есть мысли о том, правильно ли это предположение, и если да, то как исправить неустойчивость, сохранив (большинство) прирост производительности быстрой версии алгоритма?

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