O que está acontecendo aqui, quando uso perda ao quadrado na configuração de regressão logística?

Estou tentando usar a perda ao quadrado para fazer a classificação binária em um conjunto de dados de brinquedo.

Estou usando o mtcarsconjunto de dados, use milhas por galão e peso para prever o tipo de transmissão. O gráfico abaixo mostra os dois tipos de dados do tipo de transmissão em cores diferentes e o limite de decisão gerado por diferentes funções de perda. A perda ao quadrado é que é o rótulo de verdade fundamental (0 ou 1) e é a probabilidade prevista . Em outras palavras, substituo a perda logística pela perda ao quadrado na configuração de classificação, outras partes são iguais. $\sum_i (y_i-p_i)^2$ $y_i$ $p_i$ $p_i=\text{Logit}^{-1}(\beta^Tx_i)$

Para um exemplo de brinquedo com mtcarsdados, em muitos casos, recebi um modelo "semelhante" à regressão logística (veja a figura a seguir, com semente aleatória 0).

Mas em algumas situações (se o fizermos set.seed(1)), a perda ao quadrado parece não estar funcionando bem. O que esta acontecendo aqui? A otimização não converge? A perda logística é mais fácil de otimizar em comparação com a perda ao quadrado? Qualquer ajuda seria apreciada.

Código

d=mtcars[,c("am","mpg","wt")]
plot(d$mpg,d$wt,col=factor(d$am))
lg_fit=glm(am~.,d, family = binomial())
abline(-lg_fit$coefficients[1]/lg_fit$coefficients[3],
       -lg_fit$coefficients[2]/lg_fit$coefficients[3])
grid()

# sq loss
lossSqOnBinary<-function(x,y,w){
  p=plogis(x %*% w)
  return(sum((y-p)^2))
}

# ----------------------------------------------------------------
# note, this random seed is important for squared loss work
# ----------------------------------------------------------------
set.seed(0)

x0=runif(3)
x=as.matrix(cbind(1,d[,2:3]))
y=d$am
opt=optim(x0, lossSqOnBinary, method="BFGS", x=x,y=y)

abline(-opt$par[1]/opt$par[3],
       -opt$par[2]/opt$par[3], lty=2)
legend(25,5,c("logisitc loss","squared loss"), lty=c(1,2))

r machine-learning logistic optimization loss-functions Haitao Du
fonte

Talvez o valor inicial aleatório seja ruim. Por que não escolher um melhor?

whuber

A perda logística do @whuber é convexa, portanto, iniciar não importa. e quanto à perda ao quadrado em p e y? é convexo?

Haitao Du

Não consigo reproduzir o que você descreve. optimdiz que não terminou, é tudo: está convergindo. Você pode aprender muito reexecutando seu código com o argumento adicional control=list(maxit=10000), plotando seu ajuste e comparando seus coeficientes com os originais.

whuber

@amoeba obrigado por seus comentários, eu revisei a pergunta. espero que seja melhor.

Haitao Du

@amoeba Vou revisar a lenda, mas esta declaração não irá corrigir (3)? "Estou usando o conjunto de dados mtcars, use milhas por galão e peso para prever o tipo de transmissão. O gráfico abaixo mostra os dois tipos de dados do tipo de transmissão em cores diferentes e o limite de decisão gerado por diferentes funções de perda."

Haitao Du

Respostas:

Parece que você resolveu o problema em seu exemplo específico, mas acho que ainda vale a pena um estudo mais cuidadoso da diferença entre mínimos quadrados e regressão logística de máxima probabilidade.

Vamos fazer uma anotação. Seja e . Se estivermos fazendo a máxima probabilidade (ou mínima probabilidade de log negativo como estou fazendo aqui), teremos com sendo nossa função de ligação . $L_S(y_i, \hat y_i) = \frac 12(y_i - \hat y_i)^2$ $L_L(y_i, \hat y_i) = y_i \log \hat y_i + (1 - y_i) \log(1 - \hat y_i)$

{\hat{β}}_{L} := {argmin}_{b \in R^{p}} - \sum_{i = 1}^{n} y_{i} \log g^{- 1} (x_{i}^{T} b) + (1 - y_{i}) \log (1 - g^{- 1} (x_{i}^{T} b))

$\hat \beta_L := \text{argmin}_{b \in \mathbb R^p} -\sum_{i=1}^n y_i \log g^{-1}(x_i^T b) + (1-y_i)\log(1 - g^{-1}(x_i^T b))$

g

$g$

Como alternativa, temos como a solução dos mínimos quadrados. Assim, minimiza e da mesma forma para .

{\hat{β}}_{S} : = {argmin}_{b \in R^{p}} \frac{1}{2} \sum_{Eu = 1}^{n} (y_{Eu} - g^{- 1} (x_{Eu}^{T} b))^{2}

