Reconhecimento facial do zero: Eigenfaces, PCA e SVD

Este projeto implementa reconhecimento facial a partir da base Olivetti/AT&T sem usar biblioteca pronta de reconhecimento, sem OpenCV, sem face_recognition, sem scikit-learn, sem modelo pre-treinado e sem classificador escondido em API de alto nivel.

A ideia central e construir, do zero, uma geometria para faces:

1. cada imagem vira um vetor em um espaco de dimensao alta;
2. a nuvem de faces e centralizada pela face media;
3. procuramos as direcoes de maior variacao dessa nuvem;
4. projetamos cada rosto nessas direcoes;
5. classificamos uma face nova por proximidade nesse espaco reduzido.

O metodo classico resultante e conhecido como Eigenfaces.

Base de dados

A base esta em:

data/olivetti_faces.mat

Ela contem:

• 400 imagens;
• 40 pessoas;
• 10 imagens por pessoa;
• imagens em escala de cinza;
• resolucao 64 x 64 no arquivo .mat usado aqui.

Fonte historica da base Olivetti/AT&T:

https://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html

Use reconhecimento facial apenas com dados coletados e processados com consentimento.

O problema matematico

Uma imagem 64 x 64 em escala de cinza tem

d = 64 * 64 = 4096

pixels. Podemos empilhar esses pixels em um vetor coluna:

$$x_i = \begin{bmatrix} p_{i1} \\\ p_{i2} \\\ \vdots \\\ p_{id} \end{bmatrix} \in R^d$$

Aqui:

• i identifica a imagem;
• d = 4096;
• p_{ij} e a intensidade do pixel j da imagem i;
• neste projeto, normalizamos os pixels para o intervalo [0, 1].

Se temos n imagens de treino, montamos a matriz:

$$X = \begin{bmatrix} --- x_1^T --- \\\ --- x_2^T --- \\\ \vdots \\\ --- x_n^T --- \end{bmatrix} \in R^{n \times d}$$

No split padrao do projeto:

n = 280
d = 4096

Ou seja, cada rosto vive originalmente em um espaco de 4096 dimensoes. A pergunta e: precisamos mesmo de 4096 coordenadas para distinguir pessoas?

Eigenfaces responde: nao. A maior parte da variacao relevante da base vive em um subespaco muito menor.

Face media

A primeira quantidade fundamental e a media das faces de treino:

$$\mu = {1 \over n} \sum_{i=1}^{n} x_i$$

Cada imagem e entao centralizada:

$$\phi_i = x_i - \mu$$

Em matriz:

$$A = \begin{bmatrix} --- \phi_1^T --- \\\ --- \phi_2^T --- \\\ \vdots \\\ --- \phi_n^T --- \end{bmatrix} = X - 1\mu^T \in R^{n \times d}$$

onde 1 in R^n e o vetor coluna de uns.

Centralizar nao e detalhe cosmetico. Sem isso, a primeira direcao de maior energia tenderia a apontar para o brilho medio/global das imagens, nao para a variacao entre rostos. A PCA deve explicar variacao ao redor da media, nao a media em si.

O que queremos otimizar?

Queremos encontrar uma direcao unitaria w in R^d tal que, ao projetar todas as faces centralizadas nessa direcao, a variancia dos coeficientes projetados seja maxima.

A projecao escalar da face centralizada \phi_i sobre w e:

$$z_i = \phi_i^T w$$

Como \phi_i ja tem media zero ao longo da base, a variancia amostral nessa direcao e proporcional a:

$${1 \over n} \sum_{i=1}^{n} z_i^2 = {1 \over n} \sum_{i=1}^{n} (\phi_i^T w)^2$$

Reescrevendo:

$${1 \over n} \sum_{i=1}^{n} (\phi_i^T w)^2 = {1 \over n} \sum_{i=1}^{n} w^T \phi_i \phi_i^T w = w^T \left( {1 \over n} \sum_{i=1}^{n} \phi_i \phi_i^T \right) w$$

Definimos a matriz de covariancia:

$$C = {1 \over n} A^T A \in R^{d \times d}$$

Entao o problema do primeiro componente principal e:

$$\max_{w \in R^d} \quad w^T C w$$

sujeito a:

$$w^T w = 1$$

Essa restricao e obrigatoria. Sem ela, se w fosse uma solucao, 1000w teria valor objetivo um milhao de vezes maior. O problema ficaria sem sentido.

