You can use the hmean() function of scipy.stats. This function calculates the harmonic mean along the specified axis. i.e. n / (1/x1 + 1/x2 + … + 1/xn)
>>> from scipy.stats import hmean
>>> a=np.array([10,20,15])
>>> hmean(a)
13.846153846153845
>>>
You can also apply this function to a 2D array.
>>> a=np.array([[1,2],[11,12],[21,22]])
>>> a
array([[ 1, 2],
[11, 12],
[21, 22]])
>>> hmean(a,axis=0)
array([2.63498099, 4.77108434])
>>> hmean(a,axis=1)
array([ 1.33333333, 11.47826087, 21.48837209])