Эффективный пересчет взвешенной регрессии наименьших квадратов при изменении веса

Question

Я выполняю регрессию взвешенных наименьших квадратов, как описано в вики: WLS

Мне нужно решить это уравнение: $ B = (t (X)W X) ^ {- 1} * t (X)W y $

Я использую SVD, чтобы найти: $ (t (X)W X) ^ {- 1} $ и сохранить его в матрице. Кроме того, я сохраняю матрицу $ H = (t (X)W X) ^ {- 1} * t (X)W $ и просто делаю следующее для любого нового значения y: B = H y. Таким образом, я могу сэкономить на повторении SVD и матричных умножений при изменении y.

W является диагональной матрицей и обычно не изменяется. Однако иногда я меняю один или два элемента по диагонали в матрице W. В этом случае мне нужно снова сделать SVD и пересчитать H-матрицу. Это явно медленно и требует много времени.

Мой вопрос: если я знаю, что изменилось в W, и ничего не меняется в X, есть ли более эффективный метод для пересчета (t (X)W X) ^ - 1?

Или, другими словами, существует ли эффективный аналитический метод для нахождения B, учитывая, что только диагональные элементы в W могут изменяться на известную величину?

Answer 1

Существует такой метод, в случае если вычисляемое вами обратное является истинным обратным, а не обобщенным обратным (т.е. ни одно из сингулярных значений не равно 0). Однако рекомендуется соблюдать осторожность при использовании этого. Если бы вы делали свои суммы с бесконечной точностью, все было бы хорошо. С конечной точностью, и особенно с почти сингулярными проблемами - если некоторые сингулярные значения очень велики - эти формулы могут привести к потере точности.

Я назову обратное хранилище C. Если вы добавите d (который может быть положительным или отрицательным) к m-му весу, то модифицированная матрица C, скажем, C и модифицированная H, H~, может быть вычислена следующим образом: этот:

('обозначает транспонирование, а e_m - это строка вектора, в котором все 0, кроме m-го слота, равного 1)

Позволять

c = the m'th column of H, divided by the original m'th weight    
a = m'th row of the data matrix X    
f = e_m - a*H    
gamma = 1/d + a*c

(поэтому c - вектор столбцов, а a и f - векторы строк)

затем

C~ = C - c*c'/gamma
H~ = H + c*f/gamma

Если вы хотите найти новый B, скажем, для заданного y, его можно рассчитать с помощью:

r = y[m] - a*B
B~ = B + (r/gamma) * c

Вывод этих формул является простой, но утомительной алгеброй матриц. Лемма обращения матриц очень пригодится.