7
$\begingroup$

I would like to implement a LOESS regression in Mathematica. My implementation depends on LinearModelFit, which is too slow to use on large datasets. (Can't use Fit because we need weights). Here is my current version:

loess[data_, order_, smoothing_] :=
 Module[{d, weights, filtered, fData, fWeights, fit, x, weight, funcs},
  weight[z_] := Piecewise[{{(1 - Abs[z]^3)^3, Abs[z] < 1}}];
  funcs = Table[x^i, {i, 0, order}];
  Table[
   d = (data[[;; , 1]] - point[[1]])/smoothing;
   weights = weight /@ d;
   filtered = 
    Cases[Transpose[Join[Transpose[data], {weights}]], 
     a_ /; a[[3]] != 0];
   fData = filtered[[;; , {1, 2}]];
   fWeights = filtered[[;; , 3]];
   fit = LinearModelFit[fData, funcs, x, Weights -> fWeights]
   {point[[1]], fit[point[[1]]]}
   , {point, data}]
  ]

With an example dataset:

data = Table[{x, 2*x^3 - 1.5 x^2 + RandomReal[{-0.05, 0.05}]}, {x, 0, 1, 0.001}];

running loess[data, 2, 0.1] takes about 7 seconds on my machine.

Is there a way to speed up this function? The accepted answer will also allow access to the slope of the LOESS regression; for example by returning {point[[1]], fit'[point[[1]]]} in the above implementation.

(Note that my definition of the smoothing parameter is non-standard.)

$\endgroup$
  • $\begingroup$ I don't feel it's satisfying, but: Replacing LinearModelFit with NonlinearModelFit reduces the computation time on my machine from 13s to 9s. $\endgroup$ – Oscillon Sep 2 '16 at 20:16
  • $\begingroup$ Reading your motivation for implementing LOESS regression, are you interested in any other algorithms to be applied that would be fast enough? $\endgroup$ – Anton Antonov Sep 2 '16 at 21:29
  • $\begingroup$ I would certainly be interested in a completely different implementation of LOESS, if it were faster. Is that what you had in mind? If not, still curious what you are thinking of. $\endgroup$ – mdeceglie Sep 2 '16 at 21:44
  • $\begingroup$ If Fit works faster, you can simply use data * weights instead of data for weighted fit. $\endgroup$ – Alexey Popkov Sep 3 '16 at 1:13
  • 1
    $\begingroup$ You went to the trouble of creating filtered, fData and fWeights but then used data and weights in the fitting process. Typo? I get the same answer with both approaches with a factor of two speedup using the fData and fWeights. This would become more significant the larger the data set. $\endgroup$ – Jack LaVigne Sep 5 '16 at 2:51
5
$\begingroup$

As per discussion in the comments, here is an answer using Quantile Regression that produces a fitted curve for less than 0.5 seconds on my computer. (The function loess in the question took 8 seconds.)

Generate the data (as in the question):

data = Table[{x, 2*x^3 - 1.5 x^2 + RandomReal[{-0.05, 0.05}]}, {x, 0, 1, 0.001}];

Load the QuantileRegression.m package:

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]

Do quantile regression with 12 B-spline basis nodes for the quantile 0.5:

AbsoluteTiming[
 qfunc = QuantileRegression[data, 12, {0.5}][[1]];
]

(* {0.343834, Null} *)

Plot the data and the regression quantile curve:

ListPlot[{data, {#, qfunc[#]} & /@ data[[All, 1]]}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.