Logistic Regression from Scratch in C++

July 2026 · 36 min read

The linear-regression sub-series ended with a diagnostic layer for a continuous target. This sub-series keeps the same linear score $z=\beta_0+\sum_{j=1}^d\beta_j x_j$ underneath, but changes what is being predicted: a binary target modeled as a single Bernoulli draw, linked to that score through the logit and squashed into a probability by the sigmoid, then fit by minimizing log-loss instead of squared error. Standard-library C++ only, no external linear-algebra dependency, the same discipline as every earlier post in the from-scratch series. The companion Machine Learning · Logistic Regression series covers the conceptual treatment in full depth, the Bernoulli likelihood derivation, log-odds interpretation, separation theory, and calibration theory; this sub-series turns that theory into working, validated code, cross-linking to the specific conceptual post relevant to each section below.

Series Scope

Every post here compiles with a standard C++17 compiler and nothing else, continuing the linear-regression sub-series' vector<vector<double>>-only convention. The roadmap covers five posts in increasing depth: the sigmoid link, log-loss objective, and its gradient and Hessian; a second-order Newton/IRLS solver; $L_1$/$L_2$-penalized variants for collinear features and separation; a generalization to more than two classes through softmax; and the classification-specific evaluation layer, thresholds, ranking metrics, and probability calibration, that $R^2$ and Cook's distance have no direct analogue for. The remainder of this overview establishes the contract every one of those five posts implements against, so later posts change only the estimation method inside fit, never the shape of the problem.

The Bernoulli Model and Binary Targets

The shape of the problem changes in exactly one place relative to the linear-regression contract: the target. $X$ is still an $n\times d$ raw feature matrix with no bias column, but each label is $y_i\in\{0,1\}$ instead of a continuous real number, treated as one draw from a Bernoulli distribution whose success probability $p_i=P(y_i=1\mid x_i)$ is itself a function of $x_i$:

That last compact form is what makes maximum-likelihood estimation tractable: it holds for both values of $y_i$ at once, so the likelihood of the whole training set is a single product over $i=1,\dots,n$ of that expression, and its negative log is the log-loss objective every solver in this sub-series minimizes. A model that only ever produced $0$ or $1$ directly could not be fit this way; the Bernoulli parameter $p_i$ has to be a genuine probability in the open interval $(0,1)$ for the likelihood to be well defined, which is exactly the constraint the sigmoid link exists to enforce.

Log-Odds, the Logit Link, and the Sigmoid

A linear score $z_i=\beta_0+\sum_{j=1}^d\beta_j x_{ij}$ ranges over all of $\mathbb{R}$, but a probability is confined to $(0,1)$; connecting the two requires a link function. The odds of the event $y_i=1$ rescale $p_i$ onto $(0,\infty)$, and the logit, the natural log of the odds, rescales that onto all of $\mathbb{R}$, exactly matching the range of $z_i$:

Solving $\operatorname{logit}(p_i)=z_i$ for $p_i$ gives the sigmoid function, the inverse of the logit link:

Three consequences follow directly from that derivation and recur throughout the rest of this sub-series. First, $z_i=0$ corresponds exactly to $p_i=0.5$, so it is the probability-$0.5$ contour. A hard prediction threshold $t$ places its own decision boundary at $z_i=\operatorname{logit}(t)$, which equals zero only when $t=0.5$; the threshold animation below makes that distinction explicit. Second, in exact arithmetic every finite $z_i$ maps to a $p_i$ strictly inside $(0,1)$, never exactly $0$ or $1$; floating-point evaluation can round extreme finite logits to an endpoint, which is why the stable loss later works directly from $z$ rather than from the rounded probability. Perfectly separable data drives unpenalized coefficients toward infinity rather than to a finite optimum, covered in the regularization and separation section below. Third, the logit is the natural scale on which a unit change in a feature has a constant, additive effect on log-odds; the companion Probability, Log-Odds, and Interpretation post works through that coefficient-as-odds-ratio interpretation in detail.

The walkthrough below sweeps seven actual $z$ values from $-10$ to $10$, computing the odds, the sigmoid probability through the numerically stable branch-on-sign formula, the logit recovered from that probability as a round-trip check, and the log-loss the same $z$ would incur for both a true label $y=1$ and a true label $y=0$:

Probability Estimation versus Decisions

Fitting the Bernoulli model produces a probability, $\hat p=\sigma(\hat\beta_0+\sum_j\hat\beta_j x_j)$, not a label. Turning that probability into a $0/1$ decision is a second, separate step: comparing $\hat p$ against a threshold $t$, predicting $1$ when $\hat p\ge t$ and $0$ otherwise. Nothing about estimating $\hat p$ requires choosing $t$, and nothing about choosing $t$ changes $\hat p$; the two are independent by design. The default $t=0.5$ is the Bayes-optimal cutoff under $0/1$ loss with balanced misclassification costs, but it is a default, not a law, and the imbalance and cost section below works through exactly when a different $t$ is the honest choice. Keeping this separation explicit in the contract, a dedicated predictProba alongside a predict(x, threshold) that only thresholds the result, is what lets the same fitted model serve a high-recall fraud screen and a high-precision one from a single training run.

