Dimensionality Reduction

Redução de Dimensionalidade

Três métodos, três objetivos diferentes. O PCA preserva variância global e é reversível. O t-SNE preserva vizinhança local e distorce todo o resto. O UMAP tenta as duas coisas e falha de forma diferente. Saber qual distorção cada um introduz é o conteúdo inteiro desta página.¹

1. Por que a estatística empurra você a reduzir

Além dos quatro objetivos acima há uma razão estrutural: em dimensão alta, a geometria deixa de se comportar como a intuição de duas e três dimensões espera. O volume foge para a superfície, a bola inscrita praticamente some, direções aleatórias ficam quase ortogonais e — a que de fato quebra as coisas — as distâncias par a par se concentram, então "vizinho mais próximo" pouco a pouco deixa de significar algo.

Esse último efeito é a razão pela qual t-SNE e UMAP recomendam rodar PCA antes quando a entrada tem centenas ou milhares de colunas. Não é só velocidade: é a qualidade das vizinhanças sobre as quais esses métodos são construídos.

2. PCA: duas derivações, uma resposta

Centralize os dados de modo que \(\mu = 0\) e procure uma direção unitária \(\mathbf{w}\) que os resuma bem. Há duas formas de dizer "bem" e a elegância do PCA é que elas coincidem.

Tente achar o melhor ângulo antes de apertar ir para a PC1. As duas barras sempre somam o mesmo número, em qualquer ângulo: é a identidade acima, aplicada à nuvem inteira.

A segunda componente é o autovetor do segundo maior autovalor, necessariamente ortogonal ao primeiro, porque \(\Sigma\) é simétrica. Empilhando as \(k\) primeiras em \(W\): projeção \(Z = XW\) e reconstrução \(\hat X = ZW^\top\).

Escolhendo \(k\)

