Odwrócenie regresji grzbietu: biorąc pod uwagę macierz reakcji i współczynniki regresji, znajdź odpowiednie predyktory

Rozważ standardowy problem regresji OLS : Mam macierze i i chcę znaleźć aby zminimalizować L = \ | \ Y- \ X \ B \ | ^ 2. Rozwiązanie podano przez \ hat \ B = \ argmin_ \ B \ {L \} = (\ X ^ \ top \ X) ^ + \ X ^ \ top \ Y.$\newcommand{\Y}{\mathbf Y}\newcommand{\X}{\mathbf X}\newcommand{\B}{\boldsymbol\beta}\DeclareMathOperator*{argmin}{argmin}$$\Y$$\X$$\B$

$$L=\|\Y-\X\B\|^2.$$ $$\hat\B=\argmin_\B\{L\} = (\X^\top\X)^+\X^\top \Y.$$

Mogę również stanowić problem „odwrotny”: biorąc pod uwagę $\Y$ i $\B^*$ , znajdź $\hat\X$ , co da $\hat\B\approx \B^*$ , tj. Zminimalizuje $\|\argmin_\B\{L\}-\B^*\|^2$ . Innymi słowy, mam macierz odpowiedzi $\Y$ i wektor współczynnika $\B^*$ i chcę znaleźć macierz predykcyjną, która dałaby współczynniki zbliżone do $\B^*$ . Jest to oczywiście także problem regresji OLS z rozwiązaniem

$$\hat\X = \argmin_\X\Big\{\|\argmin_\B\{L\}-\B^*\|^2\Big\} = \Y\B^\top(\B\B^\top)^{+}.$$

Aktualizacja wyjaśnienia: Jak wyjaśnił @ GeoMatt22 w swojej odpowiedzi, jeśli $\Y$ jest wektorem (tj. Jeśli istnieje tylko jedna zmienna odpowiedzi), to ten $\hat \X$ będzie miał jeden stopień, a problem odwrotny jest znacznie niedookreślony. W moim przypadku $\Y$ jest właściwie macierzą (tzn. Istnieje wiele zmiennych odpowiedzi, jest to regresja wielowymiarowa ). Więc $\X$ to $n\times p$ , $\Y$ to $n\times q$ a $\B$ to $p\times q$ .

---

Jestem zainteresowany rozwiązaniem problemu „odwrotnego” regresji grzbietu. Mianowicie, moją funkcją utraty jest teraz

$$L=\|\Y-\X\B\|^2+\mu\|\B\|^2$$ a rozwiązaniem jest $$\hat\B=\argmin_\B\{L\}=(\X^\top \X+\mu\mathbf I)^{-1}\X^\top \Y.$$

Problem „wstecz” polega na znalezieniu

$$\hat\X = \argmin_\X\Big\{\|\argmin_\B\{L\}-\B^*\|^2\Big\} = \;?$$

Ponownie mam macierz odpowiedzi $\Y$ i wektor współczynnika $\B^*$ i chcę znaleźć macierz predykcyjną, która dałaby współczynniki zbliżone do $\B^*$ .

W rzeczywistości istnieją dwie powiązane formuły:

1. Znajdź $\hat\X$ dany $\Y$ i $\B^*$ i $\mu$ .
2. Znajdź $\hat\X$ i $\hat \mu$ podano $\Y$ i $\B^*$ .

Czy któryś z nich ma bezpośrednie rozwiązanie?

---

Oto krótki fragment Matlaba ilustrujący problem:

% 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(:)))

Ten kod generuje zero, jeśli mu=0nie inaczej.

Teraz, gdy pytanie skupiło się na bardziej precyzyjnym sformułowaniu interesującego nas problemu, znalazłem rozwiązanie dla przypadku 1 (znany parametr grzbietu). Powinno to również pomóc w przypadku 2 (nie do końca rozwiązanie analityczne, ale prosta formuła i niektóre ograniczenia).