Handling the Intercept

The intercept convention carries over unchanged from the linear-regression sub-series: callers pass a raw $n\times d$ matrix with no column of ones, and the model owns $\beta_0$ internally. Whether a post's solver folds an implicit ones-column into IRLS's weighted normal equations directly, or centers $X$ before optimizing and recovers $\beta_0$ afterward, the public contract stays fit(X, y) in, a fitted intercept and coef vector of length $d$ out, exactly as for every regressor in the linear series.

The Full C++ Contract

Every classifier in this sub-series follows the same struct shape, with the raw logit, the probability, and the thresholded label kept as three explicitly distinct outputs rather than one blended call:

using Matrix = vector<vector<double>>;

// X: n rows (samples) x d columns (features), raw, no bias column.
// y: n labels, each exactly 0 or 1.
struct LogisticRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;              // length nFeatures once fitted

    void fit(const Matrix& X, const vector<int>& y);

    double predictLogit(const vector<double>& x) const;   // raw z, unbounded
    double predictProba(const vector<double>& x) const;   // sigma(z) in (0,1)
    int predict(const vector<double>& x, double threshold = 0.5) const;

    vector<double> predictProba(const Matrix& X) const;    // batch overload
};

Five properties of this contract hold across all five posts and are validated explicitly, with a specific exception type and a message naming the offending value, rather than assumed or caught broadly and papered over. fit requires every row of $X$ to have exactly $d$ entries and $n$ to equal y.size(), exactly as in the linear series, plus two classification-specific checks: every entry of $X$ must be finite, and every entry of $y$ must be exactly $0$ or $1$. Failures name the offending feature coordinate or label index and value, since neither a linear score built from a non-finite feature nor a Bernoulli label outside $\{0,1\}$ has a defined likelihood. The parameters a solver hands to the model are gated the same way: setParams rejects a non-finite intercept or coefficient before it is stored, so a diverging or corrupted solver step cannot poison every later prediction. predictLogit, and everything built on it, throws std::logic_error if called before fit, rejects a feature-count mismatch or non-finite feature with std::invalid_argument, and reports arithmetic overflow in the final score with std::runtime_error. predict accepts only a finite threshold strictly between $0$ and $1$; endpoint thresholds are degenerate policies whose behavior can depend on floating-point saturation, while a non-finite threshold makes ordered comparison meaningless. Even the loss function is defensive on its own terms: stableLogLoss rejects a label outside $\{0,1\}$ itself because later solvers may call it directly rather than through validateFit. Finally, predictLogit remains a first-class member because stable loss evaluation and several ranking metrics need the raw score before it is squashed into a probability.

Train-Only Scaling

Feature scaling is not part of LogisticRegressor itself; it is a preprocessing step that happens before fit is called, using the same scaler contract as the data-preprocessing sub-series: a scaler's mean and standard deviation are estimated by calling fit on the training rows only, frozen as members, and reused unchanged by transform on both the training and test rows. None of the five posts in this sub-series recompute scaling statistics from a validation or test split, and none of them fold scaling into the model struct itself; LogisticRegressor always receives whatever $X$ it is given, already scaled or not, and treats it as the raw feature matrix for the purposes of the contract above. This separation matters more here than in ordinary least squares because the gradient-based and second-order solvers in the second post converge in a very different number of iterations depending on feature scale, while a scaler fit on test rows would leak test-distribution information into every subsequent probability the model reports, exactly the leakage the linear-regression gradient-descent post already guards against for a continuous target.

The fit / predictLogit / predictProba / predict Lifecycle

The three-stage lifecycle from the linear series carries over with the logit made an explicit intermediate stage:

• fit(X, y) validates shapes and labels, and estimates $\beta_0,\beta_{1..d}$ by minimizing log-loss, through Newton/IRLS in the second post or a penalized variant in the third.
• predictLogit(x) returns the raw unbounded score $z=\beta_0+\sum_j\beta_j x_j$, useful directly for numerically stable loss computation and for locating the decision boundary without ever forming a probability.
• predictProba(x) applies the stable sigmoid to that logit, returning a probability estimate optimized under log-loss, not yet a label. Out-of-sample calibration remains a separate property to measure.
• predict(x, threshold) thresholds that probability, defaulting to $0.5$ but accepting any caller-chosen cutoff, as covered above.
• evaluate is again free functions rather than members, and the set is larger than the regression series': accuracy, precision, recall, and ROC/PR-AUC depend on the threshold or on a ranking rather than a single residual, covered in full in the final post.

