2.1. Linear Regression

LinearRegression fits a linear map y_hat = X * w + b by least squares. RustyML’s default path is the exact closed form, like scikit-learn’s LinearRegression. The coefficients agree with Python’s results to about 1e-15. Batch gradient descent is available on demand, for cases the closed form cannot serve.

This split shapes the whole page. The default estimator has no learning rate, no iteration budget, and no convergence tolerance to set. The 3 settings that do exist belong to the iterative solver, and they travel with it. The source for this page is src/machine_learning/linear_model/linear_regression.rs and its integration tests.

2.1.1. How the model is trained

Two strategies minimize the same objective. LeastSquaresSolver picks between them.

The closed form (LeastSquaresSolver::Normal, the default) solves the ridge system in one shot. With no penalty, this is plain OLS. It runs an SVD least-squares solve on the augmented design [Xc; sqrt(lambda) * I]. This returns the minimum-norm solution even when X^T * X is singular, for example with perfectly collinear columns or more features than samples. When fit_intercept is set, the solver mean-centers the features and target first, so the intercept stays out of the penalty. It then recovers the intercept as mean(y) - mean(X) * w. The closed form has nothing to iterate. get_actual_iterations() returns Some(0) after a normal-solver fit. That 0 signals that the fit took no gradient steps.

Gradient descent (LeastSquaresSolver::GradientDescent) initializes the weight vector w to zeros and the intercept b to zero, then repeats a full-batch update. Each iteration computes the prediction and residual over the whole training set, the scalar cost, the gradients, and the parameter step.

The cost is the mean squared error, halved, plus an optional penalty. Let n be the sample count and e = X * w + b - y be the residual vector. One iteration computes cost = dot(e, e) / (2 * n) + penalty. The penalty term is 0 with no regularization. It is alpha * sum(|w_j|) for L1. It is (alpha / 2) * dot(w, w) for L2. The prediction adds the intercept b only when fit_intercept is set. Regularization never penalizes the intercept, because it shrinks slopes, not the bias.

The gradients follow directly. The weight gradient is grad_w = (X^T * e) / n, plus alpha * w for L2. The intercept gradient is grad_b = sum(e) / n, or 0 when fit_intercept is false. The update step is plain gradient descent with a fixed learning_rate. It sets w to w - learning_rate * grad_w and b to b - learning_rate * grad_b. L1 does not appear in the gradient. The solver applies it after the step, through a proximal operator (see 2.1.6). Gradient descent here has no momentum, no adaptive rate, and no line search. The step size you set on the solver is the step size used on every iteration. This is why feature scaling matters so much (see 2.1.5).

3 numerical guards run inside the loop. If the cost, any gradient component, or any updated parameter becomes NaN or infinite, fit aborts immediately with Error::NonFinite rather than returning a garbage model. A divergent learning rate trips this guard within a handful of iterations, instead of silently producing inf coefficients. The closed form has its own version of the same check on the solution it computes.

The solver computes the residual sum of squares and the intercept gradient’s sum with deterministic blocked folds. The matrix-vector products, X * w and X^T * e, run in parallel above an internal size gate. Neither solver carries any randomness. 2 runs on the same machine with the same data produce identical coefficients. The determinism test asserts bit-identical predictions across 2 independently constructed models. For the parallel-reduction machinery behind this, see 7.3. Performance Tuning and Parallelism.

2.1.2. Constructing a model

The constructor takes one argument and is infallible:

pub fn new(fit_intercept: bool) -> Self

There is nothing else to validate here, which is why new returns Self rather than Result. The iteration settings that used to sit on the estimator now live inside the solver variant that uses them. with_solver validates them instead.

LinearRegression::default() is the same as LinearRegression::new(true). It sets fit_intercept = true, selects LeastSquaresSolver::Normal, and applies no regularization. The two constructors always build the same algorithm now. default() is the analogue of Python’s LinearRegression(), and it is exact.

Two builder methods refine a constructed model. Both consume and return self, so they chain. Both return Result, because both validate what they receive:

pub fn with_solver(self, solver: LeastSquaresSolver) -> Result<Self, Error>
pub fn with_regularization(self, regularization: RegularizationType) -> Result<Self, Error>

with_solver checks the selected variant’s payload. A non-positive or non-finite learning_rate or tol, or a max_iter of 0, yields Error::InvalidParameter with the offending field’s name. LeastSquaresSolver::Normal has no payload, so selecting it can never fail. with_regularization validates the penalty coefficient alpha the same way. It rejects a negative or non-finite alpha on the spot. The regularization enum is re-exported at rustyml::machine_learning::RegularizationType:

