Inverser la régression des crêtes: compte tenu de la matrice de réponse et des coefficients de régression, trouver des prédicteurs appropriés

Considérons un problème de régression OLS standard : J'ai des matrices \ Y et \ X et je veux trouver \ B pour minimiser L = \ | \ Y- \ X \ B \ | ^ 2. La solution est donnée par \ 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.$$

Je peux aussi poser un problème "inverse": étant donné $\Y$ et $\B^*$ , trouver $\hat\X$ qui donnerait $\hat\B\approx \B^*$ , c'est-à-dire minimiserait $\|\argmin_\B\{L\}-\B^*\|^2$ . En d'autres termes, j'ai la matrice de réponse $\Y$ et le vecteur de coefficient $\B^*$ et je veux trouver la matrice de prédicteur qui donnerait des coefficients proches de $\B^*$ . Ceci est, bien sûr, également un problème de régression OLS avec la solution

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

Mise à jour de la clarification: Comme @ GeoMatt22 l'a expliqué dans sa réponse, si $\Y$ est un vecteur (c'est-à-dire s'il n'y a qu'une seule variable de réponse), alors $\hat \X$ sera de premier rang et le problème inverse est massivement sous-déterminé. Dans mon cas, $\Y$ est en fait une matrice (c'est-à-dire qu'il existe de nombreuses variables de réponse, c'est une régression multivariée ). Donc $\X$ est $n\times p$ , $\Y$ est $n\times q$ et $\B$ est $p\times q$ .

---

Je suis intéressé à résoudre un problème "inverse" pour la régression des crêtes. A savoir, ma fonction de perte est maintenant

$$L=\|\Y-\X\B\|^2+\mu\|\B\|^2$$ et la solution est $$\hat\B=\argmin_\B\{L\}=(\X^\top \X+\mu\mathbf I)^{-1}\X^\top \Y.$$

Le problème "inverse" est de trouver

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

Encore une fois, j'ai une matrice de réponse $\Y$ et un vecteur de coefficient $\B^*$ et je veux trouver une matrice prédictive qui produirait des coefficients proches de $\B^*$ .

En fait, il existe deux formulations liées:

1. Trouvez $\hat\X$ donné $\Y$ et $\B^*$ et $\mu$ .
2. Recherchez $\hat\X$ et $\hat \mu$ étant donné $\Y$ et $\B^*$ .

L'un ou l'autre a-t-il une solution directe?

---

Voici un bref extrait de Matlab pour illustrer le problème:

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

Ce code renvoie zéro si mu=0mais pas autrement.

Maintenant que la question a convergé vers une formulation plus précise du problème d'intérêt, j'ai trouvé une solution pour le cas 1 (paramètre de crête connu). Cela devrait également aider pour le cas 2 (pas exactement une solution analytique, mais une formule simple et quelques contraintes).

Résumé: Aucune des deux formulations de problèmes inverses n'a de réponse unique. Dans le cas 2 , où le paramètre de crête est inconnu, il existe une infinité de solutions , pour . Dans le cas 1, où est donné, il existe un nombre fini de solutions pour , en raison de l'ambiguïté dans le spectre de valeurs singulières.X ω ω ∈ [ 0 , ω max ] ω X $\mu\equiv\omega^2$$X_\omega$$\omega\in[0,\omega_\max]$$\omega$$X_\omega$

(La dérivation est un peu longue, donc TL, DR: il y a un code Matlab fonctionnel à la fin.)

---

Cas sous-déterminé («OLS»)

Le problème suivant est où , et . X∈ R n × p B∈ R p × q Y∈ R n ×

$$\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}$

Sur la base de la question mise à jour, nous supposons , alors est sous déterminée donnée et . Comme dans la question, nous supposerons le "défaut" (minimumB X Y L 2 B = X + Y X +$n<p<q$$B$$X$$Y$$L_2$ solution -norme) où est la pseudo - de .

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

