Question: Lets say I have a 4D array with shape `(1,2,3,3)`:

```
test = np.array([[[[11,27,33],[45,58,96],[77,85,93]],[[55,27,39],[46,51,62],[73,86,98]]]])

```

Whats the most efficient way of standardizing/calculating z scores for a 2D subset? For example, `test[0][0]` looks like this:

```
array([[11, 27, 33],
       [45, 58, 96],
       [77, 85, 93]])

```

There are 2 dimensions here, but I want to calculate the mean and standard deviation across both dimensions, and use those values to standardize each value in these 2 dimensions.

I can do it manually like this:

`(test[0][0] - np.mean(test[0][0])) / np.std(test[0][0])`

Which correctly gives:

```
array([[-1.61593336, -1.06970236, -0.86486574],
       [-0.45519249, -0.01137981,  1.2859188 ],
       [ 0.63726949,  0.91038499,  1.18350049]])

```

However, this would require me to iterate over the first 2 dimensions of the 4D array which would take too long given the size of my actual data

I see that `scipy` has a `zscore` function but that only works in 1 dimension at a time: `scipy.stats.zscore(test, axis=3)` and haven't been able to find a simple implementation that standardizes across a 2D array

Answer: 
**Approach #1 :** You could make use of using `np.mean` and `np.std` over multiple axes (in this case over the last two axes) with `axis=(2,3)` and keep their number of dims same with `keepdims=1` so that the later subtraction and division operations are broadcastable.

Thus, a vectorized implementation would be -

```
(test - test.mean(axis=(2,3),keepdims=1)) / test.std(axis=(2,3),keepdims=1)

```

**Approach #2 :** Alternative approach using the definition of `std` that would re-use the average calculations -

```
m = (test - test.mean(axis=(2,3),keepdims=1))
s = np.sqrt((np.abs(m)**2).mean(axis=(2,3),keepdims=1))
out = m/s

```

**Approach #3 :** For larger datasets, you might want to use `numexpr` module that does those summing/averaging operations quite efficiently -

```
import numexpr as ne

d0,d1 = test.shape[-2:]
m = (test - test.mean(axis=(2,3),keepdims=1))
m1 = m.reshape(-1,d0*d1)
s = np.sqrt(ne.evaluate('sum(abs(m1)**2,1)')/(d0*d1))
out = m/s[:,None,None]

```

Based on [`this post`](https://stackoverflow.com/a/40578932/3293881), we could replace those `division by s` with `1.0/s` and then multiply it with `m` for further performance boost. This would be applicable across all above mentioned three approaches.
To do normalization in numpy, just make broadcasting match.

```
def normalize_nchw(inp):
    EPS = 1e-6
    means = np.mean(inp, axis=(2,3)).expand_dims(-1).expand_dims(-1)
    inp -= means
    vars = EPS + np.mean(inp*inp, axis=(2,3)).expand_dims(-1).expand_dims(-1)
    inp *= (1./np.sqrt(vars))

```

Side note: if you are doing this for CNN, a better idea is to use batch normalization, which is built into many frameworks.