Обратная регрессия гребня: с учетом матрицы отклика и коэффициентов регрессии, найти подходящих предикторов

16

Рассмотрим стандартную задачу регрессии OLS : У меня есть матрицы Y и X и я хочу найти β чтобы минимизировать

L=YXβ2.
Решение дается
β^=argminβ{L}=(XX)+XY.

Я также могу поставить «обратную» проблему: учитывая Y и β , найдите X^ , который даст β^β , то есть сведет к минимуму \ | \ argmin_ \ B \ {L \} - \ B ^ * \ | ^ 2argminβ{L}β2 . Словом, у меня есть матрица ответов Y и вектор коэффициентов β и я хочу найти матрицу предикторов, которая выдает коэффициенты, близкие к β . Это, конечно, также проблема регрессии OLS с решением

X^=argminX{argminβ{L}β2}=Yβ(ββ)+.

Уточнение: как объяснил @ GeoMatt22 в своем ответе, если Y является вектором (т. Е. Если есть только одна переменная ответа), то этот X^ будет рангом один, а обратная задача в значительной степени недоопределена. В моем случае, Y на самом деле является матрицей (т.е. есть много переменных ответа, это многомерная регрессия). Таким образом, X равно n×p , Y равно n×q и β равно p×q .


Я заинтересован в решении "обратной" проблемы регрессии гребня. А именно, моя функция потерь теперь

L=YXβ2+μβ2
и решение
β^=argminβ{L}=(XX+μI)1XY.

«Обратной» проблемой является поиск

X^=argminX{argminβ{L}β2}=?

Опять же, у меня есть матрица ответов Y и вектор коэффициентов β и я хочу найти матрицу предикторов, которая выдает коэффициенты, близкие к β .

На самом деле есть две взаимосвязанные формулировки:

  1. Найти X^ учетом Y и β и μ .
  2. Найдите X^ и μ^ учетом Y и β .

У кого-нибудь из них есть прямое решение?


Вот краткая выдержка из Matlab для иллюстрации проблемы:

% generate some data
n = 10; % number of samples
p = 20; % number of predictors
q = 30; % number of responses
Y = rand(n,q);
X = rand(n,p);
mu = 0;
I = eye(p);

% solve the forward problem: find beta given y,X,mu
betahat = pinv(X'*X + mu*I) * X'*Y;

% backward problem: find X given y,beta,mu
% this formula works correctly only when mu=0
Xhat =  Y*betahat'*pinv(betahat*betahat');

