PERF: block-wise arithmetic for frame-with-frame
State: closed

- [ ] closes #xxxx
- [ ] tests added / passed
- [x] passes `black pandas`
- [x] passes `git diff upstream/master -u -- "*.py" | flake8 --diff`
- [x] whatsnew entry


Labels:
name: Performance
description: Memory or execution speed performance
name: Internals
description: Related to non-user accessible pandas implementation

Comments:
Did numpy recently change something that would make `np.errstate(all="ignore")` not suppress a DeprecationWarning?
also this has a linked issue IIRC? (you just commented today on it, I think)
will add comments as requested

before implementing helper functions, I'm going to wait for the where/putmask PRs (#32769,#32791, #32846) to go through, as I think that will allow us to simplify this in a more elegant way
Hello @jbrockmendel! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:























There are currently no PEP 8 issues detected in this Pull Request. Cheers! :beers: 

##### Comment last updated at 2020-05-17 21:35:14 UTC
Refactored the new function here so as to have one main execution path.  This makes the debugging assertions clunkier, but whats left is much smoother than the previous implementation.
Just pushed with a some new asvs, about 200x speedup in the best-case scenario:

```
arr = np.random.randn(10 ** 6).reshape(500, 2000).astype(np.float64)
df = pd.DataFrame(arr)
df[1000] = df[1000].astype(np.float32)
df.iloc[:, 1000:] = df.iloc[:, 1000:].astype(np.float32)

df2 = pd.DataFrame(arr)
df2[1000] = df2[1000].astype(np.int64)
df2.iloc[:, 500:1500] = df2.iloc[:, 500:1500].astype(np.int64)

df._consolidate_inplace()
df2._consolidate_inplace()

In [35]: %timeit df + df2                                                                                                                                                                         
572 ms ± 55 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # <-- master
113 ms ± 10.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)  # <-- PR

In [38]: %timeit df + df                                                                                                                                                                                
527 ms ± 69.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # <-- master
2.51 ms ± 24.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # <-- PR
```
> I assume you're happy with the existing test coverage here?

Thats a tough one; everything in tests.arithmetic is overwhelmingly single-block.  ill take a look in tests.frame.test_arithmetic to see what we have for mixed-block
likely need a whatsnew entry (can be a followon). @jorisvandenbossche if ok here.
> likely need a whatsnew entry

has one

> I would consider moving this to a new module to de-nest these functions for blockwise operations (e.g. pandas/core/ops/blockwise.py), but can be later.

yah, i think there might be a way to make this go through `mgr.apply`, will look at other options. 
So I did some timings with the example from above https://github.com/pandas-dev/pandas/pull/32779#issuecomment-607018660:

```
In [2]: %timeit df + df2
414 ms ± 43 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # master
90.8 ms ± 15.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)  # PR
28.7 ms ± 333 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)  # column-wise

In [3]: %timeit df + df 
378 ms ± 7.45 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # master
2.95 ms ± 129 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # PR
25.8 ms ± 404 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)  # column-wise
```

So for this example, column-wise is actually faster than block-wise-with-slicing-blocks (first case where the dtypes are not the same). I don't know how that conclusion generalizes, though.

What is this column-wise approach doing? Compared to master, it operates on arrays instead of series and uses the faster "iterate through column arrays" helper function. So it gives this few-line change (instead of the 75-line complex `operate_blockwise` function):
```patch
@@ -325,8 +405,18 @@ def dispatch_to_series(left, right, func, str_rep=None, axis=None):
     elif isinstance(right, ABCDataFrame):
         assert right._indexed_same(left)
 
-        def column_op(a, b):
-            return {i: func(a.iloc[:, i], b.iloc[:, i]) for i in range(len(a.columns))}
+        arrays = []
+        for l, r in zip(left._iter_column_arrays(), right._iter_column_arrays()):
+            arrays.append(func(l, r))
+        return type(left)._from_arrays(
+            arrays, left.columns, left.index, verify_integrity=False
+        )
```