Podsumowanie: Żadna z dwóch formuł odwrotnych problemów nie ma jednoznacznej odpowiedzi. W przypadku 2 , w którym parametr grzbietu jest nieznany, istnieje nieskończenie wiele rozwiązań , dla . W przypadku 1, w którym podano , istnieje skończona liczba rozwiązań dla , ze względu na niejednoznaczność w widmie o liczbie pojedynczej.$\mu\equiv\omega^2$$X_\omega$$\omega\in[0,\omega_\max]$$\omega$$X_\omega$

(Wyprowadzenie jest nieco długie, więc TL, DR: na końcu jest działający kod Matlab.)

---

Nieokreślony przypadek („OLS”)

Problem przesyłania dalej to gdzie , i

$$\min_B\|XB-Y\|^2$$ $X\in\mathbb{R}^{n\times p}$$B\in\mathbb{R}^{p\times q}$$Y\in\mathbb{R}^{n\times q}$ .

Na podstawie zaktualizowanego pytanie, będziemy zakładać, , więc jest ustalana pod podany i . Podobnie jak w kwestii, to przyjmujemy "domyślne" (minimalna -norm) Roztwór , gdzie jest pseudoinverse z .$n<p<q$$B$$X$$Y$$L_2$

$$B=X^+Y$$ $X^+$$X$

Z rozkładu wartości pojedynczej ( SVD ) , podanego przez * pseudoinwersję można obliczyć jako ** (* Pierwsze wyrażenia używają pełnego SVD, podczas gdy drugie wyrażenia używają zmniejszonego SVD. ** Dla uproszczenia zakładam, że ma pełną pozycję, tj.$X$

$$X=USV^T=US_0V_0^T$$ $$X^+=VS^+U^T=V_0S_0^{-1}U^T$$ $X$$S_0^{-1}$ .)

Tak więc problem naprzód ma rozwiązanie W przyszłości, zauważam, że , gdzie

$$B\equiv X^+Y=\left(V_0S_0^{-1}U^T\right)Y$$ $S_0=\mathrm{diag}(\sigma_0)$$\sigma_0>0$ to wektor wartości pojedynczych.

W odwrotnym problemem, daje nam i . Wiemy, że pochodzi z powyższego procesu, ale nie wiemy, . Następnie należy ustalić odpowiedni$Y$$B$$B$$X$$X$ .

Jak zaznaczono w zaktualizowanym pytanie, w tym przypadku możemy odzyskać stosując zasadniczo takie samo podejście, tzn teraz używając pseudoinverse z .$X$

$$X_0=YB^+$$ $B$

---

Przereklamowany przypadek (estymator grzbietu)

W przypadku „OLS” niedostatecznie określony problem został rozwiązany poprzez wybranie rozwiązania o minimalnej normie , tj. Nasze „unikalne” rozwiązanie zostało domyślnie uregulowane .

Zamiast wybierać rozwiązanie normy minimalnej , tutaj wprowadzamy parametr aby kontrolować „jak małą” powinna być norma, tzn. Używamy regresji grzbietu .$\omega$

W tym przypadku mamy serię problemów do przodu dla , , które są podane przez Zebranie różnych wektorów po lewej i prawej stronie do ta kolekcja problemy można sprowadzić do następującego problemu „OLS” gdzie wprowadziliśmy rozszerzone macierze $\beta_k$$k=1,\ldots,q$

$$\min_\beta\|X\beta-y_k\|^2+\omega^2\|\beta\|^2$$ $$B_{\omega}=[\beta_1,\ldots,\beta_k] \quad,\quad Y=[y_1,\ldots,y_k]$$ $$\min_B\|\mathsf{X}_\omega B-\mathsf{Y}\|^2$$ $$\mathsf{X}_\omega=\begin{bmatrix}X \\ \omega I\end{bmatrix} \quad , \quad \mathsf{Y}=\begin{bmatrix}Y \\ 0 \end{bmatrix}$$