pub enum RegularizationType {
    L1(f64), // Lasso: penalty alpha * sum(|w_j|)
    L2(f64), // Ridge: penalty (alpha / 2) * dot(w, w)
}

Every hyperparameter and fitted quantity is readable through a getter. The fitted ones return Option, None until fit has run:

There are no get_learning_rate, get_max_iterations, or get_tolerance accessors. Those numbers are fields of the solver variant. Read them back by matching on get_solver():

if let LeastSquaresSolver::GradientDescent { learning_rate, max_iter, tol } = model.get_solver() {
    println!("learning_rate = {learning_rate}, cap = {max_iter}, tol = {tol}");
}

There is deliberately no cost-history accessor either. The model stores the final iteration count, not the per-iteration cost curve. If you want to watch the cost descend live, build with the show_progress feature. Then fit renders a progress bar with the running cost and convergence counter. Otherwise get_actual_iterations() is the only window into what the optimizer did (see 2.1.7).

2.1.3. Fitting and predicting

The four data-facing methods are generic over ndarray’s storage type S: Data<Elem = f64>. This lets them accept owned arrays, views, and slices alike, without a copy:

pub fn fit<S1, S2>(&mut self, x: &ArrayBase<S1, Ix2>, y: &ArrayBase<S2, Ix1>) -> Result<&mut Self, Error>;
pub fn predict<S>(&self, x: &ArrayBase<S, Ix2>) -> Result<Array1<f64>, Error>;
pub fn fit_predict<S1, S2>(&mut self, x: &ArrayBase<S1, Ix2>, y: &ArrayBase<S2, Ix1>) -> Result<Array1<f64>, Error>;
pub fn score<S1, S2>(&self, x: &ArrayBase<S1, Ix2>, y: &ArrayBase<S2, Ix1>) -> Result<f64, Error>;

x is a feature matrix with one sample per row and one feature per column. y is the target vector. fit returns &mut Self for chaining. predict returns the prediction vector. fit_predict runs fit then predict on the same data. score returns the coefficient of determination, R^2 (1 - SS_res / SS_tot). A score of 1.0 is perfect. A score of 0.0 matches a model that always predicts the mean. A negative score is worse than the mean. The generic Fit and Predict traits expose the same operations for code that works across estimator types. Their trait methods just forward to these inherent methods.

fit rejects an empty x with Error::EmptyInput. It rejects a y whose length differs from x’s row count with Error::DimensionMismatch. It rejects any NaN or infinite entry in x with Error::NonFinite. predict and score also return Error::NotFitted when you call them before training. predict checks that the column count matches the training data, again with Error::DimensionMismatch. All of these come from the crate’s unified Error type. See 1.6. Error Handling.

A complete run on a small multivariate problem, y = 2 * x1 + 3 * x2 + 1. Nothing is configured beyond the intercept, so this is exact OLS:

use ndarray::array;
use rustyml::machine_learning::LinearRegression;

fn main() {
    // Six samples spanning the feature space.
    let x = array![
        [1.0, 1.0],
        [2.0, 1.0],
        [1.0, 2.0],
        [3.0, 2.0],
        [2.0, 3.0],
        [4.0, 1.0],
    ];
    let y = array![6.0, 8.0, 9.0, 13.0, 14.0, 12.0];

    // fit_intercept = true, closed-form solver, no regularization.
    let mut model = LinearRegression::new(true);
    model.fit(&x, &y).unwrap();

    // Coefficients land on [2.0, 3.0] and the intercept on 1.0, to machine precision.
    println!("coefficients = {:?}", model.get_coefficients().unwrap());
    println!("intercept    = {}", model.get_intercept().unwrap());
    println!("iterations   = {}", model.get_actual_iterations().unwrap()); // 0

    // Predict on unseen rows: (1,1) -> 6.0, (2,3) -> 14.0.
    let preds = model.predict(&array![[1.0, 1.0], [2.0, 3.0]]).unwrap();
    println!("predictions  = {:?}", preds);

    // R^2 on the training data, 1.0 for exactly-linear data.
    println!("R^2          = {}", model.score(&x, &y).unwrap());
}

The closed form lands on the least-squares optimum rather than approaching it. The recovered slopes match to machine precision. Gradient descent, by contrast, only converges toward that solution. The tests assert its slopes stay within 3e-3 of the true values.

2.1.4. Choosing a solver: the normal equation vs. gradient descent

LeastSquaresSolver is a payload-carrying enum. Each variant owns exactly the settings it uses. You cannot pass a learning rate to the closed form, and you cannot leave one unset for the iterative path.

pub enum LeastSquaresSolver {
    Normal, // the default
    GradientDescent { learning_rate: f64, max_iter: usize, tol: f64 },
}

LeastSquaresSolver::Normal is the right choice for small to medium dense problems. This is why it is the default. It is exact, needs no scaling, and has no hyperparameters to tune.

Prefer gradient descent in 3 cases. First, when the dataset is large enough that the closed-form factorization, roughly O(n * p^2), costs too much. Second, when you want L1, since the closed form has no L1 solution. Third, when you match a scikit-learn or Keras workflow that also uses iterative optimization.

The Normal solver supports only no regularization or L2. Pairing it with an L1 penalty makes fit return Error::InvalidInput, because Lasso has no closed form. Use gradient descent for L1.

Selecting gradient descent and configuring it are one expression, and the result is a Result, because the payload is validated:

use ndarray::array;
use rustyml::machine_learning::LinearRegression;
use rustyml::machine_learning::linear_model::LeastSquaresSolver;

fn main() {
    // y = 3*x0 - 2*x1 + 5, exactly.
    let x = array![
        [1.0, 1.0],
        [2.0, 0.0],
        [0.0, 3.0],
        [4.0, 2.0],
        [3.0, 1.0],
        [1.0, 4.0],
    ];
    let y = array![6.0, 11.0, -1.0, 13.0, 12.0, 0.0];

    // The exact solution, in one step.
    let mut exact = LinearRegression::new(true);
    exact.fit(&x, &y).unwrap();

    // The iterative solution: the settings travel with the variant that uses them.
    let mut iterative = LinearRegression::new(true)
        .with_solver(LeastSquaresSolver::GradientDescent {
            learning_rate: 0.01,
            max_iter: 10_000,
            tol: 1e-9,
        })
        .unwrap();
    iterative.fit(&x, &y).unwrap();

    println!("exact     = {:?}", exact.get_coefficients().unwrap());     // ~ [3.0, -2.0]
    println!("iterative = {:?}", iterative.get_coefficients().unwrap()); // approaches the same
    println!("iterations: exact = {:?}, iterative = {:?}",
        exact.get_actual_iterations(), iterative.get_actual_iterations());
}

The two solvers agree numerically when they solve the same objective. The gradient-descent cost divides the data term by 2 * n but scales the L2 penalty by only alpha / 2. This makes the equivalent ridge penalty on the raw sum of squares equal to lambda = n * alpha. That value is exactly what the normal solver uses. The integration tests confirm that the closed-form L2 solution matches what gradient descent converges to, for the same alpha.

2.1.5. Standardize your features first

This section is about gradient descent. The closed form sidesteps this whole concern, which is one more reason it is the default.

Gradient descent with a single global learning rate depends on how well your feature scales match. Suppose one column ranges over the thousands and another sits near unit magnitude. The loss surface then becomes a steep, narrow valley. A learning rate small enough to avoid overshooting on the steep axis crawls along the shallow one, so convergence takes far more iterations. A learning rate large enough to make progress on the shallow axis diverges on the steep one and trips the finiteness guard. Standardizing each column to zero mean and unit variance makes the curvature roughly isotropic. Then one learning rate serves every direction, and you can safely use a larger one.

StandardScaler does exactly this. It remembers the training mean and standard deviation, so it transforms the test set with the training statistics, not the test set’s own. See 4.2. Standardization and Normalization for the full train and test workflow, and for the stateless standardize free function. The coefficients you read back are then in standardized units, not raw feature units.

use ndarray::array;
use rustyml::machine_learning::LinearRegression;
use rustyml::machine_learning::linear_model::LeastSquaresSolver;
use rustyml::utils::StandardScaler;

fn main() {
    // Column 0 lives in the thousands. Column 1 is near unit scale.
    let x_raw = array![
        [1000.0, 1.0],
        [2000.0, 3.0],
        [3000.0, 2.0],
        [4000.0, 5.0],
        [5000.0, 4.0],
    ];
    let y = array![10.0, 23.0, 32.0, 45.0, 54.0];

    // Zero mean, unit variance per column. The scaler keeps the statistics for later batches.
    let mut scaler = StandardScaler::new();
    let x = scaler.fit_transform(&x_raw).unwrap();

    // Isotropic features tolerate a much larger step than the raw data would.
    let mut model = LinearRegression::new(true)
        .with_solver(LeastSquaresSolver::GradientDescent {
            learning_rate: 0.1,
            max_iter: 5_000,
            tol: 1e-9,
        })
        .unwrap();
    model.fit(&x, &y).unwrap();

    println!("iterations   = {}", model.get_actual_iterations().unwrap());
    println!("R^2          = {}", model.score(&x, &y).unwrap());
    // These slopes are in standardized units.
    println!("coefficients = {:?}", model.get_coefficients().unwrap());
}

