UFRJ Analytica
Published in

UFRJ Analytica

Uso da Decomposição em Valores Singulares na Regressão Linear

O modelo Linear

Para usar uma reta que determine a relação entre duas variáveis, assumimos que a relação entre as variáveis é a seguinte:

Assumindo que N é o número de observações que temos das variáveis, Y e e são vetores de N elementos, X é uma matriz de tamanho N×p, onde p é o número de variáveis usadas para predizer Y, além disso β é um vetor de p elementos.

Assumindo p=3, teríamos a seguinte função de regressão:

Geralmente, em uma regressão, assumimos uma das variáveis xⱼ=1, ou seja, a real equação para yᵢ quando temos p=3:

O coeficiente β₃ é chamado de termo independente.

# Gerando um dataset sintético de variáveis para explicitar a Regressão Linear em Python.

n_amostras = 100
X = np.array([np.random.normal(2, 40, n_amostras), np.random.normal(5, 5, n_amostras)]).T
e = np.random.normal(0, 50, n_amostras)

y = [(2*i[0]) + (0.5*i[1]) + 4 for i in X] + e
y_real = [(2*i[0]) + (0.5*i[1]) + 4 for i in X]
#Usando o método Linear Regression do sklearn, para realizar a Regressão Linear no dataset sintético gerado anteriormente.

reg = LinearRegression().fit(X, y_real)

beta = np.append(reg.coef_, reg.intercept_)
beta
# Calculamos o Erro Quadrático Médio da predição realizada por meio do método 'mean_squared_error' do sklearn.
mean_squared_error(y_real, reg.predict(X))
# Realizaremos um plot 3D para a regressão linear ser melhor visualizada.

fig = go.Figure(data = go.Scatter3d(x=X[:,0], y=X[:,1], z=y, mode='markers'))

x_plot,y_plot = np.meshgrid(X[:,0], X[:,1])

z = beta[0]*x_plot + beta[1]*y_plot + beta[2]

fig.add_trace(
go.Surface(
x=tuple(x_plot),
y=tuple(y_plot),
z=tuple(z),
colorscale=[[0, 'rgb(44, 160, 44)'], [1, 'rgb(44, 160, 44)']]
)
)


fig.show()
Gráfico da regressão linear

A decomposição em Valores Singulares

Dada uma matriz Xₙₓₚ, é possível decompo-la da seguinte forma:

Onde r é o posto da matriz Xₙₓₚ, Uₙₓᵣ e Vᵣₓₚ são matrizes ortogonais*, e Σᵣₓᵣ é uma matriz diagonal.

À essa decomposição é dado o nome de SVD, ou Decomposição em Valores Singulares.

Na decomposição em Valores Singulares, devemos uma especial atenção à matriz Σ, que é da seguinte forma:

O valores na diagonal de Σ, são chamados de Valores Singulares de X.

Além disso, vamos definir a matrix Σd×d, com d ≤ r, como sendo:

Ou seja, é a mesma matriz, só que substituindo todos os valores singulares maiores que d por 0.

O SVD é comumente usado para realizar a tarefa de Redução de Dimensionalidade, no contexto de Ciência de Dados, Redução de Dimensionalidade é, basicamente, reduzir o número de colunas em um dataset.

O SVD é de especial interesse nesse contexto, por conta do que é conhecido como Teorema de Eckart-Young, o teorema diz que dado uma dimensão d, a matriz ~X, definida abaixo como:

A matriz ~X, tem como propriedade posto(~X)=d, e, além disso, temos que ~X é a matriz de posto d, cuja ||X − ~X||^2 é mínima.

Ou seja, ~X é a matriz de posto d que melhor aproxima linearmente X.

*: Se uma matriz A é ortogonal, então:

O determinante de matrizes ortogonais (±1), indicam que matrizes ortogonais não alargam nem encurtam o espaço vetorial, mais especificamente, as transformações produzidas por matrizes ortogonais são exclusivamente rotações.