% verify if Xhat indeed yields betahat
betahathat = pinv(Xhat'*Xhat + mu*I)*Xhat'*Y;
max(abs(betahathat(:) - betahat(:)))

Этот код выводит ноль, если, mu=0но не иначе.

амеба говорит восстановить монику
источник
Так как и даны, они не влияют на изменения в потере. Поэтому в (1) вы все еще делаете OLS. (2) одинаково прост, потому что потерю можно сделать сколь угодно малым, если принять произвольно отрицательным в пределах любых ограничений, которые вы сравниваете, чтобы наложить на нее. Это сводит вас к случаю (1). М μBμμ^
whuber
@whuber Спасибо. Я думаю, что я не объяснил это достаточно ясно. Рассмотрим (1). и даны (назовем это ), но мне нужно найти который бы дал коэффициенты регрессии гребня, близкие к , другими словами, я хочу найти минимизирующийЯ не понимаю, почему это должно быть OLS. μ B X B X аргмин B { L r i d g e ( X , B ) } - B 2 .BμBXBX
argminB{Lridge(X,B)}B2.
говорит амеба, восстанови Монику
Как будто у меня есть и я хочу найти такой, что близко к данному . Это не то же самое, что поиск . f(v,w)vargminwf(v,w)wargminvf(v,w)
говорит амеба, восстанови Монику
Изложение в вашем посте сбивает с толку по этому вопросу, потому что, очевидно, вы на самом деле не используете в качестве функции потерь. Не могли бы вы подробно остановиться на специфике проблем (1) и (2) в посте? L
whuber
2
@ hxd1011 Многие столбцы в X обычно называют «множественной регрессией», многие столбцы в Y обычно называют «многомерной регрессией».
говорит амеба: восстанови Монику

Ответы:

11

Теперь, когда вопрос сошелся на более точную формулировку интересующей проблемы, я нашел решение для случая 1 (известный параметр гребня). Это также должно помочь в случае 2 (не аналитическое решение, а простая формула и некоторые ограничения).

Резюме: Ни одна из двух формулировок обратной задачи не имеет уникального ответа. В случае 2 , где параметр гребня неизвестен, существует бесконечно много решений для . В случае 1, где задано , существует конечное число решений для из-за неоднозначности в спектре сингулярных значений.X ω ω [ 0 , ω max ] ω X ωμω2Xωω[0,ωmax]ωXω

(Вывод немного длинный, поэтому TL, DR: в конце есть рабочий код Matlab.)


Недоопределенный случай («МНК»)

задача - это где , и . X R п × р B R р × Q Y R п × д

minBXBY2
XRn×pBRp×qYRn×q

На основе обновленного вопрос, мы будем считать , так находится под определяется с учетом и . Как и в вопросе, мы будем считать « по умолчанию» (минимум -норм) решение , где является Псевдообратным из .B X Y L 2 B = X + Y X + Xn<p<qBXYL2

B=X+Y
X+X

Из разложения по сингулярным числам ( SVD ) для , заданного * псевдообратное значение можно вычислить как ** (* Первые выражения используют полный SVD, а вторые выражения используют сокращенный SVD. ** Для простоты я предполагаю, что имеет полный ранг, т. существует.)X = U S V T = U S 0 V T 0 X + = V S + U T = V 0 S - 1 0 U TX

X=USVT=US0V0T
X+=VS+UT=V0S01UT
S - 1 0XS01

Таким образом, прямая задача имеет решение Для дальнейшего использования отмечу, что , где - это вектор сингулярных значений.S 0 = d i a

BX+Y=(V0S01UT)Y
σ 0 > 0S0=diag(σ0)σ0>0

В обратной задаче задана и . Мы знаем , что пришли из вышеописанного процесса, но мы не знаем . Задача состоит в том, чтобы определить соответствующий .YB X XBBXX

Как было отмечено в обновленном вопросе, в этом случае мы можем восстановить , используя по существу тот же самый подход, т.е. теперь используется Псевдообратный .X 0 = Y BX B

X0=YB+
В

Переопределенный случай (оценка Риджа)

В случае «OLS» недоопределенная задача была решена путем выбора решения с минимальной нормой , т.е. наше «уникальное» решение было неявно регуляризовано .

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

В этом случае у нас есть ряд прямых задач для , , которые задаются как Сбор различных левых и правых векторов в этой коллекции задачи могут быть сведены к следующей проблеме "OLS" где мы ввели расширенные матрицы k = 1 , , q min βX β - y k 2 + ω 2β 2βККзнак равно1,...,Q

минβ| |Иксβ-YК| |2+ω2| |β| |2
min BX ω B - Y 2
Вωзнак равно[β1,...,βК],Yзнак равно[Y1,...,YК]
минВ| |ИксωВ-Y| |2
Иксωзнак равно[Иксωя],Yзнак равно[Y0]

В этом переопределенном случае решение по-прежнему задается псевдообратным но теперь псевдообратное изменение изменяется, в результате чего * где новая матрица "спектра сингулярности" имеет (обратную) диагональ ** (* Несколько сжатые вычисления, необходимые для получения этого, были опущены для краткости. Это похоже на описание здесь для случая . ** Здесь записи вектор выражается через вектор , где все операции вводятся по порядку.)B ω = ( V 0 S - 2 ω U T ) Y σ 2 ω = σ

Вωзнак равноИкс+Y
Вωзнак равно(В0Sω-2UT)Y
pn
σω2знак равноσ02+ω2σ0
пNσ 0σωσ0

Теперь в этой задаче мы все еще можем формально восстановить «базовое решение» как но это больше не является верным решением.

Иксωзнак равноYВω+

Однако аналогия все еще сохраняется в том, что это «решение» имеет SVD с сингулярными значениями приведенными выше.

Иксωзнак равноUSω2В0T
σω2

Таким образом, мы можем вывести квадратное уравнение, связывающее искомые сингулярные значения с восстанавливаемыми сингулярными значениями и параметром регуляризации . Тогда решение будет σ 2 ω ωσ0σω2ω

σ0знак равноσ¯±Δσ,σ¯знак равно12σω2,Δσзнак равно(σ¯+ω)(σ¯-ω)

Демонстрация Matlab ниже (протестирована онлайн через Octave ) показывает, что этот метод решения работает как на практике, так и в теории. Последняя строка показывает, что все сингулярные значения находятся в реконструкции , но я не совсем выяснил, какой корень взять ( = против ). Для это всегда будет root. Это , как правило , кажется, держит для «малых» , в то время как для «больших» в корень , кажется, взять на себя. (Демо ниже установлено в «большой» случай в настоящее время.)Иксσ¯±Δσsgn+-ωзнак равно0+ωω-

% Matlab demo of "Reverse Ridge Regression"
n = 3; p = 5; q = 8; w = 1*sqrt(1e+1); sgn = -1;
Y = rand(n,q); X = rand(n,p);
I = eye(p); Z = zeros(p,q);
err = @(a,b)norm(a(:)-b(:),Inf);

B = pinv([X;w*I])*[Y;Z];
Xhat0 = Y*pinv(B);
dBres0 = err( pinv([Xhat0;w*I])*[Y;Z] , B )

[Uw,Sw2,Vw0] = svd(Xhat0, 'econ');

sw2 = diag(Sw2); s0mid = sw2/2;
ds0 = sqrt(max( 0 , s0mid.^2 - w^2 ));
s0 = s0mid + sgn * ds0;
Xhat = Uw*diag(s0)*Vw0';

dBres = err( pinv([Xhat;w*I])*[Y;Z] , B )
dXerr = err( Xhat , X )
sigX = svd(X)', sigHat = [s0mid+ds0,s0mid-ds0]' % all there, but which sign?

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

Что касается задачи 2 (то есть unknown), приведенное выше дает по крайней мере верхнюю границу для . Чтобы квадратичный дискриминант был неотрицательным, мы должны иметь ωω

ωωМаксимумзнак равноσ¯Nзнак равномин[12σω2]

Для неоднозначности знака с четырьмя корнями следующий фрагмент кода показывает, что независимо от знака любой даст одно и то же прямое -решение для гребня, даже если отличается от .Икс^Вσ0SВD[Икс]

Xrnd=Uw*diag(s0mid+sign(randn(n,1)).*ds0)*Vw0'; % random signs
dBrnd=err(pinv([Xrnd;w*I])*[Y;Z],B) % B is always consistent ...
dXrnd=err(Xrnd,X) % ... even when X is not
GeoMatt22
источник
1
+11. Большое спасибо за все усилия, которые вы приложили для ответа на этот вопрос, и за все обсуждения, которые у нас были. Кажется, это полностью отвечает на мой вопрос. Я чувствовал, что просто принять ваш ответ недостаточно в этом случае; это заслуживает гораздо большего, чем два отклика, что этот ответ в настоящее время имеет. Приветствия.
говорит амеба, восстанови Монику
@amoeba спасибо! Я рад, что это было полезно. Я думаю, что я опубликую комментарий на ваш ответ, который вы связываете, спрашивая, считает ли он, что это уместно и / или есть ли лучший ответ для использования. (Заметьте, что в предисловии к обсуждению SVD он указывает условие , т . пNИкс
Е. Переопределенный
@ GeoMatt22 мой комментарий на оригинальный вопрос говорит, что использование pinvне очень хорошая вещь, вы согласны?
Haitao Du
1
@ hxd1011 В общем случае вы (почти) никогда не хотите явно инвертировать матрицу численно, и это справедливо и для псевдообратной. Здесь я использовал две причины: 1) согласованность с математическими уравнениями + демонстрационный код amoeba и 2) для случая недоопределенных систем стандартные решения «слэша» Matlab могут отличаться от pinv . Почти все случаи в моем коде могут быть заменены соответствующими командами \ или /, которые обычно предпочтительнее. (Это позволяет Matlab выбирать наиболее эффективный прямой решатель.)
GeoMatt22
1
@ hxd1011, чтобы уточнить пункт 2 моего предыдущего комментария по ссылке в вашем комментарии на исходный вопрос: «Если ранг A меньше числа столбцов в A, то x = A \ B не обязательно является минимальным решение нормы. Чем дороже в вычислительном отношении x = pinv (A) * B вычисляет минимальное решение по методу наименьших квадратов. "
GeoMatt22