Regressão Linear - Batch Gradient Descent

Seguindo com a implementação da biblioteca, após ter a base das operações com as matrizes, é o momento de iniciar com o nosso primeiro modelo, as ideias apresentadas aqui são altamente fundamentas nas notas de aula do Andrew NG.

Neste post, vou passar pela teoria (na ordem em que o Andrew apresenta) e depois mostrar como cada pedaço da matemática virou código na minha implementação. No final, respondo a duas perguntas que me incomodaram durante o desenvolvimento — e que acho que vão incomodar quem está começando também.

O problema

Suponha que você tenha um conjunto de casas com área e número de quartos, e queira prever o preço. Os dados são algo como:

Área (ft²)	Quartos	Preço (1000$)
2104	3	400
1600	3	330
2400	3	369
1416	2	232
3000	4	540

Cada casa é um vetor $x^{(i)} \in \mathbb{R}^2$ (área e quartos), e o preço é o alvo $y^{(i)} \in \mathbb{R}$. O objetivo é aprender uma função $h$ (a hipótese) que, dada uma casa nova, retorne uma previsão próxima do preço real.

A hipótese

Começamos com o chute mais simples possível: assumir que $y$ é uma função linear de $x$:

\[h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2\]

Os $\theta_i$ são os parâmetros (ou pesos) que queremos aprender. Para deixar a notação mais limpa, introduzimos a convenção $x_0 = 1$ (o termo de viés, ou intercept), e aí podemos escrever tudo como um produto escalar:

\[h_\theta(x) = \sum_{i=0}^{d} \theta_i x_i = \theta^T x\]

Essa convenção de $x_0 = 1$ é a razão pela qual, na implementação, a primeira coisa que eu faço é “adicionar o viés” — ou seja, acrescentar uma coluna de 1’s na matriz de entrada.

A função custo

Com uma hipótese parametrizada por $\theta$, precisamos de uma forma de medir quão ruim ela está. A escolha clássica é a soma dos erros ao quadrado:

\[J(\theta) = \frac{1}{2} \sum_{i=1}^{n} \bigl(h_\theta(x^{(i)}) - y^{(i)}\bigr)^2\]

Essa é a função de mínimos quadrados, e o fator $\frac{1}{2}$ está ali só para a derivada ficar bonita depois (o 2 do expoente cancela).

Gradiente descendente (o algoritmo LMS)

Queremos encontrar o $\theta$ que minimiza $J(\theta)$. O algoritmo de gradiente descendente faz isso iterativamente: começa com um chute inicial para $\theta$ e vai andando na direção oposta ao gradiente (a direção em que $J$ decresce mais rápido):

\[\theta_j := \theta_j - \alpha \frac{\partial}{\partial \theta_j} J(\theta)\]

onde $\alpha$ é a taxa de aprendizado. Fazendo a derivada para um único exemplo de treino:

\[\frac{\partial}{\partial \theta_j} J(\theta) = \bigl(h_\theta(x) - y\bigr) x_j\]

O que nos dá a famosa regra de atualização LMS (Least Mean Squares), também conhecida como regra de Widrow-Hoff, para um único exemplo de treino:

\[\theta_j := \theta_j + \alpha \bigl(y^{(i)} - h_\theta(x^{(i)})\bigr) x_j^{(i)}\]

Repare no sinal positivo: como o erro é definido como $y - h$ (e não $h - y$), o sinal de menos do gradiente descendente vira um mais. Guarde esse detalhe — ele aparece no código.

Generalizando para vários exemplos com o Batch Gradient Descent:

\[\theta_j := \theta_j + \alpha \sum_{i=1}^{n} \bigl(y^{(i)} - h_\theta(x^{(i)})\bigr) x_j^{(i)} , \text{para todo } j\]

Agrupando tudo em um vetor:

\[\theta := \theta + \alpha \sum_{i=1}^{n} \bigl(y^{(i)} - h_\theta(x^{(i)})\bigr) x^{(i)}\]

Forma matricial (vetorizando tudo)