Olhemos com cuidado para a seguinte matriz, UₙₓᵣΣd×d, podemos ver que ela é uma matriz com N linhas e d colunas, tendo como imagem um hiperplano de dimensão d.

Observe agora, que a matriz ~X, definida por UₙₓᵣΣd×dVᵗᵣₓₚ, tem como imagem uma hiperplano de dimensão d em Rr.

Na verdade, Vᵗ está levando a imagem da matriz UₙₓᵣΣd×d para Rʳ, sem mudar a distância entre os pontos.

Dessa perspectiva, podemos interpretar a matriz UₙₓᵣΣd×d, como uma aproximação de dimensão d para X, originalmente de dimensão r.

À este processo damos o nome de Redução de Dimensionalidade.

# Usando o método linalg.svd do numpy na matriz de features construída anteriormente.
U,s_values,V = np.linalg.svd(X)

S = np.zeros((U.shape[0],V.shape[1]))
np.fill_diagonal(S, s_values)
X_reconstructed = U@S@V# Vemos que o SVD foi capaz de reconstruir completamente a matriz X, de maneira que o valor da norma da subtração de
#X original pelo X reconstruido pelo SVD é meramente erro numérico da linguagem Python.
np.linalg.norm(X - X_reconstructed)
#Realizando uma redução de dimensionalidade de X, com 2 valores singulares, estamos jogando X para a dimensão 2.
X_red_dim = U[:,:2]@S[:2,:2]@V[:2, :]
# Vemos que desta vez a norma não é mais 0, uma vez que a redução de dimensionalidade foi realizada e 2 valores singulares
#não são capazes de reconstruir X completamente.

np.linalg.norm(X - X_red_dim)

PCR — Principal Components Regression

O PCR (Regressão por Componentes Principais, em português), é a ideia de substituir o X na fórmula original da regressão linear pelo seu SVD. Assim, veremos que temos alguns benefícios em certos casos.

Introduzindo uma matriz αᵣₓ₁

A nova equação da regressão linear é:

A partir de tais resultados é possível chegar em uma expressão para a variância de β.
Começamos por escrever ^β em função de ^α usando a decomposição SVD de X:

Ou seja, temos para uma entrada i qualquer de β:

Ou seja, para Var[^βi]:

Desse resultado, é possível perceber que os valores singulares na diagonal da matriz Σ, são inversamente proporcionais à variância de β, ou seja, valores singulares muito pequenos farão os coeficientes da soma crescerem, atribuindo maior incerteza ao resultado de βᵢ.

O dilema do viés e variância

O dilema do viés variância é um problema inerente à qualquer problema de aprendizado supervisionado.
Tal dilema surge associado à complexidade do modelo escolhido, no caso da regressão linear a complexidade se refere ao número de coeficiente usados no modelo, ou o número de variáveis que estamos usando para predizer o target.

A ideia é que modelos de predição mais complexos tem a capacidade de aprender melhor o seu conjunto de dados, no entanto, são pouco generalizáveis, ou seja, é conferida uma grande variância à estimação das componentes que compõem aquele modelo.
Quando um modelo de aprendizado de máquina aprende o conjunto de dados de treino melhor do que deveria, perdendo a capacidade de generalização, dizemos que ocorreu um overfitting.
Além disso, modelos menos complexos possuem menos potencial de aprendizado, no entanto são muito generalizáveis, e por isso, por muitas vezes conseguem predizer situações ainda não vistas melhor do que modelos mais complexos. O problema com essa generalização é que por muitas vezes, tais modelos não são capacidades de reconhecer suficientemente os padrões dentro de conjunto de dados.
Quando um modelo de aprendizado de máquina não é capaz de aprender o conjunto de dados de treino, por conta de uma baixo nível de complexidade do modelo escolhido, dizemos que ocorreu um underfitting.

# Carregando o dataset de assentos em aviões do disponível em https://github.com/julianfaraway/faraway.

