Skip to content

Bug: axpy!/axpby! fail with MissingPrimalError (iszero check) #317

Description

@bdrhill

Description

LinearAlgebra.axpy! and LinearAlgebra.axpby! fail with TracerSparsityDetector (global sparsity) when the scalar coefficient is a tracer. This is because the generic implementation checks iszero(α) to short-circuit when the coefficient is zero.

MWE

using SparseConnectivityTracer
using LinearAlgebra

detector = TracerSparsityDetector()

f(x) = begin
    y = copy(x[1:3])
    axpy!(x[4], x[5:7], y)  # y = y + x[4] * x[5:7]
    y
end

jacobian_sparsity(f, rand(7), detector)
# ERROR: Function iszero requires primal value(s).
Stacktrace
Function iszero requires primal value(s).
A dual-number tracer for local sparsity detection can be used via `TracerLocalSparsityDetector`.
Stacktrace:
  [1] iszero(t::SparseConnectivityTracer.GradientTracer{Int64, BitSet})
    @ SparseConnectivityTracer ~/dev/SparseConnectivityTracer.jl/src/overloads/is_functions.jl:13
  [2] axpy!(α::SparseConnectivityTracer.GradientTracer{Int64, BitSet}, x::Vector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, y::Vector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}})
    @ LinearAlgebra /nix/store/.../LinearAlgebra/src/generic.jl:1636

Root Cause

LinearAlgebra.axpy! in generic.jl calls iszero(α) on the scalar coefficient to optimize the case when α=0. This fails with global tracers that don't have primal values.

The same issue affects axpby! which also checks iszero.

Workaround

Use TracerLocalSparsityDetector instead:

local_detector = TracerLocalSparsityDetector()
jacobian_sparsity(f, rand(7), local_detector)  # Works!

Potential Fix

Add overloads for axpy! and axpby! that work with tracers, avoiding the iszero check. The overload could simply perform the operation unconditionally since the sparsity pattern shouldn't depend on whether the coefficient happens to be zero.

Environment

julia> using Pkg; Pkg.status()
Project SparseConnectivityTracer v1.2.1
Status `~/dev/SparseConnectivityTracer.jl/Project.toml`
  [47edcb42] ADTypes v1.22.0
  [37e2e46d] LinearAlgebra v1.12.0
  ...

julia> versioninfo()
Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)

Activity

  1. adrhill commented on May 18, 2026

    @adrhill
    Collaborator

    For context:

      axpy!(α, x::AbstractArray, y::AbstractArray)
    
      Overwrite y with x * α + y and return y. If x and y have the same axes, it's equivalent with y .+= x .* a.
    
  2. added
    arrayFeatures regarding array overloads
    new overloadsA new method on tracers is required by a user.
    agent 🤖Agentically generated issue
    on May 18, 2026
  3. webbrain-one commented on Aug 2, 2026

    @webbrain-one

    Thanks for the detailed report! I’ve opened a PR to address this: #326. It adds custom overloads for axpy! and axpby! that skip the iszero(α) short-circuit, allowing them to work smoothly with global sparsity tracers. Please take a look and let me know if it covers your use case or if you spot any edge cases. Thanks again!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    agent 🤖Agentically generated issuearrayFeatures regarding array overloadsnew overloadsA new method on tracers is required by a user.

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions