Repository: scikit-learn/scikit-learn
Pull request: ENH Support float32 in SGDClassifier and SGDRegressor
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes: #5776

#### What does this implement/fix? Explain your changes.
- Adds float32 support for SGD classifier and regressor in addition to the current implementation for float32.
- Generates float32 and float64 variants for _plain_sgd and its related functions in _sgd_fast.pyx using Tempita.
- Makes the necessary adjustments in setup.py to generate _sgd_fast.

#### Any other comments?
- Tests to be added

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Add float32 support for SGDClassifier and SGDRegressor
Hi, 

I was wondreing why `sklearn` does not allow me to specify the `dtype` I'd like to use.

I'm working with rather large dataset, which fits into my RAM as float32, but when I'm trying to train a simple SGD on it, the model tries to copy my data into float64, causing MemoryError. 

I can change this in my local sklearn build, but I guess there is a good reason why this is not a free parameter?

Hints from issue discussion:
PR welcome ;)
It would be great to have both 32 and 64 bit everywhere, I think.
This means using fused types everywhere in Cython.

One of the reasons not to do that was the explosion in the generated C code, but I think with #5492 we need to be somewhat less careful about that.

#13243 implemented float32 support for `SequentialDataset`. This has to be propagated to `_plain_sgd` in `_sgd_fast.pyx` (and other places) in order to enable to preserve float32 `X` ndarrays in `SGDClassifier` and `SGDRegressor`.
@OmarManzoor: you might be interested in this slightly harder Cython issue.
@jjerphan Thank you for suggesting.
@jjerphan Do we need to change the _sgd_fast.pxd file in relation to this? On checking its usages it is mainly being used inside the _sag_fast module which I think already supports float32 and float64.
> @jjerphan Do we need to change the _sgd_fast.pxd file in relation to this? On checking its usages it is mainly being used inside the _sag_fast module which I think already supports float32 and float64.

On further investigation I don't think we need to change this file.

Commit message: Add preserves_dtype in _more_tags

Patch to review:
diff --git a/.gitignore b/.gitignore
index 47ec8fa2faf79..f6250b4d5f580 100644
--- a/.gitignore
+++ b/.gitignore
@@ -85,6 +85,7 @@ sklearn/utils/_seq_dataset.pxd
 sklearn/utils/_weight_vector.pyx
 sklearn/utils/_weight_vector.pxd
 sklearn/linear_model/_sag_fast.pyx
+sklearn/linear_model/_sgd_fast.pyx
 sklearn/metrics/_dist_metrics.pyx
 sklearn/metrics/_dist_metrics.pxd
 sklearn/metrics/_pairwise_distances_reduction/_argkmin.pxd
diff --git a/doc/whats_new/v1.3.rst b/doc/whats_new/v1.3.rst
index 08ebf4abc92c3..f7bb8b9950231 100644
--- a/doc/whats_new/v1.3.rst
+++ b/doc/whats_new/v1.3.rst
@@ -153,6 +153,13 @@ Changelog
 - |Enhancement| Added the parameter `fill_value` to :class:`impute.IterativeImputer`.
   :pr:`25232` by :user:`Thijs van Weezel <ValueInvestorThijs>`.
 
+:mod:`sklearn.linear_model`
+...........................
+
+- |Enhancement| :class:`SGDClassifier`, :class:`SGDRegressor` and
+  :class:`SGDOneClassSVM` now preserve dtype for `numpy.float32`.
+  :pr:`25587` by :user:`Omar Salman <OmarManzoor>`
+
 :mod:`sklearn.metrics`
 ......................
 
diff --git a/setup.py b/setup.py
index 2ac07ec44c947..882b32875caa8 100755
--- a/setup.py
+++ b/setup.py
@@ -308,7 +308,7 @@ def check_package_status(package, min_version):
     ],
     "linear_model": [
         {"sources": ["_cd_fast.pyx"], "include_np": True},
-        {"sources": ["_sgd_fast.pyx"], "include_np": True},
+        {"sources": ["_sgd_fast.pyx.tp"], "include_np": True},
         {"sources": ["_sag_fast.pyx.tp"], "include_np": True},
     ],
     "manifold": [
diff --git a/sklearn/linear_model/_sgd_fast.pyx b/sklearn/linear_model/_sgd_fast.pyx.tp
similarity index 76%
rename from sklearn/linear_model/_sgd_fast.pyx
rename to sklearn/linear_model/_sgd_fast.pyx.tp
index 84e54c40e1d28..cd10b73a74a07 100644
--- a/sklearn/linear_model/_sgd_fast.pyx
+++ b/sklearn/linear_model/_sgd_fast.pyx.tp
@@ -1,25 +1,57 @@
-# Author: Peter Prettenhofer <peter.prettenhofer@gmail.com>
-#         Mathieu Blondel (partial_fit support)
-#         Rob Zinkov (passive-aggressive)
-#         Lars Buitinck
-#
-# License: BSD 3 clause
+{{py:
 
+"""
+Template file to easily generate fused types consistent code using Tempita
+(https://github.com/cython/cython/blob/master/Cython/Tempita/_tempita.py).
 
+Generated file: _sgd_fast.pyx
+
+Each relevant function is duplicated for the dtypes float and double.
+The keywords between double braces are substituted in setup.py.
+
+Authors: Peter Prettenhofer <peter.prettenhofer@gmail.com>
+         Mathieu Blondel (partial_fit support)
+         Rob Zinkov (passive-aggressive)
+         Lars Buitinck
+
+License: BSD 3 clause
+"""
+
+# The dtypes are defined as follows (name_suffix, c_type, np_type)
+dtypes = [
+    ("64", "double", "np.float64"),
+    ("32", "float", "np.float32"),
+]
+
+}}
+
+#------------------------------------------------------------------------------
+
+"""
+SGD implementation
+WARNING: Do not edit .pyx file directly, it is generated from .pyx.tp
+"""
+
+from cython cimport floating
 import numpy as np
 from time import time
 
 from libc.math cimport exp, log, pow, fabs
 cimport numpy as cnp
 from numpy.math cimport INFINITY
-cdef extern from "_sgd_fast_helpers.h":
-    bint skl_isfinite(double) nogil
 
-from ..utils._weight_vector cimport WeightVector64 as WeightVector
-from ..utils._seq_dataset cimport SequentialDataset64 as SequentialDataset
+from ..utils._weight_vector cimport WeightVector32, WeightVector64
+from ..utils._seq_dataset cimport SequentialDataset32, SequentialDataset64
 
 cnp.import_array()
 
+{{for name_suffix, c_type, np_type in dtypes}}
+
+cdef extern from "_sgd_fast_helpers.h":
+    bint skl_isfinite{{name_suffix}}({{c_type}}) nogil
+
+{{endfor}}
+
 # Penalty constants
 DEF NO_PENALTY = 0
 DEF L1 = 1
@@ -34,6 +66,11 @@ DEF ADAPTIVE = 4
 DEF PA1 = 5
 DEF PA2 = 6
 
+# Fused type for the Dataset
+ctypedef fused DatasetType:
+    SequentialDataset32
+    SequentialDataset64
+
 
 
 # ----------------------------------------
@@ -359,37 +396,145 @@ cdef class SquaredEpsilonInsensitive(Regression):
     def __reduce__(self):
         return SquaredEpsilonInsensitive, (self.epsilon,)
 
-
-def _plain_sgd(const double[::1] weights,
-               double intercept,
-               const double[::1] average_weights,
-               double average_intercept,
-               LossFunction loss,
-               int penalty_type,
-               double alpha, double C,
-               double l1_ratio,
-               SequentialDataset dataset,
-               const unsigned char[::1] validation_mask,
-               bint early_stopping, validation_score_cb,
-               int n_iter_no_change,
-               unsigned int max_iter, double tol, int fit_intercept,
-               int verbose, bint shuffle, cnp.uint32_t seed,
-               double weight_pos, double weight_neg,
-               int learning_rate, double eta0,
-               double power_t,
-               bint one_class,
-               double t=1.0,
-               double intercept_decay=1.0,
-               int average=0):
+def _plain_sgd(
+    const floating[::1] weights,
+    double intercept,
+    const floating[::1] average_weights,
+    double average_intercept,
+    LossFunction loss,
+    int penalty_type,
+    double alpha,
+    double C,
+    double l1_ratio,
+    DatasetType dataset,
+    const unsigned char[::1] validation_mask,
+    bint early_stopping,
+    validation_score_cb,
+    int n_iter_no_change,
+    unsigned int max_iter,
+    double tol,
+    int fit_intercept,
+    int verbose,
+    bint shuffle,
+    cnp.uint32_t seed,
+    double weight_pos,
+    double weight_neg,
+    int learning_rate,
+    double eta0,
+    double power_t,
+    bint one_class,
+    double t=1.0,
+    double intercept_decay=1.0,
+    int average=0,
+):
+    """SGD for generic loss functions and penalties with optional averaging"""
+    if floating is float:
+        return _plain_sgd32(
+            weights=weights,
+            intercept=intercept,
+            average_weights=average_weights,
+            average_intercept=average_intercept,
+            loss=loss,
+            penalty_type=penalty_type,
+            alpha=alpha,
+            C=C,
+            l1_ratio=l1_ratio,
+            dataset=dataset,
+            validation_mask=validation_mask,
+            early_stopping=early_stopping,
+            validation_score_cb=validation_score_cb,
+            n_iter_no_change=n_iter_no_change,
+            max_iter=max_iter,
+            tol=tol,
+            fit_intercept=fit_intercept,
+            verbose=verbose,
+            shuffle=shuffle,
+            seed=seed,
+            weight_pos=weight_pos,
+            weight_neg=weight_neg,
+            learning_rate=learning_rate,
+            eta0=eta0,
+            power_t=power_t,
+            one_class=one_class,
+            t=t,
+            intercept_decay=intercept_decay,
+            average=average,
+        )
+    else:
+        return _plain_sgd64(
+            weights=weights,
+            intercept=intercept,
+            average_weights=average_weights,
+            average_intercept=average_intercept,
+            loss=loss,
+            penalty_type=penalty_type,
+            alpha=alpha,
+            C=C,
+            l1_ratio=l1_ratio,
+            dataset=dataset,
+            validation_mask=validation_mask,
+            early_stopping=early_stopping,
+            validation_score_cb=validation_score_cb,
+            n_iter_no_change=n_iter_no_change,
+            max_iter=max_iter,
+            tol=tol,
+            fit_intercept=fit_intercept,
+            verbose=verbose,
+            shuffle=shuffle,
+            seed=seed,
+            weight_pos=weight_pos,
+            weight_neg=weight_neg,
+            learning_rate=learning_rate,
+            eta0=eta0,
+            power_t=power_t,
+            one_class=one_class,
+            t=t,
+            intercept_decay=intercept_decay,
+            average=average,
+        )
+
+{{for name_suffix, c_type, np_type in dtypes}}
+
+def _plain_sgd{{name_suffix}}(
+    const {{c_type}}[::1] weights,
+    double intercept,
+    const {{c_type}}[::1] average_weights,
+    double average_intercept,
+    LossFunction loss,
+    int penalty_type,
+    double alpha,
+    double C,
+    double l1_ratio,
+    SequentialDataset{{name_suffix}} dataset,
+    const unsigned char[::1] validation_mask,
+    bint early_stopping,
+    validation_score_cb,
+    int n_iter_no_change,
+    unsigned int max_iter,
+    double tol,
+    int fit_intercept,
+    int verbose,
+    bint shuffle,
+    cnp.uint32_t seed,
+    double weight_pos,
+    double weight_neg,
+    int learning_rate,
+    double eta0,
+    double power_t,
+    bint one_class,
+    double t=1.0,
+    double intercept_decay=1.0,
+    int average=0,
+):
     """SGD for generic loss functions and penalties with optional averaging
 
     Parameters
     ----------
-    weights : ndarray[double, ndim=1]
+    weights : ndarray[{{c_type}}, ndim=1]
         The allocated vector of weights.
     intercept : double
         The initial intercept.
-    average_weights : ndarray[double, ndim=1]
+    average_weights : ndarray[{{c_type}}, ndim=1]
         The average weights as computed for ASGD. Should be None if average
         is 0.
     average_intercept : double
@@ -476,8 +621,8 @@ def _plain_sgd(const double[::1] weights,
     cdef Py_ssize_t n_samples = dataset.n_samples
     cdef Py_ssize_t n_features = weights.shape[0]
 
-    cdef WeightVector w = WeightVector(weights, average_weights)
-    cdef double *x_data_ptr = NULL
+    cdef WeightVector{{name_suffix}} w = WeightVector{{name_suffix}}(weights, average_weights)
+    cdef {{c_type}} *x_data_ptr = NULL
     cdef int *x_ind_ptr = NULL
 
     # helper variables
@@ -492,9 +637,9 @@ def _plain_sgd(const double[::1] weights,
     cdef double score = 0.0
     cdef double best_loss = INFINITY
     cdef double best_score = -INFINITY
-    cdef double y = 0.0
-    cdef double sample_weight
-    cdef double class_weight = 1.0
+    cdef {{c_type}} y = 0.0
+    cdef {{c_type}} sample_weight
+    cdef {{c_type}} class_weight = 1.0
     cdef unsigned int count = 0
     cdef unsigned int train_count = n_samples - np.sum(validation_mask)
     cdef unsigned int epoch = 0
@@ -507,10 +652,10 @@ def _plain_sgd(const double[::1] weights,
     cdef long long sample_index
 
     # q vector is only used for L1 regularization
-    cdef double[::1] q = None
-    cdef double * q_data_ptr = NULL
+    cdef {{c_type}}[::1] q = None
+    cdef {{c_type}} * q_data_ptr = NULL
     if penalty_type == L1 or penalty_type == ELASTICNET:
-        q = np.zeros((n_features,), dtype=np.float64, order="c")
+        q = np.zeros((n_features,), dtype={{np_type}}, order="c")
         q_data_ptr = &q[0]
     cdef double u = 0.0
 
@@ -561,12 +706,12 @@ def _plain_sgd(const double[::1] weights,
                     class_weight = weight_neg
 
                 if learning_rate == PA1:
-                    update = sqnorm(x_data_ptr, x_ind_ptr, xnnz)
+                    update = sqnorm{{name_suffix}}(x_data_ptr, x_ind_ptr, xnnz)
                     if update == 0:
                         continue
                     update = min(C, loss.loss(p, y) / update)
                 elif learning_rate == PA2:
-                    update = sqnorm(x_data_ptr, x_ind_ptr, xnnz)
+                    update = sqnorm{{name_suffix}}(x_data_ptr, x_ind_ptr, xnnz)
                     update = loss.loss(p, y) / (update + 0.5 / C)
                 else:
                     dloss = loss.dloss(p, y)
@@ -614,7 +759,7 @@ def _plain_sgd(const double[::1] weights,
 
                 if penalty_type == L1 or penalty_type == ELASTICNET:
                     u += (l1_ratio * eta * alpha)
-                    l1penalty(w, q_data_ptr, x_ind_ptr, xnnz, u)
+                    l1penalty{{name_suffix}}(w, q_data_ptr, x_ind_ptr, xnnz, u)
 
                 t += 1
                 count += 1
@@ -630,8 +775,8 @@ def _plain_sgd(const double[::1] weights,
                           % (time() - t_start))
 
             # floating-point under-/overflow check.
-            if (not skl_isfinite(intercept)
-                or any_nonfinite(&weights[0], n_features)):
+            if (not skl_isfinite{{name_suffix}}(intercept)
+                or any_nonfinite{{name_suffix}}(&weights[0], n_features)):
                 infinity = True
                 break
 
@@ -681,15 +826,27 @@ def _plain_sgd(const double[::1] weights,
         epoch + 1
     )
 
+{{endfor}}
+
+
+{{for name_suffix, c_type, np_type in dtypes}}
 
-cdef bint any_nonfinite(const double *w, int n) nogil:
+cdef bint any_nonfinite{{name_suffix}}(const {{c_type}} *w, int n) nogil:
     for i in range(n):
-        if not skl_isfinite(w[i]):
+        if not skl_isfinite{{name_suffix}}(w[i]):
             return True
     return 0
 
+{{endfor}}
 
-cdef double sqnorm(double * x_data_ptr, int * x_ind_ptr, int xnnz) nogil:
+
+{{for name_suffix, c_type, np_type in dtypes}}
+
+cdef double sqnorm{{name_suffix}}(
+    {{c_type}} * x_data_ptr,
+    int * x_ind_ptr,
+    int xnnz,
+) nogil:
     cdef double x_norm = 0.0
     cdef int j
     cdef double z
@@ -698,9 +855,18 @@ cdef double sqnorm(double * x_data_ptr, int * x_ind_ptr, int xnnz) nogil:
         x_norm += z * z
     return x_norm
 
+{{endfor}}
+
 
-cdef void l1penalty(WeightVector w, double * q_data_ptr,
-                    int *x_ind_ptr, int xnnz, double u) nogil:
+{{for name_suffix, c_type, np_type in dtypes}}
+
+cdef void l1penalty{{name_suffix}}(
+    WeightVector{{name_suffix}} w,
+    {{c_type}} * q_data_ptr,
+    int *x_ind_ptr,
+    int xnnz,
+    double u,
+) nogil:
     """Apply the L1 penalty to each updated feature
 
     This implements the truncated gradient approach by
@@ -710,7 +876,7 @@ cdef void l1penalty(WeightVector w, double * q_data_ptr,
     cdef int j = 0
     cdef int idx = 0
     cdef double wscale = w.wscale
-    cdef double *w_data_ptr = w.w_data_ptr
+    cdef {{c_type}} *w_data_ptr = w.w_data_ptr
     for j in range(xnnz):
         idx = x_ind_ptr[j]
         z = w_data_ptr[idx]
@@ -723,3 +889,5 @@ cdef void l1penalty(WeightVector w, double * q_data_ptr,
                 0.0, w_data_ptr[idx] + ((u - q_data_ptr[idx]) / wscale))
 
         q_data_ptr[idx] += wscale * (w_data_ptr[idx] - z)
+
+{{endfor}}
diff --git a/sklearn/linear_model/_stochastic_gradient.py b/sklearn/linear_model/_stochastic_gradient.py
index 6d204719aaf0e..b94583d234fb3 100644
--- a/sklearn/linear_model/_stochastic_gradient.py
+++ b/sklearn/linear_model/_stochastic_gradient.py
@@ -182,43 +182,51 @@ def _get_penalty_type(self, penalty):
         return PENALTY_TYPES[penalty]
 
     def _allocate_parameter_mem(
-        self, n_classes, n_features, coef_init=None, intercept_init=None, one_class=0
+        self,
+        n_classes,
+        n_features,
+        input_dtype,
+        coef_init=None,
+        intercept_init=None,
+        one_class=0,
     ):
         """Allocate mem for parameters; initialize if provided."""
         if n_classes > 2:
             # allocate coef_ for multi-class
             if coef_init is not None:
-                coef_init = np.asarray(coef_init, order="C")
+                coef_init = np.asarray(coef_init, dtype=input_dtype, order="C")
                 if coef_init.shape != (n_classes, n_features):
                     raise ValueError("Provided ``coef_`` does not match dataset. ")
                 self.coef_ = coef_init
             else:
                 self.coef_ = np.zeros(
-                    (n_classes, n_features), dtype=np.float64, order="C"
+                    (n_classes, n_features), dtype=input_dtype, order="C"
                 )
 
             # allocate intercept_ for multi-class
             if intercept_init is not None:
-                intercept_init = np.asarray(intercept_init, order="C")
+                intercept_init = np.asarray(
+                    intercept_init, order="C", dtype=input_dtype
+                )
                 if intercept_init.shape != (n_classes,):
                     raise ValueError("Provided intercept_init does not match dataset.")
                 self.intercept_ = intercept_init
             else:
-                self.intercept_ = np.zeros(n_classes, dtype=np.float64, order="C")
+                self.intercept_ = np.zeros(n_classes, dtype=input_dtype, order="C")
         else:
             # allocate coef_
             if coef_init is not None:
-                coef_init = np.asarray(coef_init, dtype=np.float64, order="C")
+                coef_init = np.asarray(coef_init, dtype=input_dtype, order="C")
                 coef_init = coef_init.ravel()
                 if coef_init.shape != (n_features,):
                     raise ValueError("Provided coef_init does not match dataset.")
                 self.coef_ = coef_init
             else:
-                self.coef_ = np.zeros(n_features, dtype=np.float64, order="C")
+                self.coef_ = np.zeros(n_features, dtype=input_dtype, order="C")
 
             # allocate intercept_
             if intercept_init is not None:
-                intercept_init = np.asarray(intercept_init, dtype=np.float64)
+                intercept_init = np.asarray(intercept_init, dtype=input_dtype)
                 if intercept_init.shape != (1,) and intercept_init.shape != ():
                     raise ValueError("Provided intercept_init does not match dataset.")
                 if one_class:
@@ -231,21 +239,23 @@ def _allocate_parameter_mem(
                     )
             else:
                 if one_class:
-                    self.offset_ = np.zeros(1, dtype=np.float64, order="C")
+                    self.offset_ = np.zeros(1, dtype=input_dtype, order="C")
                 else:
-                    self.intercept_ = np.zeros(1, dtype=np.float64, order="C")
+                    self.intercept_ = np.zeros(1, dtype=input_dtype, order="C")
 
         # initialize average parameters
         if self.average > 0:
             self._standard_coef = self.coef_
-            self._average_coef = np.zeros(self.coef_.shape, dtype=np.float64, order="C")
+            self._average_coef = np.zeros(
+                self.coef_.shape, dtype=input_dtype, order="C"
+            )
             if one_class:
                 self._standard_intercept = 1 - self.offset_
             else:
                 self._standard_intercept = self.intercept_
 
             self._average_intercept = np.zeros(
-                self._standard_intercept.shape, dtype=np.float64, order="C"
+                self._standard_intercept.shape, dtype=input_dtype, order="C"
             )
 
     def _make_validation_split(self, y, sample_mask):
@@ -318,12 +328,12 @@ def _make_validation_score_cb(
         )
 
 
-def _prepare_fit_binary(est, y, i):
+def _prepare_fit_binary(est, y, i, input_dtye):
     """Initialization for fit_binary.
 
     Returns y, coef, intercept, average_coef, average_intercept.
     """
-    y_i = np.ones(y.shape, dtype=np.float64, order="C")
+    y_i = np.ones(y.shape, dtype=input_dtye, order="C")
     y_i[y != est.classes_[i]] = -1.0
     average_intercept = 0
     average_coef = None
@@ -418,7 +428,7 @@ def fit_binary(
     # if average is not true, average_coef, and average_intercept will be
     # unused
     y_i, coef, intercept, average_coef, average_intercept = _prepare_fit_binary(
-        est, y, i
+        est, y, i, input_dtye=X.dtype
     )
     assert y_i.shape[0] == y.shape[0] == sample_weight.shape[0]
 
@@ -580,7 +590,7 @@ def _partial_fit(
             X,
             y,
             accept_sparse="csr",
-            dtype=np.float64,
+            dtype=[np.float64, np.float32],
             order="C",
             accept_large_sparse=False,
             reset=first_call,
@@ -596,11 +606,15 @@ def _partial_fit(
         self._expanded_class_weight = compute_class_weight(
             self.class_weight, classes=self.classes_, y=y
         )
-        sample_weight = _check_sample_weight(sample_weight, X)
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
         if getattr(self, "coef_", None) is None or coef_init is not None:
             self._allocate_parameter_mem(
-                n_classes, n_features, coef_init, intercept_init
+                n_classes=n_classes,
+                n_features=n_features,
+                input_dtype=X.dtype,
+                coef_init=coef_init,
+                intercept_init=intercept_init,
             )
         elif n_features != self.coef_.shape[-1]:
             raise ValueError(
@@ -1351,7 +1365,8 @@ def _more_tags(self):
                 "check_sample_weights_invariance": (
                     "zero sample_weight is not equivalent to removing samples"
                 ),
-            }
+            },
+            "preserves_dtype": [np.float64, np.float32],
         }
 
 
@@ -1438,22 +1453,28 @@ def _partial_fit(
             accept_sparse="csr",
             copy=False,
             order="C",
-            dtype=np.float64,
+            dtype=[np.float64, np.float32],
             accept_large_sparse=False,
             reset=first_call,
         )
-        y = y.astype(np.float64, copy=False)
+        y = y.astype(X.dtype, copy=False)
 
         n_samples, n_features = X.shape
 
-        sample_weight = _check_sample_weight(sample_weight, X)
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
         # Allocate datastructures from input arguments
         if first_call:
-            self._allocate_parameter_mem(1, n_features, coef_init, intercept_init)
+            self._allocate_parameter_mem(
+                n_classes=1,
+                n_features=n_features,
+                input_dtype=X.dtype,
+                coef_init=coef_init,
+                intercept_init=intercept_init,
+            )
         if self.average > 0 and getattr(self, "_average_coef", None) is None:
-            self._average_coef = np.zeros(n_features, dtype=np.float64, order="C")
-            self._average_intercept = np.zeros(1, dtype=np.float64, order="C")
+            self._average_coef = np.zeros(n_features, dtype=X.dtype, order="C")
+            self._average_intercept = np.zeros(1, dtype=X.dtype, order="C")
 
         self._fit_regressor(
             X, y, alpha, C, loss, learning_rate, sample_weight, max_iter
@@ -1996,7 +2017,8 @@ def _more_tags(self):
                 "check_sample_weights_invariance": (
                     "zero sample_weight is not equivalent to removing samples"
                 ),
-            }
+            },
+            "preserves_dtype": [np.float64, np.float32],
         }
 
 
@@ -2197,7 +2219,7 @@ def _fit_one_class(self, X, alpha, C, sample_weight, learning_rate, max_iter):
         # The One-Class SVM uses the SGD implementation with
         # y=np.ones(n_samples).
         n_samples = X.shape[0]
-        y = np.ones(n_samples, dtype=np.float64, order="C")
+        y = np.ones(n_samples, dtype=X.dtype, order="C")
 
         dataset, offset_decay = make_dataset(X, y, sample_weight)
 
@@ -2304,7 +2326,7 @@ def _partial_fit(
             X,
             None,
             accept_sparse="csr",
-            dtype=np.float64,
+            dtype=[np.float64, np.float32],
             order="C",
             accept_large_sparse=False,
             reset=first_call,
@@ -2313,13 +2335,20 @@ def _partial_fit(
         n_features = X.shape[1]
 
         # Allocate datastructures from input arguments
-        sample_weight = _check_sample_weight(sample_weight, X)
+        sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
 
         # We use intercept = 1 - offset where intercept is the intercept of
         # the SGD implementation and offset is the offset of the One-Class SVM
         # optimization problem.
         if getattr(self, "coef_", None) is None or coef_init is not None:
-            self._allocate_parameter_mem(1, n_features, coef_init, offset_init, 1)
+            self._allocate_parameter_mem(
+                n_classes=1,
+                n_features=n_features,
+                input_dtype=X.dtype,
+                coef_init=coef_init,
+                intercept_init=offset_init,
+                one_class=1,
+            )
         elif n_features != self.coef_.shape[-1]:
             raise ValueError(
                 "Number of features %d does not match previous data %d."
@@ -2327,8 +2356,8 @@ def _partial_fit(
             )
 
         if self.average and getattr(self, "_average_coef", None) is None:
-            self._average_coef = np.zeros(n_features, dtype=np.float64, order="C")
-            self._average_intercept = np.zeros(1, dtype=np.float64, order="C")
+            self._average_coef = np.zeros(n_features, dtype=X.dtype, order="C")
+            self._average_intercept = np.zeros(1, dtype=X.dtype, order="C")
 
         self.loss_function_ = self._get_loss_function(loss)
         if not hasattr(self, "t_"):
@@ -2543,5 +2572,6 @@ def _more_tags(self):
                 "check_sample_weights_invariance": (
                     "zero sample_weight is not equivalent to removing samples"
                 )
-            }
+            },
+            "preserves_dtype": [np.float64, np.float32],
         }
diff --git a/sklearn/linear_model/tests/test_sgd.py b/sklearn/linear_model/tests/test_sgd.py
index a3944b0497d97..58298d95367de 100644
--- a/sklearn/linear_model/tests/test_sgd.py
+++ b/sklearn/linear_model/tests/test_sgd.py
@@ -2168,3 +2168,50 @@ def test_sgd_error_on_zero_validation_weight():
 def test_sgd_verbose(Estimator):
     """non-regression test for gh #25249"""
     Estimator(verbose=1).fit(X, Y)
+
+
+@pytest.mark.parametrize(
+    "sgd_class",
+    [
+        SGDClassifier,
+        SparseSGDClassifier,
+        SGDRegressor,
+        SparseSGDRegressor,
+        SGDOneClassSVM,
+        SparseSGDOneClassSVM,
+    ],
+)
+@pytest.mark.parametrize("data_type", (np.float32, np.float64))
+def test_sgd_dtype_match(sgd_class, data_type):
+    _X = X.astype(data_type)
+    _Y = np.array(Y, dtype=data_type)
+    sgd_model = sgd_class()
+    sgd_model.fit(_X, _Y)
+    assert sgd_model.coef_.dtype == data_type
+
+
+@pytest.mark.parametrize(
+    "sgd_class",
+    [
+        SGDClassifier,
+        SparseSGDClassifier,
+        SGDRegressor,
+        SparseSGDRegressor,
+        SGDOneClassSVM,
+        SparseSGDOneClassSVM,
+    ],
+)
+def test_sgd_numerical_consistency(sgd_class):
+    X_64 = X.astype(dtype=np.float64)
+    Y_64 = np.array(Y, dtype=np.float64)
+
+    X_32 = X.astype(dtype=np.float32)
+    Y_32 = np.array(Y, dtype=np.float32)
+
+    sgd_64 = sgd_class(max_iter=20)
+    sgd_64.fit(X_64, Y_64)
+
+    sgd_32 = sgd_class(max_iter=20)
+    sgd_32.fit(X_32, Y_32)
+
+    assert_allclose(sgd_64.coef_, sgd_32.coef_)


Code review comments:
diff_hunk: @@ -359,37 +396,145 @@ cdef class SquaredEpsilonInsensitive(Regression):
     def __reduce__(self):
         return SquaredEpsilonInsensitive, (self.epsilon,)
 
-
-def _plain_sgd(const double[::1] weights,
-               double intercept,
-               const double[::1] average_weights,
-               double average_intercept,
-               LossFunction loss,
-               int penalty_type,
-               double alpha, double C,
-               double l1_ratio,
-               SequentialDataset dataset,
-               const unsigned char[::1] validation_mask,
-               bint early_stopping, validation_score_cb,
-               int n_iter_no_change,
-               unsigned int max_iter, double tol, int fit_intercept,
-               int verbose, bint shuffle, cnp.uint32_t seed,
-               double weight_pos, double weight_neg,
-               int learning_rate, double eta0,
-               double power_t,
-               bint one_class,
-               double t=1.0,
-               double intercept_decay=1.0,
-               int average=0):
+def _plain_sgd(
+    const floating[::1] weights,
+    double intercept,
+    const floating[::1] average_weights,
+    double average_intercept,
+    LossFunction loss,
+    int penalty_type,
+    double alpha,
+    double C,
+    double l1_ratio,
+    DatasetType dataset,
+    const unsigned char[::1] validation_mask,
+    bint early_stopping,
+    validation_score_cb,
+    int n_iter_no_change,
+    unsigned int max_iter,
+    double tol,
+    int fit_intercept,
+    int verbose,
+    bint shuffle,
+    cnp.uint32_t seed,
+    double weight_pos,
+    double weight_neg,
+    int learning_rate,
+    double eta0,
+    double power_t,
+    bint one_class,
+    double t=1.0,
+    double intercept_decay=1.0,
+    int average=0,
+):
+    """SGD for generic loss functions and penalties with optional averaging"""
+    if floating is float:
+        return _plain_sgd32(
+            weights=weights,
+            intercept=intercept,
+            average_weights=average_weights,
+            average_intercept=average_intercept,
+            loss=loss,
+            penalty_type=penalty_type,
+            alpha=alpha,
+            C=C,
+            l1_ratio=l1_ratio,
+            dataset=dataset,
+            validation_mask=validation_mask,
+            early_stopping=early_stopping,
+            validation_score_cb=validation_score_cb,
+            n_iter_no_change=n_iter_no_change,
+            max_iter=max_iter,
+            tol=tol,
+            fit_intercept=fit_intercept,
+            verbose=verbose,
+            shuffle=shuffle,
+            seed=seed,
+            weight_pos=weight_pos,
+            weight_neg=weight_neg,
+            learning_rate=learning_rate,
+            eta0=eta0,
+            power_t=power_t,
+            one_class=one_class,
+            t=t,
+            intercept_decay=intercept_decay,
+            average=average,
+        )
+    else:
+        return _plain_sgd64(
+            weights=weights,
+            intercept=intercept,
+            average_weights=average_weights,
+            average_intercept=average_intercept,
+            loss=loss,
+            penalty_type=penalty_type,
+            alpha=alpha,
+            C=C,
+            l1_ratio=l1_ratio,
+            dataset=dataset,
+            validation_mask=validation_mask,
+            early_stopping=early_stopping,
+            validation_score_cb=validation_score_cb,
+            n_iter_no_change=n_iter_no_change,
+            max_iter=max_iter,
+            tol=tol,
+            fit_intercept=fit_intercept,
+            verbose=verbose,
+            shuffle=shuffle,
+            seed=seed,
+            weight_pos=weight_pos,
+            weight_neg=weight_neg,
+            learning_rate=learning_rate,
+            eta0=eta0,
+            power_t=power_t,
+            one_class=one_class,
+            t=t,
+            intercept_decay=intercept_decay,
+            average=average,
+        )
path: sklearn/linear_model/_sgd_fast.pyx.tp
text: @user1:
You should be able to simplify it to assign one of the interfaces to a variable and use it as a callable rather than having two similar calls.

diff_hunk: @@ -2168,3 +2168,50 @@ def test_sgd_error_on_zero_validation_weight():
 def test_sgd_verbose(Estimator):
     """non-regression test for gh #25249"""
     Estimator(verbose=1).fit(X, Y)
+
+
+@pytest.mark.parametrize(
+    "sgd_class",
+    [
+        SGDClassifier,
+        SparseSGDClassifier,
+        SGDRegressor,
+        SparseSGDRegressor,
+        SGDOneClassSVM,
+        SparseSGDOneClassSVM,
+    ],
+)
+@pytest.mark.parametrize("data_type", (np.float32, np.float64))
+def test_sgd_dtype_match(sgd_class, data_type):
+    _X = X.astype(data_type)
+    _Y = np.array(Y, dtype=data_type)
+    sgd_model = sgd_class()
+    sgd_model.fit(_X, _Y)
+    assert sgd_model.coef_.dtype == data_type
+
+
+@pytest.mark.parametrize(
+    "sgd_class",
path: sklearn/linear_model/tests/test_sgd.py
text: @user1:
Can you replace such occurences with `SGDEstimator`?