Escrever a iteração como somas é bom para entender, mas péssimo para implementar. A beleza é que tudo isso pode ser escrito com operações matriciais. Definimos a matriz $X$ para um dataset de exemplo com 3 amostras e 2 features, o vetor $\vec{y}$ com as respostas, e o vetor $\theta$:

\[X = \begin{bmatrix} x^{(1)}_1 & x^{(1)}_2 \\ x^{(2)}_1 & x^{(2)}_2 \\ x^{(3)}_1 & x^{(3)}_2 \end{bmatrix}, \quad \vec{y} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ y^{(3)} \end{bmatrix}, \quad \vec{\theta} = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \end{bmatrix}\]

Matriz $X$ com a coluna de viés:

\[X = \begin{bmatrix} 1 & x^{(1)}_1 & x^{(1)}_2 \\ 1 & x^{(2)}_1 & x^{(2)}_2 \\ 1 & x^{(3)}_1 & x^{(3)}_2 \end{bmatrix}\]

Então a hipótese aplicada a todos os exemplos de uma vez é simplesmente $ h_\theta = X\theta$,

\[h_\theta = \begin{bmatrix} 1 & x^{(1)}_1 & x^{(1)}_2 \\ 1 & x^{(2)}_1 & x^{(2)}_2 \\ 1 & x^{(3)}_1 & x^{(3)}_2 \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \end{bmatrix} = \begin{bmatrix} h_\theta(x^{(1)}) \\ h_\theta(x^{(2)}) \\ h_\theta(x^{(3)}) \end{bmatrix}\]

e o vetor de erros é:

\[\vec{e} = \vec{y} - X\theta\]

Para o cálculo do gradiente temos o seguinte: \(X^T \vec{e}\) :

\[\frac{\partial}{\partial \theta_j} J(\theta) = \sum_{i=1}^{n} \bigl(y^{(i)} - h_\theta(x^{(i)})\bigr) x_j^{(i)}\] \[X = \begin{bmatrix} 1 & x^{(1)}_1 & x^{(1)}_2 \\ 1 & x^{(2)}_1 & x^{(2)}_2 \\ 1 & x^{(3)}_1 & x^{(3)}_2 \end{bmatrix}, \quad X^T = \begin{bmatrix} 1 & 1 & 1 \\ x^{(1)}_1 & x^{(2)}_1 & x^{(3)}_1 \\ x^{(1)}_2 & x^{(2)}_2 & x^{(3)}_2 \\ \end{bmatrix}, \quad \vec{\theta} = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \end{bmatrix}, \quad \vec{e} = \begin{bmatrix} e^{(0)} \\ e^{(1)} \\ e^{(2)} \end{bmatrix}\]

Para \(\theta_0\) temos :

\[\frac{\partial}{\partial \theta_0} J(\theta) = e^{(1)} x^{(1)}_0 + e^{(2)} x^{(2)}_0 + e^{(3)} x^{(3)}_0\]

Para \(\theta_1\) temos :

\[\frac{\partial}{\partial \theta_1} J(\theta) = e^{(1)} x^{(1)}_1 + e^{(2)} x^{(2)}_1 + e^{(3)} x^{(3)}_1\]

Assim temos:

\[X^T = \begin{bmatrix} 1 & 1 & 1 \\ x^{(1)}_1 & x^{(2)}_1 & x^{(3)}_1 \\ x^{(1)}_2 & x^{(2)}_2 & x^{(3)}_2 \\ \end{bmatrix}, \quad \vec{e} = \begin{bmatrix} e^{(0)} \\ e^{(1)} \\ e^{(2)} \end{bmatrix}\]

Onde:

• Cada linha de \(X^T\) corresponde a uma feature
• Multiplicar essa linha pelo vetor de erros é o somatório para o \(\theta_j\) específico
• \(X^T * erro\) calcula o gradiente para todos \(\theta_j\) simultaneamente

A implementação em C++

O header resume a API:

