self.matrix will be a plain tensor, not an nn.Parameter, if you wrap just single nn.Parameters in another tensor.
You can check if by calling print(rotation_matrix.matrix).
If you create the whole matrix as an nn.Parameter, the code should work:
class Rotation_Matrix(nn.Module):
    def __init__(self):
        super(Rotation_Matrix, self).__init__()
        self.a = torch.tensor(1.) # init value of Rx cosine
        self.b = torch.tensor(1.) # init value of Ry cosine
        self.c = torch.tensor(1.) # init value of Rz cosine
        self.matrix = self.get_rotation_matrix()
        
    def get_rotation_matrix(self):
        Rx = torch.tensor([[1., 0., 0.],
                          [0., self.a, -torch.sqrt(1 - torch.pow(self.a,2))],
                          [0., torch.sqrt(1 - torch.pow(self.a,2)), self.a]])

        Ry = torch.tensor([[self.b, 0., torch.sqrt(1 - torch.pow(self.b,2))],
                          [0., 1., 0.],
                          [-torch.sqrt(1 - torch.pow(self.b,2)), 0, self.b]])

        Rz = torch.tensor([[self.c, -torch.sqrt(1 - torch.pow(self.c,2)), 0.],
                          [torch.sqrt(1 - torch.pow(self.c,2)), self.c, 0.],
                          [0., 0., 1.]])

        return nn.Parameter(torch.mm(Rx, torch.mm(Ry, Rz)))
    
    def forward(self, x):
        return torch.mm(x, self.matrix)
    
model = nn.Linear(3, 1)
model.eval()

rotation_matrix = Rotation_Matrix()
optimizer = torch.optim.Adam(rotation_matrix.parameters(), lr=1.0)

rotation_matrix.train()

x = torch.randn(1, 3)
target = torch.randn(1, 1)
criterion = nn.MSELoss()

for i in range(10):
    optimizer.zero_grad()

    output = rotation_matrix(x)
    output = model(output)    
    loss = criterion(output, target)

    loss.backward()
    print(rotation_matrix.matrix.grad)
    optimizer.step()