When dtypes are equal, block-wise is of course still a lot faster. I was actually surprised column-wise being faster when the dtypes were not equal, so maybe something is not yet fully optimized in `operate_blockwise` when dtypes are not equal. But, at the same time, there is also room for improvement in the column-wise approach, since `DataFrame._from_arrays` is giving a lot of overhead finding out which block types to create. 
And also, this is of course only a single example that was timed.  

So given this, I would personally think about this approach:

- Detect the case when we have the same block-layout for left and right dataframes (so same blocks with the same mgr_locs and the same dtypes), and in this case do blockwise (assuming that the `operate_blockwise` can be much simpler if we can just iterate through the blocks and don't have to slice up the blocks to get left and right matching)
- For other cases, keep it column-wise.

That still gives the big speed-up for the simple case of all the same dtypes, but avoids the complexity of the current `operate_blockwise` otherwise.



Correction on the above: to get tests passing, I think we still need the `array_op = get_array_op(func, str_rep=str_rep)`, and call `array_op` instead the `func` on the array in my patch. 
That gives a 2x slowdown for the example, unfortunately. But, with that it is still faster than blockwise in the case of non-equal dtypes, so I think my proposal at the end of the comment above is still valid. (it also makes me wonder what `array_op` all does ..).

> it also makes me wonder what array_op all does

my guess based on your question in #33252 is that the relevant behavior in array_op is casting dt64/td64 ndarrays to DTA/TDA.

> so maybe something is not yet fully optimized in operate_blockwise when dtypes are not equal

totally plausible.  ill do some profiling.

> Detect the case when we have the same block-layout for left and right dataframes

I think this adds another layer of complexity.  I'd rather have just one path that all frame-with-frame ops go through.
Looking at the original timings, the `df + df2` (mixed-dtype) case had way more blocks than intended.  I intended for `df2` blk_locs to be `[slice(0, 500), slice(500, 1500), slice(1500, 2000)]`, but instead got 1001 blocks.  Adding `_consolidate_inplace()` calls to the construction, the timing goes from

```
In [22]: %timeit df + df2                                                                                                                                                
82 ms ± 670 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [23]: df._consolidate_inplace()                                                                                                                                       
In [24]: df2._consolidate_inplace()

In [25]: %timeit df + df2                                                                                                                                                
7.16 ms ± 95.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

Ah, that clarifies why I saw so much slicing of blocks in the timings for this case! ;)

OK, so that makes it clear that block-wise is faster (as expected) when having only a few blocks (in the example 2 blocks for 2000 columns). 
(aside: should we then consolidate before doing arithmetic ops through this path?)

There is still another case where column-wise is faster than block-wise though. And that is with extension blocks, eg the nullable ints. Similar example but now using nullable int dtype:

```
arr = np.random.randint(0, 10000, (500, 2000)) 
df = pd.DataFrame(arr, dtype="Int64")  

%timeit df + df
```

gives me about 2x slower with block-wise compared to column-wise.

So for me, caring about getting performance for nullable dtypes better, that's still a reason to think about two different code paths.
> So for me, caring about getting performance for nullable dtypes better, that's still a reason to think about two different code paths.

1) sounds like a good follow-up
2) 2D EA here we come!

> sounds like a good follow-up

Not if you don't agree we can have multiple code paths. 
> Not if you don't agree we can have multiple code paths.

Pending seeing an implementation, i'm at "shouldn't", not "can't"

We would need to come up with a criterium to decide whether to use block vs column-wise. I was thinking of "is there any extension block", but that might be too aggressive, since not only the new nullable dtypes are using ExtensionBlocks.
> thinking of "is there any extension block",

"is all extension block"? 
"nblocks > const * ncolumns"?

My main concern with having a separate path is in testing.  One way to alleviate that concern would be to have a flag like we do for USE_NUMEXPR that can make us always go down one path or the other, then run arithmetic tests both ways.

But I do think that should be a separate discussion, and shouldn't hold this up.
To clarify why I want *some* discussion about it in this PR, is because if we don't want to go down that road, my opinion on this PR would be: I would go for simply always column-wise, although slower it much simpler in code (and still gives a speed-up compared to master)
To get a bit more an idea about performance trade-offs, I did some more timing experiments with a larger dataframe and varying ratio's between number of rows and number of columns (while keeping the total number of values in the dataframe constant). Similarly to the example above, it creates a dataframe with 2 blocks (50% float64, 50% float32). 
And I added an environment variable check in my [modified version](https://github.com/jorisvandenbossche/pandas/commits/perf-arith-frame-test2) of this branch to switch between column-wise and block-wise:

```
import os
import pandas as pd
import numpy as np

res = []

for n_rows, n_cols in [(10_000_000, 10), (1_000_000, 100), (100_000, 1000), (10_000, 10_000)]:
    
    # construct dataframe with 2 blocks
    arr1 = np.random.randn(n_rows, int(n_cols/2)).astype(dtype='float64') 
    arr2 = np.random.randn(n_rows, int(n_cols/2)).astype(dtype='float32') 
    df = pd.concat([pd.DataFrame(arr1), pd.DataFrame(arr2)], axis=1, ignore_index=True)
    
    # should already be the case, but just to be sure
    df._consolidate_inplace()

    for mode in ["blockwise", "columnwise"]:
        os.environ['OPS_MODE'] = mode

        df + df
        t = %timeit -o df + df

        res.append([n_rows, n_cols, mode, t.average, t.stdev])

res = pd.DataFrame(res, columns=["n_rows", "n_cols", "mode", "average", "stdev"])
```

I ran the above 5 times (with my laptop tuned with `pyperf system tune`), and this shows the averages (with bars indicating min/max of the 5 runs):

![image](https://user-images.githubusercontent.com/1020496/78879085-b40f5780-7a53-11ea-964d-f4c11018c7b5.png)


It's of course again a single experiment (one dataframe size (around half of GB), one set of dtypes (50% foat64, 50% float32), in my dev environment, etc), but it shows *some* timings. 
And note that this is for a dataframe with 2D blocks, not for one with 1D extension blocks.

I again find it surprising that column-wise is actually slightly faster for the first cases (even up to 1000 columns for this dataframe total size), as I would expect it to be slightly slower for long dataframes, and considerably slower for wide dataframes (would be good if someone can check I didn't do anything wrong). 

If I look at the profiles of a single run for the (1_000_000, 100) size, it seems that both block-wise and column-wise spend more or less the same time in the actual computation. Block-wise additionally spend around 50% of the time in slicing the blocks (in `_slice_take_blocks_ax0`, it still seems strange this to take that much of time). On the other hand, column-wise spends around 50% of its time in stacking the resulting arrays in re-consolidation (so showing that if we would *not* consolidate, the operation could still be faster).

Of course, block-wise is faster (and always will be) for really wide dataframes. But at least this seems to indicate that the idea of *"not doing block consolidation or not operating on 2D blocks will have a huge performance impact"* doesn't have to be the case for typically sized dataframes (and the 1000 columns in the above is already a lot), at least for the arithmetic operation tested here (selecting a single row is different story ;))
Thanks for doing these timings.  If it isn't too much trouble, can you do the same for a couple of cases where we dont exactly match the numpy behavior: e.g. floordiv and basically anything dt64/td64 dtypes
Profiling `df+df2` on this branch, im seeing something like 3/8 of the runtime being in `NDFrame._init_mgr`, so would be saved by #33357.  can't speak to how the column-wise is affected
>  cases where we dont exactly match the numpy behavior: e.g. floordiv 

I quickly checked floordiv. Didn't run comprehensive timings as done above, but generally it seems to give the same trend. There is now some time spent in `mask_zero_div_zero` (the non-numpy behaviour part, I suppose). But based on a single profile run, that is for both column as blockwise around 5-6 % of the total time.
Thanks for taking a look.

Just updated to take advantage of the new constructor fastpath.
Updated benchmarks:

```
arr = np.random.randn(10 ** 6).reshape(500, 2000).astype(np.float64)
df = pd.DataFrame(arr)
df[1000] = df[1000].astype(np.float32)
df.iloc[:, 1000:] = df.iloc[:, 1000:].astype(np.float32)

df2 = pd.DataFrame(arr)
df2[1000] = df2[1000].astype(np.int64)
df2.iloc[:, 500:1500] = df2.iloc[:, 500:1500].astype(np.int64)

df._consolidate_inplace()
df2._consolidate_inplace()

In [9]: %timeit df + df2                                                                                                                                                                                            
445 ms ± 24.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # <-- master
3.29 ms ± 30.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # <-- PR

In [10]: %timeit df + df                                                                                                                                                                                            
430 ms ± 5.72 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  # <-- master
2.97 ms ± 186 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # <-- PR
```
I repeated my benchmarks, but now for some other cases as well:

Something similar as the `df + df2` example from Brock's timings (https://github.com/pandas-dev/pandas/pull/32779#issuecomment-607018660), where `df2` has a block layout that does not fully align with `df`, and has int columns:

<details>

```
res = []

for n_rows, n_cols in [(10_000_000, 10), (1_000_000, 100), (100_000, 1000), (10_000, 10_000)]: 
    
    # construct dataframe with 2 blocks (float 64 / float 32)
    arr1 = np.random.randn(n_rows, int(n_cols/2)).astype(dtype='float64') 
    arr2 = np.random.randn(n_rows, int(n_cols/2)).astype(dtype='float32') 
    df = pd.concat([pd.DataFrame(arr1), pd.DataFrame(arr2)], axis=1, ignore_index=True)

    # second dataframe with 2 blocks (float64 / int64, but not aligning with blocks of df)
    factor = int(n_cols / 10)
    arr1 = np.random.randn(n_rows, 3*factor).astype(dtype='float64') 
    arr2 = np.random.randn(n_rows, 4*factor).astype(dtype='int64') 
    arr3 = np.random.randn(n_rows, 3*factor).astype(dtype='float64') 
    df2 = pd.concat([pd.DataFrame(arr1), pd.DataFrame(arr2), pd.DataFrame(arr3)], axis=1, ignore_index=True)

    # ensure both are consolidated
    df._consolidate_inplace()
    df2._consolidate_inplace()

    for mode in ["blockwise", "columnwise"]:
        os.environ['OPS_MODE'] = mode

        df + df2
        t = %timeit -o df + df2

        res.append([n_rows, n_cols, mode, t.average, t.stdev])

res = pd.DataFrame(res, columns=["n_rows", "n_cols", "mode", "average", "stdev"])
```

</details>

![image](https://user-images.githubusercontent.com/1020496/79043172-e5c42200-7bfd-11ea-8add-40aca9b35188.png)

Then with datetime-like data. `df` has two blocks, 50% timedelta, 50% datetime64 so the dataframe can be substracted from itself):

<details>

```
res = []

for n_rows, n_cols in [(10_000_000, 10), (1_000_000, 100), (100_000, 1000), (10_000, 10_000)]:
    
    # construct dataframe with 2 blocks
    arr1 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='datetime64[ns]') 
    arr2 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='timedelta64[ns]') 
    df = pd.concat([pd.DataFrame(arr1), pd.DataFrame(arr2)], axis=1, ignore_index=True)
    
    # should already be the case, but just to be sure
    df._consolidate_inplace()

    for mode in ["blockwise", "columnwise"]:
        os.environ['OPS_MODE'] = mode

        df - df
        t = %timeit -o df - df

        res.append([n_rows, n_cols, mode, t.average, t.stdev])

