Neel Mishra
← All Posts
From Scratch · C++ · ML · Regression · Logistic Regression · Part 2

Newton’s Method and IRLS in C++

July 2026 · 42 min read

Part 1 optimized Bernoulli log-loss with first-order information. This implementation adds curvature. Each iteration forms a weighted Gram matrix, solves one linear system without explicit inversion, and accepts a damped Newton step only when stable log-loss does not increase. The classifier contract, deterministic data, and extreme-logit primitives remain the same, making the solver comparison direct.

Gradient and Hessian

For the mean objective $L(\beta)=n^{-1}\sum_i[\operatorname{softplus}(z_i)-y_i z_i]$ with $z=X_b\beta$, the gradient is

$$g=\nabla L(\beta)=\frac{1}{n}X_b^\top(p-y).$$

Since $d p_i/d z_i=p_i(1-p_i)$, another derivative gives

$$H=\nabla^2 L(\beta)=\frac{1}{n}X_b^\top W X_b,\qquad W=\operatorname{diag}(w_i),\quad w_i=p_i(1-p_i).$$

Each $w_i$ lies in $[0,1/4]$. Thus $H$ is positive semidefinite. It is positive definite when the weighted design has full column rank. Newton’s direction solves $H\Delta=g$, followed by $\beta_{\text{new}}=\beta-\Delta$. No matrix inverse is needed.

Why the Same Step Is IRLS

Define the working response

$$r_i=z_i+\frac{y_i-p_i}{w_i}.$$

Weighted least squares minimizes $\sum_i w_i(r_i-x_{bi}^\top\beta_{\text{new}})^2$. Its normal equations are

$$X_b^\top W X_b\,\beta_{\text{new}}=X_b^\top W r.$$

Substituting $r=z+W^{-1}(y-p)$ and $z=X_b\beta$ gives

$$X_b^\top W X_b(\beta_{\text{new}}-\beta)=X_b^\top(y-p),$$

which is exactly $\beta_{\text{new}}=\beta-H^{-1}g$ after the common $1/n$ scale cancels. Newton’s method and iteratively reweighted least squares are therefore two views of the same logistic update: curvature in parameter space and weighted least squares on a local working target.

Weighted Normal Equations Without Inversion

The implementation accumulates the symmetric system directly. An implicit intercept value $x_{i0}=1$ keeps the public matrix free of a bias column:

$$H_{jk}=\frac{1}{n}\sum_i w_i x_{bij}x_{bik},\qquad g_j=\frac{1}{n}\sum_i(p_i-y_i)x_{bij}.$$

A partial-pivot Gaussian solver chooses the largest available pivot magnitude in each column. Its tolerance scales with matrix magnitude and machine epsilon. A missing usable pivot reports a near-singular system instead of returning unstable coefficients. Explicit inversion would allocate and compute an entire inverse only to multiply it by one vector, with more rounding opportunities and no benefit.

Damping and Monotonic Backtracking

Pure Newton steps are locally fast but can be unreliable far from the optimum or when curvature is weak. This implementation adds a small damping value $\lambda$ to the Hessian diagonal before solving:

$$(H+\lambda I)\Delta=g.$$

Damping improves pivots and limits movement in weakly identified directions. A backtracking factor $s$ starts at one and repeatedly halves until $L(\beta-s\Delta)\le L(\beta)$. This preserves monotonic accepted losses. The method still approaches an ordinary Newton step near a well-conditioned optimum because full steps are accepted there.

Tiny Weights, Conditioning, and Separation

Extreme probabilities make $p_i(1-p_i)$ round toward zero. The program floors weights at $10^{-12}$ while building the numerical system. This does not change stable loss or reported probabilities; it prevents a row from creating an infinite IRLS working response and supplies a bounded curvature approximation. Damping then protects directions that remain weak.

Tiny weights can also signal complete or quasi-complete separation. For separated data, unpenalized log-loss decreases while coefficient norm grows without a finite minimizer. The implementation detects parameter magnitude above a configured bound and reports a separation-or-divergence error. The regularization post supplies the principled finite solution. A pivot failure can instead indicate duplicated columns or severe scaling; feature scaling, QR-based weighted least squares, or regularization is then preferable.

Explicit Convergence Contract

Data-Derived Newton and IRLS Iterations

The animation recomputes probabilities, weights, the $2\times2$ damped Hessian, and the Newton direction from the same eight one-feature rows used in both complete programs. The probability curve and minimum weight update at every accepted iteration. The label reports the IRLS interpretation of the identical step.

Newton Curvature and IRLS Reweighting

Validation and Failure Behavior

