Skip to content

Commit 8dc0bf6

Browse files
author
Python Infra CI
committed
perf: avoid a redundant full-array copy in scipy.fft hfft/hfftn adapters
1 parent a3b1cb7 commit 8dc0bf6

3 files changed

Lines changed: 177 additions & 4 deletions

File tree

‎mkl_fft/interfaces/_scipy_fft.py‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -523,8 +523,10 @@ def hfft(
523523
_check_plan(plan)
524524
x = _validate_input(x)
525525
norm = _swap_direction(norm)
526-
x = np.array(x, copy=True)
527-
np.conjugate(x, out=x)
526+
# np.conjugate(x) already allocates a fresh array, so there is no need
527+
# to copy x first and then conjugate it in place (which would touch
528+
# every element twice instead of once).
529+
x = np.conjugate(x)
528530

529531
with _Workers(workers):
530532
# Note: overwrite_x is not utilized
@@ -626,8 +628,10 @@ def hfftn(
626628
_check_plan(plan)
627629
x = _validate_input(x)
628630
norm = _swap_direction(norm)
629-
x = np.array(x, copy=True)
630-
np.conjugate(x, out=x)
631+
# np.conjugate(x) already allocates a fresh array, so there is no need
632+
# to copy x first and then conjugate it in place (which would touch
633+
# every element twice instead of once).
634+
x = np.conjugate(x)
631635
s, axes = _init_nd_shape_and_axes(x, s, axes, invreal=True)
632636

633637
with _Workers(workers):

‎pr-body.md‎

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
## Track
2+
3+
Performance (Track A), item 2 on the ranked list: an avoidable copy in
4+
`_fft_utils.py` / interface layer — here `mkl_fft/interfaces/_scipy_fft.py`.
5+
6+
## What I examined
7+
8+
I read the root `AGENTS.md`, `.github/copilot-instructions.md`, and the nested
9+
`AGENTS.md` files for `mkl_fft/`, `mkl_fft/interfaces/`, `mkl_fft/src/`, and
10+
`mkl_fft/tests/` before touching anything. I then read `_pydfti.pyx`,
11+
`_fft_utils.py`, `_mkl_fft.py`, `src/mklfft.c.src`, and the three interface
12+
adapter files (`_numpy_fft.py`, `_scipy_fft.py`, `_float_utils.py`,
13+
`_numpy_helper.py`) end to end, looking first for descriptor-handling or
14+
memory-management bugs per the bug-fix priority list, then for avoidable
15+
copies and per-call overhead per the performance priority list.
16+
17+
I specifically checked the `out=x` aliasing path in `_c2c_fft1d_impl`
18+
(`_pydfti.pyx`) against the MKL descriptor being configured
19+
`DFTI_PLACEMENT = DFTI_NOT_INPLACE` in the `*_out` C routines
20+
(`src/mklfft.c.src`), since passing identical input/output pointers to a
21+
not-in-place descriptor looked suspicious. It turned out `mkl_fft/tests/
22+
test_fft1d.py` (`test_vector5`, `test_vector6`, `test_matrix4`) already
23+
exercises exactly this case as a documented, presumably-passing feature
24+
("fft in-place is the same as fft out-of-place"), so I did not treat this as
25+
a new finding and left it alone — reopening settled, tested behaviour without
26+
being able to run the suite myself would be speculative.
27+
28+
## What I found
29+
30+
In `mkl_fft/interfaces/_scipy_fft.py`, `hfft` and `hfftn` each computed the
31+
Hermitian-conjugate of their input as two full passes over the array:
32+
33+
```python
34+
x = np.array(x, copy=True) # pass 1: copy every element
35+
np.conjugate(x, out=x) # pass 2: negate every imaginary part in place
36+
```
37+
38+
`x` at this point has already been through `_validate_input`, which calls
39+
`np.asarray(x)` (via `_supported_array_or_not_implemented`), so it is always
40+
a plain `ndarray` (subclasses are already stripped, `subok=False` is
41+
`np.asarray`'s default) — the explicit copy exists solely so that the
42+
in-place `conjugate` does not mutate the caller's own array, not for any
43+
type normalization. `np.conjugate(x)` called without `out=` already
44+
allocates a fresh output array and computes the conjugate directly from the
45+
input in a single pass, giving the identical result without the caller's
46+
array ever being touched. The two-step form is therefore a copy taken
47+
defensively where a single elementwise op already does the job — exactly the
48+
"avoidable copies" pattern called out as the second-ranked optimization
49+
target in the mandate. The sibling implementation `_numpy_fft.hfft` already
50+
uses the single-call form (`np.conjugate(x)`), so this is also a small
51+
internal inconsistency between the two interface adapters.
52+
53+
## The change
54+
55+
In both `hfft` and `hfftn` in `mkl_fft/interfaces/_scipy_fft.py`, replaced:
56+
57+
```python
58+
x = np.array(x, copy=True)
59+
np.conjugate(x, out=x)
60+
```
61+
62+
with:
63+
64+
```python
65+
x = np.conjugate(x)
66+
```
67+
68+
`ihfft`/`ihfftn` were not touched: they conjugate the *output* of `rfft`/
69+
`rfftn` in place, and that output buffer was just freshly allocated by
70+
`mkl_fft.rfft`/`rfftn` for this call, so there is no caller-owned buffer at
71+
risk and no redundant copy to remove — the one-pass in-place conjugate there
72+
is already optimal.
73+
74+
## Why this is behaviour-preserving
75+
76+
- **Correctness of the result:** `np.conjugate(a)` and
77+
`(lambda b: (np.conjugate(b, out=b), b)[1])(np.array(a, copy=True))`
78+
compute the same values element-by-element; the only difference is whether
79+
the negation of the imaginary part happens into a copy-then-mutate buffer
80+
or straight into a freshly allocated one.
81+
- **dtype:** unaffected — conjugate preserves the input dtype in both forms
82+
(real dtypes are a no-op copy, complex dtypes negate the imaginary part);
83+
`x` here is always a plain array of a dtype `_validate_input` already
84+
approved (not `float16`/`float128`/`complex256`).
85+
- **Shape:** unaffected — `conjugate` is elementwise and shape-preserving in
86+
both forms.
87+
- **Strides / memory layout:** both `np.array(x, copy=True)` and the
88+
`conjugate` ufunc default to `order='K'` (preserve layout as far as
89+
possible), so the freshly allocated array has the same layout
90+
characteristics the old two-step form produced.
91+
- **Aliasing:** the old code copied specifically so the in-place conjugate
92+
would not mutate the caller's array. `np.conjugate(x)` without `out=`
93+
never writes into `x`, so the caller's array is equally untouched under
94+
the new code — the safety property the copy existed for is preserved by a
95+
different, cheaper mechanism.
96+
- **norm/out parameters:** untouched by this diff; `hfft`/`hfftn` do not
97+
accept an `out=` argument to `mkl_fft.irfft`/`irfftn` in this file (`# Note:
98+
overwrite_x is not utilized`), so there is no interaction with an
99+
output buffer supplied by the caller.
100+
101+
## Mechanism, one sentence
102+
103+
Two full sequential passes over the array (allocate-and-copy, then negate
104+
imaginary parts in place) are replaced by one pass that allocates and
105+
negates at the same time, halving the memory traffic per element for every
106+
call to `scipy.fft.hfft`/`hfftn` routed through `mkl_fft.interfaces`.
107+
108+
## Verification
109+
110+
I have no shell access in this environment, so I could not run `pytest`,
111+
`black`, `isort`, `flake8`, `pylint`, or `cython-lint` myself. I re-read the
112+
edited file in full after the change (shown above) to confirm the two edits
113+
are the only difference from the original and that indentation/imports are
114+
unchanged; `numpy` was already imported at module scope, so no import
115+
changes were needed. I did not add a new test because this is a
116+
performance-only change with no behavioural difference to assert beyond what
117+
the existing `scipy.fft` compatibility tests already cover (`mkl_fft/tests/
118+
third_party/scipy/test_basic.py` exercises `hfft`/`hfftn` through
119+
`mkl_fft.interfaces.scipy_fft`); CI's existing suite is the check that the
120+
result is still numerically identical, and CI's ASV benchmark is the
121+
authority on whether it is faster.
122+
123+
Commands run: none (no shell in this environment). Local measurement: none
124+
taken; this report states the mechanism (one memory pass instead of two) as
125+
the basis for the hypothesis, per the mandate, and leaves confirmation to
126+
CI's benchmark job.
127+
128+
## Rejected candidates
129+
130+
- Relaxing the exact-stride-match requirement before reusing a caller's
131+
`out=` array in `_c2c_fft1d_impl` (`_pydfti.pyx`), which is marked `TODO`
132+
in the source. This would remove a real defensive copy for a wider range
133+
of strided `out=` arrays, but the comment already flags it as needing
134+
careful validation of what MKL actually tolerates for input/output stride
135+
relationships, and getting it wrong risks silently wrong FFT output. Left
136+
as a candidate for a human with the ability to run the MKL-backed test
137+
suite.
138+
- Reordering the dtype checks in `_downcast_float128_array`
139+
(`interfaces/_float_utils.py`) to test `isinstance(x, np.ndarray)` before
140+
the `longdouble`/`clongdouble` comparisons, to skip redundant work for the
141+
common plain-`ndarray` case. The saving is a couple of dtype equality
142+
checks per call — real but small enough that I could not state a clear
143+
single-sentence mechanism distinguishing it from noise, so I left it alone
144+
per "if you cannot say why a change is faster in one sentence, it is not a
145+
candidate."
146+
- Treating the `out=x` (aliasing) path through a `DFTI_NOT_INPLACE`
147+
descriptor in `_pydfti.pyx`/`src/mklfft.c.src` as a bug. Rejected because
148+
existing tests (`test_fft1d.py::Test_mklfft_vector::test_vector5/6`,
149+
`Test_mklfft_matrix::test_matrix4`) already assert this exact case is
150+
correct and treat it as a supported feature; reopening it without being
151+
able to run those tests would be guessing against settled, tested
152+
behaviour.
153+
154+
## Left to humans
155+
156+
- The stride-relaxation `TODO`s in `_pydfti.pyx` (`_c2c_fft1d_impl`,
157+
`_r2c_fft1d_impl`, `_c2r_fft1d_impl`) for reusing a caller's `out=` array
158+
under a wider set of compatible strides than exact match / both-contiguous.
159+
This is the higher-ranked "descriptor handling" and "avoidable copies"
160+
opportunity in the mandate, but validating exactly what MKL DFTI accepts
161+
for mismatched input/output strides needs to be checked against MKL
162+
documentation and run against the real test suite, which I cannot do here.
163+
- Any compiler/build flag changes (explicitly out of scope for me to
164+
implement per the mandate; I did not find a specific one to propose this
165+
run).
166+
- `CHANGELOG.md` wording: not proposed, because this run made no
167+
user-visible behaviour change (performance-only), and the mandate only
168+
asks for changelog wording for bug fixes.

‎pr-title.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
perf: avoid a redundant full-array copy in scipy.fft hfft/hfftn adapters

0 commit comments

Comments
 (0)