我在用C语言编写的逻辑回归代码中遇到了一些困难，虽然它似乎可以在较小的半随机数据集上工作，但它停止工作了（例如，分配属于类1的适当概率）在我传递43，500个观测的点附近（通过调整创建的观测的数目来确定。当创建代码中使用的150个特征时，我创建了前两个函数，作为观测值的函数，所以我不确定这是否是问题所在，尽管我使用了双精度，也许代码中的某个地方有溢出？
下面的代码应该是独立的;它将生成m= 50，000个观测值，其中包含n=150个要素。将m设置为小于43，500时，应返回“Percent class 1：0.250000”，设置为44，000或以上将返回“百分比类1：0.000000”，而不管max_iter（我们对m个观测进行采样的次数）被设置为多少。
如果类为0（前75%的观测），则第一个要素设置为1.0除以观测总数，否则设置为观测的索引除以观测总数。
第二个特征是指数除以观测总数。
所有其他特征都是随机的。
逻辑回归旨在使用随机梯度下降、随机选择观察指标、使用当前权重计算具有预测y的损失的梯度、以及利用梯度和学习率（eta）更新权重。
对Python和NumPy使用相同的初始化，即使超过50，000个观测值，我仍然可以得到正确的结果。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

// Compute z = w * x + b
double dlc( int n, double *X, double *coef, double intercept )
{
    double y_pred = intercept;
    for (int i = 0; i < n; i++)
    {
        y_pred += X[i] * coef[i];
    }
    return y_pred;
}

// Compute y_hat = 1 / (1 + e^(-z))
double sigmoid( int n, double alpha, double *X, double *coef, double beta, double intercept )
{
    double y_pred;
    y_pred = dlc(n, X, coef, intercept);
    y_pred = 1.0 / (1.0 + exp(-y_pred));

    return y_pred;
}

// Stochastic gradient descent
void sgd( int m, int n, double *X, double *y, double *coef, double *intercept, double eta, int max_iter, int fit_intercept, int random_seed )
{
    double *gradient_coef, *X_i;
    double y_i, y_pred, resid;
    int idx;

    double gradient_intercept = 0.0, alpha = 1.0, beta = 1.0;

    X_i = (double *) malloc (n * sizeof(double));
    gradient_coef = (double *) malloc (n * sizeof(double));

    for ( int i = 0; i < n; i++ )
    {
        coef[i] = 0.0;
        gradient_coef[i] = 0.0;
    }
    *intercept = 0.0;

    srand(random_seed);
    
    for ( int epoch = 0; epoch < max_iter; epoch++ )
    {
        for ( int run = 0; run < m; run++ )
        {
            // Randomly sample an observation
            idx = rand() % m;
            for ( int i = 0; i < n; i++ )
            {
                X_i[i] = X[n*idx+i];
            }
            y_i = y[idx];
            // Compute y_hat
            y_pred = sigmoid( n, alpha, X_i, coef, beta, *intercept );
            resid = -(y_i - y_pred);
            // Compute gradients and adjust weights
            for (int i = 0; i < n; i++)
            {
                gradient_coef[i] = X_i[i] * resid;
                coef[i] -= eta * gradient_coef[i];
            }
            if ( fit_intercept == 1 )
            {
                *intercept -= eta * resid;
            }
        }
    }
}

