Mostrando postagens com marcador Otimização. Mostrar todas as postagens
Mostrando postagens com marcador Otimização. Mostrar todas as postagens

sexta-feira, 5 de maio de 2017

O que é Regressão Logística?

\(\newcommand{\x}{\textbf{x}}\) \(\newcommand{\R}{\mathcal{R}}\) \(\newcommand{\x}{\textbf{x}}\) \(\newcommand{\y}{\textbf{y}}\)

Ao trabalhar com Regressão Logística nunca me lembro exatamente qual é o formato das funções de erro e se os labels são \(\{0, 1\}\) ou \(\{-1, 1\}\). Resolvi então criar este post detalhando todas as contas feitas para chegar nas diversas expressões usadas para definir a Regressão Logística. Primeiramente, supomos que temos um conjunto de treinamento \(\{(x_i, y_i)\}_{i=1}^N, x_i \in \R^m\) e iremos explorar tanto o caso em que \(y_i \in \{0, 1\}\) como o caso em que \(y_i \in \{-1, 1\}\). Denotamos \(\x\) e \(\y\) as variáveis correspondentes aos padrões de entrada \(x_i\) e aos rótulos \(y_i\).

Primeiramente, o quê é a Regressão Logística? Em outros textos do blog (veja mais em Resumo de mínimos quadrados parte 1 e parte 2) falamos um pouco sobre a Regressão Linear, um método para estimar parâmetros de uma função linear \(f:\R^m \rightarrow \R\). Os coeficientes \(w\) de \(f\) são determinados pelo seguinte problema de otimização.

\begin{equation} \min_{w \in \R^m} ||Xw - y||^2 = \sum_{i=1}^m (w^T x_i - y_i)^2, \end{equation}

Ou seja, a Regressão Linear minimiza a soma (ou média) dos erros ao quadrado (e por isso também é chamada de método dos mínimos quadrados, Least Squares em inglês). Este tipo de erro é adequado quando \(\y\) é uma variável contínua, pois o erro medido leva em conta a distância entre o valor predito \(w^T \x\) e o esperado \(\y\). Não é este o caso quando \(y_i\) é um conjunto discreto (e pequeno) de valores.

A Regressão Logística estima \(f\) para os casos em que \(\y\) é discreto usando probabilidades. Como podem ser frequentes os casos em que dois exemplos \(x_i\) e \(x_j, x_i = x_j\), possuem rótulos diferentes \(y_i \neq y_j\), vamos estimar as probabilidades \(P(\y=1 | x)\) e \(P(\y=0|x)\) (ou \(\y=-1\)). Decidimos o valor de \(f\) com base nas probabilidades calculadas.

\begin{equation} f(\x) = \begin{cases} 1 & P(\y = 1 | \x) > 0.5 \\ 0 \textrm{ ou } -1 & c.c. \end{cases} \end{equation}

Veremos que, assim como a Regressão Linear, a Regressão Logística é um modelo linear em \(\x\) e que podemos computar \(f\) sem calcular explicitamente as probabilidades \(P(\y=y|\x)\).

A motivação que costuma ser encontrada nos livros para a Regressão Logística é a de querermos modelar o log da razão entre as probabilidades como uma função linear. Eu acho isto confuso e nunca pensaria nisto, então bolei a seguinte intuição.

Gostaríamos de modelar \(P(\y=1|\x)\) como um afunção linear. A primeira opção seria usar

\begin{equation} P(\y=1|\x) = w^T \x, \end{equation}

porém isto claramente é um problema, já que \(w^T \x\) pode ser negativo. Podemos resolver isto usando exponenciação.

\begin{equation} P(\y=1|\x) = e^{w^T \x} \end{equation}

Continuamos fazendo uma combinação linear das variáveis de entrada, mas agora no expoente. Isto evita que as probabilidades sejam negativas, mas \(e^{w^T \x}\) pode ser tão grande quanto quisermos. Podemos então dividir isto por \(1+e^{w^T \x}\)! O resultado será sempre menor que 1 e maior que 0, resultando em

\begin{equation} P(\y = 1 | \x) = \frac{ e^{w^Tx} }{1 + e^{w^Tx} }. \end{equation}

Esta função tem o seguinte formato:


