A is sparse, so A-transpose is too. (Supposedly pytorch
knows this, but maybe not.)

I think you have find the culprit here: something weird happens to sparse matrices when you transpose them.
I have coded what the backpropagation algorithm is supposed to do, and I timed it, as follows:
tic = time.time()
with torch.no_grad():
    RF_pred = torch.sparse.mm(FO.A_torch, y)
    toc1 = time.time()
    print('Matrix-vector product (sec)' + str(toc1-tic))

    residual = RF.mask*(RF_pred - RF.x_data)
    toc2 = time.time()
    print('Residual calculation + masking (sec)' + str(toc2-toc1))

    At = torch.t(FO.A_torch)
    toc3 = time.time()
    print('Transpose calculation (sec)' + str(toc3 - toc2))

    grad_y = torch.sparse.mm(At, residual)
    toc4 = time.time()
    print('Matrix (A transpose)-vector product (sec)' + str(toc4 - toc3))

    y.grad = grad_y

The result I get is the following
Matrix-vector product (sec)0.3793952465057373
Residual calculation + masking (sec)0.0005414485931396484
Transpose calculation (sec)1.4321186542510986
Matrix (A transpose)-vector product (sec)24.383739471435547

First of all: taking the transpose takes waaaaaay too much time than it should.
Then, calculating the matrix product using the transposed matrix is what takes forever.
Maybe we should raise an issue on github about it?