Skip to content

Rewrite log(sqr(x)) to 2 * log(abs(x)) to avoid underflow - #2380

Merged
jessegrabowski merged 11 commits into
pymc-devs:mainfrom
jessegrabowski:fix-log-sqr-prod-stabilization
Aug 28, 2026
Merged

Rewrite log(sqr(x)) to 2 * log(abs(x)) to avoid underflow#2380
jessegrabowski merged 11 commits into
pymc-devs:mainfrom
jessegrabowski:fix-log-sqr-prod-stabilization

Conversation

@jessegrabowski

@jessegrabowski jessegrabowski commented Aug 25, 2026

Copy link
Copy Markdown
Member

Taking the log of a squared product materializes the product first, and in float32 a product of a few hundred terms underflows to zero once squared, so log(abs(sqr(prod(x)))) returns -inf. Shows up in GP code with a few hundred inducing points, on every backend.

log(sqrt(x)) and log(sqr(x)) are now one rewrite, so the square peels off and the existing log(abs(prod(x))) -> sum(log(abs(x))) case finishes the job before anything is squared. Dropping the abs is driven by a new non_negative flag on unary scalar ops instead of a hardcoded list, so abs(exp(x)) folds too. Sqrt deliberately doesn't carry the flag: sqrt(-0.0) is -0.0, and dropping the abs there would flip the sign of a downstream division.

x = pt.vector("x", dtype="float32")
f = pytensor.function([x], pt.log(pt.abs(pt.sqr(pt.prod(x)))))
f(np.full(250, 0.7, dtype="float32"))  # -inf before, -178.3 after

Closes #2379. Sharing one rewrite means sharing the constant's dtype, which is what stops the 0.5 factor truncating to 0 on integer input.

A squared product underflows to zero in float32 at a few hundred elements, sending log(abs(sqr(prod(x)))) to -inf; unwrapping the square lets the existing log(abs(prod(x))) -> sum(log(abs(x))) case fire before the product is ever materialized.
Comment thread pytensor/tensor/rewriting/math.py Outdated
@register_canonicalize
@register_specialize
@node_rewriter([pt_abs])
def local_abs_sqr(fgraph, node):

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Include any scalar op that is non negative? Also use the new syntax x.owner_op_inputs?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a non_negative flag on ScalarOp just for the occasion!

Comment thread pytensor/tensor/rewriting/math.py Outdated
@register_stabilize
@register_specialize
@node_rewriter([log])
def local_log_sqr(fgraph, node):

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

extend the sqrt rewrite instead? it's the same concept?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure why not

Comment thread tests/tensor/rewriting/test_math.py Outdated
pytest.param(lambda x: pt_abs(sqr(x)), id="abs_sqr"),
],
)
def test_sqr_rewrites_skip_complex(original_fn):

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tell your bot about RewriteTest or whatever it's called...

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

: )

I love when people talk about how "we're going to have a tool that learns over time" as if these things aren't totally blank slate every time you boot them up.

Sharing the tail means sharing the constant's dtype, which is what stops the 0.5 factor from truncating to 0 on integer input.

Closes pymc-devs#2379
Comment thread tests/tensor/rewriting/test_math.py Outdated
Comment thread tests/scalar/test_basic.py Outdated

x_test = np.array([-1e3, -2.5, -1.0, -0.0, 0.0, 1e-8, 1.0, 2.5, 1e3])
for op in flagged:
x = pt.vector("x")

@ricardoV94 ricardoV94 Aug 25, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not really a scalar test... Also do we test monotonicity? seems like over-zealous a bit?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah it's dumb

FunctionGraph.replace filters the replacement through the old variable's type, so a rewrite that changed dtype would be rejected outright and these could never fail.
The hazard it covered is still pinned by test_useless_abs_keeps_signed_zero, which fails if Sqrt is ever flagged.
def local_log_sqrt(fgraph, node):
x = node.inputs[0]
@node_rewriter([pt_abs])
def local_useless_abs(fgraph, node):

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does this supersede some other pre-existing more narrow rewrite?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no, but local_log_prod_to_sum_log was hand-rolling the check for non-negative Ops. I switched it to use _is_provably_non_negative (and by extension the new tag)

Comment thread pytensor/tensor/rewriting/math.py Outdated
return [new_out]


def _is_non_negative(var) -> bool:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we have a smarter is non negative somewhere else (used for subtensor lift rewrites maybe)? that we should repurpose/extend with this? It includes stuff like dtype or dimshuffle / ravel ... bla bla

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have _is_provably_non_negative but it doesn't have any logic about functions like exp/softplus/whatever. I could extend it but I thought that was overkill

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a case to _is_provably_non_negative that uses the new non_negative flag, then deleted the duplicate is_negative helper

The flag case is non-strict only, since abs and sqr attain zero, and is restricted to float and unsigned outputs because signed integers wrap on overflow.
Nothing in pytensor or pymc ever wrote the tag.positive hook the removed branch read, and the predicate narrows signed-integer operands, where sqr can wrap negative.
scalar_op = op.scalar_op
if not strict and isinstance(scalar_op, Cast):
return _is_provably_positive(var.owner.inputs[0], strict)
if (

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

follow up issue to move this outside of rewriting/subtensor.py

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done. #2404

@ricardoV94 ricardoV94 left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tests seem a bit overkill since the last iteration, but ok...

When a rewrite correctly does not fire, assert_eval compiles the unchanged graph twice and compares it to itself, so the int8 test now pins the actual wrapped output instead.
@jessegrabowski

Copy link
Copy Markdown
Member Author

I trimmed back the tests, will merge when CI passes again.

@jessegrabowski
jessegrabowski merged commit 597cdb7 into pymc-devs:main Aug 28, 2026
67 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

local_log_sqrt rewrites log(sqrt(x)) to zeros for integer input

2 participants