Perché il mio gradiente è sbagliato (Coursera, Logistic Regression, Julia)?

Question

Sto provando a fare la regressione logistica da Coursera a Julia, ma non funziona.

Il codice Julia per calcolare il gradiente:

sigmoid(z) = 1/(1 + e^-z) 

hypotesis(theta, x) = sigmoid(scalar(theta' * x)) 

function gradient(theta, x, y) 
    (m, n) = size(x) 
    h = [hypotesis(theta, x[i,:]') for i in 1:m] 
    g = Array(Float64, n, 1) 
    for j in 1:n 
     g[j] = sum([(h[i] - y[i]) * x[i, j] for i in 1:m]) 
    end 
    g 
end 

Se questo gradiente utilizzato produce i risultati errati. Non riesco a capire perché, il codice sembra quello giusto.

Il full Julia script. In questo script, il Theta ottimale calcolato utilizzando la mia implementazione Gradient Descent e utilizzando il pacchetto Optim incorporato, ed i risultati sono diversi.

Answer 1

ho provato confrontando gradient() nel codice del PO con derivata numerica di cost_j() (che è la funzione obiettivo di minimizzazione) usando la seguente procedura

function grad_num(theta, x, y) 
    g = zeros(3) 

    eps = 1.0e-6 
    disp = zeros(3) 

    for k = 1:3 
     disp[:] = theta[:] 

     disp[ k ]= theta[ k ] + eps 
     plus = cost_j(disp, x, y) 
     disp[ k ]= theta[ k ] - eps 
     minus = cost_j(disp, x, y) 

     g[ k ] = (plus - minus)/(2.0 * eps) 
    end 
    return g 
end 

Ma i valori del gradiente ottenuti dai due routine do non sembrano accordo molto bene (almeno per la fase iniziale di minimizzazione) ... così ho derivato manualmente il gradiente di cost_j(theta, x, y), da cui sembra che la divisione per m manca:

#/ OP's code 
# g[j] = sum([ (h[i] - y[i]) * x[i, j] for i in 1:m ]) 

#/ modified code 
    g[j] = sum([ (h[i] - y[i]) * x[i, j] for i in 1:m ])/m 

Perché non sono molto sicuro che se il codice e l'espressione di cui sopra sono veramente corretti, potresti controllarli da soli ...?

Ma in realtà, indipendentemente dal fatto che io utilizzi i gradienti originali o corretti, il programma converge allo stesso valore minimo (0.2034977016, quasi uguale a quello ottenuto da Optim), perché i due gradienti differiscono solo per un fattore moltiplicativo!Poiché la convergenza era molto lento, Ho anche modificato il passo di alpha adattivo seguendo il suggerimento di Vincent (qui usato valori più moderati per l'accelerazione/decelerazione):

function gradient_descent(x, y, theta, alpha, n_iterations) 
    ... 
    c = cost_j(theta, x, y) 

    for i = 1:n_iterations 
     c_prev = c 
     c = cost_j(theta, x, y) 

     if c - c_prev < 0.0 
      alpha *= 1.01 
     else 
      alpha /= 1.05 
     end 
     theta[:] = theta - alpha * gradient(theta, x, y) 
    end 
    ... 
end 

e chiamato questa routine come

optimal_theta = gradient_descent(x, y, [0 0 0]', 1.5e-3, 10^7)[ 1 ] 

La variazione di cost_j rispetto ai passaggi di iterazione viene riportata di seguito.

Answer 2

Il gradiente è corretto (fino a un multiplo scalare, come indicato da @roygvib). Il problema è con la discesa del gradiente.

Se si guardano i valori della funzione di costo durante la discesa del gradiente, vedrete un sacco di NaN, che probabilmente vengono dal esponenziale: abbassare la dimensione del passo (ad esempio, per 1e-5) eviterà l'overflow , ma dovrai aumentare molto il numero di iterazioni (forse a 10_000_000).

Una soluzione migliore (più veloce) sarebbe quella di far variare le dimensioni del passo. Ad esempio, si potrebbe moltiplicare la dimensione passo 1.1 se la funzione di costo migliora dopo una fase (l'optimum appare ancora lontano in questa direzione: possiamo andare più veloce), e dividerlo per 2 caso contrario (siamo andati troppo veloci e abbiamo finito oltre il minimo).

Si potrebbe anche eseguire una ricerca di linea nella direzione del gradiente per trovare la migliore dimensione del passo (ma questa operazione richiede molto tempo e può essere sostituita da approssimazioni, ad esempio la regola di Armijo).

Anche il ridimensionamento delle variabili predittive può essere d'aiuto.