Stable Sigmoid, Softplus, and Extreme Logits

Two numerical failure modes are specific to this sub-series and have no analogue in the linear-regression series. First, the naive sigmoid formula $1/(1+e^{-z})$ overflows for very negative $z$: $-z$ becomes a large positive exponent, $e^{-z}$ overflows toward infinity, and while $1/(1+\infty)$ happens to still round to the correct $0.0$ in IEEE 754, the intermediate overflow is not something to rely on across compilers or optimization levels. Every post here instead uses the branch-on-sign form, algebraically identical on both branches but never passing a positive argument to exp:

double stableSigmoid(double z) {
    if (z >= 0.0) {
        double e = exp(-z);       // -z <= 0, exp never overflows
        return 1.0 / (1.0 + e);
    }
    double e = exp(z);            // z < 0, exp never overflows
    return e / (1.0 + e);
}

Second, and more consequential, the textbook combined log-loss formula $-[y\ln p+(1-y)\ln(1-p)]$ evaluated from an already-rounded probability $p$ fails outright once $p$ saturates to exactly $0.0$ or $1.0$ in floating point, which the stable sigmoid above will do for $|z|$ around a few hundred. At $p=1.0$ and $y=0$, the zero-weighted term $(1-y)\ln(1-p)=1\cdot\ln(0)=-\infty$; the other term is $y\ln p=0\cdot\ln(1)=0\cdot0=0$, but IEEE 754 arithmetic does not simplify $0\cdot(-\infty)$ to $0$, it produces NaN, and that single NaN silently poisons a mean loss or a gradient sum computed over the whole batch. The fix is to never materialize $p$ before taking the logarithm at all, instead expressing the per-sample loss directly from the logit $z$ through the numerically stable softplus function, $\operatorname{softplus}(z)=\ln(1+e^{z})$, itself computed with the same branch-on-sign discipline together with log1p so that the argument to log1p never exceeds $1$:

The per-sample log-loss reduces to a single logaddexp-style expression in terms of softplus and the raw logit, valid for both possible labels at once:

This form stays finite for every finite $z$, including $z=\pm1000$: $\operatorname{softplus}(1000)$ rounds to exactly $1000.0$ rather than diverging, correctly reporting that a model confidently and completely wrong at that logit pays a loss proportional to how wrong it is, not an unusable infinity. The smoke test below exercises exactly this contrast, printing inf and nan from the naive formula side by side with the finite, correct values the softplus form produces at the same logits.

Deriving the Gradient and Hessian

Writing the mean log-loss over $n$ training rows using the logit form above, with $z_i$ the linear score for row $i$:

Differentiating softplus gives back the sigmoid exactly, $\operatorname{softplus}'(z)=\sigma(z)$, so the derivative of the per-sample loss with respect to its own logit is simply the residual between the predicted probability and the true label, $\partial\ell_i/\partial z_i=p_i-y_i$ where $p_i=\sigma(z_i)$. Applying the chain rule through $z_i=\beta_0+\sum_j\beta_j x_{ij}$, using $\partial z_i/\partial\beta_0=1$ and $\partial z_i/\partial\beta_j=x_{ij}$:

With $X_b$ the $n\times(d+1)$ design matrix formed by prepending a column of ones to $X$, and $p-y$ the length-$n$ residual vector, both partial derivatives combine into one matrix expression:

The residual $p_i-y_i$ is the exact same quantity that appears in ordinary least-squares gradients, $\hat y_i-y_i$; the only difference is that $p_i$ has already been passed through the sigmoid. Differentiating the gradient a second time uses $\sigma'(z)=\sigma(z)\bigl(1-\sigma(z)\bigr)$, so $\partial p_i/\partial\beta_0=p_i(1-p_i)$ and $\partial p_i/\partial\beta_j=p_i(1-p_i)\,x_{ij}$:

Collecting every entry into a $(d+1)\times(d+1)$ matrix, with $W=\operatorname{diag}\bigl(p_1(1-p_1),\dots,p_n(1-p_n)\bigr)$:

The smoke test below computes exactly this gradient vector and Hessian matrix for a tiny fixed batch, and the second post in the roadmap is built entirely around solving with this $H$ rather than approximating it.

Convexity, Newton/IRLS, and Gradient Descent

