Constraining asymptotic behavior

[MilesCranmer]

Moved from #276:

@dominik-rehse:

Dear @MilesCranmer,

Thank you for implementing this awesome feature, which I saw only now!

Would this - at least in principle - allow for symbolic constraints on the functional form (e.g., $\lim_{x\to\infty} f(x)=0$
)? If so, could you provide any hints on how to implement something like that? I am not really familiar with Julia and could not immediately find more information on symbolic operations in the docs of SymbolicRegression.jl.

Thanks,
Dominik

[MilesCranmer]

Great question. Actually the simplest and perhaps best way to constraint asymptotic behavior is to add datapoints in the direction of the asymptotes. While for other ML techniques, it wouldn't work, for SR this actually is very effective.

For example, for your constraint, I would add:

large_number = 100_000
x = np.append(x, large_number)
y = np.append(y, 0)

large_number_2 = 1_000_000
x = np.append(x, large_number2)
y = np.append(y, 0)

assuming that x is a 1D array (you would have to reshape it when passing to PySR).

You can also add some weights to enforce the constraint more.

[MilesCranmer]

One other option is to try this package: https://github.com/JuliaMath/Richardson.jl which can numerically estimate limits. This is perhaps the "proper" way to do it, but I'm not sure how effective it would be.

You can see how to install Julia packages into PySR in this example: https://astroautomata.com/PySR/examples/#7-julia-packages-and-types (or one could also use the Julia backend directly).

Here's an example of using Richardson:

using Richardson

function my_custom_objective(tree, dataset::Dataset{T,L}, options)::L where {T,L}
    prediction, flag = eval_tree_array(tree, dataset.X, options)
    if !flag
        return L(Inf)
    end
    limit_val, precision = extrapolate(1.0, x0=Inf) do x
        out, completed = eval_tree_array(tree, [x]', options)
        !completed && (out .= NaN)
        return out[1]
    end
    !isfinite(limit_val) && return L(Inf)
    prediction_loss = sum((prediction .- dataset.y) .^ 2) / dataset.n
    asymptote_loss = (limit_val - 0.0) ^ 2
    loss = prediction_loss + 100 * asymptote_loss
    return loss
end

[dominik-rehse]

Follow-up question:

Could I simply replace the second eval_tree_array here with eval_diff_tree_array in order to get $lim_{x\to\inf}f'(x)=0$ instead of $lim_{x\to\inf}f(x)=0$?

Apologies if this is stupid question but I am a Julia beginner and struggle a little with the docs.

Many thanks in advance!

[MilesCranmer]

For eval_diff_tree_array you also need to pass the direction (in this case I guess just 1). It also outputs three things rather than two: https://symbolicml.org/DynamicExpressions.jl/dev/eval/#Derivatives so use it like out, derivative, complete = …

[dominik-rehse]

Thanks!

[dominik-rehse]

Dear @MilesCranmer,

Thank you for the prompt and super-helpful replies!

-Dominik