Skip to main content

​
Overview

The MultilinearRegression class implements multivariate (multiple output) linear regression, fitting k response variables simultaneously. It provides significant performance advantages over fitting k separate models by factoring the coefficient matrix only once. Namespace: mlpp::regression Template parameter:
  • Scalar - Floating-point type (float, double, long double). Defaults to double.

​
Mathematical formulation

Solves the optimization problem: 
min_W  (1/2n) ||XW - Y||²_F  +  (λ/2) ||W||²_F
where:
  • X ∈ ℝⁿˣᵈ (feature matrix)
  • Y ∈ ℝⁿˣᵏ (response matrix)
  • W ∈ ℝᵈˣᵏ (coefficient matrix)
The Frobenius-norm objective decouples column-wise, making each response mathematically independent. However, solving them together is more efficient as the coefficient matrix A = XᵀX + nλI is factored exactly once.

​
Constructor

​
fit_intercept
bool
default:"true"
Whether to fit a bias vector b ∈ ℝᵏ (one intercept per response)
​
regularization
Scalar
default:"0"
L2 penalty λ ≥ 0 applied uniformly to all responses
​
method
SolveMethod
default:"SolveMethod::Auto"
Linear solver strategy:
  • SolveMethod::Auto - Chosen automatically (recommended)
  • SolveMethod::Cholesky - Normal equations via LDLT; O(nd² + d²k)
  • SolveMethod::SVD - Thin BDCSVD; O(nd·min(n,d))
  • SolveMethod::JacobiSVD - Full-pivoting JacobiSVD (slowest, most stable)

​
Methods

​
fit

void fit(const Matrix& X, const Matrix& Y)
Fit the model to multiple response variables.
​
X
Matrix
Feature matrix with shape (n_samples, n_features)
​
Y
Matrix
Response matrix with shape (n_samples, n_responses)

​
predict

Matrix predict(const Matrix& X) const
Predict response matrix for new samples.
​
X
Matrix
Feature matrix with shape (n_samples, n_features)
​
return
Matrix
Predicted matrix with shape (n_samples, n_responses)

​
score

Scalar score(const Matrix& X, const Matrix& Y) const
Compute mean R² across all response variables.
​
X
Matrix
Feature matrix
​
Y
Matrix
True response matrix
​
return
Scalar
Mean R² score in the range (-∞, 1]. Higher values indicate better fit

​
score_per_response

Vector score_per_response(const Matrix& X, const Matrix& Y) const
Compute per-response R² scores.
​
X
Matrix
Feature matrix
​
Y
Matrix
True response matrix
​
return
Vector
Vector of length n_responses where r2(k) = 1 - SS_res_k / SS_tot_k

​
residuals

Matrix residuals(const Matrix& X, const Matrix& Y) const
Compute residual matrix E = Y - XW - 1bᵀ.
​
X
Matrix
Feature matrix
​
Y
Matrix
True response matrix
​
return
Matrix
Residual matrix with shape (n_samples, n_responses). Requires fit()

​
gradient

Matrix gradient(const Matrix& X, const Matrix& Y) const
Compute gradient of the regularized MSE with respect to W: 
∇L = (1/n) Xᵀ(XW − Y) + λW
​
X
Matrix
Feature matrix
​
Y
Matrix
Response matrix
​
return
Matrix
Gradient matrix with shape (n_features, n_responses). Requires fit()

​
coefficients

const Matrix& coefficients() const
​
return
Matrix
Coefficient matrix W ∈ ℝᵈˣᵏ in original (unscaled) feature space

​
intercepts

const Vector& intercepts() const
​
return
Vector
Bias vector b ∈ ℝᵏ with length n_responses. Zero vector if fit_intercept is false

​
is_fitted

bool is_fitted() const noexcept
​
return
bool
True after a successful call to fit()

​
condition_number

Scalar condition_number() const
​
return
Scalar
Condition number σ_max/σ_min of the design matrix. Only available after an SVD solve

​
Usage example

#include <mlpp/regression/multilinear_regression.hpp>
#include <Eigen/Dense>

using namespace mlpp::regression;
using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;

// Create training data
Matrix X(100, 5);   // 100 samples, 5 features
Matrix Y(100, 3);   // 100 samples, 3 response variables
// ... fill X and Y with data ...

// Create and fit model
MultilinearRegression<double> model(
    true,    // fit_intercept
    0.1,     // regularization
    MultilinearRegression<double>::SolveMethod::Auto
);

model.fit(X, Y);

// Make predictions for all responses
Matrix X_test(20, 5);
Matrix predictions = model.predict(X_test);

// Evaluate overall performance
double mean_r2 = model.score(X, Y);
std::cout << "Mean R² score: " << mean_r2 << std::endl;

// Get per-response scores
Vector r2_scores = model.score_per_response(X, Y);
for (int i = 0; i < r2_scores.size(); ++i) {
    std::cout << "Response " << i << " R²: " 
              << r2_scores(i) << std::endl;
}

// Access learned parameters
const Matrix& W = model.coefficients();      // shape: (5, 3)
const Vector& b = model.intercepts();        // shape: (3,)
double cond = model.condition_number();

// Compute residuals for all responses
Matrix E = model.residuals(X, Y);

​
Performance advantages

Compared to fitting k separate LinearRegression models:
  • Cholesky solver: Factors XᵀX + nλI once, then solves k right-hand sides → O(nd² + d²k)
  • SVD solver: Computes SVD of X once, then filters k response columns → O(nd·min(n,d) + nk·min(n,d))
This is significantly faster than k × O(nd² + d³) for separate models.

​
Type aliases

using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
using Index = Eigen::Index;

Build docs developers (and LLMs) love

Get started for freeTalk to us