De la décomposition en valeurs singulières ( SVD ) deX = U S V T = U S 0 V T 0 X + = V S + U T = V 0 S - 1 0 U T X S - 1$X$ , donnée par * le pseudoinverse peut être calculé comme ** (* Les premières expressions utilisent le SVD complet, tandis que les secondes expressions utilisent le SVD réduit. ** Pour plus de simplicité, je suppose que a un rang complet, c'est-à-dire que existe.)

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

Donc, le problème avancé a la solution Pour référence future, je note que , où est le vecteur de valeurs singulières.S 0 = d i a g ( σ 0 ) σ 0 >

$$B\equiv X^+Y=\left(V_0S_0^{-1}U^T\right)Y$$ $S_0=\mathrm{diag}(\sigma_0)$$\sigma_0>0$

Dans le problème inverse, on nous donne etB B X $Y$$B$ . Nous savons que est venu du processus ci - dessus, mais nous ne savons pas . La tâche consiste alors à déterminer le approprié .$B$$X$$X$

Comme indiqué dans la question mise à jour, dans ce cas, nous pouvons récupérer X 0 = Y B +$X$ en utilisant essentiellement la même approche, à savoir en utilisant maintenant la pseudo de .

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

---

Cas surdéterminé (estimateur Ridge)

Dans le cas "OLS", le problème sous-déterminé a été résolu en choisissant la solution de norme minimale , c'est-à-dire que notre solution "unique" a été implicitement régularisée .

Plutôt que de choisir le solution de norme minimale , nous introduisons ici un paramètre pour contrôler "la taille" de la norme, c'est-à-dire que nous utilisons la régression de crête .$\omega$

Dans ce cas, nous avons une série de problèmes de pour , , qui sont donnés par Collecte des différents vecteurs gauche et droit dans cette collection de les problèmes peuvent être réduits au problème "OLS" suivant où nous avons introduit les matrices augmentées k = 1 , … , q min β ‖ X β - y k ‖ 2 + ω 2 ‖ β ‖ 2 B ω = [ β 1 , … , β k$\beta_k$$k=1,\ldots,q$

$$\min_\beta\|X\beta-y_k\|^2+\omega^2\|\beta\|^2$$ min B ‖ X ω B - Y ‖ 2 X ω = [ X ω I $$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}$$

Dans ce cas surdéterminé, la solution est toujours donnée par le pseudo-inverse mais le pseudo-inverse est maintenant changé, résultant en * où la nouvelle matrice "spectre de singularité" a une diagonale (inverse) ** (* Le calcul quelque peu compliqué requis pour dériver ceci a été omis par souci de concision. Il est similaire à l'exposé ici pour leB ω = ( V 0 S - 2 ω U T ) Y σ 2 ω = σ 2 0 + ω

$$B_\omega = \mathsf{X}^+\mathsf{Y}$$ $$B_\omega = \left(V_0S_\omega^{-2}U^T\right) Y$$ p≤nσωσ $$\sigma_\omega^2 = \frac{\sigma_0^2+\omega^2}{\sigma_0}$$ $p\leq n$ . ** Ici les entrées du vecteur est exprimé en termes de vecteur , où toutes les opérations sont en entrée.)$\sigma_\omega$$\sigma_0$

Maintenant, dans ce problème, nous pouvons toujours récupérer formellement une "solution de base" comme mais ce n'est plus une vraie solution.

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

Cependant, l'analogie tient toujours dans le fait que cette "solution" a SVD σ 2

$$X_\omega=US_\omega^2V_0^T$$ avec les valeurs singulières données ci-dessus.$\sigma_\omega^2$

Nous pouvons donc dériver une équation quadratique reliant les valeurs singulières souhaitées aux valeurs singulières récupérables et le paramètre de régularisationσ 2 ω ω σ 0 = ˉ σ ± Δ $\sigma_0$$\sigma_\omega^2$$\omega$ . La solution est alors

$$\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)}$$

---

La démonstration de Matlab ci-dessous (testée en ligne via Octave ) montre que cette méthode de solution semble fonctionner aussi bien en théorie qu'en théorie. La dernière ligne montre que toutes les valeurs singulières de sont dans la reconstruction , mais je n'ai pas complètement compris quelle racine prendre ( = vs ). Pour ce sera toujours la racine . Cela semble généralement tenir pour « petit » , alors que pour « grand » laˉ σ ± Δ σ + - ω = 0 + ω ω $X$$\bar{\sigma}\pm\Delta\sigma$sgn$+$$-$$\omega=0$$+$$\omega$$\omega$$-$ racine semble prendre le dessus. (La démo ci-dessous est actuellement réglée sur un "grand" cas.)

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

Je ne peux pas dire à quel point cette solution est robuste , car les problèmes inverses sont généralement mal posés et les solutions analytiques peuvent être très fragiles. Cependant, des expériences superficielles polluant avec du bruit gaussien (c'est-à-dire qu'il a donc un rang complet rapport à un rang réduit ) semblent indiquer que la méthode se comporte raisonnablement bien.p $B$$p$$n$

Quant au problème 2 (ie inconnu), ce qui précède donne au moins une limite supérieure sur . Pour que le discriminant quadratique soit non négatif, nous devons avoir ω ω ≤ ω max = ˉ σ n = min [ $\omega$$\omega$

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

---

Pour l'ambiguïté du signe racine quadratique, l'extrait de code suivant montre que, indépendamment du signe, tout donnera la même solution de crête avant , même lorsque diffère de .$\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