As a rule of thumb for the step size, start around 0.1 on standardized data. Back off by a factor of 10 if fit returns Error::NonFinite or the iteration count pins to max_iter. On raw, small-magnitude integer features, a rate near 0.01 is usually safe.

2.1.6. Regularization: L1 vs L2

L2 (RegularizationType::L2(alpha), ridge) adds (alpha / 2) * dot(w, w) to the cost and alpha * w to the weight gradient. It shrinks every coefficient smoothly toward zero, without forcing any coefficient to exactly zero. This suits correlated features, where you want a stable, low-variance fit. Ridge spreads weight across collinear columns instead of letting one grow too large. Both solvers support L2.

L1 (RegularizationType::L1(alpha), lasso) adds alpha * sum(|w_j|) to the cost. The solver applies it with a proximal step. After each gradient step, it soft-thresholds every weight by learning_rate * alpha. Any weight the data cannot justify lands on exactly 0.0 and stays there. This method is called ISTA, and it is what makes L1 a feature selector. An older sub-gradient form added alpha * sign(w) to the gradient instead. That form only approaches zero over time, so it never produced true sparsity. Either way, the intercept stays unpenalized. L1 needs the gradient-descent solver, because no closed form exists for it.

Setting alpha correctly matters most when you port a model from Python. Every estimator that takes RegularizationType minimizes a mean data term plus an undivided penalty:

L1:  (1 / n) * sum(loss) + alpha * ||w||_1
L2:  (1 / n) * sum(loss) + alpha * 0.5 * ||w||^2

This matches scikit-learn’s SGDRegressor and SGDClassifier objective exactly, so alpha transfers 1:1 from either one. The closed-form estimators use their own conventions, and those conventions differ from each other:

Here n is the number of training samples. The Ridge row states the lambda = n * alpha relationship from 2.1.4, read in the other direction.

This snippet uses the exact solver, so the shrinkage is unambiguous. Ridge produces a strictly smaller coefficient norm than OLS on nearly-collinear data:

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

fn main() {
    // Two nearly-collinear features make plain OLS coefficients large and unstable.
    let x = array![
        [1.0, 0.9],
        [2.0, 2.1],
        [3.0, 2.9],
        [4.0, 4.2],
        [5.0, 5.1],
    ];
    let y = array![1.0, 2.0, 3.0, 4.0, 5.0];

    // Both use the default closed-form solver.
    let mut ols = LinearRegression::new(true);
    ols.fit(&x, &y).unwrap();

    let mut ridge = LinearRegression::new(true)
        .with_regularization(RegularizationType::L2(1.0))
        .unwrap();
    ridge.fit(&x, &y).unwrap();

    let sq_norm =
        |m: &LinearRegression| m.get_coefficients().unwrap().iter().map(|c| c * c).sum::<f64>();
    println!("OLS   ||w||^2 = {}", sq_norm(&ols));
    println!("ridge ||w||^2 = {}", sq_norm(&ridge)); // strictly smaller
}

For L1, select gradient descent explicitly. The Normal solver rejects L1 with Error::InvalidInput. The crate’s own tests use 200 features, where column 0 carries the signal and the other 199 columns are pure noise. A small L1 penalty keeps the informative coefficient dominant and drives every noise coefficient to literal 0.0. Count the zeros with iter().filter(|c| **c == 0.0).count(). The tests assert exactly this behavior.

2.1.7. Convergence and diagnosing non-convergence

This section applies to LeastSquaresSolver::GradientDescent. The closed form has no notion of convergence. It either produces a solution or returns an error.

Convergence is not declared on a single small step. After each iteration, the optimizer compares the absolute change in cost against tol. A change below tol increments a counter. Training stops only after 3 consecutive iterations below tol. Any iteration whose cost change exceeds the threshold resets the counter to zero. This rule guards against a premature stop on a temporary plateau. The loop also stops unconditionally at max_iter.

After fit, compare get_actual_iterations() against the max_iter you configured. Read max_iter back by matching on get_solver(). If the actual count is strictly below the cap, the 3-in-a-row rule fired, and the model converged. If it equals the cap, the optimizer ran out of budget, and the fit may be underconverged.

The usual remedies, in order:

1. Standardize the features (see 2.1.5).
2. Raise max_iter.
3. Raise the learning rate for faster descent, but watch for divergence.
4. Loosen tol, if you do not need the last digit of precision.
5. Drop the with_solver call, and use the exact closed-form answer instead.

