Now I have the following implementation:
class Bilinear(nn.Module):
    '''
    See: torch/csrc/api/src/nn/modules/linear.cpp
    '''

    def __init__(self, input_size, hidden_size):
        super(Bilinear, self).__init__()
        stdv = input_size**-0.5

        W_a = torch.empty(input_size, input_size)
        nn.init.uniform_(W_a, -stdv, stdv)
        self.W_a = nn.Parameter(W_a)

        W_b = torch.empty(hidden_size, input_size)
        nn.init.uniform_(W_b, -stdv, stdv)
        self.W_b = nn.Parameter(W_b)
        self.b = nn.Parameter(torch.empty(input_size))

    def forward(self, x, h):
        return self.W_a.t().matmul(x) + self.b.t() + self.W_b.t().matmul(h)

which works as expected. Similar to the torch implementation of Linear but with bias being initialized with zeros.