Data checks match Part 1: nonempty rectangular finite $X$, matching target length, and labels exactly zero or one. Configuration requires positive finite tolerance, damping, weight floor, and parameter bound, plus a positive iteration limit. Prediction requires fitted state, exact feature width, finite values, and a threshold strictly inside $(0,1)$. The linear solver independently validates square finite input and a matching right-hand side. These checks isolate caller errors from numerical solver failures.

Complexity and Solver Comparison

MethodPer iterationExtra memoryBehavior
Batch gradient descent$O(nd)$$O(d)$Cheap iterations; scale-sensitive; often many iterations
Newton / IRLS normal equations$O(nd^2+d^3)$$O(d^2)$Curvature-aware; usually few iterations; condition number can be squared
IRLS with QR$O(nd^2)$$O(nd)$ or streamed variantsMore stable weighted least squares; more implementation work
Quasi-Newtontypically $O(nd+d^2)$$O(d^2)$A compromise when a full Hessian is expensive

Normal-equation IRLS is appropriate for modest, well-scaled feature counts. Gradient or limited-memory methods become attractive as $d$ grows. A production solver generally combines regularization, robust factorization, and richer line-search diagnostics.

Complete C++17 Program

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#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);
}

std::vector<double> solveLinearSystem(Matrix A, std::vector<double> b) {
    const std::size_t n = A.size();
    if (n == 0 || b.size() != n)
        throw std::invalid_argument("solveLinearSystem: incompatible dimensions");
    double scale = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        if (A[i].size() != n)
            throw std::invalid_argument("solveLinearSystem: A must be square");
        for (double value : A[i]) {
            if (!std::isfinite(value))
                throw std::invalid_argument("solveLinearSystem: non-finite matrix value");
            scale = std::max(scale, std::fabs(value));
        }
        if (!std::isfinite(b[i]))
            throw std::invalid_argument("solveLinearSystem: non-finite right side");
    }
    const double tolerance = std::numeric_limits<double>::epsilon() *
                             std::max(1.0, scale) * n * 100.0;
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        for (std::size_t i = k + 1; i < n; ++i)
            if (std::fabs(A[i][k]) > std::fabs(A[pivot][k])) pivot = i;
        if (std::fabs(A[pivot][k]) <= tolerance)
            throw std::runtime_error(
                "solveLinearSystem: singular or ill-conditioned weighted system");
        if (pivot != k) {
            std::swap(A[pivot], A[k]);
            std::swap(b[pivot], b[k]);
        }
        for (std::size_t i = k + 1; i < n; ++i) {
            const double factor = A[i][k] / A[k][k];
            A[i][k] = 0.0;
            for (std::size_t j = k + 1; j < n; ++j)
                A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    std::vector<double> x(n, 0.0);
    for (std::size_t ii = n; ii-- > 0;) {
        double value = b[ii];
        for (std::size_t j = ii + 1; j < n; ++j)
            value -= A[ii][j] * x[j];
        x[ii] = value / A[ii][ii];
        if (!std::isfinite(x[ii]))
            throw std::runtime_error("solveLinearSystem: non-finite solution");
    }
    return x;
}

struct NewtonConfig {
    int maxIterations = 100;
    double tolerance = 1e-10;
    double damping = 1e-8;
    double minimumWeight = 1e-12;
    double maxAbsParameter = 1e6;
};

struct NewtonRecord {
    int iteration;
    double loss;
    double gradientNorm;
    double stepScale;
    double minimumRawWeight;
};

struct LogisticNewton {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    std::vector<double> coef;
    std::vector<NewtonRecord> history;

    static void validateData(const Matrix& X, const std::vector<int>& y) {
        if (X.empty() || X.front().empty())
            throw std::invalid_argument("fit: X must be non-empty");
        if (X.size() != y.size())
            throw std::invalid_argument("fit: X row count must equal y size");
        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 result = total / static_cast<double>(X.size());
        if (!std::isfinite(result))
            throw std::runtime_error("meanLoss: non-finite objective");
        return result;
    }