class RegressaoLinear
{
private:
    matriz pesosTheta;
    double taxa_aprendizado;
    int epocas;
    matriz adicionarVies(const matriz& A) const;

    // normalização
    matriz minValores;
    matriz maxValores;
    matriz normalizar(const matriz& X) const;

public:
    RegressaoLinear(double taxa_aprendizado, int epocas);
    ~RegressaoLinear();

    void treinar(const matriz& dadosX, const matriz& respostaY);
    matriz prever(const matriz& dadosX);
    void calcularMinMax(const matriz& X);
};

Vou passar por cada peça.

Adicionando o viés

Essa é a implementação direta da convenção $x_0 = 1$: recebe uma matriz $X$ de tamanho $n \times d$ e devolve uma matriz $n \times (d+1)$ com uma coluna de 1’s na frente.

matriz RegressaoLinear::adicionarVies(const matriz& A) const
{
    int lin = A.getLinhas();
    int col = A.getColunas() + 1;
    matriz comV(lin, col);

    for (int i = 0; i < lin; i++)
    {
        comV.setValor(i, 0, 1.0);
        for (int j = 0; j < A.getColunas(); j++)
        {
            comV.setValor(i, j + 1, A.getValor(i, j));
        }
    }
    return comV;
}

Normalização Min-Max

Antes de treinar, cada coluna de $X$ é reescalada para o intervalo $[0, 1]$:

