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

16

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.YXβ

L=YXβ2.
β^=argminβ{L}=(XX)+XY.

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

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

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 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×p , Y é n×q e β é p×q .


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

L=YXβ2+μβ2
e a solução é
β^=argminβ{L}=(XX+μI)1XY.

O problema "reverso" é encontrar

X^=argminX{__argminβ{eu}-β__2}=?

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

Na verdade, existem duas formulações relacionadas:

  1. Encontre X^ dado Y e β e μ .
  2. Encontre e dados e .X^μ^Yβ

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.

ameba diz Restabelecer Monica
fonte
Como e são dados, eles não afetam as variações na perda. Portanto, em (1) você ainda está fazendo o OLS. (2) é igualmente simples, porque a perda pode ser arbitrariamente pequena, levando arbitrariamente negativo, dentro dos limites de quaisquer restrições que você compare para lhe impor. Isso reduz você ao caso (1). u uBμμ^
whuber
@whuber Obrigado. Eu acho que não expliquei o suficiente. Considere (1). e são dados (vamos chamá-lo de ), mas preciso encontrar que produza coeficientes de regressão de crista próximos a , ou seja, quero encontrar minimizandoNão vejo por que isso deve ser OLS. μ B * X B * X argmin B { L r i d g e ( X , B ) } - B * 2 .BμBXBX
argminB{Lridge(X,B)}B2.
Ameba diz Reinstate Monica
É como se eu tivesse e quero encontrar tal que esteja próximo a um dado . Não é o mesmo que encontrar . v argmin w f ( v , w ) w argmin v f ( v , w )f(v,W)vargminWf(v,W)Wargminvf(v,W)
Ameba diz Reinstate Monica
A exposição em sua postagem é confusa sobre esse assunto, porque evidentemente você não está realmente usando como uma função de perda. Você poderia elaborar as especificidades dos problemas (1) e (2) no post? L
whuber
2
@ hxd1011 Muitas colunas em X são geralmente chamadas de "regressão múltipla", muitas colunas em Y são geralmente chamadas de "regressão multivariada".
Ameba diz Reinstate Monica

Respostas:

11

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 ωμω2Xωω[0 0,ωmax]ωXω

(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 × q

minB__XB-Y__2
XRn×pBRp×qYRn×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 + Xn<p<qBXYeu2

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 0X

X=vocêSVT=vocêS0 0V0 0T
X+=VS+vocêT=V0 0S0 0-1vocêT
XS0 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 > 0

BX+Y=(V0 0S0 0-1vocêT)Y
S0 0=dEuumag(σ0 0)σ0 0>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 XYBBXX

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

X0 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 .ω

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 ]βkk=1,...,q

minβ__Xβ-yk__2+ω2__β__2
min B ² X ω B - Y 2 X ω = [ X ω I ]
Bω=[β1,...,βk],Y=[y1,...,yk]
minB__XωB-Y__2
Xω=[XωEu],Y=[Y0 0]

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ω=X+Y
Bω=(V0 0Sω-2vocêT)Y
σω2=σ0 02+ω2σ0 0
pnσωσ0 0

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

Xω=YBω+

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

Xω=vocêSω2V0 0T
σω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 σ0 0σω2ω

σ0 0=σ¯±Δσ,σ¯=12σω2,Δσ=(σ¯+ω)(σ¯-ω)

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σ¯±Δσsgn+-ω=0 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ã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.Bpn

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

ωωmax=σ¯n=min[12σω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 .X^Bσ0 0SVD[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
GeoMatt22
fonte
1
+11. Muito obrigado por todo o esforço que você dedicou em responder a essa pergunta e por toda a discussão que tivemos. Isso parece responder totalmente à minha pergunta. Senti que simplesmente aceitar sua resposta não é suficiente neste caso; isso merece muito mais do que dois votos positivos que esta resposta possui atualmente. Felicidades.
ameba diz Reintegrar Monica
@amoeba thanks! Fico feliz que tenha sido útil. Acho que vou postar um comentário na resposta do whuber que você vincula, perguntando se ele acha que é apropriado e / ou se há uma resposta melhor para usar. (Nota ele prefacia sua discussão SVD com a condição de , isto é, um excesso de determinados .)pnX
GeoMatt22
@ GeoMatt22 meu comentário na pergunta original diz que usar pinvnão é uma coisa boa, você concorda?
Haitao Du
1
@ hxd1011 Em geral, você (quase) nunca deseja inverter explicitamente uma matriz numericamente, e isso vale também para o pseudo-inverso. As duas razões pelas quais eu a usei aqui são 1) consistência com as equações matemáticas + código de demonstração da ameba e 2) para o caso de sistemas pouco determinados, as soluções padrão de "barra" do Matlab podem diferir das soluções pinv . Quase todos os casos no meu código podem ser substituídos pelos comandos \ ou / apropriados, que geralmente são os preferidos. (Estes permitem Matlab para decidir o solver mais eficaz direto.)
GeoMatt22
1
@ hxd1011 para esclarecer o ponto 2 do meu comentário anterior, no link do seu comentário sobre a pergunta original: "Se a classificação de A for menor que o número de colunas em A, então x = A \ B não é necessariamente o mínimo solução computacional de norma. Quanto mais caro computacionalmente, x = pinv (A) * B calcula a solução de mínimos quadrados de norma mínima. ".
GeoMatt22