Below are a few points that should help explain why the PyTorch snippet can reach around a 15% error rate (even with single‐sample updates and manually applied parameter updates), whereas the custom NumPy-based SGD code might only reach around 30%. I’ll also give some improvements you can try in your own implementation to get closer to the PyTorch results—even without batching or an official “optimizer.”

1. What “special tools” does PyTorch use internally?

Even when we manually update PyTorch parameters in a loop, PyTorch’s core autograd and nn.BCELoss (or F.binary_cross_entropy) routines offer some numerical and implementation advantages:

1. Stable log-sigmoid computations
   PyTorch’s implementation of torch.sigmoid and the internal log-based functions in nn.BCELoss use numerically stable kernels. This helps avoid exploding or vanishing values when your logits become large (positively or negatively). In NumPy, using 1 / (1 + exp(-z)) can blow up if z is large negative and can produce floating-point overflow if z is large positive.

   • Even though you are using np.clip(..., 1e-9, 1 - 1e-9) in the NumPy code, you might still be hitting less stable regions than PyTorch’s carefully implemented versions (which typically rely on $\log(1 + \exp(x))$ forms internally).

2. Consistent and correct gradient flow
   By default, PyTorch’s linear layer (nn.Linear) uses random initialization for weights (unless you explicitly zero them out). In your snippet, you did:

   python
   nn.init.constant_(self.linear.weight, 0)
   nn.init.constant_(self.linear.bias, 0)

   which indeed zeros out parameters. Although this can cause some degeneracy in certain tasks, it’s also possible that for this particular dataset, the zero initialization is not the end of the world. Still, PyTorch automatically tracks correct gradients and updates them reliably (no mismatch with the objective or indexing errors).

3. Accurate summation and in-place updates
   PyTorch’s internal usage of GPU/CPU vectorization, fused ops, etc., can lead to slight differences in the effective floating‐point arithmetic compared to custom Python loops that do single-sample updates in NumPy. Sometimes these differences help with numerical stability.

4. Loss layering
   criterion = nn.BCELoss() under the hood does:

   $$-\frac{1}{N}\sum_{i} \left[y_i \log(\hat{y}_i) + (1-y_i)\log(1-\hat{y}_i)\right]$$

   but also includes the stable $\log$ trick to handle edge cases. Even if we do not see it in a code snippet, PyTorch’s built-in functions handle a lot of those edge conditions for free.

None of these alone is “magic.” But together, they can give PyTorch an edge in both stability and speed of convergence, even if you do parameter updates manually.

2. Potential issues in your custom NumPy implementation

Below are key places where the custom code often deviates from “best practice” and might cause higher error:

1. Regularization mismatch

   • You compute $\text{L2 loss} = (l2\_penalty / 2) * \sum_{i>0} W_i^2$ (excluding the bias in the summation), which is good.
   • But in the gradient, you do:

     python
     gradient = (sigmoid - y) * X + self.l2_penalty * W

     This does penalize the bias term because you multiply the entire vector $W$ by l2_penalty. The code snippet’s comment says:

     python
     # gradient = (sigmoid - y) * X + self.l2_penalty * np.r_[[0], W[1:]]

     If you really want to exclude the bias from L2, you must similarly skip index 0 in the gradient step. Right now, you’re adding an extra penalty on the bias that the “theoretical” formula excludes. That mismatch can degrade your final accuracy because the bias might be shifting more slowly or incorrectly.

2. Initialization

   • You initialize W = np.zeros(X.shape[1]). Zero initialization for logistic regression is usually less harmful than it is for deep nets, but it can still slow you down, especially if your dataset is tricky. You might experiment with a small random normal or uniform initialization.

3. Momentum at 0

   • You’ve set momentum = 0.0. That basically disables momentum. If you want to see better gains, try a typical momentum value (e.g. 0.9). With single-sample SGD, momentum can help smooth out the updates.

4. Learning rate and number of epochs

   • Sometimes 0.001 is too small for single-sample updates (or maybe too large for your dataset—there’s no universal best). Try stepping up the epochs to 50 or 100 to see if it converges. Also test different learning rates, like 0.01 or 0.0001, to see if you get a better plateau.

5. Data Shuffling

   • You commented out np.random.shuffle(indices). In single-sample SGD, randomizing the order of samples each epoch is crucial to avoiding cyclical or stuck-in-local-minima behavior. If you do not reshuffle, you might traverse the samples in the same order each epoch, which can hamper your convergence.
   • In your code, you wrote:

     python
     indices = np.arange(X.shape[0])
     # np.random.shuffle(indices)

     Try un-commenting that shuffle line. Single-sample SGD typically benefits from a fresh shuffle every epoch.

6. Edge cases in the log and numerical stability

   • You are using:

     python
     np.clip(sigmoid, 1e-9, 1 - 1e-9)

     That’s good for preventing NaNs, but ensure that everything is consistent. Because you compute the gradient from (sigmoid - y) * X (which does not contain that same clipping logic), large logits might still cause weird curvature. In PyTorch, the underlying code for BCELoss is more gracefully handling that log and gradient.

3. How to improve your results without batching or a standard optimizer

Even if you want to keep single-sample updates and avoid PyTorch’s built-in optimizers (torch.optim.Adam, etc.), you can still do quite a bit to bring performance and accuracy closer to PyTorch’s:

1. Exclude the bias from L2 in both the loss and gradient
   Either penalize everything consistently or explicitly skip the bias. This mismatch is a common silent bug. For example:

   python
   gradient = (sigmoid - y) * X
   # skip index 0 for bias penalty
   gradient[1:] += self.l2_penalty * W[1:]

2. Initialize with small random values instead of zeros

   python
   W = 0.01 * np.random.randn(X.shape[1])

3. Shuffle your data each epoch

   python
   indices = np.arange(X.shape[0])
   np.random.shuffle(indices)
   for i in indices:
       ...

4. Increase the number of epochs
   With a single sample update (and especially with a small learning rate like 0.001), 10 epochs might not be enough. You can easily go to 50–100 epochs or more with logistic regression.

5. Try a momentum value around 0.9
   This can help your updates escape small local traps and also reduce variance in single-sample updates:

   python
   momentum = 0.9

6. Adjust learning rate
   If your error does not improve after many epochs, try a bigger or smaller $\alpha$. Start with 0.001, try 0.01, etc.

7. Check for data preprocessing
   Sometimes PyTorch’s example code normalizes or standardizes input features. If you are feeding raw data in one case and standardized data in another, that alone could explain the difference in final accuracy.

8. Verify your label conversion
   Make sure your final “0/1” labels match exactly in both the PyTorch version and the custom version. Typos in the sign flipping or in the thresholding code can creep in.

9. Test multiple seeds
   For a fair comparison, set the same random seed for everything—both PyTorch and NumPy. Even minor differences can lead to large divergences in single‐sample training.

Putting it all together

In short, there’s no single “secret sauce” in PyTorch that automatically grants ~15% error for logistic regression while your code hovers at ~30%. The main factors are:

• Numerically stable logistic/BCELoss implementation
• Slight differences in initialization
• Accurate gradient flow
• Shuffling
• Consistent regularization

By carefully aligning these practices in your custom code (particularly the gradient regularization, the initialization, the data shuffle, and enough epochs), you should see your error rate drop significantly—very likely closer to PyTorch’s results, even without batching or a more advanced optimizer.