Implement mathematical constraints

[VOLCANIC-9]

I want to implement mathematical constraints in a soft way in Pysr. I think if I apply them in a hard way I lose a lot of functions. These constraints are like I know the positiveness or negativeness of function with respect to some values, and behavior of function derivative with respect to some variables.

I find that it seems the main of your program is written in Julia. Then you wrapped it to work in Python. I tried to read "src" part of your code but it's hard for me because I am not familiar with Julia. So, I read the discussion part of the page and I saw this question: Constraints? "#304"
In answer, you said: "Update: added full_objective (PySR) and loss_function (SymbolicRegression.jl) for this purpose."

Actually, I did not find "full_objective" part in the Python code. I found and read "loss_function" in Julia's part. I also read this question: "#449" and this question:"#256"
It seems to me you are trying to say that you can define a new loss function. I want to know: can I define this new loss criterion in Python? If I define a new loss function I must download your code manipulate it and compile it again? I see your loss criterion repeated somewhere in the code, Do I need to make changes elsewhere?

[MilesCranmer]

To implement a custom loss about function constraints, the right way is to pass a Julia function as a string in Python to the loss_function parameter of PySRRegressor. You unfortunately cannot use a pure Python function as the loss function. This is technically possible but it would restrict you to serial computation only, due to Python’s “global interpreter lock”, and would also be very very slow as it couldn’t be compiled or differentiated (whereas Julia code can be).

If you describe your constraints exactly or provide an example in Python I can help get it working.

The examples that you cited are indeed the most relevant.

[VOLCANIC-9]

I tried to learn Julia. I also checked the SymbolicRegression.jl/src/LossFunctions.jl section. To clarify what is going on in my mind, I give an example:
Consider a multivariate function like y=f(x1,...,xn).
Suppose that I know that the derivative of the function is always positive concerning one of the variables (d(f(x1,...,xn))/dx1>0 for all x1,..,xn).
I think the estimated value of y and the existing value of y are present, so the function derivative is not obtained at all points.
Also, you defined a function to evaluate the complexity of mathematical expressions.
I want to define a loss function with: the distance the estimated value of y and the existing value of y, the complexity, and the positive value of the derivative concerning a variable.
What do you suggest?
Thanks

[MilesCranmer]

Suppose that I know that the derivative of the function is always positive concerning one of the variables (d(f(x1,...,xn))/dx1>0 for all x1,..,xn).

For this you could use eval_diff_tree_array: https://symbolicml.org/DynamicExpressions.jl/dev/eval/#Evaluation-and-Derivatives