By Qef (talk) - Created from scratch with gnuplot, Public Domain, https://commons.wikimedia.org/w/index.php?curid=4310325


Um vetor de pesos \(w \in \R^m\) induz a função de decisão \(f_w(\x) = 1[P(\y = 1 | \x) > 0.5]\), como visto anteriormente. Porém, note que se \(w^T x = 0\), então

\begin{equation} \frac{ e^{w^Tx} }{1 + e^{w^Tx} } = \frac{e^0}{1+e^0} = \frac{1}{2}. \end{equation}

Portanto, \(P(\y = 1 | \x) > 0.5 \Leftrightarrow w^T \x > 0\)! Não precisamos calcular explicitamente \(P(\y = 1 | \x)\). As diferenças entre as duas formulações começam a aparecer neste ponto.

Labels \(\{0, 1\}\)

Neste caso, basta definir a função de decisão como \(f_w(x) = \frac{sign(w^T x) + 1}{2}\) . Se \(sign(w^Tx) = 1\) então \(f_w(x) = 1\), se \(sign(w^Tx) = -1\) \(f_w(x) = 0\).

Podemos definir o vetor de pesos \(w \in \R^m\) como o estimador de máxima verossimilhança da amostra \(\{(x_i, y_i)\}\). Queremos, portanto, resolver o seguinte problema de otimização:

\begin{equation} \max_{w \in \R^m} \prod_{i=1}^{N} P(\y=y_i|x_i) \end{equation}

Ao invés de maximizar a verossimilhança diretamente, é mais fácil maximizar seu log. Isto transforma o produtório em um somatório, resultando na seguinte expressão.

\begin{equation} \begin{aligned} \max_{w \in \R^m} & \sum_{i=1}^{N} \log P(\y=y_i|x_i) = \\ & \sum_{i=1}^{N} y_i \log P(\y=1|x_1) + (1 - y_i) \log P(\y=0|x_i) \end{aligned} \label{eq:min-log-01} \end{equation}

Examinando cada termo \(P(\y=y_i|x_i)\) separadamente para os casos \(y_i=1\) e \(y_i=0\), temos que

\begin{equation} \begin{aligned} \log P(\y=1|x_i) = & \log \left(\frac{e^{w^T x_i}}{1+e^{w^T x_i}} \right) = \\ & \log e^{w^T x_i} - \log (1 + e^{w^T x_i}) = \\ & w^T x_i - \log (1 + e^{w^T x_i}) \end{aligned} \end{equation}

e

\begin{equation} \begin{aligned} \log P(\y=0|x_i) = & \log \left(\frac{1}{1+e^{w^T x_i}} \right) = \\ & \log 1 - \log (1 + e^{w^T x_i}) = \\ & - \log (1 + e^{w^T x_i}) \end{aligned} \end{equation}

Colocando essas expressões de volta na equação anterior obtemos