seatpos = faraway.datasets.seatpos.load()
seatpos.head()
# Matriz de correlações entre as variáveis do dataset.

seatpos.corr()
# Dividindo as variáveis entre preditoras e a variável que será predita. Em terminologia de Data Science são as variáveis
#feature e a variável target.

X = seatpos.loc[:, seatpos.columns != 'hipcenter'].values
y = seatpos.hipcenter.values
# Usaremos o método StandardScaler para padronizar as variáveis (subtrair a média e dividir pelo desvio padrão).

scaler = StandardScaler()
X = scaler.fit_transform(X)
# Podemos observar os valores singulares na matriz de features usando o método linalg.svd do numpy.
U,s_values,V = np.linalg.svd(X)

S = np.zeros((U.shape[0],V.shape[1]))
np.fill_diagonal(S, s_values)

s_values
# Calculamos o quanto cada valor singular é responsável pela variabilidade na matriz de features.
variancia_explicada = s_values/np.sum(s_values)
variancia_explicada
# Além da variância explicada pelos valores singulares, calculamos essa variância cumulativa.

variancia_explicada_cumulativa = np.cumsum(variancia_explicada)
variancia_explicada_cumulativa
# Abaixo usamos o método k-fold para treinar o nosso modelo, o método k-fold divide o dataset em k (entrada do algoritmo)
#partes iguais, onde 1 dessas k partes é separada para teste e as outras k-1 separadas para teste.
#O algoritmo vai criar k ambientes diferentes de treino e teste, em cada ambiente de treino e teste uma parte diferente
#será selecionada para treino.
# No exemplo abaixo temos dividimos o dataset em 20 partes iguais e, portanto vamos treiná-lo e testá-lo 20 vezes, além
#disso, o Kfold será aplicado à cada número de valores singulares na matriz de features, como temos 8 valores singulares,
#realizaremos 8 kfolds para o teste.

kf = KFold(n_splits=20, shuffle=True, random_state=5)
reg = LinearRegression()

MSE_train = []
MSE_test = []
for i in range(len(s_values)):
X_t = U[:,:i+1]@S[:i+1,:i+1]@V[:i+1, :]
MSE_train_t = []
MSE_test_t = []

for train_index, test_index in kf.split(X_t):
reg.fit(X_t[train_index], y[train_index])
y_pred_train = reg.predict(X_t[train_index])
y_pred_test = reg.predict(X_t[test_index])
MSE_train_t.append(mean_squared_error(y[train_index], y_pred_train))
MSE_test_t.append(mean_squared_error(y[test_index], y_pred_test))

MSE_train.append(np.mean(MSE_train_t))
MSE_test.append(np.mean(MSE_test_t))
#Podemos ver um plot do erro quadrático médio para cada número de valores singulares.

sns.set()
plt.figure(figsize=(14,7))
sns.lineplot(y=MSE_train, x=np.arange(1, 9, 1), label='Erro quadrático médio no dataset de treino')
sns.lineplot(y=MSE_test, x=np.arange(1, 9, 1), label='Erro quadrático médio no dataset de teste')
plt.xticks(ticks=np.arange(1, 9, 1))
plt.xlabel("Quantidade de valores singulares")
plt.ylabel("Erro quadrático médio da previsão gerada pela Regressão Linear")
plt.show()

A partir do gráfico acima, podemos ver que o erro quadrático médio atinge o seu menor valor com 2 valores singulares.
Realizar a redução de dimensionalidade pelo SVD no nosso dataset associou um viés à nossa previsão, uma vez que agora diminuímos o número de variáveis e o nosso modelo tem menos poder de previsão.
No entanto, o menor número de variáveis usadas para predizer o target confere uma menor variância nos resultados da Regressão Linear, o que confere um menor Erro Quadrático Médio, apesar do modelo enviesado.

--

--

Projetos desenvolvidos pela equipe de ciência de dados da maior federal do Brasil.

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store