function eval_loss(tree, dataset::Dataset{T,L}, options)::L where {T,L}
    derivative_with_respect_to = 1
    predicted, gradient, complete = eval_diff_tree_array(tree, dataset.X, derivative_with_respect_to, options)
    if !complete
        # encountered NaN/Inf, so return early
        return L(Inf)
    end
    positivity = sum(i -> gradient[i] > 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
    sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
    beta = L(1e-3)
    return (sum_square_loss + beta * positivity) / dataset.n
end

where beta is some tunable regularisation parameter.

Here, the line

gradient[i] > 0 ? L(0) : abs2(gradient[i])

says that, if the gradient of element i is above 0, then there's no penalty. Otherwise, penalise the square distance from 0.

So pass this as a string to loss_function.

[VOLCANIC-9]

Thank you for producing the sample, but I have a lot of confusion which it did not solve even by studying the link:

1. which parameters are passed to your function?
   I edited "loss_function=eval_loss()"
   it shows an error that eval_loss() needs parameters.
   I think "tree, dataset::Dataset{T,L}, options" are the parameters.
   what is the "tree"? I can't understand it even with the study link.
   what is the "T, L"? I think they are X and Y which I must change X to a new matrix with parameters in columns and Y is a column vector. Type of {T,L} is a dataset so I think I must create a structure.
   what is the "options"?
2. Can I pass your defined loss function as string PySr? I tried it but it shows multiple errors.
   I sincerely appreciate your patience in answering the questions.

[MilesCranmer]

You should do loss_function="eval_loss" in Python, or loss_function=eval_loss in Julia. But using eval_loss() is a function call. Basically you need to pass the function itself to this parameter, rather than an evaluation of it.

For these other questions you might find the docs examples useful: https://astroautomata.com/PySR/examples/

[VOLCANIC-9]

It is crazy!!! your algorithm eliminates variables instead of accepting the constraints
I have 4 variables. Two of them have constraints. All of them have units.
first I write this:

function eval_loss(tree, dataset::Dataset{T,L}, options)::L where {T,L}
    derivative_with_respect_to = 1
    predicted, gradient, complete = eval_diff_tree_array(tree, dataset.X, options, derivative_with_respect_to)
    if !complete
        # encountered NaN/Inf, so return early
        return L(Inf)
    end
	i=1
    positivity = sum(i -> gradient[i] > 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
    sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
	    
	i=2
    positivity = sum(i -> gradient[i] < 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
    sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
	
	
	
    beta = L(1e-2)
    return (sum_square_loss + beta * positivity) / dataset.n
end

code eliminates variable 1 and variable 2
Then I try to force the algorithm to use them; I write this:

function eval_loss(tree, dataset::Dataset{T,L}, options)::L where {T,L}
    derivative_with_respect_to = 1
    predicted, gradient, complete = eval_diff_tree_array(tree, dataset.X, options, derivative_with_respect_to)
    if !complete
        # encountered NaN/Inf, so return early
        return L(Inf)
    end
	i=1
    positivity = sum(i -> gradient[i] > 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
    sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
	
	contains_x1 = any(
        node -> (
            node.degree == 0
            && !(node.constant)
            && node.feature == 1
        ),
        tree
    )
	if !(contains_x1)
        # penalty term
        sum_square_loss += L(1e10)
    end
	    
	i=2
    positivity = sum(i -> gradient[i] < 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
    sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
	
	contains_x3 = any(
        node -> (
            node.degree == 0
            && !(node.constant)
            && node.feature == 1
        ),
        tree
    )
	if !(contains_x1)
        # penalty term
        sum_square_loss += L(1e10)
    end
	
    beta = L(1e-2)
    return (sum_square_loss + beta * positivity) / dataset.n
end

It doesn't work. I increased the regularisation parameter but the algorithm can't converge.
Also, algorithm did not use variables that constant for some cycles in data.

[MilesCranmer]

For

contains_x3 = any(
        node -> (
            node.degree == 0
            && !(node.constant)
            && node.feature == 1
        ),
        tree
    )

I’m assuming the node.feature == 1 was supposed to be == 3?

[MilesCranmer]

Also note you have if !(contains_x1) twice

[VOLCANIC-9]

I just tried to produce a sample; in the main code, all is correct. the interesting part is the algorithm keeps two parameters that do not match with y unit.

[MilesCranmer]

FYI I think the constraint you wrote only enforces that the model has x1

[MilesCranmer]

Also note that when you write

i = 1
i -> gradient[i]

This does not set gradient[1]. It is a new definition of i.

[VOLCANIC-9]

so you say that this part
" i=1
positivity = sum(i -> gradient[i] > 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
"
doesn't mean: the derivative with respect to variable 1 must be greater than zero
and this part:
" i=2
positivity = sum(i -> gradient[i] > 0 ? L(0) : abs2(gradient[i]), eachindex(gradient))
sum_square_loss = sum(i -> abs2(predicted[i] - dataset.y[i]), eachindex(predicted, dataset.y))
"
doesn't mean: the derivative with respect to variable 2 must be greater than zero ???
So how I say derive with respect to variable 1?

" FYI I think the constraint you wrote only enforces that the model has x1"


Ok, but the model just uses variables x3 and x4

[MilesCranmer]

gradient is a 2D array. Also You are using the symbol i for two different things which is why it’s no longer equal to 1