Every diagonal entry of $W$ is $p_i(1-p_i)$, a product of two numbers in $(0,1)$, so $W$ is positive semidefinite for any finite $\beta$. For any direction vector $v$, the quadratic form $v^\top H v=\frac1n\sum_i p_i(1-p_i)\bigl((X_bv)_i\bigr)^2\ge0$, which makes $J(\beta)$ convex everywhere, not merely near a particular optimum. This is not automatic: composing squared error directly with a sigmoid output, rather than using log-loss, produces a non-convex objective with spurious local minima, which is exactly why log-loss, the Bernoulli negative log-likelihood, is the loss paired with this link function rather than an arbitrary alternative. Full column rank makes the finite-data objective strictly convex while all weights remain positive, so at most one finite minimizer exists. Convexity characterizes that optimum but does not guarantee that an arbitrary Newton or gradient step reaches it: gradient descent still needs a stable learning rate, and Newton/IRLS needs a nonsingular Hessian plus damping or line search when a full step is unsafe. Separation can remove the finite minimizer entirely, which is the failure mode covered next.

Newton's method updates $\beta\leftarrow\beta-H^{-1}\nabla J$. Substituting the matrix forms above:

That expression is exactly a weighted least-squares normal-equation solve against a working response $z^{*}=X_b\beta^{(t)}+W^{-1}(y-p)$, which is why the algorithm is called Iteratively Reweighted Least Squares: each iteration reuses the same weighted-normal-equations solver the linear-regression series already built for Huber and quantile regression, only with a different weight formula and a different working response, recomputed after every step as $\beta$ (and therefore $p$ and $W$) changes. Gradient descent instead updates $\beta\leftarrow\beta-\alpha\nabla J(\beta)$ directly, at $O(nd)$ cost per step with no matrix solve, against Newton/IRLS's $O(nd^2+d^3)$ per step; gradient descent needs many more iterations to reach the same tolerance since it uses no curvature information, while Newton's method converges quadratically once close enough to the optimum, at the cost of forming and solving against $H$ every single iteration. The second post in the roadmap covers exactly this trade-off in the context of choosing $d$-dependent solver cost against convergence speed.

Regularization and Separation

When the training classes are linearly separable, or close enough to it that only a handful of rows sit on the wrong side of some hyperplane, the unpenalized maximum-likelihood fit has no finite optimum. Because the sigmoid asymptotes toward $0$ and $1$ without ever reaching them, log-loss can always be driven closer to zero by scaling $\beta$ up further in the separating direction, so the coefficient magnitudes diverge toward infinity across Newton/IRLS iterations, and the Hessian $H=\frac1n X_b^\top WX_b$ becomes numerically singular in that same direction as every $p_i(1-p_i)$ weight collapses toward zero. Adding an $L_2$ penalty $\tfrac{\lambda}{2}\lVert\beta\rVert_2^2$ (excluding the intercept, by convention) adds $\lambda I$ to the slope block of $H$. With both classes represented, this supplies a finite, unique slope solution under separation while the likelihood anchors the unpenalized intercept; a one-class training set still needs an explicit validation policy because its intercept remains unidentifiable without an intercept penalty or prior. An $L_1$ penalty instead induces exact sparsity in $\beta$, at the cost of the objective no longer being differentiable at $\beta_j=0$, so a plain Newton step no longer applies and coordinate descent or a proximal method is needed instead. The third post in the roadmap implements both, and the companion Regularization, Identifiability, and Perfect Separation post works through the separation failure mode and its diagnostics in full.

Class Imbalance, Thresholds, and Cost

Class imbalance in the training labels does not, by itself, require any change to fit: as long as the training rows are a genuine random sample of the population the model will be applied to, the Bernoulli likelihood already accounts for the true class proportions, and $\hat p$ remains a meaningful probability estimate at any imbalance ratio. What imbalance does change is which threshold and which metric are trustworthy to report. On a $99{:}1$ imbalanced problem, a classifier that always predicts the majority class reaches $99\%$ accuracy at the default $t=0.5$ while never once identifying the minority class, which is why accuracy at a fixed threshold is a poor summary once the classes are far from balanced, and why the final post covers precision, recall, and threshold-free ranking metrics such as ROC-AUC and PR-AUC alongside it. Separately, when the two kinds of mistake carry different real costs, a false positive costing $c_{FP}$ and a false negative costing $c_{FN}$, the cost-minimizing threshold shifts away from $0.5$ to $t^{*}=c_{FP}/(c_{FP}+c_{FN})$, derived by setting the expected cost of predicting $1$ equal to the expected cost of predicting $0$ at the decision boundary. Changing $t$ this way only moves the boundary in feature space; it never touches $\hat\beta$ or $\hat p$, which is the same probability-versus-decision separation established earlier and demonstrated concretely below. A distinct technique, per-sample class weighting inside fit itself, rescales each row's contribution to the gradient and Hessian sums and does change $\hat\beta$; the companion Decision Thresholds, Class Imbalance, and Error Costs post separates the two techniques in detail.