int main(void)
{
    double *X, *y, *coef, *y_pred;
    double intercept;
    double eta = 0.05;
    double alpha = 1.0, beta = 1.0;
    long m = 50000;
    long n = 150;
    int max_iter = 20;

    long class_0 = (long)(3.0 / 4.0 * (double)m);
    double pct_class_1 = 0.0;

    clock_t test_start;
    clock_t test_end;
    double test_time;

    printf("Constructing variables...\n");
    X = (double *) malloc (m * n * sizeof(double));
    y = (double *) malloc (m * sizeof(double));
    y_pred = (double *) malloc (m * sizeof(double));
    coef = (double *) malloc (n * sizeof(double));

    // Initialize classes
    for (int i = 0; i < m; i++)
    {
        if (i < class_0)
        {
            y[i] = 0.0;
        }
        else
        {
            y[i] = 1.0;
        }
    }

    // Initialize observation features
    for (int i = 0; i < m; i++)
    {
        if (i < class_0)
        {
            X[n*i] = 1.0 / (double)m;
        }
        else
        {
            X[n*i] = (double)i / (double)m;
        }
        X[n*i + 1] = (double)i / (double)m;
        for (int j = 2; j < n; j++)
        {
            X[n*i + j] = (double)(rand() % 100) / 100.0;
        }
    }

    // Fit weights
    printf("Running SGD...\n");
    test_start = clock();
    sgd( m, n, X, y, coef, &intercept, eta, max_iter, 1, 42 );
    test_end = clock();
    test_time = (double)(test_end - test_start) / CLOCKS_PER_SEC;
    printf("Time taken: %f\n", test_time);

    // Compute y_hat and share of observations predicted as class 1
    printf("Making predictions...\n");
    for ( int i = 0; i < m; i++ )
    {
        y_pred[i] = sigmoid( n, alpha, &X[i*n], coef, beta, intercept );
    }

    printf("Printing results...\n");
    for ( int i = 0; i < m; i++ )
    {
        //printf("%f\n", y_pred[i]);
        if (y_pred[i] > 0.5)
        {
            pct_class_1 += 1.0;
        }
        // Troubleshooting print
        if (i < 10 || i > m - 10)
        {
            printf("%g\n", y_pred[i]);
        }
    }
    printf("Percent class 1: %f", pct_class_1 / (double)m);

    return 0;
}

作为参考，下面是我的（大概）等价Python代码，它在超过50，000次观察时返回正确的类百分比：

import numpy as np
import time

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

class LogisticRegressor:
    def __init__(self, eta, init_runs, fit_intercept=True):
        self.eta = eta
        self.init_runs = init_runs
        self.fit_intercept = fit_intercept
    
    def fit(self, x, y):
        m, n = x.shape
        self.coef = np.zeros((n, 1))
        self.intercept = np.zeros((1, 1))
        
        for epoch in range(self.init_runs):
            for run in range(m):
                idx = np.random.randint(0, m)
                x_i = x[idx:idx+1, :]
                y_i = y[idx]
                y_pred_i = sigmoid(x_i.dot(self.coef) + self.intercept)
                gradient_w = -(x_i.T * (y_i - y_pred_i))
                self.coef -= self.eta * gradient_w
                if self.fit_intercept:
                    gradient_b = -(y_i - y_pred_i)
                    self.intercept -= self.eta * gradient_b
        
    def predict_proba(self, x):
        m, n = x.shape
        y_pred = np.ones((m, 2))
        y_pred[:,1:2] = sigmoid(x.dot(self.coef) + self.intercept)
        y_pred[:,0:1] -= y_pred[:,1:2]
        return y_pred
    
    def predict(self, x):
        return np.round(sigmoid(x.dot(self.coef) + self.intercept))
    

m = 50000
n = 150
class1 = int(3.0 / 4.0 * m)

X = np.random.rand(m, n)
y = np.zeros((m, 1))

for obs in range(m):
    if obs < class1:
        continue
    else:
        y[obs,0] = 1

for obs in range(m):
    if obs < class1:
        X[obs, 0] = 1.0 / float(m)
    else:
        X[obs, 0] = float(obs) / float(m)
    X[obs, 1] = float(obs) / float(m)

logit = LogisticRegressor(0.05, 20)
start_time = time.time()
logit.fit(X, y)
end_time = time.time()
print(round(end_time - start_time, 2))
y_pred = logit.predict(X)
print("Percent:", y_pred.sum() / len(y_pred))