Skip to content

Python scalar operands in binary elementwise functions synchronize the queue #3148

Open

Description

@cakedev0

A binary elementwise function with a Python or NumPy scalar operand (2 * x, x == 0, x / n) blocks the host until all work already queued on the device is done. Array API code is full of these scalars, so every such op becomes a full sync and the host can no longer run ahead of the device.

On an Arc B580, 2 * x issued right after a ~80 ms matmul returns only when the matmul is done:

call issued after big @ big host time
x * x 0.05 ms
2 * x 78 ms
x == 0 78 ms
x / 3 78 ms
x * np.float32(2) 78 ms

In a loop with one sync at the end, the scalar version costs ~40 us more per call on a small array and ~115 us more on a (65536, 128) array. torch XPU does not have this gap:

op, float32 dpnp, x.shape=(16,) dpnp, x.shape=(65536, 128) torch XPU, x.shape=(65536, 128)
x * x 19 us 190 us 169 us
2 * x 60 us 305 us 170 us

In #2087, hoisting the scalar out of the loop was suggested as a workaround. That works in user code, but not in libraries written against the array API standard, such as scikit-learn, where 2 * x and x == 0 are written inline and the namespace is not known in advance.

Impact

scikit-learn LogisticRegression (lbfgs, array_api_dispatch=True) runs ~14 binary ops with a scalar operand per iteration. With these ops made non-blocking (a local patch: same values, same iterations), the fit is 10% to 16% faster on four of the six shapes below, 46% faster on the wide one, and unchanged on the largest one. Fit times on make_classification data, median of 3 fits, scikit-learn 1.9.1, dpnp 0.20.0, torch 2.14.0+xpu:

n_samples x n_features, dtype iterations dpnp dpnp, scalar ops not blocking torch XPU
36000 x 126, float64 17 70 ms 60 ms 14 ms
200000 x 100, float32 18 99 ms 83 ms 18 ms
20000 x 1000, float64 20 92 ms 78 ms 29 ms
3200 x 25300, float32 12 142 ms 76 ms 37 ms
2000000 x 20, float32 13 221 ms 199 ms 41 ms
8000000 x 20, float64 11 1111 ms 1105 ms 237 ms

This sync is only part of the gap to torch XPU. On the real datasets of our benchmarks, LogisticRegression is 1.4x to 5.8x slower on dpnp than on torch XPU on this GPU. I'll report the other causes separately once I've found them.