The animation below fixes one small two-feature model and dataset, $z=-5+x_1-x_2$, six labeled points, and sweeps the threshold $t$ through five values from $0.1$ to $0.9$. Each probability $\hat p_i=\sigma(z_i)$ is computed once and never recomputed; only the decision boundary, derived from $\operatorname{logit}(t)$ as the offset added to the fixed linear score, slides through the fixed point cloud, changing which points fall on the predicted-$1$ side of the line without changing a single probability:

Generalizing to Multiclass Softmax

With $K$ classes instead of two, each row gets $K$ linear scores instead of one, $z_{ik}=\beta_{0k}+\sum_j\beta_{jk}x_{ij}$ for $k=1,\dots,K$, and the sigmoid is replaced by softmax, which normalizes all $K$ scores into a probability distribution at once instead of squashing one score into a single probability:

The implementation never evaluates that expression with raw logits when their magnitude can be large. With $m_i=\max_k z_{ik}$, the algebraically identical stable form $p_{ik}=e^{z_{ik}-m_i}/\sum_l e^{z_{il}-m_i}$ makes every exponent non-positive and leaves the normalized probabilities unchanged.

Log-loss generalizes to categorical cross-entropy, $-\sum_k y_{ik}\ln p_{ik}$ with $y_{ik}$ a one-hot indicator, and the gradient with respect to each $z_{ik}$ reduces to the same residual form found above, $p_{ik}-y_{ik}$, generalizing the binary case rather than replacing it; setting $K=2$ and fixing one class's scores to zero recovers the sigmoid exactly. The fourth post implements this directly alongside a one-vs-rest baseline built from $K$ independent binary fits, comparing the two on the same dataset; the companion Multiclass Logistic and Softmax Regression post covers the identifiability constraint that makes the all-classes parameterization above overcomplete by one class's worth of parameters.

Metrics and Calibration

A fitted classifier needs to be judged by more than training log-loss before it is trusted for a threshold-based decision. Accuracy, precision, recall, and the full confusion matrix all depend on whichever threshold is in force, while ROC-AUC and PR-AUC rank predictions without committing to any single cutoff, summarizing how well $\hat p$ orders positives ahead of negatives regardless of where a threshold eventually lands. A separate question, calibration, asks whether $\hat p$ means what it claims to mean: among all rows where the model predicts $\hat p\approx0.7$, do roughly $70\%$ actually carry $y=1$. Minimizing log-loss encourages calibrated probabilities on the training distribution but does not guarantee them out-of-sample, particularly once regularization or model misspecification is present, which is why reliability diagrams, the Brier score, and post-hoc corrections such as Platt scaling or isotonic regression exist as a separate diagnostic layer. The final post in the roadmap implements all of the above; the companion Confusion Matrix, ROC-AUC, PR-AUC, and Log-Loss and Probability Calibration and Reliability posts cover the statistical reasoning behind each metric in full.

Series Roadmap

Every post below shares the same Matrix alias, the same stableSigmoid/softplus free functions, the same LogisticRegressor struct with its predictLogit/predictProba/predict trio, and a gradient-and-Hessian routine shaped like evaluateBatch in the smoke test below; only the estimation method inside fit, and in the fourth post the shape of the output itself, changes from post to post.

Testing Invariants and Edge Cases

Every post in this sub-series is expected to hold the following invariants, checked directly rather than assumed, and the smoke test below exercises each one against the shared LogisticRegressor contract:

• Shape validation. fit rejects a row-count mismatch between $X$ and $y$, an empty $X$, and a ragged $X$ where rows disagree on column count, each with std::invalid_argument naming the offending row.
• Label validation. fit rejects any label outside $\{0,1\}$, since the Bernoulli likelihood above has no meaning for any other value; stableLogLoss rejects the same condition again on its own terms, rather than trusting a caller to have already validated upstream.
• Non-finite input rejection. fit rejects a non-finite entry anywhere in $X$, setParams rejects a non-finite intercept or coefficient, and predictLogit rejects a non-finite feature value, each with std::invalid_argument naming the offending row/column or index, since a single NaN or Inf silently poisons every arithmetic expression built from it rather than failing where it was introduced.
• Fitted-state validation. predictLogit, and everything built on it, throws std::logic_error if called before fit, never silently returning a zero or default score.
• Threshold validation. predict rejects a threshold that is not a finite value strictly inside the open interval $(0,1)$, since either endpoint would make the decision independent of the model and a NaN threshold would fail every ordered comparison silently in the caller's favor.
• Extreme-logit stability. stableSigmoid, softplus, and every loss/gradient/Hessian entry built from them stay finite for $z=\pm1000$ and beyond, verified with std::isfinite rather than by inspection.
• Property-style finiteness. for a spread of logits from $-1000$ to $1000$ and both possible labels, the mean loss, the full gradient vector, and the full Hessian matrix are all finite, checked in a loop rather than at a single hand-picked value; this sweep is unaffected by the new non-finite-input checks above, since every logit in it is a deliberately extreme but still finite double.

