Case Study 1: Credit Scoring Optimization at Meridian Financial

Deriving and Implementing Gradient Descent for Regularized Logistic Regression

The Problem

Meridian Financial, a consumer lending institution operating under ECOA and FCRA regulatory constraints, processes approximately 2.3 million credit applications annually. Their risk modeling team maintains a logistic regression model for initial credit scoring — a deliberate choice, since regulators require model interpretability and the ability to provide applicants with specific reasons for adverse decisions.

The current model uses 47 features derived from credit bureau data (payment history, credit utilization, account age, inquiry counts) and demographic-adjacent variables (zip code median income, employment sector). The modeling team is rebuilding the optimizer pipeline to replace their legacy solver with a modern gradient-based approach. They want to understand — precisely, mathematically — what happens during optimization, because regulators may ask.

The task: derive the gradient of the regularized logistic regression loss from first principles, implement multiple optimizers, and compare their convergence behavior on Meridian's credit data.

Deriving the Gradient

The loss function for $N$ applicants with features $\mathbf{x}_i \in \mathbb{R}^d$ and binary labels $y_i \in \{0, 1\}$ (1 = default):

$$ \mathcal{L}(\boldsymbol{\theta}) = -\frac{1}{N}\sum_{i=1}^{N}\left[y_i \log \sigma(\boldsymbol{\theta}^T \mathbf{x}_i) + (1 - y_i)\log(1 - \sigma(\boldsymbol{\theta}^T \mathbf{x}_i))\right] + \frac{\lambda}{2}\|\boldsymbol{\theta}\|^2 $$

We derive the gradient in three steps.

Step 1: Sigmoid derivative. Let $z_i = \boldsymbol{\theta}^T \mathbf{x}_i$ and $\hat{y}_i = \sigma(z_i)$. The sigmoid derivative is:

$$ \frac{d\sigma}{dz} = \sigma(z)(1 - \sigma(z)) = \hat{y}_i(1 - \hat{y}_i) $$

Step 2: Per-sample gradient. By the chain rule:

$$ \frac{\partial \ell_i}{\partial \boldsymbol{\theta}} = \left[-\frac{y_i}{\hat{y}_i} \cdot \hat{y}_i(1 - \hat{y}_i) + \frac{1 - y_i}{1 - \hat{y}_i} \cdot \hat{y}_i(1 - \hat{y}_i)\right]\mathbf{x}_i $$

Simplifying: $-y_i(1 - \hat{y}_i)\mathbf{x}_i + (1 - y_i)\hat{y}_i \mathbf{x}_i = (\hat{y}_i - y_i)\mathbf{x}_i$.

Step 3: Full gradient with regularization.

$$ \nabla_{\boldsymbol{\theta}}\mathcal{L} = \frac{1}{N}\sum_{i=1}^{N}(\hat{y}_i - y_i)\mathbf{x}_i + \lambda \boldsymbol{\theta} $$

This expression is elegant: the gradient is a weighted average of feature vectors, weighted by prediction error. The model concentrates its updates on applicants it misclassifies.

Verifying Convexity

Before optimizing, the team verifies that the problem is convex — guaranteeing that any optimizer will find the global minimum.

The Hessian:

$$ \mathbf{H} = \frac{1}{N}\sum_{i=1}^{N}\hat{y}_i(1 - \hat{y}_i)\mathbf{x}_i \mathbf{x}_i^T + \lambda \mathbf{I} $$

Since $\hat{y}_i(1 - \hat{y}_i) > 0$ for all $\hat{y}_i \in (0, 1)$, each term $\hat{y}_i(1 - \hat{y}_i)\mathbf{x}_i\mathbf{x}_i^T$ is positive semidefinite. The regularization term $\lambda \mathbf{I}$ with $\lambda > 0$ makes the full Hessian positive definite. The loss is therefore strictly convex with a unique global minimum.

import numpy as np
import matplotlib.pyplot as plt

