16 votos

Revertir la regresión ridge: dada la respuesta de la matriz y de los coeficientes de regresión, encontrar predictores adecuados

Considere la posibilidad de un estándar de regresión OLS problema$\newcommand{\Y}{\mathbf Y}\newcommand{\X}{\mathbf X}\newcommand{\B}{\boldsymbol\beta}\DeclareMathOperator*{argmin}{argmin}$: he matrices $\Y$ $\X$ y quiero encontrar a $\B$ a minimizar $$L=\|\Y-\X\B\|^2.$$ La solución está dada por $$\hat\B=\argmin_\B\{L\} = (\X^\top\X)^+\X^\top \Y.$$

I can also pose a "reverse" problem: given $\Y$ and $\B^*$, find $\hat\X$ that would yield $\hat\B\approx \B^*$, i.e. would minimize $\|\argmin_\B\{L\}-\B^*\|^2$. In words, I have the response matrix $\Y$ and the coefficient vector $\B^*$ and I want to find the predictor matrix that would yield coefficients close to $\B^*$. This is, of course, also an OLS regression problem with solution $$\hat\X = \argmin_\X\Big\{\|\argmin_\B\{L\}-\B^*\|^2\Big\} = \Y\B^\top(\B\B^\top)^{+}.$$

Clarification update: As @GeoMatt22 explained in his answer, if $\Y$ is a vector (i.e. if there is only one response variable), then this $\hat \X$ will be rank one, and the reverse problem is massively underdetermined. In my case, $\Y$ is actually a matrix (i.e. there are many response variables, it is a multivariate regression). So $\X$ is $n\times p$, $\Y$ is $n\times p$ and $\B$ is $p\veces q$.


I am interested in solving a "reverse" problem for ridge regression. Namely, my loss function now is $$L=\|\Y-\X\B\|^2+\mu\|\B\|^2$$ y la solución es $$\hat\B=\argmin_\B\{L\}=(\X^\top \X+\mu\mathbf I)^{-1}\X^\top \Y.$$

The "reverse" problem is to find $$\hat\X = \argmin_\X\Big\{\|\argmin_\B\{L\}-\B^*\|^2\Big\} = \;?$$

Again, I have a response matrix $\Y$ and a coefficient vector $\B^*$ and I want to find a predictor matrix that would yield coefficients close to $\B^*$.

Actually there are two related formulations:

  1. Find $\hat\X$ given $\Y$ and $\B^*$ and $\mu$.
  2. Encontrar$\hat\X$$\hat \mu$$\Y$$\B^*$.

¿Cualquiera de ellos tiene una solución directa?


Aquí está un breve Matlab extracto para ilustrar el problema:

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

Este código devuelve cero si mu=0 , pero no lo contrario.

11voto

GeoMatt22 Puntos 1290

Ahora que la cuestión ha convergido en una más precisa formulación del problema de interés, he encontrado una solución para el caso 1 (conocido ridge parámetro). Esto también debería ayudar para el caso 2 (no una solución analítica exacta, pero en una fórmula simple y algunas restricciones).

Resumen: Ninguno de los dos problema inverso formulaciones tiene una única respuesta. En el caso 2, donde la cresta del parámetro $\mu\equiv\omega^2$ es desconocida, hay infinitamente muchas soluciones $X_\omega$$\omega\in[0,\omega_\max]$. En el caso 1, donde $\omega$ es dado, hay un número finito de soluciones para $X_\omega$, debido a la ambigüedad en el singular valor de espectro.

(La derivación es un poco largo, así que TL,DR: hay un código de Matlab en el final).


En virtud de determinado Caso ("OLS")

El problema hacia adelante es $$\min_B\|XB-Y\|^2$$ donde $X\in\mathbb{R}^{n\times p}$, $B\in\mathbb{R}^{p\times q}$, y $Y\in\mathbb{R}^{n\times q}$.

Basado en la actualización de la pregunta, vamos a suponer $n<p<q$, lo $B$ es bajo determinadas determinado$X$$Y$. Como en la pregunta, vamos a suponer que el "default" (mínimo $L_2$-norma) solución $$B=X^+Y$$ donde $X^+$ es el pseudoinverse de $X$.

A partir de la descomposición de valor singular (SVD) de $X$, dado por* $$X=USV^T=US_0V_0^T$$ el pseudoinverse puede ser calculada como** $$X^+=VS^+U^T=V_0S_0^{-1}U^T$$ (*Las primeras expresiones de utilizar la totalidad de enfermedad vesicular porcina, mientras que el segundo expresiones de uso de la reducción de la enfermedad vesicular porcina. **Por simplicidad asumo $X$ tiene rango completo, es decir, $S_0^{-1}$ existe).

Para que se transmita problema tiene solución $$B\equiv X^+Y=\left(V_0S_0^{-1}U^T\right)Y$$ Para referencia en el futuro, tengo en cuenta que $S_0=\mathrm{diag}(\sigma_0)$ donde $\sigma_0>0$ es el vector de valores singulares.

En el problema inverso, se nos da $Y$$B$. Sabemos que $B$ provenían del proceso anterior, pero no sabemos $X$. La tarea es, entonces, determinar el $X$.

Como se señaló en la actualización de la pregunta, en este caso podemos recuperar $X$ utilizando esencialmente el mismo enfoque, es decir, $$X_0=YB^+$$ ahora, usando el pseudoinverse de $B$.


La sobre-determinado Caso (estimador Ridge)

En el "OLS", el bajo-determinado problema fue resuelto por la elección de la mínima norma de solución, es decir, nuestra "única" solución fue implícitamente regularizado.

En lugar de elegir la mínima norma de la solución, aquí se introduce un parámetro de $\omega$ a "control de cómo los pequeños" la norma debe ser, es decir, utilizamos la regresión ridge.

En este caso, tenemos una serie de problemas hacia adelante para $\beta_k$, $k=1,\ldots,q$, que se dan por
$$\min_\beta\|X\beta-y_k\|^2+\omega^2\|\beta\|^2$$ La recolección de los diferentes derecha y a la izquierda vectores en $$B_{\omega}=[\beta_1,\ldots,\beta_k] \quad,\quad Y=[y_1,\ldots,y_k]$$ esta colección de problemas pueden reducirse a los siguientes "OLS" problema $$\min_B\|\mathsf{X}_\omega B-\mathsf{Y}\|^2$$ donde hemos introducido la aumentada matrices $$\mathsf{X}_\omega=\begin{bmatrix}X \\ \omega I\end{bmatrix} \quad \quad \mathsf{Y}=\begin{bmatrix}Y \\ 0 \end{bmatrix}$$

En este caso determinado, la solución todavía está dado por la pseudo-inversa $$B_\omega = \mathsf{X}^+\mathsf{Y}$$ pero la pseudo-inversa ahora ha cambiado, lo que resulta en* $$B_\omega = \left(V_0S_\omega^{-2}U^T\right) Y$$ cuando el nuevo "singularidad espectro" de la matriz (inversa) de diagonal** $$ \sigma_\omega^2 = \frac{\sigma_0^2+\omega^2}{\sigma_0} $$ (*Algo que participan de cálculo necesario para obtener esa ha sido omitido en aras de la brevedad. Es similar a la exposición aquí para la $p\leq n$ de los casos. **Aquí las entradas de la $\sigma_\omega$ vector se expresa en términos de la $\sigma_0$ vector, donde todas las operaciones de entrada-sabio.)

Ahora en este problema se puede todavía formalmente recuperar una "solución en base a" como $$X_\omega=YB_\omega^+$$ pero esto no es una verdadera solución ya.

Sin embargo, la analogía todavía tiene en que esta "solución" se ha SVD $$X_\omega=US_\omega^2V_0^T$$ con los valores singulares de a $\sigma_\omega^2$ dado anteriormente.

Así que podemos derivar una ecuación cuadrática relativa deseada de valores singulares de a $\sigma_0$ a la cuantía de valores singulares de a $\sigma_\omega^2$ y el parámetro de regularización $\omega$. La solución es entonces $$\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)}$$


The Matlab demo below (tested online via Octave) shows that this solution method appears to work in practice as well as theory. The last line shows that all the singular values of $X$ are in the reconstruction $\bar{\sigma}\pm\Delta\sigma$, but I have not completely figured out which root to take (sgn = $+$ vs. $-$). For $\omega=0$ it will always be the $+$ root. This generally seems to hold for "small" $\omega$, whereas for "large" $\omega$ the $-$ root seems to take over. (Demo below is set to "large" case currently.)

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

I cannot say how robust this solution is, as inverse problems are generally ill-posed, and analytical solutions can be very fragile. However cursory experiments polluting $B$ with Gaussian noise (i.e. so it has full rank $p$ vs. reduced rank $n$) seem to indicate the method is reasonably well behaved.

As for problem 2 (i.e. $\omega$ unknown), the above gives at least an upper bound on $\omega$. Para la cuadrática discriminante para ser no negativo, debemos tener $$\omega \leq \omega_{\max} = \bar{\sigma}_n = \min[\tfrac{1}{2}\sigma_\omega^2]$$


For the quadratic-root sign ambiguity, the following code snippet shows that independent of sign, any $\hat{X}$ will give the same forward $B$ ridge-solution, even when $\sigma_0$ differs from $\mathrm{CAMPO}[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

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X