Standard-Library Smoke Test

The complete program below compiles with g++ -std=c++17 -Wall -Wextra -Werror -O2 and nothing else: standard headers only, no external dependency, no randomness anywhere, so its output is exactly reproducible. It walks through extreme logits and the stable-versus-naive log-loss contrast described above, evaluates loss, gradient, and Hessian on a fixed six-row two-feature batch with a fixed illustrative (not yet fitted) coefficient vector, exercises every validation path listed above, and finishes with the property-style finiteness sweep across nine logits from $-1000$ to $1000$:

#include <cmath>
#include <cstdio>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;

using Matrix = vector<vector<double>>;

// Numerically stable logistic sigmoid: branch on the sign of z so the
// argument passed to exp() never exceeds zero, whichever branch runs.
double stableSigmoid(double z) {
    if (z >= 0.0) {
        double e = exp(-z);
        return 1.0 / (1.0 + e);
    }
    double e = exp(z);
    return e / (1.0 + e);
}

// Numerically stable softplus, log(1 + e^z), via log1p and the same
// branch-on-sign discipline as stableSigmoid.
double softplus(double z) {
    if (z >= 0.0) return z + log1p(exp(-z));
    return log1p(exp(z));
}

// Per-sample log-loss expressed directly from the logit z rather than from
// an already-rounded probability: y=1 costs softplus(-z), y=0 costs
// softplus(z). Never divides, never takes log(0). Rejects any label
// outside {0,1} explicitly rather than silently treating it as y=0, since
// this function is called directly by evaluateBatch and could be called
// directly by later posts too, not only after validateFit has already run.
double stableLogLoss(double z, int y) {
    if (y != 0 && y != 1)
        throw invalid_argument("stableLogLoss: y must be exactly 0 or 1");
    return (y == 1) ? softplus(-z) : softplus(z);
}

// The textbook combined formula -[y*log(p) + (1-y)*log(1-p)] evaluated from
// an already-computed probability p. Included only to demonstrate the
// failure mode stableLogLoss avoids: once p rounds to exactly 0.0 or 1.0,
// the zero-weighted term becomes 0 * (-inf), which is NaN in IEEE 754.
double naiveLogLoss(double p, int y) {
    return -(y * log(p) + (1 - y) * log(1.0 - p));
}

struct LogisticRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;

    // Rejects a shape mismatch, a ragged row, a non-finite feature value
    // anywhere in X, and any label outside {0,1}, each with a specific
    // invalid_argument naming the offending row/column so a caller never
    // has to guess which entry failed.
    static void validateFit(const Matrix& X, const vector<int>& y) {
        if (X.size() != y.size())
            throw invalid_argument("fit: X and y must have the same number of rows");
        if (X.empty())
            throw invalid_argument("fit: X must have at least one row");
        size_t d = X[0].size();
        if (d == 0)
            throw invalid_argument("fit: X must have at least one column");
        for (size_t i = 0; i < X.size(); ++i) {
            if (X[i].size() != d)
                throw invalid_argument("fit: row " + to_string(i) + " has " +
                                        to_string(X[i].size()) + " columns, expected " + to_string(d));
            for (size_t j = 0; j < X[i].size(); ++j)
                if (!isfinite(X[i][j]))
                    throw invalid_argument("fit: X[" + to_string(i) + "][" + to_string(j) +
                                            "] is not finite (NaN or Inf)");
        }
        for (size_t i = 0; i < y.size(); ++i)
            if (y[i] != 0 && y[i] != 1)
                throw invalid_argument("fit: y[" + to_string(i) +
                                       "] must be 0 or 1, got " + to_string(y[i]));
    }

    // Sets already-estimated parameters directly; later posts replace this
    // with an actual solver (Newton/IRLS, gradient descent, penalized fits),
    // but every solver still has to hand its result through this same gate:
    // a non-finite intercept or coefficient is rejected outright rather than
    // stored, since every downstream prediction would otherwise silently
    // inherit a NaN or Inf logit.
    void setParams(double b0, const vector<double>& b) {
        if (!isfinite(b0))
            throw invalid_argument("setParams: intercept must be finite");
        for (size_t j = 0; j < b.size(); ++j)
            if (!isfinite(b[j]))
                throw invalid_argument("setParams: coef[" + to_string(j) + "] must be finite");
        intercept = b0;
        coef = b;
        nFeatures = (int)b.size();
        fitted = true;
    }

    void requireFitted() const {
        if (!fitted) throw logic_error("predict called before fit");
    }

    // Requires a fitted model, an exact feature-count match, and every
    // feature to be finite; a NaN or Inf feature would silently poison the
    // dot product below rather than fail loudly at the boundary.
    double predictLogit(const vector<double>& x) const {
        requireFitted();
        if ((int)x.size() != nFeatures)
            throw invalid_argument("predictLogit: expected " + to_string(nFeatures) +
                                    " features, got " + to_string(x.size()));
        for (size_t j = 0; j < x.size(); ++j)
            if (!isfinite(x[j]))
                throw invalid_argument("predictLogit: feature " + to_string(j) +
                                        " is not finite (NaN or Inf)");
        double z = intercept;
        for (int j = 0; j < nFeatures; ++j) z += coef[j] * x[j];
        if (!isfinite(z))
            throw runtime_error("predictLogit: score became non-finite");
        return z;
    }

    double predictProba(const vector<double>& x) const { return stableSigmoid(predictLogit(x)); }

    // A threshold of exactly 0 or 1 would make the decision independent of
    // the fitted model, and a non-finite threshold (NaN in particular,
    // which fails every ordered comparison) is rejected explicitly rather
    // than relying on the comparison below to fail silently in the caller's
    // favor.
    int predict(const vector<double>& x, double threshold = 0.5) const {
        if (!(isfinite(threshold) && threshold > 0.0 && threshold < 1.0))
            throw invalid_argument("predict: threshold must be a finite value strictly between 0 and 1");
        return predictProba(x) >= threshold ? 1 : 0;
    }
};