Divergence is the other failure mode. A learning rate that is too large sends the cost, and then the parameters, to infinity. The in-loop finiteness guard converts that into Error::NonFinite within a few iterations, instead of returning coefficients filled with inf. A NonFinite error from fit, on clean and finite data, almost always means the learning rate is too high, or the features need scaling. The model keeps no per-iteration cost history. To watch the descent directly, build with the show_progress feature. Then fit renders the running cost and the k/3 convergence counter on each iteration.

2.1.8. Inspecting the fitted parameters

The getters let you read the fitted model without running inference again. get_coefficients() returns Option<&Array1<f64>>, a borrow, because the model owns its weights. get_intercept() and get_actual_iterations() return owned Option values instead. The solver’s settings come back inside the enum. Read them with a match or an if let.

use ndarray::array;
use rustyml::machine_learning::LinearRegression;
use rustyml::machine_learning::linear_model::LeastSquaresSolver;

fn main() {
    let x = array![[1.0], [2.0], [3.0], [4.0], [5.0]];
    let y = array![3.0, 5.0, 7.0, 9.0, 11.0]; // y = 2x + 1

    let mut model = LinearRegression::new(true)
        .with_solver(LeastSquaresSolver::GradientDescent {
            learning_rate: 0.01,
            max_iter: 10_000,
            tol: 1e-10,
        })
        .unwrap();
    model.fit(&x, &y).unwrap();

    // Weights come back by reference, so iterate without moving them out of the model.
    let coefs = model.get_coefficients().unwrap();
    for (j, w) in coefs.iter().enumerate() {
        println!("w[{j}] = {w:.6}");
    }

    // Intercept and iteration count are owned values.
    println!("intercept       = {:.6}", model.get_intercept().unwrap());
    let ran = model.get_actual_iterations().unwrap();
    println!("iterations run  = {ran}");

    // The iteration settings live in the solver variant.
    match model.get_solver() {
        LeastSquaresSolver::Normal => println!("closed form, nothing to tune"),
        LeastSquaresSolver::GradientDescent { learning_rate, max_iter, tol } => {
            println!("learning_rate   = {learning_rate}");
            println!("tolerance       = {tol}");
            println!("max_iter        = {max_iter}");
            println!("converged early = {}", ran < max_iter);
        }
    }

    // Regularization is readable whether or not the model has been fitted.
    println!("regularization  = {:?}", model.get_regularization_type());
}

When fit_intercept is false, get_intercept() returns Some(0.0) after fitting, by contract. The intercept is fixed at zero, not merely small. Reading it back as exactly 0.0 confirms that no bias term was learned.

2.1.9. Saving and loading a model

A fitted model serializes to a compact postcard binary blob. The blob includes the coefficients, the intercept, the hyperparameters, and the training metadata. save_to_path(&self, path: &str) writes it. load_from_path(path: &str) -> Result<Self, Error> reads it back. Both methods report I/O and serialization failures as Error::Io. A round trip reproduces predictions exactly, because the raw f64 coefficients stay bit for bit the same.

use ndarray::array;
use rustyml::machine_learning::LinearRegression;

fn main() {
    let x = array![[1.0], [2.0], [3.0], [4.0], [5.0]];
    let y = array![3.0, 5.0, 7.0, 9.0, 11.0];

    let mut model = LinearRegression::new(true);
    model.fit(&x, &y).unwrap();
    let before = model.predict(&array![[6.0]]).unwrap();

    // Persist to disk, then reload into a fresh instance.
    model.save_to_path("lr_model.bin").unwrap();
    let loaded = LinearRegression::load_from_path("lr_model.bin").unwrap();
    let after = loaded.predict(&array![[6.0]]).unwrap();

    println!("before = {:?}, after = {:?}", before, after); // identical

    std::fs::remove_file("lr_model.bin").unwrap();
}

The serialized layout changed when the iteration settings moved into the solver variant. What used to be 3 loose fields is now 1 payload. A blob written by an older version will not load. Re-fit and re-save any persisted models.

The file extension carries no meaning. The format is always postcard binary, whether you name the file .bin, .dat, or anything else. LinearRegression also derives Clone and Debug. An in-process copy is a plain model.clone(), and {:?} prints the full struct. For full details on cross-model persistence, including neural-network weights, see 7.2. Model Persistence in Depth.

For the classification counterpart built on the same gradient-descent machinery, continue to 2.2. Logistic Regression. For more regression scores you can compute on a fitted model’s predictions, see 5.1. Regression Metrics.