I am trying to write a custom function that computes the dominant eigenvector and its derivative of a symmetric matrix using Eq. (68) in the matrix cookbook 41 and numpy.
Here is my code:
import torch
from torch.autograd import Function, Variable, gradcheck
import numpy as np


def to_tensor(x):
    return torch.from_numpy(x).float()

class TopEigen(Function):

    @staticmethod
    def forward(ctx, matrix):
        cov = matrix.numpy()
        lams, vecs = np.linalg.eigh(cov)
        ctx.intrm = (lams, vecs, cov)
        return to_tensor(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) - cov)
        grad_np = np.dot(np.dot(output.T, pinv).T, vecs[None, :, -1])
        return Variable(to_tensor(grad_np))


# Testing
topeig = TopEigen.apply

p, q = 5, 3
in_tensor = Variable(torch.rand(p, q), requires_grad=True)
cov_in = torch.mm(in_tensor.t(), in_tensor)

out = topeig(cov_in).mean()
out.backward(retain_graph=True)

test = gradcheck(topeig, (cov_in, ), eps=1e-6, atol=1e-4)
print(test)

When I run this script, I get the following error message:
RuntimeError: for output no. 0,
 numerical:(
-0.0680  0.0404  0.0532
 0.0000  0.0000  0.0000
 0.0000  0.0000  0.0000
-0.0015 -0.1016  0.0791
 0.0249 -0.0779  0.0285
 0.0000  0.0000  0.0000
-0.0019  0.0790 -0.0578
 0.0655 -0.0740 -0.0246
 0.0431  0.0375 -0.0817
[torch.FloatTensor of size 9x3]
,)
analytical:(
    0     0     0
    0     0     0
    0     0     0
    0     0     0
    0     0     0
    0     0     0
    0     0     0
    0     0     0
    0     0     0
[torch.FloatTensor of size 9x3]
,)

Why is the analytical output all 0?  I check the code with double() tensors but the result was the same.
Can anyone help me with spotting the error in my code?