// Mean log-loss, gradient (length 1+d, index 0 = d/d(beta0)), and Hessian
// ((1+d) x (1+d)) for a batch, computed from stable logits/probabilities.
struct BatchResult {
    double meanLoss = 0.0;
    vector<double> grad;
    Matrix hessian;
};

BatchResult evaluateBatch(const LogisticRegressor& model, const Matrix& X, const vector<int>& y) {
    LogisticRegressor::validateFit(X, y);
    if ((int)X[0].size() != model.nFeatures)
        throw invalid_argument("evaluateBatch: X column count does not match model.nFeatures");

    int n = (int)X.size();
    int p1 = model.nFeatures + 1;
    BatchResult r;
    r.grad.assign(p1, 0.0);
    r.hessian.assign(p1, vector<double>(p1, 0.0));

    for (int i = 0; i < n; ++i) {
        double z = model.predictLogit(X[i]);
        double p = stableSigmoid(z);
        r.meanLoss += stableLogLoss(z, y[i]);

        vector<double> xAug(p1);
        xAug[0] = 1.0;
        for (int j = 0; j < model.nFeatures; ++j) xAug[j + 1] = X[i][j];

        double residual = p - y[i];
        double w = p * (1.0 - p);
        for (int a = 0; a < p1; ++a) {
            r.grad[a] += residual * xAug[a];
            for (int b = 0; b < p1; ++b) r.hessian[a][b] += w * xAug[a] * xAug[b];
        }
    }
    r.meanLoss /= n;
    for (double& g : r.grad) g /= n;
    for (auto& row : r.hessian)
        for (double& h : row) h /= n;
    return r;
}

bool allFinite(const vector<double>& v) {
    for (double x : v) if (!isfinite(x)) return false;
    return true;
}
bool allFinite(const Matrix& m) {
    for (const auto& row : m) if (!allFinite(row)) return false;
    return true;
}

