Checking for monotonicity through derivatives

[gm89uk]

Hi all,

I'm struggling to get the automatic derivative functionality to work:

I set enable_autodiff=True in the python code. When trying to asses the gradient with respect to p (X[3,:]) to (y) (negative correlated, so negative grad expected) I'm getting a "JuliaError: TaskFailedException".

EDIT: a more detailed error message: MethodError: no method matching eval_grad_tree_array(::Node{Float32}, ::Matrix{Float32}, ::Int64, ::Options{SymbolicRegression.CoreModule.OptionsStructModule.ComplexityMapping{Int64, Int64}, DynamicExpressions.OperatorEnumModule.OperatorEnum, Node, false, true, nothing, StatsBase.Weights{Float64, Float64, Vector{Float64}}})

So far I've checked out:
#625 (comment)
#548 (comment)
#447 (comment)
#299 (comment)
#667 (reply in thread)

I've also checked out the API but I have struggled to apply it.

jl.seval("using Zygote")
elementwise_loss = """

function my_loss_function(tree, dataset::Dataset{T,L}, options, idx)::L where {T,L}
    X = idx === nothing ? dataset.X : view(dataset.X, :, idx)
    y = idx === nothing ? dataset.y : view(dataset.y, idx)
    weights = idx === nothing ? dataset.weights : view(dataset.weights, idx)
    derivative_with_respect_to = 3
    prediction, grad, complete = eval_grad_tree_array(tree, X, derivative_with_respect_to, options)
    if !complete
        return L(Inf)
    end
    is_monotonic = all(grad .<= 0) || all(grad .>= 0)
    if !is_monotonic
        return L(Inf)
    end

Thanks in advance

edit: fixed line:

prediction, grad, complete = eval_diff_tree_array(tree, X, options,derivative_with_respect_to)

[gm89uk]

Pysr has really exploded various aspects of my research that is resulting in several manuscripts being written. I am extremely grateful to the team not just for Pysr, but for the contribution you offer the community in your spare time @MilesCranmer and I'm sorry to bump.

Essentially, at lower complexity, there is monotonicity between a feature in x and y; pysr has found great expressions for my problem. At higher complexities, there is overfitting and monotonicity is violated. I'm trying to see if there are slightly more complex expressions it can find if I enforce a monotonicity constraint through the custom loss function.

This is and one other aspect (which I edited out of this request to keep this post clean and I can post it separately), are the final pieces of the puzzle for me for one that is almost ready. Please let me know if anything is unclear or you require any further information, thank you again.

p.s. I'm sorry for the multiple edits, I had updated the code as I made some progress :)

[MilesCranmer]

Hey sorry for the late reply! I was off on vacation last week.

Awesome to hear your experience is going well, that makes me very happy to hear 🙂

So the issue here is that eval_grad_tree_array gives you the whole gradient, and thus only takes 3 arguments. The other function is eval_diff_tree_array (sorry for the weird naming) which only differentiates with respect to a single feature – that function is where you provide it with a 4th argument (as you have done) giving the feature index.

These API calls are documented here – https://symbolicml.org/DynamicExpressions.jl/dev/eval/#Derivatives. (Note that you can pass options instead of the operators::OperatorEnum which is just the lower-level call)

[gm89uk]

Thank you so much, Miles! I hope you had a lovely holiday.

I’m amazed that the solution was so simple—I should have read the API more closely!

The custom loss feature is truly the secret sauce of PySR. It adds another level of versatility.

A breakthrough for me was using the custom loss feature to implement a 'pseudo' multilevel regression approach. I knew that differences between certain groups had a linear difference effect on y, and this approach allowed me to get the best fit without being thrown off by constant differences between the groups.

    unique_group = unique(weights)
    mean_residuals_by_group = Dict{L, T}()
    for member in unique_group
        group_residuals = residuals[weights .== member]
        mean_residuals_by_group[member] = sum(group_residuals) / length(group_residuals)
    end
    adjusted_predictions = copy(prediction)
    for i in eachindex(y)
        adjusted_predictions[i] -= mean_residuals_by_group[weights[i]]
    end
    mse = sum(i -> (adjusted_predictions[i] - y[i])^2, eachindex(y)) / length(y)
    return mse

Thanks again for your help!

[MilesCranmer]

Awesome to hear it worked!! Always happy to hear about your latest use-cases :)