res = pd.DataFrame(res, columns=["n_rows", "n_cols", "mode", "average", "stdev"])
```
</details>

![image](https://user-images.githubusercontent.com/1020496/79043304-9b8f7080-7bfe-11ea-9e34-bc8ee8c50ade.png)

And something similar, but also comparing `datetime64[ns]` vs `datetime64[ns, UTC]`, which gives basically the same arithmetic operation, but where the first is a 2D consolidated block, while the second is a 1D extensionblock

<details>

```
res = []

for n_rows, n_cols in [(10_000_000, 10), (1_000_000, 100), (100_000, 1000), (10_000, 10_000)]: 
  
    for dtype in ["datetime64[ns]", "datetime64[ns, UTC]"]:
        # construct dataframe with 2 blocks dt+td
        arr1 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='datetime64[ns]') 
        arr2 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='timedelta64[ns]') 
        df = pd.concat([pd.DataFrame(arr1).astype(dtype), pd.DataFrame(arr2)], axis=1, ignore_index=True)
        
        # construct dataframe with 2 blocks td+dt (switched order, so we can do +)
        arr1 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='datetime64[ns]') 
        arr2 = np.random.randint(0, 1000000, (n_rows, int(n_cols/2))).astype(dtype='timedelta64[ns]') 
        df2 = pd.concat([pd.DataFrame(arr2), pd.DataFrame(arr1).astype(dtype)], axis=1, ignore_index=True)
        
        # ensure they are consolidated
        df._consolidate_inplace()
        df2._consolidate_inplace()

        for mode in ["blockwise", "columnwise"]:
            os.environ['OPS_MODE'] = mode

            df + df2
            t = %timeit -o -r 1 -n 1 df + df2

            res.append([n_rows, n_cols, dtype, mode, t.average, t.stdev])