$\hat \beta_S := \text{argmin}_{b \in \mathbb R^p} \frac 12 \sum_{i=1}^n (y_i - g^{-1}(x_i^T b))^2$

{\hat{β}}_{S}

$\hat \beta_S$

L_{S}

$L_S$

L_{L}

$L_L$

Seja e as funções objetivas correspondentes à minimização de e respectivamente, como é feito para e . Finalmente, deixe so . Observe que, se estivermos usando o link canônico, temos $f_S$ $f_L$ $L_S$ $L_L$ $\hat \beta_S$ $\hat \beta_L$ $h = g^{-1}$ $\hat y_i = h(x_i^T b)$

h (z) = \frac{1}{1 + e^{- z}} ⟹ h^{'} (z) = h (z) (1 - h (z)) .

$h(z) = \frac{1}{1+e^{-z}} \implies h'(z) = h(z) (1 - h(z)).$

Para regressão logística regular, temos Usando , podemos simplificar isso para então

\frac{\partial f_{L}}{\partial b_{j}} = - \sum_{i = 1}^{n} h^{'} (x_{i}^{T} b) x_{i j} (\frac{y_{i}}{h (x_{i}^{T} b)} - \frac{1 - y_{i}}{1 - h (x_{i}^{T} b)}) .

$\frac{\partial f_L}{\partial b_j} = -\sum_{i=1}^n h'(x_i^T b)x_{ij} \left( \frac{y_i}{h(x_i^T b)} - \frac{1-y_i}{1 - h(x_i^T b)}\right).$

h^{'} = h \cdot (1 - h)

$h' = h \cdot (1 - h)$

\frac{\partial f_{L}}{\partial b_{j}} = - \sum_{i = 1}^{n} x_{i j} (y_{i} (1 - {\hat{y}}_{i}) - (1 - y_{i}) {\hat{y}}_{i}) = - \sum_{i = 1}^{n} x_{i j} (y_{i} - {\hat{y}}_{i})

$\frac{\partial f_L}{\partial b_j} = -\sum_{i=1}^n x_{ij} \left( y_i(1 - \hat y_i) - (1-y_i)\hat y_i\right) = -\sum_{i=1}^n x_{ij}(y_i - \hat y_i)$

\nabla f_{L} (b) = - X^{T} (Y - \hat{Y}) .

$\nabla f_L(b) = -X^T (Y - \hat Y).$

Em seguida, vamos fazer as derivadas secundárias. The Hessian

H_{L} := \frac{\partial^{2} f_{L}}{\partial b_{j} \partial b_{k}} = \sum_{i = 1}^{n} x_{i j} x_{i k} {\hat{y}}_{i} (1 - {\hat{y}}_{i}) .

$H_L:= \frac{\partial^2 f_L}{\partial b_j \partial b_k} = \sum_{i=1}^n x_{ij} x_{ik} \hat y_i (1 - \hat y_i).$ Isso significa que onde . depende dos valores atuais ajustados mas caiu, e é PSD. Portanto, nosso problema de otimização é convexo em .

H_{L} = X^{T} A X

$H_L = X^T A X$

A = diag (\hat{Y} (1 - \hat{Y}))

$A = \text{diag} \left(\hat Y (1 - \hat Y)\right)$

H_{L}

$H_L$

\hat{Y}

$\hat Y$

Y

$Y$

H_{L}

$H_L$

b

$b$

Vamos comparar isso com mínimos quadrados.

\frac{\partial f_{S}}{\partial b_{j}} = - \sum_{Eu = 1}^{n} (y_{Eu} - {\hat{y}}_{Eu}) h^{'} (x_{Eu}^{T} b) x_{Eu j} .

$\frac{\partial f_S}{\partial b_j} = - \sum_{i=1}^n (y_i - \hat y_i) h'(x^T_i b)x_{ij}.$

Isso significa que temos Este é um ponto vital: o gradiente é quase o mesmo, exceto por todos os portanto, basicamente estamos nivelando o gradiente em relação a . Isso tornará a convergência mais lenta.

\nabla f_{S} (b) = - X^{T} UMA (Y - \hat{Y}) .

$\nabla f_S(b) = -X^T A (Y - \hat Y).$

i

$i$

{\hat{y}}_{i} (1 - {\hat{y}}_{i}) \in (0, 1)

$\hat y_i (1 - \hat y_i) \in (0,1)$

\nabla f_{L}

$\nabla f_L$

Para o Hessiano, primeiro podemos escrever

\frac{\partial f_{S}}{\partial b_{j}} = - \sum_{i = 1}^{n} x_{i j} (y_{i} - {\hat{y}}_{i}) {\hat{y}}_{i} (1 - {\hat{y}}_{i}) = - \sum_{i = 1}^{n} x_{i j} (y_{i} {\hat{y}}_{i} - (1 + y_{i}) {\hat{y}}_{i}^{2} + {\hat{y}}_{i}^{3}) .

$\frac{\partial f_S}{\partial b_j} = - \sum_{i=1}^n x_{ij}(y_i - \hat y_i) \hat y_i (1 - \hat y_i) = - \sum_{i=1}^n x_{ij}\left( y_i \hat y_i - (1+y_i)\hat y_i^2 + \hat y_i^3\right).$

Isso nos leva a

H_{S} := \frac{\partial^{2} f_{S}}{\partial b_{j} \partial b_{k}} = - \sum_{i = 1}^{n} x_{i j} x_{i k} h^{'} (x_{i}^{T} b) (y_{i} - 2 (1 + y_{i}) {\hat{y}}_{i} + 3 {\hat{y}}_{i}^{2}) .

$H_S:=\frac{\partial^2 f_S}{\partial b_j \partial b_k} = - \sum_{i=1}^n x_{ij} x_{ik} h'(x_i^T b) \left( y_i - 2(1+y_i)\hat y_i + 3 \hat y_i^2 \right).$

Seja . Agora temos $B = \text{diag} \left( y_i - 2(1+y_i)\hat y_i + 3 \hat y_i ^2 \right)$

H_{S} = - X^{T} A B X .

$H_S = -X^T A B X.$

Infelizmente para nós, não é garantido que os pesos em sejam negativos: se então que é positivo se . Da mesma forma, se então que é positivo quando (também é positivo para mas isso não é possível). Isso significa que não é necessariamente PSD, então não apenas estamos esmagando nossos gradientes, o que tornará o aprendizado mais difícil, mas também confundimos a convexidade do nosso problema. $B$ $y_i = 0$ $y_i - 2(1+y_i)\hat y_i + 3 \hat y_i ^2 = \hat y_i (3 \hat y_i - 2)$ $\hat y_i > \frac 23$ $y_i = 1$ $y_i - 2(1+y_i)\hat y_i + 3 \hat y_i ^2 = 1-4 \hat y_i + 3 \hat y_i^2$ $\hat y_i < \frac 13$ $\hat y_i > 1$ $H_S$

Em suma, não é nenhuma surpresa que a regressão logística dos mínimos quadrados às vezes lute, e no seu exemplo você tem valores ajustados suficientes próximos de ou para que possa ser bem pequeno e, portanto, o gradiente é bastante achatado. $0$ $1$ $\hat y_i (1 - \hat y_i)$

Conectando isso a redes neurais, mesmo que seja apenas uma humilde regressão logística, acho que com a perda ao quadrado, você está experimentando algo como o que Goodfellow, Bengio e Courville estão se referindo em seu livro Deep Learning quando escrevem o seguinte:

Um tema recorrente no design de redes neurais é que o gradiente da função de custo deve ser grande e previsível o suficiente para servir como um bom guia para o algoritmo de aprendizado. As funções que saturam (tornam-se muito planas) minam esse objetivo porque tornam o gradiente muito pequeno. Em muitos casos, isso acontece porque as funções de ativação usadas para produzir a saída das unidades ocultas ou as unidades de saída saturam. A probabilidade de log negativa ajuda a evitar esse problema em muitos modelos. Muitas unidades de saída envolvem uma função exp que pode saturar quando seu argumento é muito negativo. A função de log na função de custo de probabilidade de log negativo desfaz a exp de algumas unidades de saída. Discutiremos a interação entre a função de custo e a escolha da unidade de saída no segundo. 6.2.2

e, em 6.2.2,

Infelizmente, o erro quadrado médio e o erro absoluto médio geralmente levam a resultados ruins quando usados com a otimização baseada em gradiente. Algumas unidades de saída que saturam produzem gradientes muito pequenos quando combinadas com essas funções de custo. Esse é um dos motivos pelos quais a função de custo da entropia cruzada é mais popular que o erro médio quadrático ou o erro absoluto médio, mesmo quando não é necessário estimar uma distribuição inteira . $p(y|x)$

(ambos os trechos são do capítulo 6).

jld
fonte

Eu realmente gosto que você me ajudou a derivar o derivado e o hessian. Vou verificar com mais cuidado amanhã.

Haitao Du 02/02

@ hxd1011, de nada, e obrigado pelo link dessa sua pergunta mais antiga! Eu tenho realmente o que significa que passar por isso com mais cuidado, que foi uma ótima desculpa :)

JLD

Eu li atentamente a matemática e verifiquei com o código. Eu achei Hessian para perda ao quadrado não corresponde à aproximação numérica. Pode verificar? Fico feliz em mostrar o código, se você quiser.

Haitao Du

@ hxd1011 Acabei de passar pela derivação novamente e acho que há um erro de sinal: para , acho que em todos os lugares que tenho , deve ser . Você poderia verificar novamente e me dizer se isso corrige? Muito obrigado pela correção.

H_{S}

$H_S$

y_{i} - 2 (1 - y_{i}) {\hat{y}}_{i} + 3 {\hat{y}}_{i}^{2}

$y_i - 2(1-y_i)\hat y_i + 3 \hat y_i^2$

y_{i} - 2 (\underset{⏟}{1 + y_{i}}) {\hat{y}}_{i} + 3 {\hat{y}}_{i}^{2}

$y_i - 2(\underbrace{1+y_i})\hat y_i + 3 \hat y_i^2$

JLD

@ hxd1011 feliz que corrigiu isso! obrigado novamente para descobrir que

JLD

Agradeço a @whuber e @Chaconne pela ajuda. Especialmente @Chaconne, essa derivação é o que eu queria ter há anos.

O problema está na parte de otimização. Se definirmos a semente aleatória como 1, o BFGS padrão não funcionará. Mas se alterarmos o algoritmo e alterar o número máximo de iterações, ele funcionará novamente.

Como o @Chaconne mencionou, o problema é que a perda ao quadrado da classificação não é convexa e é mais difícil de otimizar. Para adicionar à matemática de @ Chaconne, gostaria de apresentar algumas visualizações sobre a perda logística e a perda ao quadrado.

Mudaremos os dados de demonstração de mtcars, uma vez que o exemplo original de brinquedo possui coeficientes, incluindo a interceptação. Usaremos outro conjunto de dados de brinquedo gerado a partir de , nesse conjunto de dados, parâmetros, o que é melhor para visualização. $3$ mlbench $2$

Aqui está a demo

Os dados são mostrados na figura da esquerda: temos duas classes em duas cores. x, y são dois recursos para os dados. Além disso, usamos a linha vermelha para representar o classificador linear da perda logística, e a linha azul representa o classificador linear da perda ao quadrado.
A figura do meio e a figura à direita mostram o contorno da perda logística (vermelho) e da perda ao quadrado (azul). x, y são dois parâmetros que estamos ajustando. O ponto é o ponto ideal encontrado pelo BFGS.

A partir do contorno, podemos ver facilmente como é mais difícil otimizar a perda ao quadrado: como Chaconne mencionou, é não convexo.

Aqui está mais uma visão do persp3d.

Código

set.seed(0)
d=mlbench::mlbench.2dnormals(50,2,r=1)
x=d$x
y=ifelse(d$classes==1,1,0)

lg_loss <- function(w){
  p=plogis(x %*% w)
  L=-y*log(p)-(1-y)*log(1-p)
  return(sum(L))
}
sq_loss <- function(w){
  p=plogis(x %*% w)
  L=sum((y-p)^2)
  return(L)
}

w_grid_v=seq(-15,15,0.1)
w_grid=expand.grid(w_grid_v,w_grid_v)

opt1=optimx::optimx(c(1,1),fn=lg_loss ,method="BFGS")
z1=matrix(apply(w_grid,1,lg_loss),ncol=length(w_grid_v))

opt2=optimx::optimx(c(1,1),fn=sq_loss ,method="BFGS")
z2=matrix(apply(w_grid,1,sq_loss),ncol=length(w_grid_v))

par(mfrow=c(1,3))
plot(d,xlim=c(-3,3),ylim=c(-3,3))
abline(0,-opt1$p2/opt1$p1,col='darkred',lwd=2)
abline(0,-opt2$p2/opt2$p1,col='blue',lwd=2)
grid()
contour(w_grid_v,w_grid_v,z1,col='darkred',lwd=2, nlevels = 8)
points(opt1$p1,opt1$p2,col='darkred',pch=19)
grid()
contour(w_grid_v,w_grid_v,z2,col='blue',lwd=2, nlevels = 8)
points(opt2$p1,opt2$p2,col='blue',pch=19)
grid()


# library(rgl)
# persp3d(w_grid_v,w_grid_v,z1,col='darkred')

Haitao Du
fonte

Eu não vejo nenhuma não-convexidade na terceira subtrama da sua primeira figura ...

ameba diz Reinstate Monica

@amoeba Eu pensei que o contorno convexo é mais parecido com elipse, duas curvas em forma de U de costas para trás não são convexas, certo?

Haitao Du

Não por que? Talvez seja parte de um contorno maior como uma elipse? Quero dizer, pode muito bem não ser convexo, estou apenas dizendo que não o vejo nesta figura em particular.

Ameba diz Reinstate Monica