In case someone is trying to look for a solution:
def cov(m):
    # m = m.type(torch.double)  # uncomment this line if desired
    fact = 1.0 / (m.shape[-1] - 1)  # 1 / N
    m -= torch.mean(m, dim=(1, 2), keepdim=True)
    mt = torch.transpose(m, 1, 2)  # if complex: mt = m.t().conj()
    return fact * m.matmul(mt).squeeze()


def compute_rank_correlation(x, y):
    x, y = rankmin(x), rankmin(y)
    return corrcoef(x, y)


def corrcoef(x, y):
    batch_size = x.shape[0]
    x = torch.stack((x, y), 1)
    # calculate covariance matrix of rows
    c = cov(x)
    # normalize covariance matrix
    d = torch.diagonal(c, dim1=1, dim2=2)
    stddev = torch.pow(d, 0.5)
    stddev = stddev.repeat(1, 2).view(batch_size, 2, 2)
    c = c.div(stddev)
    c = c.div(torch.transpose(stddev, 1, 2))
    return c[:, 1, 0]

There is small difference between the above implementation and scipy’s spearmanr implementation:

This operate on batches, that is, it expects input to be of dimension (N, x) where N is the batch size for both parameters x and y.
scipy uses rank average while this uses rank min to compute the ranking before corrcoef. There is 10^-3 disagreement and although I suspect that this is the cause, this could possibly be.

Please feel free to share your improvements on this as it could be useful for other people and me