<p>For generalized linear models (i.e. logistic regression, ridge regression, poisson regression),
you can efficiently tune many regularization hyperparameters
using exact derivatives and approximate leave-one cross-validation.</p>
<p>But don't stop at just the gradient, compute the full hessian and use a second-order optimizer -- it's
both more efficient and robust.</p>
<p>sklearn doesn't currently have this functionality, but there are other tools available that can do it.</p>
<p>For example, here's how you can use the python package <a href="https://buildingblock.ai/" rel="nofollow noreferrer">bbai</a> to fit the
hyperparameter for ridge regularized logistic regression to maximize the log likelihood of the
approximate leave-one-out cross-validation of the training data set for the Wisconsin Breast Cancer Data Set.</p>
<p><strong>Load the data set</strong></p>
<pre><code>from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler

data = load_breast_cancer()
X = data['data']
X = StandardScaler().fit_transform(X)
y = data['target']
</code></pre>
<p><strong>Fit the model</strong></p>
<pre><code>import bbai.glm

model = bbai.glm.LogisticRegression()
# Note: it automatically fits the C parameter to minimize the error on
# the approximate leave-one-out cross-validation.
model.fit(X, y)
</code></pre>
<p>Because it uses use both the gradient and hessian with efficient exact formulas
(no automatic differentiation), it can dial into an exact hyperparameter quickly with only a few
evaluations.</p>
<p>YMMV, but when I compare it to sklearn's LogisticRegressionCV with default parameters, it runs
in a fraction of the time.</p>
<pre><code>t1 = time.time()
model = bbai.glm.LogisticRegression()
model.fit(X, y)
t2 = time.time()
print('***** approximate leave-one-out optimization')
print('C = ', model.C_)
print('time = ', (t2 - t1))

from sklearn.linear_model import LogisticRegressionCV
print('***** sklearn.LogisticRegressionCV')
t1 = time.time()
model = LogisticRegressionCV(scoring='neg_log_loss', random_state=0)
model.fit(X, y)
t2 = time.time()
print('C = ', model.C_[0])
print('time = ', (t2 - t1))
</code></pre>
<p>Prints</p>
<pre><code>***** approximate leave-one-out optimization
C =  0.6655139682151275
time =  0.03996014595031738
***** sklearn.LogisticRegressionCV
C =  0.3593813663804626
time =  0.2602980136871338
</code></pre>
<h2>How it works</h2>
<p>Approximate leave-one-out cross-validation (ALOOCV) is a close approimxation to leave-one-out
cross-validation that's much more efficient to evaluate for generalized linear models.</p>
<p>It first fits the regularized model. Then uses a single step of Newton's algorithm to approximate what
the model weights would be when we leave a single data point out. If the regularized cost function for
the generalized linear model is represented as</p>
<img src="https://i.stack.imgur.com/zZ0zp.png" width="300"/>
<p>Then the ALOOCV can be computed as</p>
<img src="https://i.stack.imgur.com/JICyA.png" width="350"/>
<p>where</p>
<img src="https://i.stack.imgur.com/FLMS7.png" width="300"/>
<p>(Note: H represents the hessian of the cost function at the optimal weights)</p>
<p>For more background on ALOOCV, you can check out this <a href="https://buildingblock.ai/logistic-regression-guide#approximate-leave-one-out-cross-validation" rel="nofollow noreferrer">guide</a>.</p>
<p>It's also possible to compute exact derivatives for ALOOCV which makes it efficient to optimize.</p>
<p>I won't put the derivative formulas here as they are quite involved, but see the paper
<a href="https://arxiv.org/abs/2011.10218" rel="nofollow noreferrer">Optimizing Approximate Leave-one-out Cross-validation</a>.</p>
<p>If we plot out ALOOCV and compare to leave-one-out cross-validation for the example data set,
you can see that it tracks it very closely and the ALOOCV optimum is nearly the same as the
LOOCV optimum.</p>
<p><strong>Compute Leave-one-out Cross-validation</strong></p>
<pre><code>import numpy as np

def compute_loocv(X, y, C):
    model = bbai.glm.LogisticRegression(C=C)
    n = len(y)
    loo_likelihoods = []
    for i in range(n):
        train_indexes = [i_p for i_p in range(n) if i_p != i]
        test_indexes = [i]
        X_train, X_test = X[train_indexes], X[test_indexes]
        y_train, y_test = y[train_indexes], y[test_indexes]
        model.fit(X_train, y_train)
        pred = model.predict_proba(X_test)[0]
        loo_likelihoods.append(pred[y_test[0]])
    return sum(np.log(loo_likelihoods))
</code></pre>
<p><strong>Compute Approximate Leave-one-out Cross-validation</strong></p>
<pre><code>import scipy

def fit_logistic_regression(X, y, C):
    model = bbai.glm.LogisticRegression(C=C)
    model.fit(X, y)
    return np.array(list(model.coef_[0]) + list(model.intercept_))

def compute_hessian(p_vector, X, alpha):
    n, k = X.shape
    a_vector = np.sqrt((1 - p_vector)*p_vector)
    R = scipy.linalg.qr(a_vector.reshape((n, 1))*X, mode='r')[0]
    H = np.dot(R.T, R)
    for i in range(k-1):
        H[i, i] += alpha
    return H

def compute_alo(X, y, C):
    alpha = 1.0 / C
    w = fit_logistic_regression(X, y, C)
    X = np.hstack((X, np.ones((X.shape[0], 1))))
    n = X.shape[0]
    y = 2*y - 1
    u_vector = np.dot(X, w)
    p_vector = scipy.special.expit(u_vector*y)
    H = compute_hessian(p_vector, X, alpha)
    L = np.linalg.cholesky(H)
    T = scipy.linalg.solve_triangular(L, X.T, lower=True)
    h_vector = np.array([np.dot(ti, ti) for pi, ti in zip(p_vector, T.T)])
    loo_u_vector = u_vector - \
        y * (1 - p_vector)*h_vector / (1 - p_vector*(1 - p_vector)*h_vector)
    loo_likelihoods = scipy.special.expit(y*loo_u_vector)
    return sum(np.log(loo_likelihoods))
</code></pre>
<p><strong>Plot out the results (along with the ALOOCV optimum)</strong></p>
<pre><code>import matplotlib.pyplot as plt

Cs = np.arange(0.1, 2.0, 0.1)
loocvs = [compute_loocv(X, y, C) for C in Cs]
alos = [compute_alo(X, y, C) for C in Cs]

fig, ax = plt.subplots()
ax.plot(Cs, loocvs, label='LOOCV', marker='o')
ax.plot(Cs, alos, label='ALO', marker='x')
ax.axvline(model.C_, color='tab:green', label='C_opt')
ax.set_xlabel('C')
ax.set_ylabel('Log-Likelihood')
ax.set_title(&quot;Breast Cancer Dataset&quot;)
ax.legend()
</code></pre>
<p>Displays</p>
<p><a href="https://i.stack.imgur.com/oMQsw.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/oMQsw.png" alt="enter image description here" /></a></p>
sklearn: Hyperparameter tuning by gradient descent?
<p>Is there a way to perform hyperparameter tuning in scikit-learn by gradient descent? While a formula for the gradient of hyperparameters might be difficult to compute, numerical computation of the hyperparameter gradient by evaluating two close points in hyperparameter space should be pretty easy. Is there an existing implementation of this approach? Why is or isn't this approach a good idea?</p>