Stable sigmoid function for binary logistic regression and effect on derivatives?

Question

I'm trying to implement binary logistic regression in python. When x variables are closer to 0, the model performs fantastic! But when x variables are very large positive or negative numbers, that's when troubles arise.

Numpy can only store so many decimal places. When is raised to one of the above numbers, it returns 0 or infinity, respectively. Because of this, my gradients explode or vanish and my model's cost function returns NaN.

I've been searching for a numerically stable form of the sigmoid function and I've found several variations that aren't all that similar. Some suggested using a normalizing constant. Others suggested normalization or standardization to pull x-values closer to 0. But even if I were to implement a change to my sigmoid function (defined as "forward" below) I'm not sure if that would affect my derivatives/gradients?

What's the preferred solution?

Please see my python code below:

class LogReg():
"""
Implementation of binary logistic regression
"""

def __init__(self,dataframe,y_column):
    """
    Initialize LogReg class with a pandas dataframe and the name of the dependent variable column, y.
    All other variable will be assumed to be independent variables, x.
    """
    import pandas as pd
    import numpy as np
    self.y = dataframe[y_column].to_numpy()
    self.x = dataframe.drop(y_column,axis=1).to_numpy()
    self.weights = np.random.uniform(low=0.001, high=0.1, size=(1,self.x.shape[1]))[0]
    self.bias = np.random.uniform(low=0.001, high=0.1, size=(1))[0]

def forward(self,row):
    """
    Produce a predicted y-value given a row (vector of x values).
    The prediction is the dot product of weights and vector x added to the bias.
    """
    x_vars = self.x[row]
    z = np.dot(x_vars,self.weights) + self.bias
    y_pred = 1/(1 + np.exp(-z))
    return y_pred

def loss(self,row):
    """
    Standard logistic regression loss function. 

    When y is 1, the loss function evaluates to the natural log of y_pred.
        When y_pred is near 1, ln(1) will assign a loss near 0.
        When y_pred is near 0, ln(0) will assign a penalty >0.

    When y is 0, the loss function evaluates to natural log of (1-y_pred)
        When y_pred is near 0, ln(1-0) will assign a penalty near 0.
        When y_pred is near 1, ln(1-0) will assign a penalty >0.    

    """
    y_pred = self.forward(row)
    y_true = self.y[row]
    loss = -1*(y_true*np.log(y_pred) + (1-y_true)*np.log(1-y_pred))
    return loss

def cost(self):
    """
    Cost function finds average loss function applied to all rows.
    """
    return sum([self.loss(row) for row in range(self.x.shape[0])])/self.x.shape[0]

def gradient_loss_y_pred(self,row):
    """
    Function that returns the gradient of the loss function with respect to y_pred
    """
    y_true = self.y[row]
    y_pred = self.forward(row)
    return -1*y_true/y_pred + (1-y_true)/(1-y_pred)

def gradient_y_pred_z(self,row):
    """
    Function that returns the gradient of y_pred with respect to z (dot product(weights.X) + bias)
    """
    y_pred = self.forward(row)
    return y_pred*(1-y_pred)

def gradient_z_w(self,row):
    """
    Function that returns the gradient of z with respect to weights. Note, this is not the dot product;
    each weight will be multiplied by its corresponding x-value, but not summed.
    """
    return self.weights * self.x[row]

def gradient_weights(self,row):
    """
    Function that returs the gradient of loss with respect to weights by stringing together the above 3 gradients
    by observing the chain rule.
    """
    return self.gradient_loss_y_pred(row)*self.gradient_y_pred_z(row)*self.gradient_z_w(row)

def gradient_bias(self,row):
    """
    Function that returns the gradient of loss with respect to bias by observing the chain rule, similar to above.
    Note, the gradient of z with respect to bias evaluates to 1; thus there are only 2 functions to determine the
    this gradient
    """
    return self.gradient_loss_y_pred(row)*self.gradient_y_pred_z(row)

def train_model(self, iterations, lr):
    """
    Function that trains the LogReg model via gradient descent. Given a number of iterations and a learning rate,
    this function will begin an epoch (iteration), initialize weight and bias gradients to zero, iterate through
    rows in training data. Once all rows have been operated on, the gradients will be applied to the weights
    and biases according to the learning rate specified.
    """
    for epoch in range(iterations):
        d_w = np.zeros_like(self.weights)
        d_b = 0.0

        for i in range(self.x.shape[0]):
            d_w += self.gradient_weights(i)
            d_b += self.gradient_bias(i)
        d_w = d_w/self.x.shape[0]
        d_b = d_b/self.x.shape[0]

        self.weights -= lr*d_w
        self.bias -= lr*d_b
        if epoch % 50 == 0:
            print(self.cost())

def predict(self,row):
    z = np.dot(row,self.weights) + self.bias
    y_pred = 1/(1 + np.exp(-z))
    if y_pred >= 0.5:
        return 1
    else:
        return 0