Problem solved. In Method2, when do
z2 = (torch.sum(x2 ** 2, dim=1) + torch.sum(y2 ** 2, dim=0)) - 2 * cos_sim_matrix

The broadcast is triggered by -2 * cos_sim_matrix, which is wrong. The correct broadcast should be
z2 = (torch.sum(x2 ** 2, dim=1, keepdim=True) + torch.sum(y2 ** 2, dim=0, keepdim=True)) - 2 * cos_sim_matrix