KMeans. An array API KMeans is 2x to 3.5x slower on dpnp than on torch XPU (ogrisel/cursor-sandbox#17 (comment)). This sync is one of several causes there. I'll report the others separately.

Reproducer
import time

import dpnp
import numpy as np

try:  # optional reference column
    import torch
except ImportError:
    torch = None


def bench(fn, sync, reps=200):
    """Best of 5, microseconds per call, one sync per `reps` calls."""
    fn()
    sync()
    best = float("inf")
    for _ in range(5):
        t0 = time.perf_counter()
        for _ in range(reps):
            fn()
        sync()
        best = min(best, (time.perf_counter() - t0) / reps)
    return 1e6 * best


rng = np.random.default_rng(0)
for shape in [(16,), (65536, 128)]:
    a = rng.standard_normal(shape, dtype=np.float32)
    x = dpnp.asarray(a)
    if torch is not None:
        t = torch.as_tensor(a, device="xpu")
    print(f"\ndpnp {dpnp.__version__}, x.shape={shape}, float32\n")
    print("| op | dpnp (us) | torch XPU (us) |\n| --- | ---: | ---: |")
    for name, fd, ft in [
        ("x * x", lambda: x * x, lambda: t * t),
        ("2 * x", lambda: 2 * x, lambda: 2 * t),
    ]:
        b = f"{bench(ft, torch.xpu.synchronize):.0f}" if torch is not None else ""
        print(f"| `{name}` | {bench(fd, x.sycl_queue.wait):.0f} | {b} |")

# Does a scalar operand wait for work already queued on the device?
x = dpnp.asarray(rng.standard_normal(16, dtype=np.float32))
big = dpnp.asarray(rng.standard_normal((8192, 8192), dtype=np.float32))
q = x.sycl_queue
print("\nHost time of one call issued right after an ~80 ms matmul:\n")
for name, f in [
    ("x * x", lambda: x * x),
    ("2 * x", lambda: 2 * x),
    ("x == 0", lambda: x == 0),
    ("x / 3", lambda: x / 3),
    ("x * np.float32(2)", lambda: x * np.float32(2)),
]:
    times = []
    for _ in range(3):
        big @ big
        q.wait()
        t0 = time.perf_counter()
        big @ big
        t1 = time.perf_counter()
        f()
        t2 = time.perf_counter()
        q.wait()
        times.append(t2 - t1)
    print(f"- `{name}`: {1e3 * min(times):.2f} ms")

Same numbers with dpnp 0.21.0 (dpctl 0.23.0), and master (98155b1) has no change in this area since 0.21.0. Environment: Intel Arc B580, compute runtime 26.31.39395 (Level Zero 1.17.39395), Ubuntu 26.04, kernel 7.0 (xe driver), Python 3.12 (microbenchmarks) and 3.14 (scikit-learn benchmarks), torch 2.14 XPU from PyPI.

Found with Claude Code while profiling the benchmarks above.

Activity

  1. ndgrigorian commented on Nov 3, 2026

    @ndgrigorian
    Collaborator

    @cakedev0
    Yes, this is due to bringing the scalar to the device as a 0D array/array scalar. We have seen issues with performance as a result of this before.

    The only real solution given the use case, that I can immediately think of, is the creation of kernels that operate with a scalar brought to the device without allocation (i.e., as an argument to the kernel, relying on device-copyable semantics). But previously, we found this to be of limited benefit performance-wise. May have changed with newer compiler versions.

  2. cakedev0 commented on Nov 3, 2026

    @cakedev0
    Author

    Thanks @ndgrigorian! I think the sync comes from how the 0-D array is filled rather than from the 0-D array itself, so there may be a cheap fix before going to scalar kernel arguments.

    The scalar goes through dpt.asarray (_elementwise_common.py#L655-L659, also L603 and L953), which ends in _copy_from_numpy_into. There, _copy_numpy_ndarray_into_usm_ndarray is a synchronizing call that depends on all submitted events (_copy_utils.py#L135-L138), so the host waits for the whole queue, matmul included.

    dpt.full((), value, dtype=..., sycl_queue=exec_q) also makes a 0-D device array, but it passes the value to _full_usm_ndarray and does not wait. The "dpnp, scalar ops not blocking" column in the LogisticRegression table above comes from exactly that local patch: Python and NumPy scalars go through dpt.full at these three places, everything else stays on dpt.asarray. It gives the same values and the same iteration counts.

    This only removes the sync. Each scalar operand still costs one small allocation and one fill launch, which is what passing the scalar as a kernel argument would remove. I haven't measured how much of the 60 us vs 19 us gap on small arrays that part accounts for.

  3. added 2 commits that reference this issue on Nov 3, 2026
  4. ndgrigorian commented on Nov 3, 2026

    @ndgrigorian
    Collaborator
    `dpt.full((), value, dtype=..., sycl_queue=exec_q)` also makes a 0-D device array, but it passes the value to `_full_usm_ndarray` and does not wait. The "dpnp, scalar ops not blocking" column in the LogisticRegression table above comes from exactly that local patch: Python and NumPy scalars go through `dpt.full` at these three places, everything else stays on `dpt.asarray`. It gives the same values and the same iteration counts.
    This only removes the sync. Each scalar operand still costs one small allocation and one fill launch, which is what passing the scalar as a kernel argument would remove. I haven't measured how much of the 60 us vs 19 us gap on small arrays that part accounts for.
    

    Very good idea! For some reason, it didn't cross my mind. I'll think about how this might look (maybe np.isscalar helper could work) and we can evaluate after a refactor.

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

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions