Coursera ML - A escolha do algoritmo de otimização afeta a precisão da regressão logística multiclasse?

Recentemente, completei o exercício 3 do Machine Learning de Andrew Ng no Coursera usando Python .

Ao concluir inicialmente as partes 1.4 a 1.4.1 do exercício, tive dificuldades para garantir que meu modelo treinado tenha a precisão que corresponde aos 94,9% esperados. Mesmo após a depuração e a garantia de que minhas funções de custo e gradiente estavam livres de erros e que meu código de previsão estava funcionando corretamente, eu ainda estava obtendo apenas 90,3% de precisão. Eu estava usando o algoritmo de gradiente conjugado (CG) em scipy.optimize.minimize.

Por curiosidade, decidi tentar outro algoritmo e usei Broyden – Fletcher – Goldfarb – Shannon (BFGS). Para minha surpresa, a precisão melhorou drasticamente para 96,5% e, portanto, excedeu a expectativa. A comparação desses dois resultados diferentes entre CG e BFGS pode ser vista no meu notebook, sob o cabeçalho Diferença de precisão devido a diferentes algoritmos de otimização .

O motivo dessa diferença de precisão é devido à escolha diferente do algoritmo de otimização? Se sim, alguém poderia explicar o porquê?

Além disso, eu apreciaria muito qualquer revisão do meu código apenas para garantir que não haja um bug em nenhuma das minhas funções que está causando isso.

Obrigado.

EDIT: Aqui abaixo, adicionei o código envolvido na pergunta, a pedido dos comentários que faço nesta página, em vez de encaminhar os leitores para os links para meus cadernos Jupyter.

Funções de custo do modelo:

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def compute_cost_regularized(theta, X, y, lda):
    reg =lda/(2*len(y)) * np.sum(theta[1:]**2) 
    return 1/len(y) * np.sum(-y @ np.log(sigmoid(X@theta)) 
                             - (1-y) @ np.log(1-sigmoid(X@theta))) + reg

def compute_gradient_regularized(theta, X, y, lda):
    gradient = np.zeros(len(theta))
    XT = X.T
    beta = sigmoid(X@theta) - y
    regterm = lda/len(y) * theta
    # theta_0 does not get regularized, so a 0 is substituted in its place
    regterm[0] = 0 
    gradient = (1/len(y) * XT@beta).T + regterm
    return gradient

Função que implementa o treinamento de classificação one-vs-all:

from scipy.optimize import minimize

def train_one_vs_all(X, y, opt_method):
    theta_all = np.zeros((y.max()-y.min()+1, X.shape[1]))
    for k in range(y.min(),y.max()+1):
        grdtruth = np.where(y==k, 1,0)
        results = minimize(compute_cost_regularized, theta_all[k-1,:], 
                           args = (X,grdtruth,0.1),
                           method = opt_method, 
                           jac = compute_gradient_regularized)
        # optimized parameters are accessible through the x attribute
        theta_optimized = results.x
        # Assign thetheta_optimized vector to the appropriate row in the 
        # theta_all matrix
        theta_all[k-1,:] = theta_optimized
    return theta_all

Chamada de função para treinar o modelo com diferentes métodos de otimização:

theta_all_optimized_cg = train_one_vs_all(X_bias, y, 'CG')  # Optimization performed using Conjugate Gradient
theta_all_optimized_bfgs = train_one_vs_all(X_bias, y, 'BFGS') # optimization performed using Broyden–Fletcher–Goldfarb–Shanno

Vemos que os resultados das previsões diferem com base no algoritmo usado:

def predict_one_vs_all(X, theta):
    return np.mean(np.argmax(sigmoid(X@theta.T), axis=1)+1 == y)*100

In[16]: predict_one_vs_all(X_bias, theta_all_optimized_cg)
Out[16]: 90.319999999999993

In[17]: predict_one_vs_all(X_bias, theta_all_optimized_bfgs)
Out[17]: 96.480000000000004

Para quem quiser obter dados para experimentar o código, eles podem ser encontrados no meu Github, conforme vinculado neste post.

Limites de precisão e estabilidade numéricas estão causando dificuldades nas rotinas de otimização.

Você pode ver isso mais facilmente alterando o termo de regularização para 0,0 - não há razão para que isso não funcione em princípio e você não está usando nenhuma engenharia de recursos que precise particularmente. Com a regularização definida como 0,0, você verá os limites de precisão atingidos e tenta obter log de 0 ao calcular a função de custo. As duas rotinas de otimização diferentes são afetadas de maneira diferente, devido ao menor número de pontos de amostra na rota.

Eu acho que com o termo de regularização definido alto, você remove a instabilidade numérica, mas à custa de não ver o que realmente está acontecendo com os cálculos - na verdade, os termos de regularização se tornam dominantes para os difíceis exemplos de treinamento.

Você pode compensar alguns dos problemas de precisão modificando a função de custo:

def compute_cost_regularized(theta, X, y, lda):
    reg =lda/(2*len(y)) * np.sum(theta[1:]**2) 
    return reg - 1/len(y) * np.sum(
      y @ np.log( np.maximum(sigmoid(X@theta), 1e-10) ) 
      + (1-y) @ np.log( np.maximum(1-sigmoid(X@theta), 1e-10) ) )

Também para obter algum feedback durante o treinamento, você pode adicionar

                       options = {
                           'disp': True
                       }

Para a chamada para minimize.

Com essa alteração, você pode tentar com o termo de regularização definido como zero. Quando faço isso, recebo:

predict_one_vs_all(X_bias, theta_all_optimized_cg)
Out[156]:
94.760000000000005
In [157]:

predict_one_vs_all(X_bias, theta_all_optimized_bfgs)
/usr/local/lib/python3.6/site-packages/ipykernel/__main__.py:2: RuntimeWarning: overflow encountered in exp
  from ipykernel import kernelapp as app
Out[157]:
98.839999999999989

O valor de CG de 94,76 parece corresponder muito bem ao resultado esperado - então eu me pergunto se isso foi feito sem regularização. O valor BFGS ainda é "melhor", embora não tenha certeza do quanto confio, dadas as mensagens de aviso durante o treinamento e a avaliação. Para saber se esse resultado aparentemente melhor do treinamento realmente se traduz em uma melhor detecção de dígitos, você precisará medir os resultados em um conjunto de testes de espera.

O CG não converge para o mínimo, assim como o BFGS

Se eu também puder adicionar uma resposta aqui à minha própria pergunta, créditos concedidos a um bom amigo que se ofereceu para examinar meu código. Ele não está na troca de pilha de ciência de dados e não sentiu a necessidade de criar uma conta apenas para postar a resposta, então passou a chance de postar para mim.

Eu também faria referência a @Neil Slater, pois há chances de que sua análise sobre a questão da estabilidade numérica possa explicar isso.

Portanto, a principal premissa por trás da minha solução é:

Sabemos que a função de custo é convexa, o que significa que não possui locais e apenas um mínimo global. Como a previsão usando parâmetros treinados com BFGS é melhor do que aqueles treinados com GC, isso implica que o BFGS convergiu mais próximo do mínimo que o GC. Quer o BFGS tenha convergido ou não para o mínimo global, não podemos dizer com certeza, mas podemos dizer com certeza que é mais próximo do que o CG.

Portanto, se pegarmos os parâmetros que foram treinados usando CG e os passarmos pela rotina de otimização usando BFGS, veremos que esses parâmetros serão otimizados ainda mais, pois o BFGS aproxima tudo ao mínimo. Isso deve melhorar a precisão da previsão e aproximá-la da obtida com o treinamento simples de BFGS.

Aqui abaixo está o código que verifica isso, os nomes das variáveis ​​seguem o mesmo da pergunta:

# Copy the old array over, else only a reference is copied, and the 
# original vector gets modified
theta_all_optimized_bfgs_from_cg = np.copy(theta_all_optimized_cg)

for k in range(y.min(),y.max()+1):
    grdtruth = np.where(y==k, 1,0)
    results = minimize(compute_cost_regularized,theta_all_optimized_bfgs_from_cg[k-1,:], 
                       args = (X_bias,grdtruth,0.1),
                       method = "BFGS", 
                       jac = compute_gradient_regularized, options={"disp":True})
    # optimized parameters are accessible through the x attribute
    theta_optimized = results.x
    # Assign thetheta_optimized vector to the appropriate row in the 
    # theta_all matrix
    theta_all_optimized_bfgs_from_cg[k-1,:] = theta_optimized

Durante a execução do loop, apenas uma das iterações produziu uma mensagem que mostrou um número diferente de zero de iterações de rotina de otimização, o que significa que uma otimização adicional foi executada:

Optimization terminated successfully.
         Current function value: 0.078457
         Iterations: 453
         Function evaluations: 455
         Gradient evaluations: 455

E os resultados foram aprimorados:

In[19]:  predict_one_vs_all(X_bias, theta_all_optimized_bfgs_from_cg)
Out[19]:  96.439999999999998

Ao treinar ainda mais os parâmetros, que foram inicialmente obtidos do CG, por meio de uma execução adicional de BFGS, nós os otimizamos ainda mais para fornecer uma precisão de previsão 96.44%muito próxima da 96.48%obtida diretamente usando apenas BFGS!

Atualizei meu notebook com esta explicação.

É claro que isso levanta mais questões, como por que o CG não funcionou tão bem quanto o BFGS nessa função de custo, mas acho que essas são perguntas destinadas a outro post.