Pandas groupby throws an error when using sum()

I am trying to calculate the between cluster scatter matrix. In order to do that, for each cluster (named "group" in the example below), I need to perform an operation which results in a matrix and subsequently perform an element-wise addition of the matrices from each cluster.

To do this I try the following:

import pandas as pd
import numpy as np

df = pd.DataFrame({'group': [1, 2, 1, 0, 0, 0, 1, 2],
                   'A': [1.5, 0.5, 2.5, 0.5, 1.5, 0.5, 1.5, 0.5],
                   'B': [3.5, 2.5, 3.5, 2.5, 3.5, 2.5, 3.5, 2.5]})
features = list(df.columns)
features.remove('group')


def g(x, mu):
    y = np.array([np.mean(x) - mu])
    print((y.T @ y)*len(x))
    print("")
    return (y.T @ y)*len(x)


m = len(df.index)
mu = df.groupby('group')[features].apply(lambda x: (np.multiply(x.count(), np.mean(x)))/m).sum()
print("mu:")
print(mu)

Sb = df.groupby('group')[features].apply(g, mu=(mu)).sum()

This example throws a TypeError: Series.name must be a hashable type error on the last line. The print statement in the g function shows the result as expected, see below, so I believe the error is due to the .sum() operation.

[[0.2553 0.1458]
 [0.1458 0.0833]]

[[1.5052 1.0625]
 [1.0625 0.75  ]]

[[0.7912 0.625 ]
 [0.625 0.5    ]]

The result I was expecting by adding the .sum() operation was the element-wise addition of the three matrices above.

The expected output is:

[[2.5416 1.8333]
 [1.8333 1.3333]]

Any ideas why this is giving me an error and what I can do to correct it?

Update 1:
Using:

Sb = df.groupby('group').apply(g, mu=(mu)).sum()

instead of

Sb = df.groupby('group')[features].apply(g, mu=(mu)).sum()

gives the correct matrix, padded with nans. Why does features cause an error?

Have you tried this ?

sb=df.groupby('group').apply(g, mu=(mu)).sum()

it gives the following result:

[[2.54166667 1.83333333        nan]
 [1.83333333 1.33333333        nan]
 [       nan        nan        nan]]

Is it what you want ?

You still have to deal with the nans though

Edit to answer your comments:

To answer you problem in the comments you could change your function as below:

def g(x, mu):
    x=x[["A","B"]] #or x=x[features]
    y = np.array([np.mean(x) - mu])
    print((y.T @ y)*(len(x)))
    print("")
    return (y.T @ y)*(len(x))

and then:

sb=df.groupby(['group']).apply(g, mu=(mu)).sum()
print(sb)

which gives:

[[2.54166667 1.83333333]
 [1.83333333 1.33333333]]