Sigmoid, Log-Loss, and Gradients in C++
This first implementation turns the logistic-regression contract into a complete batch gradient-descent classifier. The public surface distinguishes logits, probabilities, and decisions; the objective is evaluated from logits; every accepted update is checked for finite, non-increasing loss; and a numerical gradient check tests the derivation before fitting begins.
From Bernoulli Likelihood to Log-Loss
For a binary label $y_i\in\{0,1\}$ and probability $p_i=P(y_i=1\mid x_i)$, one Bernoulli probability mass is
$$P(y_i\mid x_i)=p_i^{y_i}(1-p_i)^{1-y_i}.$$Conditional independence turns the dataset likelihood into a product. Products of many probabilities underflow, so optimization uses the negative log and divides by $n$:
$$L(\beta)=-\frac{1}{n}\sum_{i=1}^{n}\left[y_i\log p_i+(1-y_i)\log(1-p_i)\right].$$The division changes the scale, not the minimizer. It also keeps a learning rate comparable as the number of rows changes.
Stable Sigmoid and Extreme Logits
The linear score $z_i=\beta_0+x_i^\top\beta$ is unrestricted. The sigmoid maps it to a probability. A branch on the sign avoids evaluating an exponential with a large positive argument:
$$\sigma(z)= \begin{cases} 1/(1+e^{-z}), & z\ge 0,\\ e^z/(1+e^z), & z<0. \end{cases}$$For $z=1000$, the positive branch evaluates $e^{-1000}$ and returns a value rounded to one. For $z=-1000$, the negative branch evaluates $e^{-1000}$ and returns a value rounded to zero. Both results are finite. Rounded endpoints are acceptable probability outputs, but taking their logarithms is not a stable loss computation.
Softplus Loss Directly from Logits
Substituting $p=\sigma(z)$ into Bernoulli log-loss gives an equivalent per-row form:
$$\ell(z,y)=\log(1+e^z)-yz=\operatorname{softplus}(z)-yz.$$A stable softplus identity keeps the exponential bounded:
$$\operatorname{softplus}(z)=\max(z,0)+\log(1+e^{-|z|}).$$The implementation uses log1p for accuracy near zero. The equivalent branch form $\ell(z,1)=\operatorname{softplus}(-z)$ and $\ell(z,0)=\operatorname{softplus}(z)$ avoids a potentially cancelling subtraction for a large positive $z$ with $y=1$. Consequently, all four checks $\ell(\pm1000,0)$ and $\ell(\pm1000,1)$ remain finite; correct extreme predictions have loss near zero and confidently wrong predictions have loss near $1000$.
Batch Mean Gradient
The derivative identities $\sigma'(z)=p(1-p)$ and $\partial\ell/\partial z=p-y$ collapse the chain rule. With $X_b$ denoting $X$ with an implicit leading column of ones,
$$\nabla_\beta L=\frac{1}{n}X_b^\top(p-y).$$The intercept component is the mean residual probability, $\frac{1}{n}\sum_i(p_i-y_i)$. Feature component $j$ is $\frac{1}{n}\sum_i(p_i-y_i)x_{ij}$. A batch update subtracts $\alpha\nabla L$.
Differentiating once more gives the curvature used in Part 2:
$$H=\frac{1}{n}X_b^\top W X_b,\qquad W=\operatorname{diag}(p_i(1-p_i)).$$Since every diagonal weight is nonnegative, the Hessian is positive semidefinite and log-loss is convex. A full-rank, non-separated design gives a unique finite optimum.
Central Finite-Difference Check
For each parameter $k$, a central difference independently approximates the derivative:
$$g_k^{FD}=\frac{L(\beta+\varepsilon e_k)-L(\beta-\varepsilon e_k)}{2\varepsilon}.$$The program uses $\varepsilon=10^{-6}$ at a nonzero test parameter and requires the largest absolute discrepancy to stay below $10^{-6}$. Central differences have $O(\varepsilon^2)$ truncation error, while extremely small $\varepsilon$ amplifies floating-point cancellation. The chosen value balances those effects for this dataset.
Training and Convergence Safeguards
- The full batch gradient and mean loss are finite before an update is considered.
- Backtracking halves the trial learning rate until candidate loss is no larger than current loss, up to 30 trials.
- Convergence occurs when gradient norm is below tolerance or relative loss change is below tolerance.
- Parameter magnitude is bounded as a separation and divergence guard.
- History records the initial state, every 100th accepted iteration, and the final state.
- Failure to find a finite descent step or converge within the configured limit raises
std::runtime_error.
Gradient Descent from the Underlying Rows
The plot below uses the same eight rows as the C++ program. JavaScript recomputes stable sigmoid probabilities, mean softplus loss, and full-batch gradients, then stores selected optimization states. The curve and loss label therefore arise from the displayed data and objective.
Batch Gradient Descent on Bernoulli Log-Loss
Validation and Error Behavior
| Condition | Result |
|---|---|
| Empty or ragged $X$; length mismatch | std::invalid_argument |
| Non-finite feature or label outside $\{0,1\}$ | std::invalid_argument naming the index |
| Nonpositive learning rate, tolerance, or iteration limit | std::invalid_argument |
| Prediction before fitting | std::logic_error |
| Wrong prediction width, non-finite feature, or threshold outside $(0,1)$ | std::invalid_argument |
| Non-finite score, stalled line search, divergence, or iteration exhaustion | std::runtime_error |
Complexity
Each loss and gradient pass costs $O(nd)$ time. Backtracking can add more loss passes, bounded here by 30, while accepted iterations remain $O(nd)$ for a fixed limit. Parameters and gradient use $O(d)$ auxiliary memory; history uses $O(h)$ records. The caller-owned $X$ remains $O(nd)$. Gradient descent avoids the $O(d^2)$ Hessian storage and $O(d^3)$ solve used by Newton’s method, trading cheaper iterations for more of them.
Complete C++17 Program
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>
using Matrix = std::vector<std::vector<double>>;
double stableSigmoid(double z) {
if (!std::isfinite(z))
throw std::invalid_argument("stableSigmoid: z must be finite");
if (z >= 0.0) {
const double e = std::exp(-z);
return 1.0 / (1.0 + e);
}
const double e = std::exp(z);
return e / (1.0 + e);
}
double softplus(double z) {
if (!std::isfinite(z))
throw std::invalid_argument("softplus: z must be finite");
return std::max(z, 0.0) + std::log1p(std::exp(-std::fabs(z)));
}
double logLossFromLogit(double z, int y) {
if (y != 0 && y != 1)
throw std::invalid_argument("logLossFromLogit: y must be 0 or 1");
return y == 1 ? softplus(-z) : softplus(z);
}
struct IterationRecord {
int iteration;
double loss;
double gradientNorm;
};
struct GDConfig {
double learningRate = 0.5;
int maxIterations = 20000;
double tolerance = 1e-10;
double maxAbsParameter = 1e6;
};
struct LogisticGD {
bool fitted = false;
int nFeatures = 0;
double intercept = 0.0;
std::vector<double> coef;
std::vector<IterationRecord> history;
static void validateData(const Matrix& X, const std::vector<int>& y) {
if (X.empty())
throw std::invalid_argument("fit: X must contain at least one row");
if (X.size() != y.size())
throw std::invalid_argument("fit: X row count must equal y size");
if (X.front().empty())
throw std::invalid_argument("fit: X must contain at least one feature");
const std::size_t d = X.front().size();
for (std::size_t i = 0; i < X.size(); ++i) {
if (X[i].size() != d)
throw std::invalid_argument("fit: ragged X at row " +
std::to_string(i));
for (std::size_t j = 0; j < d; ++j) {
if (!std::isfinite(X[i][j]))
throw std::invalid_argument("fit: non-finite feature at row " +
std::to_string(i) + ", column " +
std::to_string(j));
}
if (y[i] != 0 && y[i] != 1)
throw std::invalid_argument("fit: label at index " +
std::to_string(i) + " must be 0 or 1");
}
}
static double score(const std::vector<double>& beta,
const std::vector<double>& x) {
double z = beta[0];
for (std::size_t j = 0; j < x.size(); ++j)
z += beta[j + 1] * x[j];
if (!std::isfinite(z))
throw std::runtime_error("score: non-finite linear predictor");
return z;
}
static double meanLoss(const Matrix& X, const std::vector<int>& y,
const std::vector<double>& beta) {
double total = 0.0;
for (std::size_t i = 0; i < X.size(); ++i)
total += logLossFromLogit(score(beta, X[i]), y[i]);
const double value = total / static_cast<double>(X.size());
if (!std::isfinite(value))
throw std::runtime_error("meanLoss: non-finite objective");
return value;
}
static std::vector<double> gradient(
const Matrix& X, const std::vector<int>& y,
const std::vector<double>& beta) {
std::vector<double> g(beta.size(), 0.0);
for (std::size_t i = 0; i < X.size(); ++i) {
const double residual = stableSigmoid(score(beta, X[i])) - y[i];
g[0] += residual;
for (std::size_t j = 0; j < X[i].size(); ++j)
g[j + 1] += residual * X[i][j];
}
for (double& value : g)
value /= static_cast<double>(X.size());
return g;
}
void fit(const Matrix& X, const std::vector<int>& y,
const GDConfig& config = {}) {
validateData(X, y);
if (!std::isfinite(config.learningRate) || config.learningRate <= 0.0)
throw std::invalid_argument("fit: learningRate must be finite and positive");
if (config.maxIterations <= 0)
throw std::invalid_argument("fit: maxIterations must be positive");
if (!std::isfinite(config.tolerance) || config.tolerance <= 0.0)
throw std::invalid_argument("fit: tolerance must be finite and positive");
if (!std::isfinite(config.maxAbsParameter) ||
config.maxAbsParameter <= 0.0)
throw std::invalid_argument("fit: maxAbsParameter must be positive");
std::vector<double> beta(X.front().size() + 1, 0.0);
std::vector<IterationRecord> nextHistory;
double currentLoss = meanLoss(X, y, beta);
bool converged = false;
for (int iteration = 0; iteration < config.maxIterations; ++iteration) {
const std::vector<double> g = gradient(X, y, beta);
double squaredNorm = 0.0;
for (double value : g) squaredNorm += value * value;
const double gradientNorm = std::sqrt(squaredNorm);
if (!std::isfinite(gradientNorm))
throw std::runtime_error("fit: non-finite gradient");
if (iteration == 0 || iteration % 100 == 0)
nextHistory.push_back({iteration, currentLoss, gradientNorm});
if (gradientNorm <= config.tolerance) {
converged = true;
break;
}
double step = config.learningRate;
std::vector<double> candidate(beta.size());
double candidateLoss = std::numeric_limits<double>::infinity();
bool accepted = false;
for (int trial = 0; trial < 30; ++trial) {
bool bounded = true;
for (std::size_t k = 0; k < beta.size(); ++k) {
candidate[k] = beta[k] - step * g[k];
bounded = bounded && std::isfinite(candidate[k]) &&
std::fabs(candidate[k]) <= config.maxAbsParameter;
}
if (bounded) {
candidateLoss = meanLoss(X, y, candidate);
if (candidateLoss <= currentLoss) {
accepted = true;
break;
}
}
step *= 0.5;
}
if (!accepted)
throw std::runtime_error("fit: line search found no finite descent step");
const double relativeChange =
std::fabs(currentLoss - candidateLoss) / std::max(1.0, currentLoss);
beta = candidate;
currentLoss = candidateLoss;
if (relativeChange <= config.tolerance) {
const std::vector<double> finalGradient = gradient(X, y, beta);
double finalSquaredNorm = 0.0;
for (double value : finalGradient)
finalSquaredNorm += value * value;
nextHistory.push_back(
{iteration + 1, currentLoss, std::sqrt(finalSquaredNorm)});
converged = true;
break;
}
}
if (!converged)
throw std::runtime_error("fit: maximum iterations reached");
nFeatures = static_cast<int>(X.front().size());
intercept = beta[0];
coef.assign(beta.begin() + 1, beta.end());
history = std::move(nextHistory);
fitted = true;
}
double predictLogit(const std::vector<double>& x) const {
if (!fitted)
throw std::logic_error("predictLogit: model is not fitted");
if (x.size() != static_cast<std::size_t>(nFeatures))
throw std::invalid_argument("predictLogit: feature count mismatch");
for (std::size_t j = 0; j < x.size(); ++j)
if (!std::isfinite(x[j]))
throw std::invalid_argument("predictLogit: non-finite feature at column " +
std::to_string(j));
std::vector<double> beta(1, intercept);
beta.insert(beta.end(), coef.begin(), coef.end());
return score(beta, x);
}
double predictProba(const std::vector<double>& x) const {
return stableSigmoid(predictLogit(x));
}
int predict(const std::vector<double>& x,
double threshold = 0.5) const {
if (!std::isfinite(threshold) || threshold <= 0.0 || threshold >= 1.0)
throw std::invalid_argument("predict: threshold must be inside (0, 1)");
return predictProba(x) >= threshold ? 1 : 0;
}
};
double gradientCheck(const Matrix& X, const std::vector<int>& y) {
const std::vector<double> beta = {0.2, -0.35};
const std::vector<double> analytic = LogisticGD::gradient(X, y, beta);
const double epsilon = 1e-6;
double largest = 0.0;
for (std::size_t k = 0; k < beta.size(); ++k) {
std::vector<double> plus = beta;
std::vector<double> minus = beta;
plus[k] += epsilon;
minus[k] -= epsilon;
const double finiteDifference =
(LogisticGD::meanLoss(X, y, plus) -
LogisticGD::meanLoss(X, y, minus)) / (2.0 * epsilon);
largest = std::max(largest, std::fabs(analytic[k] - finiteDifference));
}
return largest;
}
int main() {
const Matrix X = {{-3.0}, {-2.0}, {-1.0}, {0.0},
{1.0}, {2.0}, {3.0}, {4.0}};
const std::vector<int> y = {0, 0, 0, 1, 0, 1, 1, 1};
const double gradientError = gradientCheck(X, y);
if (gradientError > 1e-6)
throw std::runtime_error("main: finite-difference gradient check failed");
LogisticGD model;
model.fit(X, y);
int correct = 0;
double fittedLoss = 0.0;
for (std::size_t i = 0; i < X.size(); ++i) {
correct += model.predict(X[i]) == y[i] ? 1 : 0;
fittedLoss += logLossFromLogit(model.predictLogit(X[i]), y[i]);
}
fittedLoss /= static_cast<double>(X.size());
for (std::size_t i = 1; i < model.history.size(); ++i)
if (model.history[i].loss > model.history[i - 1].loss + 1e-12)
throw std::runtime_error("main: accepted-loss monotonicity failed");
const double low = stableSigmoid(-1000.0);
const double high = stableSigmoid(1000.0);
if (!(low >= 0.0 && low <= 1.0 &&
high >= 0.0 && high <= 1.0))
throw std::runtime_error("main: sigmoid range property failed");
for (int label : {0, 1}) {
if (!std::isfinite(logLossFromLogit(-1000.0, label)) ||
!std::isfinite(logLossFromLogit(1000.0, label)))
throw std::runtime_error("main: extreme-logit loss property failed");
}
std::cout << std::fixed << std::setprecision(6);
std::cout << "gradient max error: " << gradientError << '\n';
std::cout << "intercept: " << model.intercept << '\n';
std::cout << "coefficient: " << model.coef.at(0) << '\n';
std::cout << "mean log-loss: " << fittedLoss << '\n';
std::cout << "training accuracy: "
<< static_cast<double>(correct) / X.size() << '\n';
std::cout << "P(y=1 | x=0.5): "
<< model.predictProba({0.5}) << '\n';
std::cout << "history records: " << model.history.size() << '\n';
std::cout << "property checks: passed\n";
}Expected Output
The exact values below correspond to the stated compiler-independent deterministic rows. The final digits can vary slightly with the standard library’s elementary functions.
gradient max error: 0.000000 intercept: -0.640986 coefficient: 1.282227 mean log-loss: 0.312881 training accuracy: 0.750000 P(y=1 | x=0.5): 0.500032 history records: 4 property checks: passed
From First-Order to Curvature
Gradient descent sees only slope. Newton’s method and IRLS add the Hessian, solve a weighted least-squares system, and use damping to preserve the same monotonic-loss guarantee. The public classifier contract and stable numerical primitives remain unchanged.