One can compute the full jacobian matrix by doing backward multiple times through the graph. For example, given a  f  that maps R^n to R^m, then we can compute the jacobian via the following:
def unit_vectors(length):
    result = []
    for i in range(0, length):
        x = torch.zeros(length)
        x[i] = 1
        result.append(x)
    return result

x = torch.randn(3, requires_grad=True)
y = f(x)
result = [torch.autograd.grad(outputs=[y], inputs=[x], grad_outputs=[unit], retain_graph=True)[0] for unit in unit_vectors(y.size(0))]
jacobian = torch.stack(result, dim=0)

Following up from the questions in https://github.com/pytorch/pytorch/issues/10223 6,

@zou3519 also, where do you use the  unit_vectors()  function?

unit_vectors is used in the second-to-last line.

Why are multiple grad calls necessary? Shouldn’t one be sufficient?

Given a scalar output f_1 and an input vector x, autograd gives the vector (\frac{df}{dx_1}, \frac{df}{dx_2}, ..., \frac{df}{dx_n}).
For the full matrix Jacobian we have an output (f_1, f_2, ..., f_m) and we want all \frac{df_i}{dx_j}. One way to get this is to run autograd once on each f_i, which is what is being done in the example above.
.
There might be a way to compress all of this computation into one backward pass; were you thinking of something specific?