    void fit(const Matrix& X, const std::vector<int>& y,
             const NewtonConfig& config = {}) {
        validateData(X, y);
        if (config.maxIterations <= 0)
            throw std::invalid_argument("fit: maxIterations must be positive");
        if (!std::isfinite(config.tolerance) || config.tolerance <= 0.0 ||
            !std::isfinite(config.damping) || config.damping <= 0.0 ||
            !std::isfinite(config.minimumWeight) || config.minimumWeight <= 0.0 ||
            config.minimumWeight > 0.25 ||
            !std::isfinite(config.maxAbsParameter) ||
            config.maxAbsParameter <= 0.0)
            throw std::invalid_argument("fit: numerical settings are invalid");

        const std::size_t p = X.front().size() + 1;
        std::vector<double> beta(p, 0.0);
        std::vector<NewtonRecord> nextHistory;
        double currentLoss = meanLoss(X, y, beta);
        bool converged = false;

        for (int iteration = 0; iteration < config.maxIterations; ++iteration) {
            Matrix hessian(p, std::vector<double>(p, 0.0));
            std::vector<double> gradient(p, 0.0);
            double minimumRawWeight = 0.25;
            for (std::size_t i = 0; i < X.size(); ++i) {
                std::vector<double> row(1, 1.0);
                row.insert(row.end(), X[i].begin(), X[i].end());
                const double probability = stableSigmoid(score(beta, X[i]));
                const double residual = probability - y[i];
                const double rawWeight = probability * (1.0 - probability);
                minimumRawWeight = std::min(minimumRawWeight, rawWeight);
                const double weight = std::max(rawWeight, config.minimumWeight);
                for (std::size_t j = 0; j < p; ++j) {
                    gradient[j] += residual * row[j];
                    for (std::size_t k = j; k < p; ++k)
                        hessian[j][k] += weight * row[j] * row[k];
                }
            }
            const double n = static_cast<double>(X.size());
            double gradientSquaredNorm = 0.0;
            for (std::size_t j = 0; j < p; ++j) {
                gradient[j] /= n;
                gradientSquaredNorm += gradient[j] * gradient[j];
                for (std::size_t k = j; k < p; ++k) {
                    hessian[j][k] /= n;
                    hessian[k][j] = hessian[j][k];
                }
                hessian[j][j] += config.damping;
            }
            const double gradientNorm = std::sqrt(gradientSquaredNorm);
            if (!std::isfinite(gradientNorm))
                throw std::runtime_error("fit: non-finite gradient");
            if (gradientNorm <= config.tolerance) {
                nextHistory.push_back(
                    {iteration, currentLoss, gradientNorm, 0.0, minimumRawWeight});
                converged = true;
                break;
            }

            const std::vector<double> direction =
                solveLinearSystem(hessian, gradient);
            double stepScale = 1.0;
            std::vector<double> candidate(p);
            double candidateLoss = std::numeric_limits<double>::infinity();
            bool accepted = false;
            for (int trial = 0; trial < 30; ++trial) {
                bool bounded = true;
                for (std::size_t j = 0; j < p; ++j) {
                    candidate[j] = beta[j] - stepScale * direction[j];
                    bounded = bounded && std::isfinite(candidate[j]) &&
                              std::fabs(candidate[j]) <= config.maxAbsParameter;
                }
                if (bounded) {
                    candidateLoss = meanLoss(X, y, candidate);
                    if (candidateLoss <= currentLoss) {
                        accepted = true;
                        break;
                    }
                }
                stepScale *= 0.5;
            }
            if (!accepted)
                throw std::runtime_error(
                    "fit: no monotonic Newton step; possible separation");

            nextHistory.push_back(
                {iteration, currentLoss, gradientNorm, stepScale, minimumRawWeight});
            const double relativeChange =
                std::fabs(currentLoss - candidateLoss) / std::max(1.0, currentLoss);
            beta = candidate;
            currentLoss = candidateLoss;
            if (relativeChange <= config.tolerance) {
                converged = true;
                break;
            }
        }
        if (!converged)
            throw std::runtime_error(
                "fit: maximum iterations reached; scaling or separation may be responsible");

        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;
    }
};

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};

    LogisticNewton model;
    model.fit(X, y);

    int correct = 0;
    double fittedLoss = 0.0;
    for (std::size_t i = 0; i < X.size(); ++i) {
        const double probability = model.predictProba(X[i]);
        if (!(probability >= 0.0 && probability <= 1.0))
            throw std::runtime_error("main: probability range property failed");
        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");
    if (!std::isfinite(logLossFromLogit(-1000.0, 1)) ||
        !std::isfinite(logLossFromLogit(1000.0, 0)))
        throw std::runtime_error("main: extreme-logit loss property failed");

    std::cout << std::fixed << std::setprecision(6);
    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 << "Newton iterations: " << model.history.size() << '\n';
    std::cout << "property checks: passed\n";
}

Expected Output

Newton/IRLS and gradient descent minimize the same unpenalized objective, so their fitted parameters agree within convergence tolerance. Minor final-digit differences can follow from elementary-function implementations.

intercept: -0.641147
coefficient: 1.282293
mean log-loss: 0.312881
training accuracy: 0.750000
P(y=1 | x=0.5): 0.500000
Newton iterations: 6
property checks: passed

Beyond the Unpenalized Solver

The next post adds $L_1$ and $L_2$ penalties. Regularization controls variance, improves conditioning, and creates a finite target when separation would otherwise send coefficients toward infinity.