<jupyter_start><jupyter_text>Implicit SGD procedure Below we repeat the expression of the implicit Stochastic Gradient Descent (SGD) update.$$\theta_{n}=\theta_{n-1}+\gamma_{n}C_{n}\nabla\texttt{log}f(y_{n};x_{n},\theta_{n})$$<jupyter_code>def sgd_iteration(sgd_obj, D: data_set, m: model, averaging: bool=True) -> tuple:
    N, p = D._X.shape
    n_passes = sgd_obj.get_value_of("n_passes")
    good_gradient = True
    good_validity = True
    theta_old = np.ones(p) / 100.0
    theta_old_ave = theta_old
    max_iters = N * n_passes
    converged = False

    t = 1
    while True:
        theta_new = sgd_obj.update(t, theta_old, D, m, good_gradient)
        if not averaging:
            sgd_obj.sync_members(theta_new)
            converged = sgd_obj.convergence(theta_new, theta_old)
        else:
            theta_new_ave = 0.5 * theta_old_ave + 0.5 * theta_new
            sgd_obj.sync_members(theta_new_ave)
            converged = sgd_obj.convergence(theta_new_ave, theta_old_ave)
            theta_old_ave = theta_new_ave
        if converged: break
        theta_old = theta_new
        if t == max_iters: break
        t += 1
    if averaging: theta_new = theta_new_ave
    return theta_new, converged, sgd_obj
def simulate_normal_data(theta: np.ndarray, size: tuple=(1500, 20)) -> data_set:
    N, p = size
    # Create the covariance matrix with values between 0.5 and 5
    S = np.eye(p) * np.random.uniform(0.5, 5, p)
    # Generate the design matrix
    X = np.random.multivariate_normal(mean=np.zeros(p), cov=S, size=N)
    # Generate the response vector
    Y = np.random.normal(X @ theta, 1, size=N)
    return data_set(X, Y)<jupyter_output><empty_output><jupyter_text>Asymptotic variance  Normal linear model First-order SGD is used for simplicity, namely $C_{n}=I$<jupyter_code>N, p = 1500, 20
theta = np.ones(p)
D = simulate_normal_data(theta=theta)
m = glm("gaussian", "identity")
timer = time
# n: number of observations
# p: number of features
N, p = D._X.shape
gammas = np.linspace(1.2, 10, 25)

details = {
    "lr": "one-dim",
    "lr_controls": {
        "scale": None,
        "alpha": 1.0,
        "gamma": 1.0,
        "c": 1.0
    },
    "reltol": 5e-4,
    "npasses": 20,
    "size": 10,
    "check": True,
    "truth": theta
}

result = np.empty((25, 150, p))
for j, gamma in enumerate(gammas):
    for i in range(150):
        details["lr_controls"]["scale"] = gamma
        # tester = ExplicitSGD(n, p, timer, **details)
        tester = ImplicitSGD(N, p, timer, **details)
        result[j, i, :], _, _ = sgd_iteration(sgd_obj=tester, D=D, m=m, averaging=False)<jupyter_output><empty_output><jupyter_text>Asymptotic normality  Running time comparison<jupyter_code>N, p = int(1e5), int(1e2)
theta = np.ones(p)
D = simulate_normal_data(theta=theta, size=(N, p))
m = glm("gaussian", "identity")
familia = sm.families.Gaussian(link=sm.families.links.identity())

glm_bnch = sm.GLM(D._Y, D._X, family = familia)
%time bnch_est = glm_bnch.fit()
timer = time
# n: number of observations
# p: number of features
N, p = D._X.shape

details = {
    "lr": "adagrad",
    "lr_controls": {
        "eta": 1.0,
        "eps": 1e-6
    },
    "reltol": 5e-4,
    "npasses": 20,
    "size": 10,
    "check": True,
    "truth": theta
}

tester = ImplicitSGD(N, p, timer, **details)
%time sgd_est, converged, _ = sgd_iteration(sgd_obj=tester, D=D, m=m, averaging=False)
sgd_est<jupyter_output><empty_output><jupyter_text>Poisson regression<jupyter_code># X = np.empty((n, p), dtype=np.float64)
X = np.random.normal(0, 1, size=(n, p))
# X = (X - np.mean(X, axis=0)) / np.std(X, axis=0)
X = np.hstack((np.ones((n, 1)), X)) ## add intercept
theta = np.ones(p+1) * 0.4
## Gaussian
eps = np.random.normal(0, 1, size=n)
# Y = X @ theta + eps
## Poisson
eta = np.exp(X @ theta)
poisson = np.vectorize(lambda x: np.random.poisson(x, 1))
Y = poisson(eta) + eps
D = data_set(X, Y)
familia = sm.families.Poisson(link=sm.families.links.log())

glm_bnch = sm.GLM(D._Y, D._X, family = familia)
true_est = glm_bnch.fit()
# true_est.params
true_est.params
m = glm("poisson", "exponential")
timer = time
# n: number of observations
# p: number of features
n, p = D._X.shape

details = {
    "lr": "adagrad",
    "lr_controls": {
        "scale": 1.0,
        "alpha": 1.0,
        "gamma": 0.6,
        "c": 0.5,
        "eta": 1.0,
        "eps": 1e-6
    },
    "reltol": 5e-4,
    "npasses": 20,
    "size": 10,
    "pass": True,
    "check": True,
    "truth": theta
}

# tester = ExplicitSGD(n, p, timer, **details)
tester = ImplicitSGD(n, p, timer, **details)
n_passes = tester.get_value_of("n_passes")
good_gradient = True
good_validity = True
averaging = False
theta_old = np.ones(p) / 100.0
theta_old_ave = theta_old
max_iters = n * n_passes
converged = False

t = 1

while True:
    theta_new = tester.update(t, theta_old, D, m, good_gradient)
    if not averaging:
        tester.sync_members(theta_new)
        converged = tester.convergence(theta_new, theta_old)
    else:
        theta_new_ave = 0.5 * theta_old_ave + 0.5 * theta_new
        tester.sync_members(theta_new_ave)
        converged = tester.convergence(theta_new_ave, theta_old_ave)
        theta_old_ave = theta_new_ave
    if converged: break
    theta_old = theta_new
    if t == max_iters: break
    t += 1
theta_new
converged<jupyter_output><empty_output><jupyter_text>Data from file<jupyter_code>#data = sm.datasets.scotland.load()
#data.exog = sm.add_constant(data.exog, prepend=False)
#D = data_set(data.exog, data.endog)
#data = sm.datasets.scotland.load(as_pandas=True)
#data.exog = sm.add_constant(data.exog, prepend=False)
#data.exog["RESPONSE"] = data.endog
#data.exog.to_csv(os.path.join("/path/to/test_data.csv"), index=False, encoding='UTF-8')
# from numpy import genfromtxt
# data = genfromtxt('test_data.csv', delimiter=',', encoding='UTF-8', skip_header=1)
# D = data_set(data[:,1:], data[:,0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14982284460941456, 'fraction_numerical': 0.0195714526742028, 'mean_word_length': 3.4909090909090907, 'pattern_counts': {'":': 34, '<': 21, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '551258', 'license': 'MIT', 'n_tokens_mistral': 2467, 'n_tokens_neox': 2316, 'n_words': 503, 'path': 'sgd_test.ipynb', 'repo_name': 'georgios-vassos1/fastsgd', 'size': '5927'}