\[x'_j = \frac{x_j - \min_j}{\max_j - \min_j}\]

matriz RegressaoLinear::normalizar(const matriz& X) const
{
    matriz retorno(X.getLinhas(), X.getColunas());

    for (int i = 0; i < X.getLinhas(); i++)
    {
        for (int j = 0; j < X.getColunas(); j++)
        {
            double numerador   = X.getValor(i, j) - this->minValores.getValor(0, j);
            double denominador = this->maxValores.getValor(0, j) - this->minValores.getValor(0, j);
            double valor = (denominador == 0) ? 0 : numerador / denominador;
            retorno.setValor(i, j, valor);
        }
    }
    return retorno;
}

Os mínimos e máximos de cada coluna são calculados apenas uma vez, durante o treinamento, e ficam guardados na classe para serem reaproveitados no prever().

Atenção: as estatísticas de normalização (minValores e maxValores) devem ser calculadas apenas nos dados de treino e então aplicadas aos dados de teste. Calcular de novo no prever() é um erro clássico de data leakage — o modelo veria informação do conjunto de teste que não deveria ver.

O treinamento: o coração do algoritmo

É aqui que a matemática da seção anterior aparece quase literalmente:

void RegressaoLinear::treinar(const matriz& dadosX, const matriz& respostaY)
{
    this->calcularMinMax(dadosX);
    matriz X_norm = this->normalizar(dadosX);
    matriz X_vies = this->adicionarVies(X_norm);

    int m = X_vies.getLinhas();
    this->pesosTheta = matriz(X_vies.getColunas(), 1);
    this->pesosTheta.preencherAleatorio();

    for (int i = 0; i < this->epocas; i++)
    {
        matriz h_theta = X_vies * this->pesosTheta;           // Xθ
        matriz erro    = respostaY - h_theta;                 // y - Xθ
        matriz custo   = erro.transposta() * erro;            // eᵀe

        if (i % 50 == 0)
        {
            std::cout << "Epoca: " << i
                      << "  Função Custo(J): "
                      << custo.getValor(0, 0) / (2.0 * m)
                      << std::endl;
        }

        matriz gradiente = X_vies.transposta() * erro;        // Xᵀ(y - Xθ)
        double passo     = this->taxa_aprendizado / m;
        this->pesosTheta = this->pesosTheta + (gradiente * passo);
    }
}

Compare linha por linha com as fórmulas:

Fórmula	Código
$h_\theta = X\theta$	X_vies * this->pesosTheta
$\vec{e} = \vec{y} - X\theta$	respostaY - h_theta
$J(\theta) = \frac{1}{2m}\vec{e}^{\,T}\vec{e}$	custo.getValor(0,0) / (2.0*m)
$X^T \vec{e}$	X_vies.transposta() * erro
$\theta := \theta + \alpha X^T \vec{e}$	pesosTheta = pesosTheta + (gradiente * passo)

Note que eu divido por m (o número de exemplos) dentro do passo. Isso transforma a soma numa média — em vez de soma dos erros ao quadrado, estamos usando o erro quadrático médio. É uma escolha comum, porque deixa o custo menos sensível ao tamanho do dataset e torna a mesma taxa de aprendizado útil para conjuntos de tamanhos diferentes.

A previsão

Com $\theta$ aprendido, prever é só aplicar a hipótese em dados novos — lembrando de normalizar e adicionar o viés antes (exatamente como no treino):

matriz RegressaoLinear::prever(const matriz& dadosX)
{
    matriz X_norm   = this->normalizar(dadosX);
    matriz comVies  = this->adicionarVies(X_norm);

    std::cout << "\033[31m" <<  R"(
        __,---.
       /__|o\  )
        `-\ / /
          ,) (,
         //   \\
        {(     )}
  =======""===""===============
          |||||
           |||
            |
    )" << "\033[0m" << std::endl;
    std::cout << "Uirapuru-laranja" << std::endl;

    return (comVies * this->pesosTheta);
}

Dois pontos que me deixaram curioso

Durante o desenvolvimento, deixei um comentário no código com duas coisas para revisar. Já que a resposta delas ilumina bastante o algoritmo, vale a pena destrinchar aqui.

Por que normalizar?

Sem normalização, cada coluna de $X$ pode viver em uma escala completamente diferente. No exemplo das casas, a área está na faixa de 1000 a 3000, e o número de quartos está na faixa de 2 a 4. Isso significa que o gradiente em relação ao $\theta$ da área vai ser muito maior do que o gradiente em relação ao $\theta$ dos quartos — simplesmente porque os valores de $x$ são maiores.

O efeito prático: a função custo vira uma “calha” alongada, e o gradiente descendente fica zigue-zagueando de um lado para o outro em vez de descer direto. Para compensar, você precisaria usar uma taxa de aprendizado minúscula, e o treinamento ficaria lento.

Colocando todas as colunas na mesma escala $[0, 1]$, as curvas de nível de $J$ ficam mais esféricas, o gradiente aponta mais diretamente para o mínimo, e a convergência é bem mais rápida e estável.

Por que $\vec{e}^{\,T}\vec{e}$ é a soma dos erros ao quadrado?

Essa era a pergunta que me incomodava mais. Escrevendo o erro como um vetor coluna:

\[\vec{e} = \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{bmatrix}\]

A transposta $\vec{e}^{\,T}$ é um vetor linha $1 \times n$. Multiplicar $\vec{e}^{\,T}$ (tamanho $1 \times n$) por $\vec{e}$ (tamanho $n \times 1$) dá uma matriz $1 \times 1$ — um escalar:

\[\vec{e}^{\,T} \vec{e} = \begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{bmatrix} = e_1^2 + e_2^2 + \cdots + e_n^2 = \sum_{i=1}^{n} e_i^2\]

É exatamente a soma dos erros ao quadrado! Por isso, no código, erro.transposta() * erro me devolve uma matriz $1 \times 1$, e eu pego getValor(0, 0) para acessar o escalar lá dentro.

Por que isso importa na prática? Porque agora toda a função custo pode ser computada com uma única multiplicação de matrizes, sem loops explícitos. É o tipo de vetorização que deixa implementações de ML rápidas.

Fechamento

O que mais me surpreendeu nessa jornada foi o quão curta a implementação acaba sendo. O loop principal do treino cabe em cinco ou seis linhas de código, e cada uma delas mapeia quase direto para uma linha da matemática . Toda a mágica de “aprender com dados” se resume a: calcule o erro, projete-o de volta no espaço dos parâmetros via $X^T$, e dê um pequeno passo nessa direção. Repita até ficar bom.

A próxima parada é o SGD, Stochastic Gradient Descent