Logistic Regression from scratch

Algorithm and implementation
#

Logistic Regression formula is the same as the Linear Regression’s one passed to a sigmoid function to obtain a probability score from logits and can be replicated on numpy with:
sigmoid(X @ w + b) (or sigmoid(np.dot(X,w)).
The implementation of gradient descent as an optimizer is the same of Linear Regression from scratch, what changes is the loss and the error since we deal with probabilities:

errors = probabilities - y
The loss is binary cross entropy:
$$BCE(y,p)=-[y\log(p)+(1-y)\log(1-p)], \qquad L=mean(BCE).$$ And it can be read as two cases:
- when y=1 we have only -log(p)
- when y=0 we have only -log(1-p) The loss is small when the correct outcome receives high probability.
  The derivative is then: $$dw=X^T errors/N, \qquad db=mean(errors).$$
weights w and biases are both initialized with zero values. After choosing a learning rate, at further steps w and b will be updated as:
- w-=dw * learning_rate
- b-=db * learning_rate

Class structure
#

main attributes:
- learning_rate
- n_iterations
- weights w_ and bias b_ initialized as None
fit method
predict_proba method
predict method
_sigmoid static method (it uses the @staticmethod decorator, that simply allows us to not include the self , since this function does not depend on any self. parameter)
_binary_cross_entropy method

import numpy as np

class LogisticRegression:
	def __init__(self, learning_rate=0.05, n_iterations=1000):
		self.learning_rate = learning_rate
		self.n_iterations = n_iterations
		self.w_ = None
		self.b_ = None
		
	def fit(self, X, y):
		pass
		
	def predict_proba(self, X):
		pass
	
	@staticmethod
	def _sigmoid(score):
		pass
		
	@staticmethod
	def _binary_cross_entropy(probabilities, y, epsilon):
		pass
		
	def predict(self, X):
		pass

Sigmoid
#

To pass from logits to probabilities, we need the sigmoid function to transform the linear score.
Its base implementation is just as follows:

@staticmethod
def _sigmoid(scores):
	return 1/(1+np.exp(-scores))

This formula is correct for every real number, but a computer has a largest representable floating-point value that can lead to overflow.
If z = -1000, the direct formula asks NumPy to compute exp(1000), which is too large.
The final sigmoid should simply be very close to zero, so an intermediate overflow is unnecessary.

We choose the formula according to the sign of z so the value passed to exp is always zero or negative:

for non-negative score (z >= 0) we keep the standard formula. Here -z <= 0, so exp(-z) is safe. $$sigmoid(z)=\frac{1}{1+e^{-z}}.$$
for negative score (z < 0): rearrange the same formula by multiplying its numerator and denominator by exp(z): $$\frac{1}{1+e^{-z}}\times\frac{e^z}{e^z}=\frac{e^z}{1+e^z}.$$Now, if z = -1000, we compute exp(-1000), which safely becomes approximately zero. The probability also becomes approximately zero. We never compute the dangerous exp(1000).

The two branches are not two different sigmoid functions. They are two algebraically equivalent ways to calculate the same function, selected only to keep the intermediate exponential safe.

@staticmethod
def _sigmoid(scores):
	#create an array with the desired shape filled with zeros
	scores = np.asarray(scores, dtype=float)
	probabilities = np.zeros_like(scores)
	
	#create the branch filters
	nonnegative = scores >= 0
    negative = scores < 0
	
	#fill the probabilities array from the positive branch
	probabilities[nonnegative] = (1.0/ (1.0 + np.exp(-scores[nonnegative])))
	
	#fill from the negative branch
	exp_scores = np.exp(scores[negative])
	probabilities[negative] = exp_scores / (1.0 + exp_scores)
	
	return probabilities

Binary cross-entropy
#

Binary cross-entropy clips probabilities, calculates the positive-label and negative-label terms, adds them and averages.

@staticmethod
def _binary_cross_entropy(probabilities, y, epsilon=1e-12):
	clipped_probabilities = np.clip(probabilities, epsilon, 1.0 - epsilon)
	
	positive_label_loss = -y * np.log(clipped_probabilities)
	negative_label_loss = -(1.0 - y) * np.log(1.0 - clipped_probabilities)
	
	per_example_loss = positive_label_loss + negative_label_loss
	loss = np.mean(per_example_loss)
	return loss

Training
#

Similar to Linear Regression from scratch, the fit method can be implemented following these steps:

input shapes setup with X of shape (N,D) and y of shape (N,)
weights and bias initialization
“SPELGU loop”:
- Score: score = X @ self.w_ + self.b_ (shape (N,))
- Probabilities: self._sigmoid(score) (shape (N,))
- Errors: errors = probabilities - y (shape (N,))
- Loss: loss = binary_cross_entropy(probabilities, y) (scalar).
- Gradient:
  - dw = X.T @ errors / n_samples (shape (D,))
  - db = errors.mean() (scalar)
- Update:
  - self.w_ = self.learning_rate * dw (shape (D,))
  - self.b_ = self.learning_rate * db (scalar)

def fit(self, X, y):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = 0.0
        self.loss_history_ = []

        for _ in range(self.n_iterations):
            score = X @ self.w_ + self.b_       # S: (N,)
            probabilities = self._sigmoid(scores)       # P: (N,)
            errors = probabilities - y                  # E: (N,)
            loss = self._binary_crossentropy(probabilities, y) # L: scalar

            dw = X.T @ errors / n_samples              # G: (D,)
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
           
        return self

Inference (probability)
#

The prediction of the probability can be implemented as:

def predict_proba(self, X):
	if self.w_ is None:
		raise RuntimeError("fit must be called before predict")
	X = np.asarray(X, dtype=float)
	scores = X @ self.w_ + self.b_
	probabilities = self._sigmoid(scores)
	return probabilities

Inference (label)
#

In case we need a label prediction and not a probability prediction, we use the predict method setting up a threshold:

def predict(self, X, threshold=0.5):
	probabilities = self.predict_proba(X)
	labels = (probabilities >= threshold).astype(int)
	return labels

Features scaling
#

To make learning more robust, Standardization can be used.

Adding regularization
#

Similar to Linear Regression, adding regularization is a common step to tackle situations in which a model can fit training noise with unnecessarily large weights (overfitting), especially when features are correlated or data is limited.

We can penalize the model for large weights adding a sort of “penalty” to the loss we saw in the base implementation:
regularization_loss = 0.5 * l2 * sum(w ** 2) Its gradient is regularization_dw = l2 * w, which pulls every weight toward zero, with l2 controlling the “strength of the pull”.

This affects the total loss and the weights gradient dw. Bias b does not need to be regularized because it does not measure the dependence on a feature.
This step usually comes after scaling.

The fit method of the LogisticRegression class would look like:

def fit(self, X, y, l2):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = 0.0
        self.loss_history_ = []

        for _ in range(self.n_iterations):
            score = X @ self.w_ + self.b_       # S: (N,)
            probabilities = self._sigmoid(scores)       # P: (N,)
            errors = probabilities - y                  # E: (N,)
            
            data_loss = self._binary_crossentropy(probabilities, y) # L: scalar
            regularization_loss = 0.5 * l2 * np.sum(self.w_ ** 2)
            loss = data_loss + regularization_loss
            
            data_dw = X.T @ errors / n_samples     # G: (D,)
            regularization_dw = l2 * self.w_
			dw = data_dw + regularization_dw
                         
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
           
        return self

Adding early-stopping
#

When an algorithm converges, further iterations barely change the parameters. A fixed iterations count is safe but it may waste work. The idea behind early stopping is to stop when the largest absolute update of w and b is below tolerance.
If, on the contrary, learning rate diverges we must stop the training as we are likely finding nulls in the loss.

The implementation is the same as the one used in Linear Regression from scratch.

Algorithm and implementation#

Class structure#

Sigmoid#

Binary cross-entropy#

Training#

Inference (probability)#

Inference (label)#

Features scaling#

Adding regularization#

Adding early-stopping#