dml_core/learner/
logistic.rs

1use ndarray::{concatenate, Array1, Array2, ArrayView1, ArrayView2, Axis};
2
3use super::Classifier;
4use crate::error::{DmlError, Result};
5use crate::linalg::cholesky_solve;
6
7/// L2-penalised logistic regression with an unpenalised intercept, fitted by
8/// Newton's method with backtracking line search.
9///
10/// Minimises `C · Σ logloss + ½‖w‖²`, the same objective as scikit-learn's
11/// `LogisticRegression(C=c)`; the default `C = 1` matches scikit-learn's default.
12#[derive(Debug, Clone)]
13pub struct LogisticRegression {
14    c: f64,
15    tol: f64,
16    max_iter: usize,
17    coef: Option<Array1<f64>>,
18    intercept: f64,
19}
20
21impl Default for LogisticRegression {
22    fn default() -> Self {
23        Self::new(1.0)
24    }
25}
26
27impl LogisticRegression {
28    /// Inverse regularisation strength `c` (larger = weaker penalty).
29    pub fn new(c: f64) -> Self {
30        Self {
31            c,
32            tol: 1e-10,
33            max_iter: 100,
34            coef: None,
35            intercept: 0.0,
36        }
37    }
38
39    /// No penalty. Fails to converge if the classes are perfectly separable.
40    pub fn unpenalized() -> Self {
41        Self::new(f64::INFINITY)
42    }
43
44    pub fn tol(mut self, tol: f64) -> Self {
45        self.tol = tol;
46        self
47    }
48
49    pub fn max_iter(mut self, max_iter: usize) -> Self {
50        self.max_iter = max_iter;
51        self
52    }
53
54    pub fn coef(&self) -> Option<&Array1<f64>> {
55        self.coef.as_ref()
56    }
57
58    pub fn intercept(&self) -> f64 {
59        self.intercept
60    }
61}
62
63/// `log(1 + e^η)` without overflow.
64fn softplus(eta: f64) -> f64 {
65    if eta > 0.0 {
66        eta + (-eta).exp().ln_1p()
67    } else {
68        eta.exp().ln_1p()
69    }
70}
71
72fn sigmoid(eta: f64) -> f64 {
73    if eta >= 0.0 {
74        1.0 / (1.0 + (-eta).exp())
75    } else {
76        let e = eta.exp();
77        e / (1.0 + e)
78    }
79}
80
81impl LogisticRegression {
82    /// `[intercept, coef…]` of the fitted model, to warm-start another fit.
83    pub(super) fn params(&self) -> Option<Array1<f64>> {
84        let coef = self.coef.as_ref()?;
85        let mut beta = Array1::zeros(coef.len() + 1);
86        beta[0] = self.intercept;
87        beta.slice_mut(ndarray::s![1..]).assign(coef);
88        Some(beta)
89    }
90
91    /// Newton's method started from `init = [intercept, coef…]` if given,
92    /// else from the intercept-only fit.
93    pub(super) fn fit_from(
94        &mut self,
95        x: ArrayView2<'_, f64>,
96        y: ArrayView1<'_, f64>,
97        init: Option<&Array1<f64>>,
98    ) -> Result<()> {
99        self.fit_impl(x, y, None, init)
100    }
101
102    /// Newton's method on `Σ wᵢ·loglossᵢ + ‖coef‖² / (2C)` (unit weights if
103    /// `weights` is `None`), the minimiser scikit-learn's `sample_weight` gives.
104    fn fit_impl(
105        &mut self,
106        x: ArrayView2<'_, f64>,
107        y: ArrayView1<'_, f64>,
108        weights: Option<ArrayView1<'_, f64>>,
109        init: Option<&Array1<f64>>,
110    ) -> Result<()> {
111        let (n, p) = x.dim();
112        let sw = match weights {
113            Some(w) if w.len() != n => {
114                return Err(DmlError::DimensionMismatch(
115                    "LogisticRegression::fit_weighted weights length".into(),
116                ))
117            }
118            Some(w) if !w.iter().all(|&v| v.is_finite() && v >= 0.0) || w.sum() <= 0.0 => {
119                return Err(DmlError::InvalidArgument(
120                    "sample weights must be finite, non-negative and not all zero".into(),
121                ))
122            }
123            Some(w) => w.to_owned(),
124            None => Array1::ones(n),
125        };
126        if y.len() != n {
127            return Err(DmlError::DimensionMismatch(
128                "LogisticRegression::fit x/y rows".into(),
129            ));
130        }
131        if !y.iter().all(|&v| v == 0.0 || v == 1.0) {
132            return Err(DmlError::InvalidArgument(
133                "LogisticRegression labels must be 0 or 1".into(),
134            ));
135        }
136        let n1 = y.sum();
137        if n1 == 0.0 || n1 == n as f64 {
138            return Err(DmlError::InvalidArgument(
139                "LogisticRegression needs both classes in the training data".into(),
140            ));
141        }
142        let lambda = if self.c.is_finite() {
143            1.0 / self.c
144        } else {
145            0.0
146        };
147        // Design with a leading column of ones; beta = [intercept, w].
148        let xa = concatenate![Axis(1), Array2::ones((n, 1)), x];
149        let objective = |beta: &Array1<f64>| -> f64 {
150            let eta = xa.dot(beta);
151            let loss: f64 = eta
152                .iter()
153                .zip(y.iter())
154                .zip(sw.iter())
155                .map(|((&e, &t), &w)| w * (softplus(e) - t * e))
156                .sum();
157            loss + 0.5
158                * lambda
159                * beta
160                    .slice(ndarray::s![1..])
161                    .dot(&beta.slice(ndarray::s![1..]))
162        };
163
164        let mut beta = match init {
165            Some(b) if b.len() == p + 1 => b.clone(),
166            _ => {
167                let mean = sw.dot(&y) / sw.sum();
168                let mut b = Array1::<f64>::zeros(p + 1);
169                b[0] = (mean / (1.0 - mean)).ln();
170                b
171            }
172        };
173        let mut obj = objective(&beta);
174        let mut converged = false;
175        for _ in 0..self.max_iter {
176            let prob = xa.dot(&beta).mapv(sigmoid);
177            let mut grad = xa.t().dot(&(&sw * &(&prob - &y)));
178            let w = &sw * &prob.mapv(|q| q * (1.0 - q));
179            let mut hess = xa.t().dot(&(&xa * &w.insert_axis(Axis(1))));
180            for j in 1..=p {
181                grad[j] += lambda * beta[j];
182                hess[[j, j]] += lambda;
183            }
184            let step = cholesky_solve(&hess, &grad)?;
185            // Backtracking: halve the Newton step until the objective decreases.
186            let mut t = 1.0;
187            let mut candidate = &beta - &step;
188            let mut cand_obj = objective(&candidate);
189            while cand_obj > obj && t > 1e-10 {
190                t *= 0.5;
191                candidate = &beta - &(t * &step);
192                cand_obj = objective(&candidate);
193            }
194            let change = (t * &step).iter().fold(0.0_f64, |m, v| m.max(v.abs()));
195            beta = candidate;
196            obj = cand_obj;
197            if change <= self.tol * beta.iter().fold(1.0_f64, |m, v| m.max(v.abs())) {
198                converged = true;
199                break;
200            }
201        }
202        if !converged || !beta.iter().all(|v| v.is_finite()) {
203            return Err(DmlError::Convergence(format!(
204                "LogisticRegression did not converge in {} Newton steps (separable classes?)",
205                self.max_iter
206            )));
207        }
208        self.intercept = beta[0];
209        self.coef = Some(beta.slice(ndarray::s![1..]).to_owned());
210        Ok(())
211    }
212}
213
214impl Classifier for LogisticRegression {
215    fn fit(&mut self, x: ArrayView2<'_, f64>, y: ArrayView1<'_, f64>) -> Result<()> {
216        self.fit_from(x, y, None)
217    }
218
219    fn fit_weighted(
220        &mut self,
221        x: ArrayView2<'_, f64>,
222        y: ArrayView1<'_, f64>,
223        w: ArrayView1<'_, f64>,
224    ) -> Result<()> {
225        self.fit_impl(x, y, Some(w), None)
226    }
227
228    fn predict_proba(&self, x: ArrayView2<'_, f64>) -> Result<Array1<f64>> {
229        let coef = self.coef.as_ref().ok_or(DmlError::NotFitted)?;
230        if x.ncols() != coef.len() {
231            return Err(DmlError::DimensionMismatch(
232                "LogisticRegression::predict_proba feature count".into(),
233            ));
234        }
235        Ok((x.dot(coef) + self.intercept).mapv(sigmoid))
236    }
237
238    fn fresh(&self) -> Box<dyn Classifier> {
239        Box::new(Self::new(self.c).tol(self.tol).max_iter(self.max_iter))
240    }
241}
242
243#[cfg(test)]
244mod tests {
245    use super::*;
246    use ndarray::array;
247
248    #[test]
249    fn integer_weights_equal_repeated_rows() {
250        let x = array![
251            [0.0, 1.0],
252            [1.0, 0.5],
253            [2.0, -1.0],
254            [3.0, 0.0],
255            [1.5, 2.0],
256            [0.5, -0.5]
257        ];
258        let y = array![0.0, 0.0, 1.0, 1.0, 0.0, 1.0];
259        let w = array![1.0, 2.0, 1.0, 3.0, 1.0, 2.0];
260        let mut weighted = LogisticRegression::new(2.0).tol(1e-12);
261        weighted.fit_weighted(x.view(), y.view(), w.view()).unwrap();
262        let rows: Vec<usize> = (0..6)
263            .flat_map(|i| std::iter::repeat(i).take(w[i] as usize))
264            .collect();
265        let mut repeated = LogisticRegression::new(2.0).tol(1e-12);
266        repeated
267            .fit(
268                x.select(Axis(0), &rows).view(),
269                y.select(Axis(0), &rows).view(),
270            )
271            .unwrap();
272        let (a, b) = (
273            weighted.predict_proba(x.view()).unwrap(),
274            repeated.predict_proba(x.view()).unwrap(),
275        );
276        assert!(a.iter().zip(&b).all(|(p, q)| (p - q).abs() < 1e-10));
277        assert!(weighted
278            .fit_weighted(x.view(), y.view(), (-&w).view())
279            .is_err());
280    }
281
282    #[test]
283    fn gradient_is_zero_at_optimum() {
284        let x = array![
285            [0.0, 1.0],
286            [1.0, 0.5],
287            [2.0, -1.0],
288            [3.0, 0.0],
289            [1.5, 2.0],
290            [0.5, -0.5]
291        ];
292        let y = array![0.0, 0.0, 1.0, 1.0, 0.0, 1.0];
293        let mut m = LogisticRegression::new(2.0);
294        m.fit(x.view(), y.view()).unwrap();
295        let p = m.predict_proba(x.view()).unwrap();
296        let r = &p - &y;
297        // Intercept gradient: Σ (p − y) = 0; slopes: C·Xᵀ(p − y) + w = 0.
298        assert!(r.sum().abs() < 1e-9);
299        let g = 2.0 * x.t().dot(&r) + m.coef().unwrap();
300        assert!(g.iter().all(|v| v.abs() < 1e-9));
301    }
302
303    #[test]
304    fn rejects_single_class_and_bad_labels() {
305        let x = array![[0.0], [1.0]];
306        let mut m = LogisticRegression::default();
307        assert!(m.fit(x.view(), array![1.0, 1.0].view()).is_err());
308        assert!(m.fit(x.view(), array![0.0, 2.0].view()).is_err());
309    }
310
311    #[test]
312    fn separable_unpenalized_reports_non_convergence() {
313        let x = array![[0.0], [1.0], [2.0], [3.0]];
314        let y = array![0.0, 0.0, 1.0, 1.0];
315        let err = LogisticRegression::unpenalized()
316            .fit(x.view(), y.view())
317            .unwrap_err();
318        assert!(matches!(err, DmlError::Convergence(_) | DmlError::Singular));
319    }
320}
dml_core/learner/logistic.rs

dml_core/learner/
logistic.rs