\begin{equation} \begin{aligned} \max_{w \in \R^m} & \sum_{i=1}^{N} y_i \log P(\y=1|x_1) + (1 - y_i) \log P(\y=0|x_i) = \\ & \sum_{i=1}^{N} y_i (w^T x_i - \log (1+e^{w^T x_i} ) + (1 - y_i) (-\log (1+e^{w^T x_i})) = \\ & \sum_{i=1}^{N} y_i w^T x_i - y_i \log (1+e^{w^T x_i}) - \log(1+e^{w^T x_i}) + y_i\log (1+e^{w^T x_i}) = \\ & \sum_{i=1}^{N} y_i w^T x_i - \log (1+e^{w^T x+i}) \end{aligned} \end{equation}

Em Aprendizado de Máquina a maioria dos problemas de otimização são de minimização. Logo, a forma mais comum da Regresssão logística é

\begin{equation} \min_{w \in \R^m} \sum_{i=1}^{N} -y_i w^T x_i + \log(1+e^{w^T x_i}) \label{eq:regr-log-01} \end{equation}

Labels \(\{-1, 1\}\)

Dependendo da fonte estudada pode-se encontrar a Regressão Logística escrita usando labels \(\{-1, 1\}\) ao invés de \(\{0, 1\}\). Isto facilita na hora de definir a a função de decisão \(f(x) = sign(w^T x)\) e não torna o problema necessariamente mais complicado. As definições relativas às funções de probabilidade são iguais:

\begin{equation} P(\y=1|\x) = \frac{ e^{w^T \x} }{1 + e^{w^T \x} }. \end{equation} \begin{equation} P(\y=-1|\x) = \frac{1}{1 + e^{w^T \x} }. \end{equation}

Porém, note que

\begin{equation} \begin{aligned} P(\y=1|\x) = \frac{e^{w^T \x}}{1 + e^{w^T \x}} \\ \frac{e^{-w^T \x}}{e^{-w^T \x}} \frac{e^{w^T \x}}{1 + e^{w^T \x}} = \\ \frac{1}{e^{-w^T \x}(1+ e^{w^T \x})} = \\ \frac{1}{1 + e^{-w^T \x}} \end{aligned} \end{equation}

Esta expressão é idêntica a de \(P(\y=-1|x)\) a menos do sinal nas exponenciações. Como \(y \in \{-1, 1\}\), podemos definir então uma só expressão para as probabilidades.

\begin{equation} \begin{aligned} P(\y=y|\x) = & \frac{1}{1 + e^{-\y w^T \x}} \end{aligned}\end{equation}

Usamos o mesmo método (Máxima Verossimilhança) para encontrar um \(w\) que torne um conjunto de treinamento \(\{(x_i, y_i)\}\) o mais provável possível, mas desta vez usando a expressão acima para as probabilidades. As passagens abaixo iniciam já na log-Máxima Verossimilhança.

\begin{equation}\begin{aligned} \min_{w \in \R^m} & -\sum_{i=1}^N \log P(\y=y_i | x_i) = \\ & -\sum_{i=1}^N \log \frac{1}{1 + e^{-y_i w^T x_i}} = \\ & -\sum_{i=1}^N - \log (1 + e^{-y_i w^T x_i}) = \\ & \sum_{i=1}^N \log (1 + e^{-y_i w^T x_i}) \end{aligned}\end{equation}

Portanto, mostramos acima as duas formas mais comuns da Regressão Logística e vimos as funções de custo usadas em cada um dos casos. Espero que o post tenha sido útil. Comentários e dúvidas são muito benvindos nas seção de comentários abaixo!

sexta-feira, 20 de novembro de 2015

Resumo de mínimos quadrados - parte 2

Em diversos textos anteriores abordei alguns temas sobre otimização e um texto que recebeu muita atenção foi o de Regressão Linear usando Mínimos Quadrados. No primeiro texto expliquei a motivação inicial da minimização de normas $l_2$ e apresentei alguns exemplos de regressões lineares usando o método de Mínimos quadrados. Neste texto irei um pouco mais a fundo no assunto e apresentarei a formulação básica do problema de otimização, sua solução direta e uma extensão deste método que visa controlar o erro de generalização.

Relembrando a Regressão Linear

O objetivo de uma regressão linear é estimar um vetor $y \in \mathbb{R}^n$ a partir das entradas de uma matrix $X \in \mathbb{R}^{n \times p}, p < n$ utilizando uma função linear $g(x) = x^t w + b$. Tipicamente, $X$ contém em suas linhas observações de um vetor aleatório $\mathbf{X}$ e $y$ contém observações de uma variável aleatória $\mathbf{y}$ que queremos estimar a partir de $\mathbf{X}$. Neste cenário podemos obter facilmente novas observações de $\mathbf{X}$, mas não de $\mathbf{y}$.

Como vimos anteriormente, estimamos o vetor de pesos $w \in \mathbb{R}^p$ minimizando a soma dos erros ao quadrado (Eq. \eqref{custo}).

\begin{equation} \min_w ||Xw - y||_2^2 \label{custo} \end{equation}

Primeiramente,

$$ ||Xw - y||_2^2 = (Xw - y)^T (Xw - y) = \\ (Xw)^TXw - (Xw)^T y - y^T Xw + y^T y = \\ w^T X^T X w - 2 w^T X^T y + ||y||_2^2. $$

Como $||y||_2^2$ é um termo constante (não depende de $w$), podemos ignorá-lo durante a minimização. Sabemos que, no ponto ótimo $w^*$, o gradiente do custo é 0. Logo,

$$ 2X^T X w^* - 2X^T y = 0 \\ 2X^T X w^* = 2X^T y \\ w^* = (X^T X)^{-1} X^T y $$

Portanto, a solução ótima pode ser determinada resolvendo um sistema de equações lineares se $(X^TX)$ for invertível, ou seja, se todas as colunas (ou linhas) desta matrix forem linearmente independentes (LI). Isto ocorre se as colunas/linhas de $X$ forem LI.


Prova:

Suponha, primeiramente, que as colunas de $X$ são LI mas $G = (X^TX)$ possui uma (ou mais) colunas LD e seja $j$ o índice desta coluna. Logo, o elemento $i, j$ de $G$ pode ser escrito como

$$ G_{i,j} = \sum_{k \neq j} \alpha_k G_{i, k} = \sum_{k \neq j} \alpha_k X_i^T X_k = X_i^T (\sum_{k \neq j} \alpha_k X_k) $$

Como $G_{i, j} = X_i^T X_j$, isto implica que $X_j = \sum_{k \neq j} \alpha_k X_k$, uma contradição pois as colunas de $X$ são LI. Logo, se as colunas de $X$ são LI entõ as colunas de $(X^T X)$ também devem ser LI.

Por outro lado, suponha que as colunas de $X$ são LD. Então existe um $j$ tal que $X_j = \sum_{k \neq j} \alpha_k X_k$. Vamos então analisar a coluna $G_j$

$$ G_j = \left( \begin{array}{c} X_1^T X_j \\ \dots \\ X_i^T X_j \\ \dots \end{array} \right) = \left( \begin{array}{c} X_1^T (\sum_{k \neq j} \alpha_k X_k) \\ \dots \\ X_i^T (\sum_{k \neq j} \alpha_k X_k) \\ \dots \end{array} \right) = \left( \begin{array}{c} \sum_{k \neq j} \alpha_k X_1^T X_k \\ \dots \\ \sum_{k \neq j} \alpha_k X_i^T X_k \\ \dots \end{array} \right) = \sum_{k \neq j} \alpha_k G_k $$

Logo, se a $X$ não possui rank completo, $X^T X$ também não. Portanto, uma regressão linear usando mínimos quadrados só possui solução única se a matrix $X$ possuir rank completo.


Apesar de este método poder ser usado para calcular a solução, outros métodos (como a decomposição QR) podem ser mais rápidos na prática e também numericamente mais estáveis, o que significa que eles cometem erros de aproximação menores.

Controlando overfitting

A solução tradicional usando OLS possui duas desvantagens que são especialmente custosas quando consideramos problemas com muitas variáveis.

O primeiro problema é que soluções OLS tendem a privilegiar o encaixe em pontos com erro grande. Dado que um aumento no número de variáveis implica em um aumento nas distâncias entre os pontos, um único outlier pode literalmente estragar o modelo predito. A segunda desvantagem é que a solução OLS nem sempre existe ou é única. Isto ocorre principalmente quando $p > n$, ou seja, quando temos mais variáveis que exemplos.

Estes dois pontos negativos podem ser mitigados adicionando um termo de regularização da solução. O termo mais comum (e matematicamente conveniente) de ser adicionado é a norma $\mathcal{l}_2$ ao quadrado (ou seja, $w^T w$). Este tipo de regularização é conhecido como Ridge Regression e é representado pela equação abaixo (Eq. \eqref{rigde}). Iremos mostrar que este problema sempre tem solução e argumentar que esta formulação privilegia soluções com coeficientes menores, evitando situações como a mostrada anteriormente.

\begin{equation} \min ||Xw - y||_2^2 + \lambda||w||_2^2 \label{ridge} \end{equation}

Como anteriormente, calculamos o gradiente da função de custo em um mínimo $w^*$:

$$ \nabla(w^*) ||Xw^* - y||_2^2 + \lambda||w^*||_2^2 = 0\\ 2X^T X w^* - 2X^T y + 2\lambda w^* = 0 (X^T X + \lambda I)w^* = X^T y. $$

Logo, \begin{equation} w^* = (X^T X + \lambda I)^{-1} X^T y \label{ridge-sol} \end{equation}

A matriz $(X^T X + \lambda I)$ é invertível se ela não possui nenhum autovalor $0$. Primeiramente, sabemos que $X^T X$ é positiva semi-definida, o que implica que todos seus autovalores são não negativos. Logo, tomemos um par de autovalor/vetor $(\alpha, w)$ de $X^T X$. Pela definição,

$$ (X^T X)w = \alpha w \\ (X^T X)w + \lambda w = \alpha w + \lambda w \\ (X^T X + \lambda I)w = (\alpha + \lambda) w $$

Logo, o par $(\alpha + \lambda, w)$ é autovalor/vetor de $X^T X + \lambda I$. Como $X^T X$ é PSD, $\alpha \geq 0$, logo $\alpha + \lambda > 0$, pois $\lambda >0$. Portanto, $X^T X + \lambda I$ é invertível.


O termo $||w||_2^2$ penaliza soluções que possuem componentes com grande magnitude. Isto significa que soluções com custo muito baixo, mas com um ou mais componentes muito maiores que os outros irão ter seu custo aumentado significativamente, enquanto soluções ligeiramente piores mas com magnitudes baixas não irão ter grande aumento no custo. O parâmetro $\lambda$ controla justamente quais tipos de soluções queremos privilegiar. Uma $\lambda$ grande privilegia soluções com baixa magnitude, enquanto um $\lambda$ pequeno interfere menos nos pesos estimados.

Vemos abaixo um exemplo

In [14]:
%matplotlib inline
import numpy as np
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import numpy.linalg

def poly2(inter):
    return np.array([np.ones(len(inter)), [x**2 for x in inter]]).T

def lstsq(X, y):
    inv1 = np.linalg.inv(np.dot(X.T, X))
    return np.dot(inv1, np.dot(X.T, y))

def ridge(X, y, lamb):
    inv1 = np.linalg.inv(np.dot(X.T, X) + lamb)
    return np.dot(inv1, np.dot(X.T, y))

def gen_data():
    inter = np.random.uniform(-50, 50, 21)
    X = poly2(inter)
    y = np.array([x**2 + 56 for x in inter])
    
    inter2 = inter.copy()
    y2 = y.copy()
    inter2[-1] = -15
    y2[-1] = -2000
    X2 = poly2(inter2)
    return inter, X, y, inter2, X2, y2

def plot_curve(inter, alpha, y):
    plt.plot(inter, y, 'o')
    ic = range(-50, 50)
    y_curve = np.dot(poly2(ic), alpha)
    plt.plot(ic, y_curve, 'r-')

inter, X, y, inter2, X2, y2 = gen_data()

plt.figure(figsize=(15,6))
plt.suptitle('Regressão linear')
plt.subplot(121)
alpha = lstsq(X, y)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = lstsq(X2, y2)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)

