@albanD, @tom : I am still not able to get this function right.  Here is my new code:
import torch
from torch.autograd import Function, gradcheck
from torch import from_numpy
import numpy as np


class TopEigen(Function):
    """Copmpute the top eigenvector of a matrix."""

    @staticmethod
    def forward(ctx, matrix):
        cov = matrix.detach().numpy()
        lams, vecs = np.linalg.eigh(cov)
        ctx.intrm = lams, vecs
        return from_numpy(vecs[:, -1, None])

    @staticmethod
    def backward(ctx, grad_output):
        lams, vecs, cov = ctx.intrm
        output = grad_output.data.numpy()
        pinv = np.linalg.pinv(lams[-1]*np.eye(cov.shape[0]) - cov)
        return from_numpy(np.dot(np.dot(output.T, pinv).T, vecs[None, :, -1]))


# Create a random symmetric matrix
p, q = 5, 3
torch.manual_seed(0)
in_tensor = torch.rand(p, q, dtype=torch.float64, requires_grad=True)
cov_in = torch.mm(in_tensor.t(), in_tensor)

# Testing
topeig = TopEigen.apply
out = topeig(cov_in).mean()
test = gradcheck(topeig, (cov_in, ), eps=1e-10, atol=1e-3)
print(test)

Here is the output in pytorch 0.4.1:
Traceback (most recent call last):
  File "eigen.py", line 41, in <module>
    test = gradcheck(topeig, (cov_in, ), eps=1e-10, atol=1e-3)
  File ".../anaconda3/lib/python3.6/site-packages/torch/autograd/gradcheck.py", line 214, in gradcheck
    'numerical:%s\nanalytical:%s\n' % (i, j, n, a))
  File ".../anaconda3/lib/python3.6/site-packages/torch/autograd/gradcheck.py", line 194, in fail_test
    raise RuntimeError(msg)
RuntimeError: Jacobian mismatch for output 0 with respect to input 0,
numerical:tensor([[-0.0822,  0.0487,  0.0460],
        [ 0.0000,  0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000],
        [-0.0256, -0.0462,  0.0814],
        [ 0.0441, -0.0817,  0.0360],
        [ 0.0000,  0.0000,  0.0000],
        [-0.0219,  0.0801, -0.0611],
        [ 0.0819, -0.0387, -0.0568],
        [ 0.0381,  0.0329, -0.0820]], dtype=torch.float64)
analytical:tensor([[-0.0822,  0.0487,  0.0460],
        [-0.0743,  0.0441,  0.0416],
        [-0.0680,  0.0403,  0.0381],
        [ 0.0487, -0.0903,  0.0398],
        [ 0.0441, -0.0817,  0.0360],
        [ 0.0403, -0.0747,  0.0329],
        [ 0.0460,  0.0398, -0.0992],
        [ 0.0416,  0.0360, -0.0897],
        [ 0.0381,  0.0329, -0.0820]], dtype=torch.float64)

Interestingly, some rows are similar.  Would you please spot my mistake?