int main() {
    printf("=== 1. Extreme logits: stable sigmoid and softplus ===\n");
    double extremeZ[2] = {1000.0, -1000.0};
    for (double z : extremeZ) {
        double p = stableSigmoid(z);
        printf("z=%.1f   sigma(z)=%.6f   softplus(z)=%.6f   softplus(-z)=%.6f   finite=%s\n",
               z, p, softplus(z), softplus(-z), (isfinite(p) ? "yes" : "no"));
    }

    printf("\n=== 2. Stable vs naive log-loss at saturation ===\n");
    for (double z : extremeZ) {
        double p = stableSigmoid(z);
        for (int y = 0; y <= 1; ++y) {
            double stable = stableLogLoss(z, y);
            double naive = naiveLogLoss(p, y);
            printf("z=%.1f y=%d   stable=%.6f   naive=%.6f\n", z, y, stable, naive);
        }
    }

    printf("\n=== 3. Tiny-batch loss, gradient, and Hessian ===\n");
    LogisticRegressor model;
    model.setParams(-0.5, {0.8, -0.3});
    Matrix X = {{1.0, 3.0}, {2.0, 1.0}, {3.0, 2.0}, {4.0, 0.0}, {1.5, 4.0}, {5.0, 1.0}};
    vector<int> y = {0, 0, 1, 1, 0, 1};
    for (size_t i = 0; i < X.size(); ++i) {
        double z = model.predictLogit(X[i]);
        double p = model.predictProba(X[i]);
        printf("row %zu: x=(%.1f,%.1f) y=%d   z=%.6f   p=%.6f   loss=%.6f\n",
               i, X[i][0], X[i][1], y[i], z, p, stableLogLoss(z, y[i]));
    }
    BatchResult batch = evaluateBatch(model, X, y);
    printf("mean loss = %.6f\n", batch.meanLoss);
    printf("gradient  = [%.6f, %.6f, %.6f]\n", batch.grad[0], batch.grad[1], batch.grad[2]);
    for (int a = 0; a < 3; ++a)
        printf("hessian row %d = [%.6f, %.6f, %.6f]\n", a, batch.hessian[a][0], batch.hessian[a][1], batch.hessian[a][2]);

    printf("\n=== 4. Explicit validation ===\n");
    try {
        Matrix badX = {{1.0, 2.0}, {3.0}};
        vector<int> badY = {0, 1};
        LogisticRegressor::validateFit(badX, badY);
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        Matrix okX = {{1.0}, {2.0}};
        vector<int> badLabels = {0, 2};
        LogisticRegressor::validateFit(okX, badLabels);
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        Matrix nanX = {{1.0, numeric_limits<double>::infinity()}, {2.0, 3.0}};
        vector<int> okLabels = {0, 1};
        LogisticRegressor::validateFit(nanX, okLabels);
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        LogisticRegressor badParams;
        badParams.setParams(numeric_limits<double>::quiet_NaN(), {1.0, 2.0});
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        LogisticRegressor unfitted;
        unfitted.predictProba({1.0});
    } catch (const logic_error& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        model.predictProba({numeric_limits<double>::quiet_NaN(), 1.0});
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        model.predict({1.0, 1.0}, 1.5);
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }
    try {
        stableLogLoss(0.0, 2);
    } catch (const invalid_argument& e) {
        printf("caught: %s\n", e.what());
    }

    printf("\n=== 5. Property tests: finite gradient/Hessian under extreme logits ===\n");
    double testZ[9] = {-1000.0, -700.0, -50.0, -1.0, 0.0, 1.0, 50.0, 700.0, 1000.0};
    bool allPass = true;
    for (double z : testZ) {
        LogisticRegressor probe;
        probe.setParams(z, {0.0});
        for (int yy = 0; yy <= 1; ++yy) {
            BatchResult r = evaluateBatch(probe, Matrix{{1.0}}, vector<int>{yy});
            bool pass = isfinite(r.meanLoss) && allFinite(r.grad) && allFinite(r.hessian);
            allPass = allPass && pass;
            printf("z=%7.1f y=%d   loss=%12.6f   grad finite=%s   hessian finite=%s   %s\n",
                   z, yy, r.meanLoss, (allFinite(r.grad) ? "yes" : "no"),
                   (allFinite(r.hessian) ? "yes" : "no"), (pass ? "PASS" : "FAIL"));
        }
    }
    printf("overall: %s\n", (allPass ? "ALL PROPERTY TESTS PASS" : "PROPERTY TEST FAILURE"));
    return 0;
}

Compiled with g++ -std=c++17 -Wall -Wextra -Werror -O2, this prints:

Section 1 confirms the stable sigmoid stays finite at $z=\pm1000$, both branches rounding to exactly $1.0$ and $0.0$ rather than overflowing. Section 2 is the concrete version of the numerical-stability argument above: at $z=1000$ with $y=0$, the naive combined formula evaluates to inf because $\ln(1-p)$ diverges once $p$ rounds to exactly $1.0$, and at $z=1000$ with $y=1$ it evaluates to nan because the zero-weighted term becomes $0\cdot(-\infty)$; the stable softplus form reports $1000.000000$ and $0.000000$ respectively in both cases, finite and directly interpretable. Section 3 matches the gradient and Hessian derivation above exactly: the printed gradient is $\frac1nX_b^\top(p-y)$ and the printed Hessian is $\frac1nX_b^\top WX_b$ for the fixed six-row batch and fixed coefficients. Section 4 exercises every validation path from the invariants list, eight in total: a ragged row, a bad label, a non-finite entry inside $X$, a non-finite parameter passed to setParams, a call before fit, a non-finite feature passed to predictProba, an out-of-range threshold, and a bad label passed straight to stableLogLoss itself, each caught by its specific exception type, std::invalid_argument or std::logic_error, never a catch-all. Section 5 is the property-style sweep: eighteen combinations of nine logits from $-1000$ to $1000$ and both labels, every one still a finite, legitimately extreme double rather than a NaN or Inf, so none of it is rejected by the new checks above; every combination reports a finite loss, a finite gradient, and a finite Hessian, ending in the single deterministic line overall: ALL PROPERTY TESTS PASS.

Why Each Post Matters