PWL on LOESS smoothers?

[JockLawrie]

Hi there,

How much work would it be to enable piecewise linear approximations of functions that aren't supplied as expressions?

In particular, I have a 2-D scatter plot which I convert into a smooth continuous function using Loess.jl. Since the resulting loess object is made up of cubic splines, at any point the derivatives are easily computable.

When I tried putting this through pwl, I get an error because the differentiate function failed. What needs to happen to get this working? Implement a differentiate function for the loess object?

Cheers, Jock

[Codsilla]

Hi Jock,

The package works with both expressions and native Julia functions.

As of now, the package assumes differentiability (it would be possible to have a derivative-free version of Linearize but it would be much slower).

Specifically, your function needs to be differentiable by the package ForwardDiff.jl (or Calculus.jl if you're using an expression instead of a native Julia function). If you make ForwardDiff.gradient and Calculus.derivative work on your function, then LinA will work.

PS: Consider giving a Minimal reproducible example.

[JockLawrie]

Thanks for responding.

Below is a MWE of the smooth curve constructed from a scatterplot using the LOESS package. Note that value of the curve (f(x)) and its derivatives (df(x) and ddf(x)) are easily found for any point x. I'd like to be able to run the function f through Linearize.

Since the derivatives are already known exactly, is ForwardDiff.gradient unnecessary?

using Loess

# Data
xs = 6 .* rand(1000);
ys = sin.(xs) .+ 0.5 .* randn(1000);

# Smooth fitted curve
model = loess(xs, ys; degree=3, span=0.9)

# Univariate function using the Loess curve
f(x) = Loess.predict(model, x)

# Derivative functions df and ddf
first_deriviative(b, x)  = b[2] + 2.0*b[3]*x + 3.0*b[4]*x*x
second_deriviative(b, x) = 2.0*b[3] + 6.0*b[4]*x

function deriviative(model, x, dfunc)
    adjacent_vertices = Loess.traverse(model.kdtree, [x])
    knot1, knot2 = adjacent_vertices[1][1], adjacent_vertices[2][1]
    i1, i2 = model.verts[[knot1]], model.verts[[knot2]]
    b1, b2 = view(model.bs, i1, :), view(model.bs, i2, :)
    if x == knot1
        return dfunc(b1, x)
    elseif x == knot2
        return dfunc(b2, x)
    else
        p = (x - knot1) / (knot2 - knot1)
        return p * dfunc(b1, x) + (1.0 - p) * dfunc(b2, x)
    end
end

df(x)  = deriviative(model, x, first_deriviative)   # First  derivative
ddf(x) = deriviative(model, x, second_deriviative)  # Second derivative

# Usage
f(4.0)
df(4.0)
ddf(4.0)