# What is Mathematica's ML Linear Regression Method?

When I run p = Predict[trainingSet, Method -> "LinearRegression"];, how is that different from just running ordinary least squares?

It's quite slow on a training set of 10000 examples, each with six features, so I presume something a bit more sophisticated than OLS must be going on? It's also then quite slow to apply the predictor function: predictions = p /@ trainingSet[[;;, 1]];.

As noted in the documentation for Predict[], it claims to do $L^2$ regularization, so its results for Method -> "LinearRegression" are not equivalent to what you'll get for e.g. LinearModelFit[].

To get results equivalent to LinearModelFit[], you need to set "L2Regularization" -> 0.

train = {1 -> 1.3, 2 -> 2.4, 3 -> 4.4, 4 -> 5.1, 6 -> 7.3};

p1 = Predict[train, Method -> {"LinearRegression", "L2Regularization" -> 0}];

p2 = LinearModelFit[List @@@ train, {1, x}, x];

{p1[5/2], p2[5/2]}
{3.25338, 3.25338}


However, there is something off in their description for $L^2$ (Tikhonov) regularization. If it is true that Tikhonov is being done behind the scenes, then a manual reimplementation ought to give the same results:

λ = 1; (* regularization parameter *)

p1b = Predict[trainingset,
Method -> {"LinearRegression", "L1Regularization" -> 0,
"L2Regularization" -> λ,
"OptimizationMethod" -> "NormalEquation"}];

{dm, rv} = Through[{DesignMatrix[#, {1, x}, x] &, #[[All, -1]] &}[List @@@ train]];

(* minimize the loss function *)
{ar, br} = FindArgMin[#.# &[rv - dm.{a, b}] + λ {a, b}.{a, b}, {a, b}]
{0.256849, 1.18493}

(* directly solve the modified normal equations *)
{ar, br} = LinearSolve[Transpose[dm].dm + λ IdentityMatrix[2], Transpose[dm].rv]
{0.256849, 1.18493}


but

{p1b[5/2], ar + br 5/2}
{3.39448, 3.21918}


so there is something additional/different being done behind the scenes.