Deduzindo os autovetores

Para resolver

$$\max_w \quad w^T C w \quad \text{sujeito a} \quad w^T w = 1$$

usamos multiplicadores de Lagrange.

Defina:

$$L(w, \lambda) = w^T C w - \lambda (w^T w - 1)$$

Derivando em relacao a w:

$$\nabla_w L = 2Cw - 2\lambda w$$

No ponto otimo:

$$2Cw - 2\lambda w = 0$$

logo:

$$Cw = \lambda w$$

Ou seja: a melhor direcao w e um autovetor da matriz de covariancia. O valor da variancia nessa direcao e o autovalor correspondente:

$$w^T C w = w^T \lambda w = \lambda w^T w = \lambda$$

Portanto:

• o primeiro componente principal e o autovetor de maior autovalor;
• o segundo e o autovetor de segundo maior autovalor;
• e assim por diante.

De uma direcao para um subespaco

Uma unica direcao nao basta. Queremos um subespaco de dimensao k. Escreva:

$$W = \begin{bmatrix} --- w_1^T --- \\\ --- w_2^T --- \\\ \vdots \\\ --- w_k^T --- \end{bmatrix} \in R^{k \times d}$$

com:

$$W W^T = I_k$$

A projecao de uma face centralizada e:

$$y_i = W\phi_i$$

A energia/variancia preservada pelo subespaco e:

$$\sum_{j=1}^{k} w_j^T C w_j$$

Usando traco:

$$\sum_{j=1}^{k} w_j^T C w_j = \mathrm{tr}(W C W^T)$$

O problema passa a ser:

$$\max_W \quad \mathrm{tr}(W C W^T) \quad \text{sujeito a} \quad W W^T = I_k$$

Como C e simetrica semidefinida positiva, ela admite decomposicao espectral:

$$C = Q \Lambda Q^T$$

onde:

$$\Lambda = \mathrm{diag}(\lambda_1, \lambda_2, ..., \lambda_d)$$

e assumimos:

$$\lambda_1 \ge \lambda_2 \ge ... \ge \lambda_d \ge 0$$

Pelo principio de Rayleigh-Ritz, o maximo da variancia preservada por qualquer subespaco k-dimensional e:

$$\lambda_1 + \lambda_2 + ... + \lambda_k$$

e e atingido escolhendo w_1, ..., w_k como os autovetores associados aos k maiores autovalores. Assim, PCA nao escolhe um subespaco plausivel; ela escolhe o subespaco otimo para preservar variancia sob restricao de ortonormalidade.

Por que nao calcular C diretamente?

No nosso caso:

C in R^{4096 x 4096}

Isso significa mais de 16 milhoes de entradas. Ainda e possivel em uma maquina moderna, mas e desnecessario e piora conforme a imagem cresce.

O script usa SVD diretamente na matriz centralizada:

$$A = U \Sigma V^T$$

onde:

• A in R^{n x d};
• U in R^{n x r} tem colunas ortonormais;
• Sigma in R^{r x r} contem valores singulares;
• V in R^{d x r} tem colunas ortonormais;
• r = rank(A) na SVD fina.

Agora observe:

$$C = {1 \over n} A^T A$$

Substituindo a SVD:

$$A^T A = (U \Sigma V^T)^T (U \Sigma V^T)$$ $$= V \Sigma^T U^T U \Sigma V^T$$

Como U^T U = I:

$$A^T A = V \Sigma^2 V^T$$

Logo:

$$C = V \left({\Sigma^2 \over n}\right) V^T$$

Isso mostra que as colunas de V sao autovetores da covariancia C, e os autovalores sao:

$$\lambda_j = {\sigma_j^2 \over n}$$

onde \sigma_j e o j-esimo valor singular.

No codigo:

_, _, vt = np.linalg.svd(centered, full_matrices=False)
components = vt[:n_components]

vt e V^T. Portanto, as primeiras linhas de vt sao exatamente as direcoes principais transpostas. Cada linha pode ser remodelada como uma imagem 64 x 64; visualmente, essas direcoes sao as eigenfaces.

Reconstrucao e perda de informacao

Se guardamos apenas k componentes, representamos a face por:

$$y = W(x - \mu)$$

Podemos reconstruir uma aproximacao:

$$\hat{x} = \mu + W^T y$$