plt.figure(figsize=(15, 6))
plt.suptitle('Ridge com lambda=300')
plt.subplot(121)
alpha = ridge(X, y, 300)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = ridge(X2, y2, 300)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)


plt.figure(figsize=(15, 6))
plt.suptitle('Ridge com lambda=10')
plt.subplot(121)
alpha = ridge(X, y, 10)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = ridge(X2, y2, 10)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)

No caso sem outliers, a regressão linear "direta" obtém os melhores resultados, tanto visualmente quanto em termos do erro. No caso com outliers, o erro da regressão direta é menor, mas a curva está, visualmente, menos ajustada aos pontos do que as duas versões regularizadas. A versão com $\lambda$ grande "toca" todos os pontos da curva (menos o outlier) apesar de ter um erro quadrático médio pior e mesmo a versão com $\lambda$ pequeno é claramente melhor que a sem regularização. Este exemplo ilustra ambas as propriedades que comentamos antes e motiva a utilização de alguma regularização ao estimar os coeficientes de uma regressão linear. É recomendado, ao usar Ridge Regression, manter todas as variáveis em uma mesma escala (ou ao menos em uma escala parecida), já que a otimização dará preferência a soluções com norma menor. Um procedimento comum é, para cada variável $x_i$, aplicar a transformação $x_i = \frac{x_i - \overline{x_i}}{\sigma_{x_i}}$ para centralizar as variáveis e deixá-las com desvio padrão 1.

Isto conclui este novo post sobre regressão linear. Se houverem detalhes que não ficaram claros ou sugestões de outros assuntos (ou erros, por que não?), não hesite em deixar um comentário ;)

Até a próxima!