res = pd.DataFrame(res, columns=["n_rows", "n_cols", "dtype", "mode", "average", "stdev"])
```
</details>

![image](https://user-images.githubusercontent.com/1020496/79043335-dabdc180-7bfe-11ea-9f82-848e205a7c93.png)


So the conclusions are all quite similar as for the initial timing I posted above (https://github.com/pandas-dev/pandas/pull/32779#issuecomment-611427183). It depends a bit on the exact dataframe layout and operation whether column-wise is slightly faster or not for the cases up to 1000 columns (with the given dataframe size), and column-wise becomes slower for larger number of columns.  
Also interesting to note is that tz-aware (1D block) is in some cases faster as tz-naive (2D block) datetimes (but, the difference is not very big), except for the widest dataframe (where again the overhead of doing it column-wise becomes significant). 
Just looking at tz-aware (a dataframe with only 1D blocks), column-wise is always slightly faster as block-wise (as expected I think, as the blocks give some overhead).

So this gives a different picture as the timing that Brock posted just above. Of course, it's not comparing the same (Brock is comparing against master, I am comparing with an optimized column-wise version). 
But an important factor that also explains the difference is that I used larger dataframes (and different n_rows/n_cols ratio). With the size of dataframe (570-760MB) and the time the operation takes (around 0.5-1.5s for most cases), the dominant factor seems to be the actual arithmetic operation (up to 1000 columns). And as long as the columns have a certain length, it doesn't seem to matter much if this operation is done on the 2D array or per 1D array. It's only when going to many columns (or more columns than rows), that the overhead of doing it per column becomes significant (an overhead for which there is also still a lot of room to improve this).

To do the same timing as Brock this with column-wise:

```
In [8]: %timeit df + df
2.83 ms ± 59.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # block-wise
55.3 ms ± 532 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)  # column-wise

In [10]: %timeit df + df2
3.79 ms ± 57.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  # block-wise
58.8 ms ± 337 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)  # column-wise
```

So column-wise is significantly slower, although the difference is not that big as compared to master. 
But speaking in absolute terms, for this case column-wise is around 50ms slower, while in one of the cases above column-wise is a few 100ms faster. It's of course all comparing apples with oranges, with many factors that influence it, but it does pose the question which case we want to optimize for. 

@jorisvandenbossche are there any cases you've identified where block-wise is slower than master?
I am not comparing with master. All the comparisons I show above are with the optimized version of column-wise (iterating over the columns as arrays instead of as series, [branch](https://github.com/jorisvandenbossche/pandas/commits/perf-arith-frame-test)).
> I am not comparing with master

That's what I thought.  If this is unamiguously better than master, we should move forward, as discussed on the call.