# 2.2. Logistic Regression
`LogisticRegression` is RustyML's linear classifier for **binary** problems. It shares most of its machinery with [`LinearRegression`](./2.1._Linear_Regression.md): 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](#226-standardization-class-imbalance-and-separable-data).
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](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md) 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.
```rust
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`](../Chapter-01/1.6._Error_Handling.md), instead of letting it fail later at fit time.
| Parameter | Type | Default | Constraint |
| --- | --- | --- | --- |
| `fit_intercept` | `bool` | `true` | none |
| `learning_rate` | `f64` | `0.01` | strictly positive and finite |
| `max_iterations` | `usize` | `100` | at least 1 |
| `tolerance` | `f64` | `1e-4` | strictly positive and finite |
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:
| Getter | Returns |
| --- | --- |
| `get_fit_intercept()` | `bool` |
| `get_learning_rate()` | `f64` |
| `get_max_iterations()` | `usize` |
| `get_tolerance()` | `f64` |
| `get_regularization_type()` | `Option<RegularizationType>` |
| `get_actual_iterations()` | `Option<usize>` (iterations actually run, `None` before fit) |
| `get_weights()` | `Option<&Array1<f64>>` (`None` before fit) |
`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`](../Chapter-01/1.6._Error_Handling.md). This matches the label convention used by [`SVC`](./2.5._Support_Vector_Machines.md) and [`LinearSVC`](./2.5._Support_Vector_Machines.md), 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](../Chapter-04/4.3._Label_Encoding.md). Pick which class is positive, the `1`, on purpose. That choice defines what precision, recall, and the probability output mean.
```rust
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.
```rust
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`:
```rust
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)`](./2.1._Linear_Regression.md#216-regularization-l1-vs-l2) (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.1._Linear_Regression.md#216-regularization-l1-vs-l2).
## 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](../Chapter-04/4.2._Standardization_and_Normalization.md). 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](../Chapter-05/5.2._Classification_Metrics.md), 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**.
```rust
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](./2.5._Support_Vector_Machines.md), 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](../Chapter-05/5.0._Model_Evaluation.md). `predict` already returns hard `{0.0, 1.0}` labels as an `Array1<f64>`. That is exactly what [`ConfusionMatrix::new`](../Chapter-05/5.2._Classification_Metrics.md) 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.
```rust
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:
```text
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.
```rust
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](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md) for where seeds do matter. See [7.2. Model Persistence in Depth](../Chapter-07/7.2._Model_Persistence_in_Depth.md) 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.