I checked your code and compared the gradients as you did here 21 and here 6 and they are the same up to the floating precision.
Unfortunately, I don’t know, what fminunc 8 is doing exactly.
Are you sure it’s SGD? Maybe it calls some iterations internally on its own.
I replaced the SGD optimizer with LBFGS and it you’ll get a loss of ~0.2 in approx. 200 epochs.
model = LogisticRegression(input_size, num_class)
criterion = nn.BCELoss(reduce=True)
optimizer = torch.optim.LBFGS(model.parameters(), lr=1e-3)

loss_history = []

def closure():
    optimizer.zero_grad()
    output = model(X)
    loss = criterion(output, Y)
    loss.backward()
    return loss


for epoch in range(400):
    loss = optimizer.step(closure)

    if (epoch + 1) % 10 == 0:
        print('Epoch: [%d/%d], Loss: %.4f' %(epoch+1, total_epoch, loss))
    loss_history.append(loss)

EDIT:
Also I initialized the weights to zeros to compare both codes:
class LogisticRegression(nn.Module):
    def __init__(self, input_size, num_class):
        super(LogisticRegression, self).__init__()
        self.linear = nn.Linear(input_size, num_class)
        self.initialize()
        
    def initialize(self):
        self.linear.weight.data.zero_()
        self.linear.bias.data.zero_()

    def forward(self, x):
        out = self.linear(x)
        return F.sigmoid(out)