"""PCA on Iris: what each component costs, and what it is worth.

Four standardized measurements, four components. `explained` is the share of
total variance each one carries, `recon_mse` is what is lost by keeping only
the first k and projecting back, `trustworthiness` is how many of each point's
twelve nearest neighbours in four dimensions are still its neighbours in k, and
`knn_acc` is a 5-NN classifier trained on the k components.

The last column is the one to look at twice. Two components carry 95.8% of the
variance and reconstruct almost perfectly — and they classify *worse* than
three. PCA is unsupervised: it maximizes variance, and nothing guarantees that
the directions with the most variance are the directions that separate the
classes.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.manifold import trustworthiness
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_iris(return_X_y=True)
X_std = StandardScaler().fit_transform(X)        # PCA maximizes variance, and variance has units
cv = StratifiedKFold(5, shuffle=True, random_state=0)

ratios = PCA().fit(X_std).explained_variance_ratio_

print("| `k` | `explained` | `cumulative` | `recon_mse` | `trustworthiness` | `knn_acc` |")
print("|---:|---:|---:|---:|---:|---:|")
for k in (1, 2, 3, 4):
    pca = PCA(k).fit(X_std)
    scores = pca.transform(X_std)
    recon = np.mean((X_std - pca.inverse_transform(scores)) ** 2)
    # the accuracy is cross-validated with PCA *inside* the pipeline, or it would leak
    acc = cross_val_score(
        make_pipeline(StandardScaler(), PCA(k), KNeighborsClassifier(5)), X, y, cv=cv
    ).mean()
    print(f"| {k} | {ratios[k - 1]:.4f} | {ratios[:k].sum():.4f} | {recon:.4f} "
          f"| {trustworthiness(X_std, scores, n_neighbors=12):.3f} | **{acc:.3f}** |")

Iris, padronizado, quatro componentes. As duas primeiras carregam 95,8% da variância e reconstroem com erro de 0,042 contra 0,270 de uma só. Pelos critérios de sempre — "guarde 95% da variância", "procure o cotovelo no gráfico de escarpa" — a resposta é \(k = 2\) e você para por aí.

Agora leia a última coluna. \(k = 2\) classifica a 0,900 e \(k = 3\) a 0,960. A terceira componente carrega 3,7% da variância e vale seis pontos de acurácia.

from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

pipe = make_pipeline(StandardScaler(), PCA(n_components=0.95), MLPClassifier())
# n_components como float = "guarde componentes suficientes para esta fração da variância"
# e, dentro do pipeline, ele é reajustado a cada dobra de treino

Fazendo à mão, uma vez

Vale escrever a decomposição em autovalores uma vez, para o PCA(n_components=2) deixar de ser caixa-preta. Os dois scripts produzem a mesma projeção, a menos do sinal arbitrário de cada componente.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from io import StringIO
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

# Loading Iris dataset
iris = load_iris()

# Transform in dataframe
df = pd.DataFrame(
    data=iris.data,
    columns=['sepal_l', 'sepal_w', 'petal_l', 'petal_w']
)
df['class'] = iris.target_names[iris.target]

X = df.iloc[:,0:4].values
y = df.iloc[:,4].values

# Standardizing
X_std = StandardScaler().fit_transform(X)

# Covariance
cov_mat = np.cov(X_std.T)

# Calculate autovalues and autovectors
eig_vals, eig_vecs = np.linalg.eig(cov_mat)

print('Eigenvectors \n%s' %eig_vecs)
print('\nEigenvalues \n%s' %eig_vals)

# Make a list of (eigenvalue, eigenvector) tuples
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low
eig_pairs.sort(key=lambda x: x[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues
print('Eigenvalues in descending order:')
for i in eig_pairs: print(i[0])

# Sum the cummulative of each eigen value
tot = sum(eig_vals)
var_exp = [(i / tot)*100 for i in sorted(eig_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

n_eigen = [1, 2, 3, 4]

# Plot the cumulative for each eign value
plt.figure(figsize=(6, 4))
plt.bar(n_eigen, var_exp, alpha=0.5, align='center',
    label='individual explained variance')
plt.step(n_eigen, cum_var_exp, where='mid',
    label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()

# Take the only the two firsts eigen values
matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1),
                      eig_pairs[1][1].reshape(4,1)))

print('*' * 10)
print('Reduced to 2-D')
print('Matrix W:\n', matrix_w)

# Calculate the new Y for all samples
Y = X_std.dot(matrix_w)

# Plot the data for the 2 firsts principal components
plt.figure(figsize=(6, 4))
for lab, col in zip(('setosa', 'versicolor', 'virginica'), ('blue', 'red', 'green')):
    plt.scatter(Y[y==lab, 0],
                Y[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend(loc='lower center')
plt.tight_layout()

# Para imprimir na página HTML
buffer = StringIO()
plt.savefig(buffer, format="svg", transparent=True)
print(buffer.getvalue())
plt.close()

import matplotlib.pyplot as plt
import pandas as pd
from io import StringIO
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA as pca
from sklearn.preprocessing import StandardScaler

# Loading Iris dataset
iris = load_iris()

# Transform in dataframe
df = pd.DataFrame(
    data=iris.data,
    columns=['sepal_l', 'sepal_w', 'petal_l', 'petal_w']
)
df['class'] = iris.target_names[iris.target]

X = df.iloc[:,0:4].values
y = df.iloc[:,4].values

# Standardizing
X_std = StandardScaler().fit_transform(X)

sklearn_pca = pca(n_components=2)
Y = sklearn_pca.fit_transform(X_std)

# Plot the data for the 2 firsts principal components
plt.figure(figsize=(6, 4))
for lab, col in zip(('setosa', 'versicolor', 'virginica'), ('blue', 'red', 'green')):
    plt.scatter(Y[y==lab, 0],
                Y[y==lab, 1],
                label=lab,
                c=col)
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend(loc='lower center')
plt.tight_layout()

# Para imprimir na página HTML
buffer = StringIO()
plt.savefig(buffer, format="svg", transparent=True)
print(buffer.getvalue())
plt.close()

3. Onde o PCA falha

O PCA encontra um subespaço linear. Quando os dados moram sobre uma variedade curva — uma espiral, um rocambole suíço, o conjunto de todas as imagens de um objeto girando — plano algum os resume e as primeiras componentes vão alegremente descrever uma direção que atravessa a curva em linha reta.

O número sob cada projeção conta os vizinhos falsos: dos dez pontos mais próximos no mapa, quantos estão a mais de um décimo da curva de distância ao longo dela. O PCA marca 72% — três em cada quatro dos vizinhos aparentes de um ponto estão, na verdade, em outro lugar. O par destacado é o caso extremo: quatro voltas de distância na mola, 0,007 de distância depois de projetar.

Duas coisas merecem atenção antes de seguir.

Primeiro, essa falha é invisível para a trustworthiness, a métrica que a seção 7 recomenda. Ela dá 0,956 ao PCA aqui, porque compara o mapa contra distâncias no espaço ambiente de 3-D — e em 3-D pontos de voltas vizinhas de fato estão perto. O que a projeção destrói é a estrutura intrínseca — a distância ao longo da curva — e isso só se enxerga medindo contra ela. A métrica escolhida decide o que conta como falha.

Segundo, veja o que o t-SNE fez. Ele acerta as vizinhanças quase perfeitamente — 1% de falsos — e rasga a curva em arcos e os espalha. Isso não é um defeito a consertar; é a troca. As vizinhanças estão certas e a disposição dos pedaços não significa nada, que é exatamente por que existem as sete afirmações proibidas.

É essa a abertura para os métodos não lineares e eles mudam a pergunta. O PCA pergunta "quais direções concentram variância?" — pergunta global, resposta linear e reversível. O t-SNE e o UMAP perguntam "quem está perto de quem?" — pergunta local, respondida sacrificando deliberadamente a geometria global.

4. t-SNE: casar distribuições de vizinhança²

Em uma frase: transformar distâncias em probabilidades de vizinhança nos dois espaços e mover os pontos do mapa 2D até que as duas distribuições coincidam.

Passo 1 — vizinhança no espaço original. Para cada ponto \(i\), um núcleo gaussiano sobre os vizinhos dele:

Cada ponto ganha seu próprio \(\sigma_i\), achado por busca binária até a perplexidade \(2^{H(P_i)}\) atingir o valor pedido. Perplexidade ≈ número efetivo de vizinhos. Regiões densas ganham \(\sigma\) pequeno; regiões esparsas, \(\sigma\) grande — uma adaptação automática de escala. O resultado é simetrizado, \(p_{ij} = (p_{j|i} + p_{i|j})/2n\), para que todo ponto contribua para a perda, inclusive pontos isolados.

Passo 2 — vizinhança no mapa, com cauda pesada.

Uma \(t\) de Student com um grau de liberdade (uma Cauchy) e essa é a diferença essencial em relação ao SNE original. Por que a cauda pesada — o problema do amontoamento: o volume de uma bola cresce como \(r^d\), então um ponto em dimensão alta tem espaço para muitos vizinhos a distância moderada e em 2D esse espaço não existe. Com gaussiana dos dois lados, todo vizinho moderado seria forçado ao centro e o mapa colapsaria. A cauda pesada permite que pares moderadamente distantes fiquem muito distantes no mapa a um custo baixo, o que libera espaço para separar os grupos.

O efeito colateral é o que se precisa lembrar: as distâncias entre grupos na figura deixam de ser interpretáveis.

Passo 3 — minimizar a divergência.

A KL é assimétrica e a assimetria é o comportamento inteiro do método. \(p\) alto com \(q\) baixo — vizinhos reais afastados — é punido caro. \(p\) baixo com \(q\) alto — pontos distantes desenhados perto — é quase de graça. Daí a regra:

O que está junto num mapa t-SNE pode não estar junto de verdade.

Veja acontecer. O painel roda o algoritmo de verdade — a busca da perplexidade, o exagero inicial, o gradiente descendente — algumas iterações por quadro.

O conjunto em que vale passar mais tempo é o último. São doze colunas de ruído uniforme, sem nada a encontrar — e o t-SNE lhe entrega grupos separados e arrumados assim mesmo.

5. UMAP: um grafo difuso, desenhado

O UMAP chega a um resultado parecido por um caminho diferente. Em vez de distribuições sobre todos os pares, constrói explicitamente um grafo de \(k\) vizinhos com pesos difusos e depois o desenha.

Passo 1 — grafo local com conectividade garantida.

\(\rho_i\) é a distância ao vizinho mais próximo e \(\sigma_i\) é calibrado para que \(\sum_j w_{i \to j} = \log_2 k\). Subtrair \(\rho_i\) garante que todo ponto tenha ao menos uma aresta de peso 1, então nada fica isolado. As duas visões direcionadas são combinadas por uma união difusa — o análogo da simetrização do t-SNE.

Passo 2 — desenhar o grafo. A perda é uma entropia cruzada binária, não uma KL:

O segundo termo é o que o t-SNE não tem: uma penalidade explícita por aproximar o que está longe. É a origem da preservação global um pouco melhor do UMAP. A otimização usa SGD com amostragem negativa — uma aresta é sorteada e os extremos dela se atraem, alguns pontos aleatórios são sorteados como negativos e se repelem — o que evita a soma sobre todos os pares e é a razão de o UMAP ser mais rápido.

# não é dependência do curso: pip install umap-learn
import umap
Z = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=0).fit_transform(X_std)

O n_neighbors troca detalhe local por estrutura global; o min_dist controla o quanto os pontos podem se empacotar no mapa e é puramente estético — muda a figura sem mudar o que foi aprendido.

O painel abaixo constrói o grafo e depois o desenha, com os mesmos quatro conjuntos. Ligue grafo para ver aquilo que está sendo otimizado: o layout é esse grafo, desembaraçado.

Rode a curva enrolada nos dois painéis e a diferença entre as perdas aparece direto — o segundo termo do UMAP, a penalidade explícita por aproximar o que está longe, mantém a curva conectada por bem mais tempo.

6. Escolhendo — e o que se pode dizer depois

flowchart TD
    A["para que serve a redução?"] --> B{"um modelo vai consumir<br/>a saída?"}
    B -->|sim| P["<b>PCA</b><br/><small>fit/transform honesto, reversível</small>"]
    B -->|"não — é uma figura"| C{"quantas linhas?"}
    C -->|"até ~10 mil"| T["<b>t-SNE</b><br/><small>melhor fidelidade local</small>"]
    C -->|"mais, ou chegam linhas novas"| U["<b>UMAP</b><br/><small>mais rápido, transform aproximado</small>"]
    T --> V["e sempre: PCA para ~50 antes"]
    U --> V

    classDef ok fill:#e6f4ea,stroke:#3fb950,color:#14532d
    classDef q  fill:#eef2f7,stroke:#8b949e,color:#1f2937
    class P,T,U,V ok
    class A,B,C q

A última caixa não é opcional. Com centenas ou milhares de colunas de entrada, rode PCA para cerca de 50 antes do t-SNE ou do UMAP — não só por velocidade, mas porque as vizinhanças que esses métodos consomem são calculadas justamente com as distâncias de que o A3 acima trata.

7. Como avaliar uma redução

Um gráfico não é resultado. Se a redução vai sustentar uma afirmação, ela precisa de métrica. As três usadas:

• Erro de reconstrução — só faz sentido para métodos reversíveis (PCA, autoencoder). Mede o que foi perdido.
• Confiabilidade (trustworthiness) — dos \(k\) vizinhos mais próximos de cada ponto no mapa, quantos eram de fato vizinhos próximos dele no espaço original. Aplicável a qualquer método, inclusive t-SNE e UMAP. Está em sklearn.manifold.trustworthiness.
• Desempenho a jusante — a acurácia de um kNN treinado no espaço reduzido contra a do espaço original. É o teste mais honesto quando a redução é pré-processamento.

"""What t-SNE is good at, and the thing it will do to you if you let it.

