PySR Hyperparameter setup for large dataset and Dimension Consistancy

[TrailblazerH]

Hi awesome PySR developers,
Thanks for building such a great tool! I'm using PySR on a relatively large dataset (~100,000 data points) and have tuned the hyperparameters following the suggested workflow. I'm using three input features: Feature 1 and Feature 2 each have 40 unique values and Feature 3 has 7 unique values.
The input data consists of two variables with units of radians and one dimensionless variable. Based on previous runs without dimensional constraints, symbolic regression typically generates expressions containing trigonometric functions like sine, cosine, or tangent, which is the expected behavior. However, I've observed cases where the radian-valued variables appear outside of trigonometric functions, resulting in dimensionally inconsistent expressions. To enforce dimensional correctness, I'm implementing dimensional constraints that should restrict radian-valued variables to only appear within trigonometric functions.
However, I am running into an issue with dimensional constraints. As you can see below, I add dimensional constraints, but the weird thing happens after I add the rad units here:
model.fit(X, y, X_units=["", "rad", "rad"], y_units="")
The expression it found how does not include cos/sin/tan function anymore? Do you have any idea what could be the cause?
If possible, could you please take a quick look and let me know if the whole setup seems reasonable for this large dataset symbolic regression?
Here’s the hyperparameter setup I’m using:

model = PySRRegressor(procs=32,
                      niterations=5000, 
                      populations=96,
                      maxsize=50, 
                      maxdepth=10, 
                      population_size=50,                      
                      ncycles_per_iteration=5000,
                      binary_operators=["+", "-", "*", "/"],
                      unary_operators=["safe_sin(x) = isfinite(x) ? sin(x) : convert(typeof(x), NaN)", 
                                       "safe_cos(x) = isfinite(x) ? cos(x) : convert(typeof(x), NaN)",
                                       "safe_tan(x) = isfinite(x) ? tan(x) : convert(typeof(x), NaN)",
                                       "log", "sqrt"],
                      constraints={'safe_sin': 3, 'safe_cos': 3,'safe_tan': 3, 'log': 3, 'sqrt': 3},
                      nested_constraints={"safe_sin": {"safe_sin": 0, "safe_cos": 0, "safe_tan": 0, "log": 0},
                                          "safe_cos": {"safe_sin": 0, "safe_cos": 0, "safe_tan": 0, "log": 0},
                                          "safe_tan": {"safe_sin": 0, "safe_cos": 0, "safe_tan": 0, "log": 0},
                                          "log": {"safe_sin": 0, "safe_cos": 0, "safe_tan": 0, "log": 0},
                                          "sqrt": {"sqrt": 0, "log": 0}},
                      extra_sympy_mappings = {
                          "safe_sin": lambda x: sympy.sin(x),
                          "safe_cos": lambda x: sympy.cos(x),
                          "safe_tan": lambda x: sympy.tan(x),},
                      dimensional_constraint_penalty=1e2,
                      dimensionless_constants_only=True,
                      parsimony=1e-6,
                      early_stop_condition=("stop_if(loss,complexity) = loss<1e-6 && complexity<10"),
                      logger_spec=logger_spec,
                      progress=True,
                      batching=True,
                      weight_optimize=0.001,
                      turbo=True)
                      

model.fit(X, y, X_units=["", "rad", "rad"], y_units="")

Any feedback or suggestions for improvement would be much appreciated!
Thanks again! 🙌

[MilesCranmer]

Hi @TrailblazerH,

In SI units (the dimensional analysis system uses SI), note that 1 rad = 1, so a radian as a unit has no effect on dimensional analysis. cos and sin are allowed to take unitless input, since there is no difference to radian.

If you would like the units to be cancelled out, one option is to switch to a different unit instead (like distance)

Also be sure to read the tuning page: https://ai.damtp.cam.ac.uk/pysr/tuning/

Cheers,
Miles

[MilesCranmer]

Also, you can use regular cos, sin; no need to write your own safe version (they get mapped internally to safe versions!)

[TrailblazerH]

Thanks for your quick response. I didn’t quite understand what you meant by “switching a unit.” Did you mean converting from radians to degrees?
I’ve found this part a bit tricky. I’m trying to define my_loss so that sin/cos only apply to variables with radian units, as you mentioned earlier in #927 . However, even with that setup, I still get expressions where the first feature (which is unitless) appears inside sin/cos/tan, and the third feature (which is in radians) doesn’t always show up as expected.
Do you have any idea why this might be happening?

objective = """
function my_loss(tree, dataset::Dataset{T,L}, options) where {T,L}
    feature1_idx = 1    # 1st feature no unit
    feature3_idx = 3   # 3rd feature unit = rad
    
    bad_yaw_count = 0
    bad_r_count = 0
    
    for node in tree
        if node.degree == 0 && !node.constant
            if node.feature == feature3_idx
                parent_is_trig = false
                current = node
                while hasproperty(current, :parent) && current.parent !== nothing
                    parent = current.parent
                    if parent.degree == 1 && (parent.op in [
                        options.operators.unary[:safe_sin],
                        options.operators.unary[:safe_cos],
                        options.operators.unary[:safe_tan]
                    ])
                        parent_is_trig = true
                        break
                    end
                    current = parent
                end
                if !parent_is_trig
                    bad_yaw_count += 1
                end
            end
            
            if node.feature == feature1_idx
                current = node
                while hasproperty(current, :parent) && current.parent !== nothing
                    parent = current.parent
                    if parent.degree == 1 && (parent.op in [
                        options.operators.unary[:safe_sin],
                        options.operators.unary[:safe_cos],
                        options.operators.unary[:safe_tan]
                    ])
                        bad_r_count += 1
                        break
                    end
                    current = parent
                end
            end
        end
    end
    
    prediction, valid = eval_tree_array(tree, dataset.X, options)
    if !valid
        return L(Inf)
    end
    y = dataset.y
    mse = sum(i -> abs2(prediction[i] - y[i]), eachindex(prediction, y)) / length(y)
    
    total_penalty = 0.0
    if bad_yaw_count > 0
        total_penalty += 100.0 * bad_yaw_count * mse
    end
    if bad_r_count > 0
        total_penalty += 100.0 * bad_r_count * mse  # Same penalty weight for both violations
    end
    
    return mse + total_penalty
end

"""

[MilesCranmer]

I recommend testing this loss on the Julia side first if possible, since it is fairly complex. e.g., options.operators.unary[:safe_sin] is an invalid call.

[MilesCranmer]

Note that hasproperty(current, :parent) is always false. There is no such property. The node type is detailed here: https://ai.damtp.cam.ac.uk/symbolicregression/dev/types/

[MilesCranmer]

You also don’t need to pass X_units since it doesn’t appear to be used anywhere in the loss.