2.2. Logistic Regression

LogisticRegression is RustyML’s linear classifier for binary problems. It shares most of its machinery with LinearRegression: a weight vector, an optional intercept, gradient descent, and optional L1 or L2 penalties. It replaces the squared-error objective with the logistic loss. It also passes the linear score through a sigmoid. The output then reads as a probability. This model corresponds to scikit-learn’s sklearn.linear_model.LogisticRegression, stripped to the essentials. RustyML supports one binary class boundary and plain full-batch gradient descent. It has no built-in multiclass support. It also expresses regularization strength directly, not as the inverse C.

2.2.1. What the model actually optimizes

Each sample gets a linear score z = w * x (plus a bias when the model fits an intercept). The sigmoid sigmoid(z) = 1 / (1 + e^-z) maps that score into (0, 1). RustyML reads that value as the probability of the positive class. Training minimizes the mean binary cross-entropy between those probabilities and the labels. The optimizer is plain full-batch gradient descent. Every iteration computes the gradient over the entire training set, (1/n) * X^T * (sigmoid(X * w) - y). It then takes one step of size learning_rate against that gradient.

This choice has 2 consequences in practice. First, there is no stochastic sampling and no random initialization. Weights start at exactly zero, and the update is deterministic. 2 fits on the same data and hyperparameters produce bit-identical weights (the test suite verifies this). Second, full-batch gradient descent is a first-order method with one global step size. It is far more sensitive to feature scaling and to learning_rate than the quasi-Newton solvers (lbfgs, liblinear) that scikit-learn defaults to. This is the main behavioral difference from those solvers. It drives the practical advice in 2.2.6.

RustyML evaluates the loss in the numerically stable log-sum-exp form, instead of taking the log of a sigmoid. The formula is max(z, 0) - z * y + ln(1 + e^-|z|). Large-magnitude logits therefore do not overflow the loss computation. The weights can still overflow, which section 2.2.6 covers as a real failure mode. The per-iteration logits, gradient, and loss run through RustyML’s parallel GEMV and deterministic reduction primitives above their size gates. Results are therefore reproducible on a given machine. See 7.3. Performance Tuning and Parallelism for more.

2.2.2. Constructing the model

2 constructors exist. LogisticRegression::default() gives a reasonable starting point. LogisticRegression::new(...) sets every hyperparameter and validates each one up front.

use rustyml::machine_learning::LogisticRegression;

fn main() {
    // Defaults: fit_intercept = true, lr = 0.01, max_iter = 100, tol = 1e-4, no penalty
    let _a = LogisticRegression::default();

    // Explicit: new(fit_intercept, learning_rate, max_iterations, tolerance)
    let _b = LogisticRegression::new(true, 0.1, 1000, 1e-6).unwrap();
}

new returns Result<Self, Error>. It rejects an invalid hyperparameter immediately with Error::InvalidParameter, instead of letting it fail later at fit time.

The defaults are deliberately conservative. learning_rate = 0.01 with only 100 iterations rarely converges, except on trivial, well-scaled data. Treat default() as a smoke test, not a production configuration. Most real fits need a larger learning_rate (after standardizing the data) and a max_iterations value in the thousands.

Every stored value has a getter. 2 of them are the model’s after-fit diagnostics:

get_actual_iterations() gives the convergence check. Training stops when the loss change falls below tolerance between iterations, or when it reaches max_iterations. If the returned count equals max_iterations, assume the model did not converge to the tolerance. Raise max_iterations, raise learning_rate, or standardize the inputs instead of trusting that boundary.

2.2.3. The label contract: strictly 0 and 1

fit takes a feature matrix x (rows are samples, columns are features) and a target vector y. Both must be f64 arrays with the same storage type. The label domain is exact. Every entry of y must be 0.0 or 1.0. RustyML rejects anything else, such as 0.5, 2.0, or -1.0, with Error::InvalidInput. This matches the label convention used by SVC and LinearSVC, which also predict 0.0 and 1.0. If your labels use a different encoding, map them to {0, 1} first. For string or categorical labels, encode them first with 4.3. Label Encoding. Pick which class is positive, the 1, on purpose. That choice defines what precision, recall, and the probability output mean.

use ndarray::array;
use rustyml::machine_learning::LogisticRegression;

fn main() {
    // Logical AND, encoded as {0.0, 1.0}
    let x_train = array![[0.0, 0.0], [0.0, 1.0], [1.0, 0.0], [1.0, 1.0]];
    let y_train = array![0.0, 0.0, 0.0, 1.0];

    let mut model = LogisticRegression::new(true, 0.5, 500, 1e-7).unwrap();
    model.fit(&x_train, &y_train).unwrap();

    let preds = model.predict(&x_train).unwrap(); // Array1<f64> with values in {0.0, 1.0}
    println!("predictions: {:?}", preds);
    println!("iterations:  {:?}", model.get_actual_iterations());
}

When fit_intercept is true, the model prepends the bias as weight index 0. get_weights() then returns n_features + 1 values. When fit_intercept is false, it returns exactly n_features values. This indexing matters when you inspect the weights, because the regularizer treats the intercept specially (see below).

2.2.4. Predicting: hard labels versus probabilities

RustyML has 3 prediction entry points. The distinction between hard labels and probabilities is the part most people get wrong.

• predict(&x) -> Result<Array1<f64>, Error> returns hard class labels, 0.0 or 1.0. It thresholds the positive-class probability at 0.5.
• predict_proba(&x) -> Result<Array1<f64>, Error> returns the raw positive-class probability in (0, 1). It gives one value per sample, the output of the sigmoid.
• fit_predict(&mut self, x, y) runs fit, then predict on the same x, for a quick check on the training set.

predict is exactly predict_proba with a fixed >= 0.5 cutoff. The model has no threshold parameter. That fixed cutoff works for balanced problems, but the 0.5 boundary is a modeling choice, not a rule. When a false positive costs more or less than a false negative, or the classes are imbalanced, call predict_proba and set your own threshold.

use ndarray::array;
use rustyml::machine_learning::LogisticRegression;

fn main() {
    let x_train = array![[-3.0], [-2.0], [-1.0], [1.0], [2.0], [3.0]];
    let y_train = array![0.0, 0.0, 0.0, 1.0, 1.0, 1.0];

    let mut model = LogisticRegression::new(true, 0.3, 500, 1e-7).unwrap();
    model.fit(&x_train, &y_train).unwrap();

    let x_test = array![[-0.5], [0.5]];
    let proba = model.predict_proba(&x_test).unwrap(); // positive-class probability
    let default = model.predict(&x_test).unwrap();     // 0.5 threshold -> Array1<f64> of {0.0, 1.0}

    // A stricter operating point: only call positive when p >= 0.8
    let strict: Vec<i32> = proba.iter().map(|&p| if p >= 0.8 { 1 } else { 0 }).collect();

    println!("proba:   {:?}", proba);
    println!("default: {:?}", default);
    println!("strict:  {:?}", strict);
}

All 3 methods validate the input against the trained feature count, excluding the implicit bias column. Calling an unfitted model returns Error::NotFitted. A wrong number of columns returns Error::DimensionMismatch. A NaN or infinite entry returns Error::NonFinite. Pass features without a manual bias column. The model adds and removes that column internally to match how it trained.

2.2.5. Regularization

By default there is no penalty. This departs from scikit-learn, whose default is L2. Add a penalty with the builder method with_regularization. It consumes and returns the model, so it chains directly off new:

use ndarray::array;
use rustyml::machine_learning::{LogisticRegression, RegularizationType};

fn main() {
    let x = array![
        [-4.0, -3.0], [-3.0, -4.0], [-2.0, -1.0],
        [ 2.0,  1.0], [ 3.0,  4.0], [ 4.0,  3.0],
    ];
    let y = array![0.0, 0.0, 0.0, 1.0, 1.0, 1.0];

    let mut plain = LogisticRegression::new(true, 0.1, 2000, 1e-8).unwrap();
    plain.fit(&x, &y).unwrap();

    let mut ridge = LogisticRegression::new(true, 0.1, 2000, 1e-8)
        .unwrap()
        .with_regularization(RegularizationType::L2(5.0))
        .unwrap();
    ridge.fit(&x, &y).unwrap();

    // L2 norm of the feature weights, skipping the (unpenalized) bias at index 0
    let feature_norm = |m: &LogisticRegression| {
        m.get_weights().unwrap().iter().skip(1).map(|w| w * w).sum::<f64>()
    };
    println!("no penalty: {:.4}", feature_norm(&plain));
    println!("L2(5.0):    {:.4}", feature_norm(&ridge));
}

The f64 inside each variant is the penalty strength alpha. It must be non-negative and finite (alpha = 0 is accepted and equals no penalty). RegularizationType::L2(alpha) (ridge) adds alpha * 0.5 * ||w||^2 to the loss. It shrinks weights smoothly toward zero. RegularizationType::L1(alpha) (lasso) adds alpha * ||w||_1. It drives individual weights to exactly zero, giving a sparse model useful for feature selection.

That “exactly zero” result is literal. RustyML does not fold L1 into the gradient as alpha * sign(w). A sub-gradient step only approaches zero, so that approach gives no sparsity, however long it runs. Instead, the optimizer takes its ordinary gradient step and then applies a proximal step. Every feature weight is soft-thresholded by learning_rate * alpha, so a weight the data cannot support lands on 0.0 and stays there. This method is ISTA. It is what makes L1 an actual feature selector. Count the zeros with w.iter().skip(1).filter(|v| **v == 0.0).count().

3 implementation details change how you pick alpha:

• The intercept is never penalized. When the model fits an intercept, the penalty gradient starts at feature index 1, so the bias stays free to move. This is the standard, correct choice. The regularizer should not fight the model’s ability to shift its decision boundary. Skip index 0 when you compare weight norms, as the example does.
• The penalty is not divided by the sample count. The data term is the mean log-loss, but RustyML adds the penalty as the absolute value alpha * R(w). So alpha keeps a fixed meaning no matter how many rows the data has. Replicating a dataset leaves the regularized optimum unchanged (the test suite checks this invariance). Here, a larger alpha means more regularization, the opposite of scikit-learn’s inverse C.
• alpha transfers 1:1 from scikit-learn’s SGD estimators. It needs conversion from the rest. The objective above, mean data term plus undivided penalty, matches SGDClassifier exactly. An alpha tuned there moves over untouched. From LogisticRegression(C=c), use alpha = 1 / (c * n), where n is the number of training samples. The full conversion table lives on RegularizationType, shared with LinearRegression.

2.2.6. Standardization, class imbalance, and separable data

3 failure modes are common enough to name here.

Standardize your features. The optimizer is full-batch gradient descent with one global learning_rate. Features on very different scales therefore converge at very different rates. The step that suits a feature ranging over [0, 1] is far too small for one ranging over [0, 10000]. Training then crawls, and the loss plateaus before it reaches tolerance within max_iterations. Center and scale the data first with the tools in 4.2. Standardization and Normalization. This step gives the largest gain in both convergence speed and numerical stability. It also lets you use a useful learning_rate, such as 0.1 to 1.0, instead of the cautious default.

There is no class weighting. The model has no class_weight parameter, and predict fixes its threshold at 0.5. On imbalanced data, it can learn to predict only the majority class and still report a high accuracy. Use 2 defenses. Threshold predict_proba yourself, at an operating point chosen from a precision-recall or ROC analysis. Also evaluate with metrics that imbalance does not fool, such as balanced_accuracy, mcc, or roc_auc from 5.2. Classification Metrics, instead of raw accuracy.

Perfectly separable data makes the unregularized maximum likelihood estimate diverge. When a hyperplane separates the classes cleanly, the likelihood maximizes by pushing the weight norm toward infinity, since the probabilities saturate at 0 and 1. Unregularized gradient descent then keeps growing the weights the longer it runs. In practice, max_iterations and tolerance cut training off, and the classification stays correct. The weights, and so the probabilities, still become arbitrary in magnitude and poorly calibrated. A large learning_rate on large-magnitude separable inputs can also overflow the weight update completely. RustyML catches this with an in-loop Error::NonFinite guard, instead of silently returning NaN. Adding even a small L2 penalty bounds the optimum, keeps the probabilities meaningful, and removes the overflow risk. This is the main reason to keep a penalty on by default in production.

2.2.7. Polynomial features for non-linear boundaries

The decision boundary is linear in the feature space you give it. Non-linear problems therefore need a richer space. generate_polynomial_features(&x, degree) expands each row into all monomials up to degree. For 2 features at degree 2, this gives the 5 columns [x1, x2, x1^2, x1*x2, x2^2], with no constant column (the intercept already supplies that). Fit on the expansion, and predict on the same expansion.

use ndarray::array;
use rustyml::machine_learning::{LogisticRegression, generate_polynomial_features};

fn main() {
    // Inner ring = class 0, outer ring = class 1: not linearly separable in (x1, x2)
    let x = array![
        [ 1.0, 0.0], [0.0,  1.0], [-1.0, 0.0], [0.0, -1.0],
        [ 5.0, 0.0], [0.0,  5.0], [-5.0, 0.0], [0.0, -5.0],
    ];
    let y = array![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0];

    // 2 features, degree 2 -> [x1, x2, x1^2, x1*x2, x2^2]
    let x_poly = generate_polynomial_features(&x, 2);
    assert_eq!(x_poly.ncols(), 5);

    let mut model = LogisticRegression::new(true, 0.01, 3000, 1e-7).unwrap();
    model.fit(&x_poly, &y).unwrap();

    // The x1^2 + x2^2 term makes the two rings linearly separable
    let preds = model.predict(&x_poly).unwrap();
    println!("{:?}", preds);
}

The column count grows combinatorially with both feature count and degree. 3 features at degree 3 already yields 19 columns. Use this tool only for a handful of features and a low degree. Reach for a kernel method, such as SVC, when the expansion gets large.

2.2.8. A worked example with evaluation

This example fits on a small 2-feature set, then evaluates it with the classification metrics from Chapter 5. predict already returns hard {0.0, 1.0} labels as an Array1<f64>. That is exactly what ConfusionMatrix::new requires. It thresholds nothing itself, so it rejects a probability vector instead of binarizing it silently. roc_auc works the other way. It wants a boolean truth vector alongside the continuous predict_proba scores, because ranking is the whole point.

use ndarray::{array, Array1};
use rustyml::machine_learning::LogisticRegression;
use rustyml::metrics::{ConfusionMatrix, accuracy, roc_auc};

fn main() {
    let x_train = array![
        [-2.0, -1.5], [-1.5, -2.0], [-1.0, -0.5], [-2.5, -1.0], [-0.5, -1.0],
        [ 2.0,  1.5], [ 1.5,  2.0], [ 1.0,  0.5], [ 2.5,  1.0], [ 0.5,  1.0],
    ];
    let y_train = array![0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0];

    let mut model = LogisticRegression::new(true, 0.5, 500, 1e-7).unwrap();
    model.fit(&x_train, &y_train).unwrap();

    // Already hard {0.0, 1.0} labels, which is what ConfusionMatrix::new requires
    let preds = model.predict(&x_train).unwrap();

    let cm = ConfusionMatrix::new(&y_train, &preds);
    println!("{}", cm.summary());
    println!("accuracy:  {:.3}", accuracy(&y_train, &preds));

    // ROC AUC ranks the probabilities, so it needs the raw scores, not the 0/1 labels
    let truth: Array1<bool> = y_train.mapv(|v| v > 0.5);
    let scores = model.predict_proba(&x_train).unwrap();
    println!("ROC AUC:   {:.3}", roc_auc(&truth, &scores));
}

ConfusionMatrix::summary() prints the counts alongside accuracy, balanced accuracy, precision, recall, specificity, F1, and MCC in one table. It gives a fast way to check a binary classifier. Feeding predict_proba into roc_auc, instead of the thresholded labels, makes the AUC a threshold-independent measure of ranking quality. It answers how well the model orders positives above negatives, regardless of where you later set the cutoff. The output above depends on the data. Here is the shape to expect:

Confusion Matrix:
+-----------------+--------------------+--------------------+
| ...             | Predicted Positive | Predicted Negative |
...
Performance Metrics:
- Accuracy:          <0.0..1.0>
- ...
accuracy:  <0.0..1.0>
ROC AUC:   <0.0..1.0>

2.2.9. Persistence and reproducibility

A trained model serializes to a compact postcard binary through save_to_path and load_from_path. This binary carries the weights, hyperparameters, and iteration count. A round-trip is byte-exact, so a loaded model’s predictions match the original’s predictions exactly.

use ndarray::array;
use rustyml::machine_learning::LogisticRegression;

fn main() {
    let x = array![[-2.0, 1.0], [-1.0, -1.0], [1.0, 1.0], [2.0, -1.0]];
    let y = array![0.0, 0.0, 1.0, 1.0];

    let mut model = LogisticRegression::new(true, 0.3, 500, 1e-7).unwrap();
    model.fit(&x, &y).unwrap();

    let path = "logreg_model.bin";
    model.save_to_path(path).unwrap();
    let loaded = LogisticRegression::load_from_path(path).unwrap();

    assert_eq!(model.predict(&x).unwrap(), loaded.predict(&x).unwrap());

    std::fs::remove_file(path).unwrap();
    println!("round-trip OK");
}

Fitting has no randomness anywhere, so reproducibility comes free. There is no seed to set, unlike the sampling-based models in this chapter. 2 fits on identical data and hyperparameters yield identical weights. This is what makes the persistence round-trip exact. See 7.1. Reproducibility and Random Seeds for where seeds do matter. See 7.2. Model Persistence in Depth for the serialization format and its cross-version limits. If you build with the show_progress feature, fit also renders a live progress bar with the running loss. This bar is a convenient way to watch for the non-convergence and divergence behaviors described above.