Two measurements, each answering one half of the question.

First: does a two-dimensional map keep the neighbourhoods of the original
space? `trustworthiness` counts, for every point, how many of its twelve
nearest neighbours in the map were genuinely among its nearest neighbours in
64 dimensions. This is what t-SNE optimizes and PCA does not, and the gap is
the reason t-SNE exists.

Second: what does t-SNE do to data with no structure at all? The input is
uniform noise in twenty dimensions — no clusters, by construction. A 5-means
silhouette near zero is the correct answer, and the original space gives it.
The t-SNE map does not: the optimization has to put the points somewhere, and
where it puts them looks like groups.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, trustworthiness
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler

NEIGHBOURS, PERPLEXITY = 12, 30

digits = StandardScaler().fit_transform(load_digits(return_X_y=True)[0])
embed = lambda X: TSNE(n_components=2, init="pca", perplexity=PERPLEXITY,
                       random_state=0).fit_transform(X)

print(f"| `digits_1797x64 -> 2d` | `trustworthiness` |")
print("|---|---:|")
for name, Z in (("PCA", PCA(2, random_state=0).fit_transform(digits)),
                ("t-SNE", embed(digits))):
    print(f"| `{name}` | **{trustworthiness(digits, Z, n_neighbors=NEIGHBOURS):.3f}** |")

print()
noise = np.random.default_rng(0).random((600, 20))       # no structure whatsoever
print("| `uniform_noise_600x20` | `silhouette_of_5_means` |")
print("|---|---:|")
for name, Z in (("original_20d", noise), ("t-SNE_map_2d", embed(noise))):
    labels = KMeans(5, n_init=10, random_state=0).fit_predict(Z)
    print(f"| `{name}` | **{silhouette_score(Z, labels):.3f}** |")