ENH: Return a frame with a MultiIndex when doing math on inner indices
#### Is your feature request related to a problem?

Given a datafame with a `MultiIndex` axis, it would be nice to be able to perform a mathematical operation based on the outer axis (or axes) and obtain a result that maintained the number of dimensions in the index. For example:

```python
cols = [(x, y) for p in ["r", "s"] for n in ["a","b","c"]]
dct = {col: [0, 1, 2] for col in cols}
df = pd.DataFrame(dct)
#   r        s      
#   a  b  c  a  b  c
#0  0  0  0  0  0  0
#1  1  1  1  1  1  1
#2  2  2  2  2  2  2

type(df.columns)
# <class 'pandas.core.indexes.multi.MultiIndex'>

diff = df["r"] - df["s"]
type(diff.columns)
# <class 'pandas.core.indexes.base.Index'>
```
This would allow for expressions like the following:
```python
df["diff"] = df["r"] - df["s"]  # THIS DOES NOT WORK
```
At the moment, the only way to add a new set of columns for `"diff"` is to create them one at a time.

#### Describe the solution you'd like

I know it's tricky, but it would be nice if "inner" data frames were aware that they were part of a `MultiIndex`.

#### API breaking implications

Note sure...

#### Describe alternatives you've considered

As stated before, the goal can be achieved by iterating over columns.

#### Additional context

I was not able to find duplicate issues, but feel free to close if there is one.