Como W tem linhas ortonormais, W^T W e o projetor ortogonal no subespaco gerado pelas k eigenfaces:

$$\hat{x} = \mu + W^T W (x - \mu)$$

O erro de reconstrucao e:

$$e = x - \hat{x}$$

e sua norma quadratica:

$$||e||_2^2 = ||x - \mu - W^T W(x - \mu)||_2^2$$

A PCA escolhe o subespaco de dimensao k que minimiza a soma dos erros quadraticos de reconstrucao sobre o conjunto de treino. Equivalentemente, ela maximiza a variancia preservada nas projecoes.

Esse e um fato profundo: maximizar variancia projetada e minimizar erro quadratico de reconstrucao linear sao duas faces do mesmo problema.

Por que isso reconhece rostos?

Rostos humanos alinhados e iluminados de forma parecida nao ocupam arbitrariamente R^4096. Eles formam uma estrutura muito mais concentrada.

Variacoes como:

• formato geral do rosto;
• posicao dos olhos;
• sombra;
• presenca de oculos;
• expressao;
• orientacao leve da cabeca;

criam direcoes dominantes de variacao. Eigenfaces descobre essas direcoes sem receber semantica alguma. O algoritmo nao sabe o que e um olho, nariz ou boca. Ele so ve vetores e variancia.

Depois da projecao, cada imagem vira um ponto em R^k. Se duas imagens sao da mesma pessoa, esperamos que seus pontos fiquem relativamente proximos nesse espaco. Se sao de pessoas diferentes, esperamos distancia maior.

Classificacao por vizinho mais proximo

Depois do treino, para cada imagem de treino x_i, guardamos:

$$y_i = W(x_i - \mu)$$

e seu rotulo c_i, que indica a pessoa.

Para uma imagem nova x, calculamos:

$$y = W(x - \mu)$$

Entao procuramos o exemplo de treino mais perto:

$$i^* = \arg\min_i ||y - y_i||_2$$

e predizemos:

$$\hat{c} = c_{i^*}$$

No codigo, isso aparece como:

distances = np.linalg.norm(
    projected[:, None, :] - model.train_projection[None, :, :],
    axis=2,
)
nearest = distances.argmin(axis=1)
return model.train_labels[nearest], distances[np.arange(len(images)), nearest]

Isso e deliberadamente simples. Nao ha rede neural, nao ha detector facial, nao ha descritor pre-treinado. A geometria e construida a partir da propria base.

Classificacao por centroide

Tambem existe uma alternativa:

python3 face_recognition_from_scratch.py eval --classifier centroid

Para cada classe c, calculamos o centroide no espaco projetado:

$$m_c = {1 \over |I_c|} \sum_{i \in I_c} y_i$$

onde:

$$I_c = \{i : c_i = c\}$$

Para uma face nova:

$$\hat{c} = \arg\min_c ||y - m_c||_2$$

O centroide e mais compacto, pois guarda um vetor por pessoa. O vizinho mais proximo guarda todos os exemplos de treino. Na base atual, o vizinho mais proximo funciona melhor porque cada pessoa tem poucas imagens e as variacoes internas de expressao, oculos e iluminacao podem ser grandes demais para um unico centroide representar bem.

Escolha de k

O parametro k e --components.

Se k e pequeno demais:

• descartamos informacao discriminativa;
• rostos diferentes podem colapsar para pontos parecidos;
• a acuracia cai por subajuste.

Se k e grande demais:

• preservamos ruido e detalhes acidentais;
• a distancia pode passar a depender de iluminacao ou pixels especificos;
• a acuracia pode cair por sobreajuste.

O padrao atual e:

k = 80

Com split fixo de 7 imagens por pessoa para treino e 3 para teste, o resultado observado foi:

Imagens de treino: 280
Imagens de teste: 120
Componentes PCA: 80
Classificador: nearest-neighbor
Acuracia: 89.17%

Esse numero nao deve ser tratado como verdade universal. Ele depende do split, do pre-processamento, do numero de componentes e da metrica de distancia.

Ligacao com algebra linear geometrica

PCA pode ser vista como escolher uma base ortonormal adaptada aos dados.

A base canonica de R^4096 olha para pixels individuais:

$$e_1, e_2, ..., e_{4096}$$

Cada e_j pergunta: "qual e a intensidade do pixel j?"

Eigenfaces troca essa base por direcoes globais:

$$w_1, w_2, ..., w_k$$

Cada w_j pergunta: "quanto esta face se parece com este padrao global de variacao?"

Por isso uma coordenada PCA nao e um pixel. Ela e um coeficiente de semelhanca com uma direcao visual aprendida.

Se w_j for remodelado para 64 x 64, ele pode parecer uma face fantasmagorica com regioes claras e escuras. Matematicamente, isso e apenas um autovetor da covariancia. Visualmente, virou uma eigenface.

Normalizacao e escalas

Os pixels originais sao inteiros. O codigo converte para float64 e divide por 255:

faces = mat["faces"].T.astype(np.float64) / 255.0

Isso coloca todas as intensidades em [0, 1].

Como a distancia euclidiana e sensivel a escala, essa normalizacao e importante. Se os pixels estivessem em [0, 255], as distancias seriam 255 vezes maiores. A classificacao por argmin nao mudaria necessariamente se todos os vetores fossem escalados igualmente, mas a escala numerica dos limiares e distancias ficaria menos interpretavel.

O algoritmo completo

Entrada:

X: matriz de imagens de treino
c: rotulos das imagens de treino
k: numero de componentes

Treino:

1. Calcule mu = media das linhas de X.
2. Calcule A = X - mu.
3. Calcule A = U Sigma V^T.
4. Guarde W = primeiras k linhas de V^T.
5. Projete Y = A W^T.
6. Guarde Y e os rotulos c.

Predicao:

1. Receba uma imagem x.
2. Centralize: phi = x - mu.
3. Projete: y = W phi.
4. Ache i* = argmin_i ||y - y_i||_2.
5. Retorne c_i*.

Mapeamento direto para o codigo

Carregamento:

mat = loadmat(path)
faces = mat["faces"].T.astype(np.float64) / 255.0
labels = np.repeat(np.arange(faces.shape[0] // IMAGES_PER_PERSON), IMAGES_PER_PERSON)

Media e centralizacao:

mean_face = x_train.mean(axis=0)
centered = x_train - mean_face

SVD:

_, _, vt = np.linalg.svd(centered, full_matrices=False)
components = vt[:n_components]

Projecao:

train_projection = centered @ components.T

Predicao:

projected = (images - model.mean_face) @ model.components.T

Distancia:

distances = np.linalg.norm(
    projected[:, None, :] - model.train_projection[None, :, :],
    axis=2,
)

Comandos

Treinar e avaliar:

python3 face_recognition_from_scratch.py eval

Mudar o numero de componentes:

python3 face_recognition_from_scratch.py eval --components 40
python3 face_recognition_from_scratch.py eval --components 120

Comparar com classificador por centroide:

python3 face_recognition_from_scratch.py eval --classifier centroid

Treinar e salvar o modelo:

python3 face_recognition_from_scratch.py train --model models/eigenfaces.npz

Predizer uma imagem da propria base:

python3 face_recognition_from_scratch.py predict --model models/eigenfaces.npz --index 0

Salvar uma imagem da base como PNG:

python3 face_recognition_from_scratch.py export-image --index 0 --out outputs/face_0.png

O que foi usado

• numpy: arrays, SVD e operacoes vetoriais.
• scipy: apenas para ler o arquivo .mat.
• matplotlib: apenas para exportar imagens de exemplo.

O reconhecimento em si e implementado no arquivo:

face_recognition_from_scratch.py

Licenca

Este projeto esta sob a licenca MIT. Veja LICENSE.

Limites do metodo

Eigenfaces e historicamente importante e matematicamente elegante, mas tem limites claros:

• depende de faces razoavelmente alinhadas;
• sofre com mudancas fortes de iluminacao;
• sofre com pose muito diferente;
• nao detecta rosto na imagem;
• nao lida naturalmente com multiplas faces;
• nao aprende invariancias complexas como redes neurais modernas.

Ainda assim, como primeiro sistema completo de reconhecimento, ele e excelente: mostra como uma ideia limpa de algebra linear transforma pixels crus em uma representacao discriminativa.

Resumo em uma linha

Eigenfaces encontra uma base ortonormal que captura as direcoes de maior variancia das faces de treino; cada rosto vira um vetor curto de coeficientes nessa base, e o reconhecimento acontece por proximidade geometrica nesse espaco reduzido.