Reversão de regressão de crista: dada matriz de resposta e coeficientes de regressão, encontre preditores adequados

Considere um problema de regressão OLS padrão : Eu tenho matrizes e \ X e quero encontrar \ B para minimizar L = \ | \ Y- \ X \ B \ | ^ 2. A solução é dada por \ 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.$$

Também posso apresentar um problema "inverso": dado $\Y$ e $\B^*$ , localize $\hat\X$ que produziria $\hat\B\approx \B^*$ , ou seja, minimizaria $\|\argmin_\B\{L\}-\B^*\|^2$ . Em palavras, eu tenho a matriz de resposta $\Y$ e o vetor de coeficiente $\B^*$ e quero encontrar a matriz preditora que produziria coeficientes próximos a $\B^*$ . Obviamente, esse também é um problema de regressão do OLS com a solução

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

Atualização de esclarecimentos: Como @ GeoMatt22 explicou em sua resposta, se $\Y$ é um vetor (ou seja, se houver apenas uma variável de resposta), então esse $\hat \X$ será o primeiro, e o problema inverso será maciçamente indeterminado. No meu caso, $\Y$ é realmente uma matriz (ou seja, existem muitas variáveis ​​de resposta, é uma regressão multivariada ). Então $\X$ é $n\times p$ , $\Y$ é $n\times q$ e $\B$ é $p\times q$ .

---

Estou interessado em resolver um problema "reverso" para regressão de crista. Nomeadamente, minha função de perda agora é

$$L=\|\Y-\X\B\|^2+\mu\|\B\|^2$$ e a solução é $$\hat\B=\argmin_\B\{L\}=(\X^\top \X+\mu\mathbf I)^{-1}\X^\top \Y.$$

O problema "reverso" é encontrar

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

Novamente, eu tenho uma matriz de resposta $\Y$ e um vetor de coeficiente $\B^*$ e quero encontrar uma matriz preditiva que produza coeficientes próximos a $\B^*$ .

Na verdade, existem duas formulações relacionadas:

1. Encontre $\hat\X$ dado $\Y$ e $\B^*$ e $\mu$ .
2. Encontre e dados e .$\hat\X$$\hat \mu$$\Y$$\B^*$

Algum deles tem uma solução direta?

---

Aqui está um breve trecho do Matlab para ilustrar o 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(:)))

Esse código gera zero se mu=0não for o contrário.

Agora que a pergunta convergiu para uma formulação mais precisa do problema de interesse, encontrei uma solução para o caso 1 (parâmetro conhecido da crista). Isso também deve ajudar no caso 2 (não exatamente uma solução analítica, mas uma fórmula simples e algumas restrições).

Resumo: Nenhuma das duas formulações inversas de problemas tem uma resposta única. No caso 2 , onde o parâmetro cume é desconhecido, existem infinitas soluções , para . No caso 1, onde é fornecido, há um número finito de soluções para , devido à ambiguidade no espectro de valor singular.X ω ω ∈ [ 0 , ω max ] ω X $\mu\equiv\omega^2$$X_\omega$$\omega\in[0,\omega_\max]$$\omega$$X_\omega$

(A derivação é um pouco demorada, portanto, TL, DR: existe um código Matlab em funcionamento no final.)

---

Caso sub-determinado ("OLS")

O problema de encaminhamento é que , e . X∈ R n × p B∈ R p x 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}$

Com base na pergunta actualizado, vamos assumir , então é sob determinadas dada e . Como no questão, vamos supor que o "padrão" (mínimo -norm) solução onde é a pseudo-inversa de .B X Y L 2 B = X + Y X +$n<p<q$$B$$X$$Y$$L_2$

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

A partir do valor singular decomposição ( SVD ) de , dado por * o pseudoinverse pode ser calculado como ** (* As primeiras expressões usam o SVD completo, enquanto as segundas expressões usam o SVD reduzido. ** Por uma questão de simplicidade, presumo que tenha uma classificação completa, ou seja, existe.)X = 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$

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

Portanto, o problema de encaminhamento tem a solução Para referência futura, observe que , em que é o vetor de valores singulares.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$

No problema inverso, que são dados e . Sabemos que veio do processo acima, mas não sabemos . A tarefa é determinar o apropriado .B B X $Y$$B$$B$$X$$X$

Como observado na pergunta actualizado, neste caso, pode recuperar utilizando essencialmente a mesma abordagem, isto é, agora usando o pseudoinverse de .X 0 = Y B +$X$

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

---

Caso sobredeterminado (estimador de Ridge)

No caso "OLS", o problema sub-determinado foi resolvido escolhendo a solução de norma mínima , ou seja, nossa solução "única" foi implicitamente regularizada .

Em vez de escolher a solução de norma mínima , aqui apresentamos um parâmetro para controlar "quão pequena" a norma deve ser, ou seja, usamos regressão de crista .$\omega$

Nesse caso, temos uma série de problemas avançados para , , dados por Coletando os diferentes vetores do lado esquerdo e direito em nesta coleção de os problemas podem ser reduzidos para o seguinte problema "OLS" onde introduzimos as matrizes aumentadas 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}$$

Nesse caso sobredeterminado, a solução ainda é dada pelo pseudo-inverso mas o pseudo-inverso agora é alterado, resultando em * onde a nova matriz "espectro de singularidade" possui diagonal (inversa) ** (* O cálculo um tanto envolvido necessário para derivar isso foi omitido por uma questão de brevidade. É semelhante à exposição aqui para o caso . ** Aqui as entradas do vetor é expresso em termos do vetor , onde todas as operações são iniciantes.)

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

Agora, nesse problema, ainda podemos recuperar formalmente uma "solução base" como mas essa não é mais uma solução verdadeira.

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

No entanto, a analogia ainda se mantém em que essa "solução" possui SVD com os valores singulares dados acima.

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

Portanto, podemos derivar uma equação quadrática relacionando os valores singulares desejados aos valores singulares recuperáveis e o parâmetro de regularização . A solução é então $\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)}$$

---

A demonstração do Matlab abaixo (testada on-line via Octave ) mostra que esse método de solução parece funcionar tanto na prática quanto na teoria. A última linha mostra que todos os valores singulares de estão na reconstrução , mas ainda não descobri completamente qual raiz usar ( = vs. ). Para , será sempre a raiz . Isso geralmente é consistente com as "pequenas" , enquanto que para "grande" da raiz parece assumir. (A demonstração abaixo está atualmente definida como "grande".$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?

Não posso dizer quão robusta é essa solução, pois os problemas inversos geralmente são mal colocados e as soluções analíticas podem ser muito frágeis. Contudo, experimentos superficiais poluindo com ruído gaussiano (isto é, para que ele tenha uma classificação completa vs. uma classificação reduzida ) parecem indicar que o método é razoavelmente bem comportado.$B$$p$$n$

Quanto ao problema 2 (isto é, desconhecido), o acima fornece pelo menos um limite superior em . Para que o discriminante quadrático seja não negativo, precisamos ter $\omega$$\omega$

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

---

Para a ambiguidade do sinal da raiz quadrática, o seguinte trecho de código mostra que, independentemente do sinal, qualquer fornecerá a mesma solução cume direta , mesmo quando diferente 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