W tym przereklamowanym przypadku rozwiązanie jest nadal podawane przez pseudo-odwrotność ale pseudo-odwrotność jest teraz zmieniana, w wyniku czego * gdzie nowa macierz „spektrum osobliwości” ma (odwrotną) przekątną ** (* Nieco obliczenia wymagane do uzyskania tego zostały pominięte ze względu na zwięzłość. Jest to podobne do opisu tutaj dla przypadku . ** Tutaj wpisy wektor jest wyrażany w kategoriach wektora , gdzie wszystkie operacje są wprowadzane).

$$B_\omega = \mathsf{X}^+\mathsf{Y}$$ $$B_\omega = \left(V_0S_\omega^{-2}U^T\right) Y$$ $$\sigma_\omega^2 = \frac{\sigma_0^2+\omega^2}{\sigma_0}$$ $p\leq n$$\sigma_\omega$$\sigma_0$

Teraz w tym problemie nadal możemy formalnie odzyskać „rozwiązanie podstawowe” jako ale to już nie jest prawdziwe rozwiązanie.

$$X_\omega=YB_\omega^+$$

Jednak analogia nadal utrzymuje się w tym, że to „rozwiązanie” ma SVD z pojedynczymi wartościami podanymi powyżej.

$$X_\omega=US_\omega^2V_0^T$$ $\sigma_\omega^2$

Możemy więc wyprowadzić równanie kwadratowe dotyczące pożądanych wartości pojedynczych z do odzyskania wartościami pojedynczymi i parametrem regularyzacji . Rozwiązaniem jest zatem $\sigma_0$$\sigma_\omega^2$$\omega$

$$\sigma_0=\bar{\sigma} \pm \Delta\sigma \quad , \quad \bar{\sigma} = \tfrac{1}{2}\sigma_\omega^2 \quad , \quad \Delta\sigma = \sqrt{\left(\bar{\sigma}+\omega\right)\left(\bar{\sigma}-\omega\right)}$$

---

Poniższe demo Matlaba (przetestowane online przez Octave ) pokazuje, że ta metoda rozwiązania wydaje się działać zarówno w praktyce, jak i teorii. Ostatni wiersz pokazuje, że wszystkie pojedyncze wartości znajdują się w rekonstrukcji , ale nie do końca ustaliłem, który root wybrać ( = vs. ). Dla zawsze będzie to root. Wydaje się, że ogólnie dotyczy to „małego” , podczas gdy dla „dużego” wydaje się, że root. (Poniższa wersja demonstracyjna jest obecnie ustawiona na „dużą” wielkość).$X$$\bar{\sigma}\pm\Delta\sigma$sgn$+$$-$$\omega=0$$+$$\omega$$\omega$$-$

% 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?

Nie mogę powiedzieć, jak solidne jest to rozwiązanie, ponieważ odwrotne problemy są na ogół źle postawione, a rozwiązania analityczne mogą być bardzo kruche. Jednak pobieżne eksperymenty zanieczyszczające szumem Gaussa (tj. Mają pełną rangę porównaniu ze zmniejszoną rangą ) wydają się wskazywać, że metoda jest właściwie zachowana.$B$$p$$n$

Co do problemu 2 (tj nieznany), powyżej daje przynajmniej górna granica w . Aby kwadratowy dyskryminator był nieujemny, musimy mieć $\omega$$\omega$

$$\omega \leq \omega_{\max} = \bar{\sigma}_n = \min[\tfrac{1}{2}\sigma_\omega^2]$$

---

W przypadku dwuznaczności znaku kwadratowego-root poniższy fragment kodu pokazuje, że niezależnie od znaku każdy da to samo rozwiązanie grzbietu przód , nawet jeśli różni się od .$\hat{X}$$B$$\sigma_0$$\mathrm{SVD}[X]$

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