def generate_credit_data(
    n_samples: int = 5000,
    n_features: int = 47,
    default_rate: float = 0.08,
    seed: int = 42
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Generate synthetic credit scoring data mimicking Meridian Financial's
    application characteristics.

    Returns:
        X_train, y_train, X_test, y_test
    """
    rng = np.random.default_rng(seed)

    # True coefficients: most features are weak, a few are strong
    theta_true = rng.normal(0, 0.1, n_features)
    theta_true[:5] = np.array([1.2, -0.8, 0.6, -0.4, 0.9])  # strong predictors

    # Generate features: correlated (as real credit features are)
    cov = np.eye(n_features)
    for i in range(n_features):
        for j in range(i + 1, min(i + 5, n_features)):
            cov[i, j] = cov[j, i] = 0.3 ** (j - i)

    X = rng.multivariate_normal(np.zeros(n_features), cov, n_samples)

    # Add intercept bias to shift default rate
    logits = X @ theta_true - np.log((1 - default_rate) / default_rate)
    probs = 1 / (1 + np.exp(-logits))
    y = rng.binomial(1, probs)

    # Train/test split (80/20)
    split = int(0.8 * n_samples)
    return X[:split], y[:split], X[split:], y[split:]

def sigmoid(z: np.ndarray) -> np.ndarray:
    """Numerically stable sigmoid."""
    return np.where(z >= 0,
                    1 / (1 + np.exp(-z)),
                    np.exp(z) / (1 + np.exp(z)))

def logistic_loss_and_gradient(
    theta: np.ndarray,
    X: np.ndarray,
    y: np.ndarray,
    lam: float = 0.01
) -> tuple[float, np.ndarray]:
    """Regularized logistic regression loss and gradient."""
    N = X.shape[0]
    z = X @ theta
    y_hat = sigmoid(z)

    loss = -np.mean(
        y * np.log(y_hat + 1e-12) + (1 - y) * np.log(1 - y_hat + 1e-12)
    ) + 0.5 * lam * np.dot(theta, theta)

    grad = (1 / N) * X.T @ (y_hat - y) + lam * theta
    return loss, grad

# Generate Meridian Financial's credit data
X_train, y_train, X_test, y_test = generate_credit_data()
print(f"Training: {X_train.shape[0]} applications, {y_train.mean():.1%} default rate")
print(f"Testing:  {X_test.shape[0]} applications, {y_test.mean():.1%} default rate")

Comparing Optimizers

The team implements three optimizers and compares convergence on the credit scoring loss.

def sgd(loss_grad_fn, theta_init, lr=0.1, n_steps=500):
    """Vanilla gradient descent."""
    theta = theta_init.copy()
    history = []
    for _ in range(n_steps):
        loss, grad = loss_grad_fn(theta)
        history.append(loss)
        theta -= lr * grad
    return theta, history

def sgd_momentum(loss_grad_fn, theta_init, lr=0.1, beta=0.9, n_steps=500):
    """SGD with momentum."""
    theta = theta_init.copy()
    v = np.zeros_like(theta)
    history = []
    for _ in range(n_steps):
        loss, grad = loss_grad_fn(theta)
        history.append(loss)
        v = beta * v + grad
        theta -= lr * v
    return theta, history

def adam(loss_grad_fn, theta_init, lr=0.001, beta1=0.9, beta2=0.999,
         eps=1e-8, n_steps=500):
    """Adam optimizer."""
    theta = theta_init.copy()
    m = np.zeros_like(theta)
    v = np.zeros_like(theta)
    history = []
    for t in range(1, n_steps + 1):
        loss, grad = loss_grad_fn(theta)
        history.append(loss)
        m = beta1 * m + (1 - beta1) * grad
        v = beta2 * v + (1 - beta2) * grad**2
        m_hat = m / (1 - beta1**t)
        v_hat = v / (1 - beta2**t)
        theta -= lr * m_hat / (np.sqrt(v_hat) + eps)
    return theta, history

# Optimization setup
lam = 0.01
theta_init = np.zeros(X_train.shape[1])
loss_fn = lambda theta: logistic_loss_and_gradient(theta, X_train, y_train, lam)

# Run optimizers
theta_sgd, hist_sgd = sgd(loss_fn, theta_init, lr=0.5, n_steps=500)
theta_mom, hist_mom = sgd_momentum(loss_fn, theta_init, lr=0.05, beta=0.9, n_steps=500)
theta_adam, hist_adam = adam(loss_fn, theta_init, lr=0.01, n_steps=500)

# Plot convergence
plt.figure(figsize=(10, 6))
plt.semilogy(hist_sgd, label=f"SGD (lr=0.5), final={hist_sgd[-1]:.6f}", linewidth=2)
plt.semilogy(hist_mom, label=f"Momentum (lr=0.05, β=0.9), final={hist_mom[-1]:.6f}", linewidth=2)
plt.semilogy(hist_adam, label=f"Adam (lr=0.01), final={hist_adam[-1]:.6f}", linewidth=2)
plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss (log scale)", fontsize=14)
plt.title("Meridian Financial: Optimizer Convergence on Credit Scoring Loss", fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("credit_scoring_convergence.png", dpi=150)
plt.show()

Results and Analysis

Convergence behavior on convex loss:

On this convex problem, all three optimizers converge to the same minimum (as guaranteed by theory). The differences are in convergence speed:

1. SGD converges fastest in terms of iteration count when the learning rate is well-tuned, but is highly sensitive to the learning rate choice. At $\eta = 1.0$, it diverges. At $\eta = 0.1$, it is slower than momentum.

2. Momentum provides the most consistent convergence across learning rate choices. The velocity accumulation smooths the trajectory and dampens oscillations in the correlated feature space (credit features are notoriously correlated — utilization ratio and balance are mechanically linked).

3. Adam converges more slowly on this well-conditioned convex problem but is the most robust to hyperparameter choice. Its advantage is more pronounced on ill-conditioned or non-convex problems.

Visualizing the Loss Landscape

To build geometric intuition, the team projects the 47-dimensional loss landscape onto two dimensions using the dominant eigenvectors of the Hessian at the optimum.

def compute_hessian_at_optimum(
    theta_star: np.ndarray,
    X: np.ndarray,
    lam: float
) -> np.ndarray:
    """Compute the Hessian of regularized logistic loss at theta_star."""
    y_hat = sigmoid(X @ theta_star)
    W = np.diag(y_hat * (1 - y_hat))  # diagonal weight matrix
    H = (1 / X.shape[0]) * X.T @ W @ X + lam * np.eye(X.shape[1])
    return H

H_star = compute_hessian_at_optimum(theta_adam, X_train, lam)
eigvals, eigvecs = np.linalg.eigh(H_star)
print(f"Condition number: {eigvals[-1] / eigvals[0]:.1f}")
print(f"Smallest eigenvalue: {eigvals[0]:.6f}")
print(f"Largest eigenvalue: {eigvals[-1]:.6f}")

# Project onto top-2 eigenvectors for visualization
d1, d2 = eigvecs[:, -1], eigvecs[:, -2]  # directions of highest curvature

Regulatory Implications

The derivation above answers a specific regulatory question: "Can you demonstrate that your model's parameters are uniquely determined by the training data and the regularization parameter?" Because the loss is strictly convex ($\lambda > 0$), the answer is yes — there is exactly one global minimum, and any convergence optimizer will find it. This determinism is important for FCRA compliance: two independent runs of the training pipeline should produce identical model coefficients (up to numerical precision), yielding identical decisions for identical applicants.

Production Reality: In practice, Meridian's team uses L-BFGS (via scipy.optimize.minimize(method='L-BFGS-B')) for production logistic regression, not Adam. For convex problems of moderate dimension ($d < 10{,}000$), L-BFGS converges in far fewer iterations because it uses approximate second-order information. Adam's advantage emerges only for non-convex problems or when $d$ is very large. The choice of optimizer should be driven by the problem structure, not by familiarity.

Key Lessons

1. Convexity provides guarantees. For convex problems like logistic regression, all optimizers reach the same solution. The question is speed, not quality.
2. The gradient has a natural interpretation. Each data point contributes proportionally to its prediction error — the model focuses on its mistakes.
3. Condition number determines difficulty. Feature correlation in credit data increases the condition number, slowing convergence. L2 regularization improves conditioning (one of its less appreciated benefits).
4. Choose the optimizer for the problem. L-BFGS for convex + moderate dimension. Adam for non-convex or very high dimension. SGD with momentum when generalization matters more than convergence speed.