argmin/solver/landweber/
mod.rs

1// Copyright 2018-2024 argmin developers
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
5// http://opensource.org/licenses/MIT>, at your option. This file may not be
6// copied, modified, or distributed except according to those terms.
7
8//! # Landweber iteration
9//!
10//! The Landweber iteration is a solver for ill-posed linear inverse problems.
11//! See [`Landweber`] for details.
12//!
13//! ## References
14//!
15//! Landweber, L. (1951): An iteration formula for Fredholm integral equations of the first
16//! kind. Amer. J. Math. 73, 615–624
17//!
18//! <https://en.wikipedia.org/wiki/Landweber_iteration>
19
20use crate::core::{ArgminFloat, Error, Gradient, IterState, Problem, Solver, KV};
21use argmin_math::ArgminScaledSub;
22#[cfg(feature = "serde1")]
23use serde::{Deserialize, Serialize};
24
25/// # Landweber iteration
26///
27/// The Landweber iteration is a solver for ill-posed linear inverse problems.
28///
29/// In iteration `k`, the new parameter vector `x_{k+1}` is calculated from the previous parameter
30/// vector `x_k` and the gradient at `x_k` according to the following update rule:
31///
32/// `x_{k+1} = x_k - omega * \nabla f(x_k)`
33///
34/// ## Requirements on the optimization problem
35///
36/// The optimization problem is required to implement [`Gradient`].
37///
38/// ## References
39///
40/// Landweber, L. (1951): An iteration formula for Fredholm integral equations of the first
41/// kind. Amer. J. Math. 73, 615–624
42///
43/// <https://en.wikipedia.org/wiki/Landweber_iteration>
44#[derive(Clone, Copy)]
45#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
46pub struct Landweber<F> {
47    /// omega
48    omega: F,
49}
50
51impl<F> Landweber<F> {
52    /// Construct a new instance of [`Landweber`]
53    ///
54    /// # Example
55    ///
56    /// ```
57    /// # use argmin::solver::landweber::Landweber;
58    /// let omega: f64 = 0.5;
59    /// let landweber = Landweber::new(omega);
60    /// ```
61    pub fn new(omega: F) -> Self {
62        Landweber { omega }
63    }
64}
65
66impl<O, F, P, G> Solver<O, IterState<P, G, (), (), (), F>> for Landweber<F>
67where
68    O: Gradient<Param = P, Gradient = G>,
69    P: Clone + ArgminScaledSub<G, F, P>,
70    F: ArgminFloat,
71{
72    fn name(&self) -> &str {
73        "Landweber"
74    }
75
76    fn next_iter(
77        &mut self,
78        problem: &mut Problem<O>,
79        mut state: IterState<P, G, (), (), (), F>,
80    ) -> Result<(IterState<P, G, (), (), (), F>, Option<KV>), Error> {
81        let param = state.take_param().ok_or_else(argmin_error_closure!(
82            NotInitialized,
83            concat!(
84                "`Landweber` requires an initial parameter vector. ",
85                "Please provide an initial guess via `Executor`s `configure` method."
86            )
87        ))?;
88        let grad = problem.gradient(&param)?;
89        let new_param = param.scaled_sub(&self.omega, &grad);
90        Ok((state.param(new_param), None))
91    }
92}
93
94#[cfg(test)]
95mod tests {
96    use super::*;
97    use crate::core::{test_utils::TestProblem, ArgminError, State};
98    use approx::assert_relative_eq;
99
100    test_trait_impl!(landweber, Landweber<f64>);
101
102    #[test]
103    fn test_new() {
104        let omega_in: f64 = 0.5;
105        let Landweber { omega } = Landweber::new(omega_in);
106        assert_eq!(omega.to_ne_bytes(), omega_in.to_ne_bytes());
107    }
108
109    #[test]
110    fn test_next_iter_param_not_initialized() {
111        let omega: f64 = 0.5;
112        let mut landweber = Landweber::new(omega);
113        let state = IterState::new();
114        let res = landweber.next_iter(&mut Problem::new(TestProblem::new()), state);
115        assert_error!(
116            res,
117            ArgminError,
118            concat!(
119                "Not initialized: \"`Landweber` requires an initial parameter vector. ",
120                "Please provide an initial guess via `Executor`s `configure` method.\""
121            )
122        );
123    }
124
125    #[test]
126    fn test_next_iter() {
127        let omega: f64 = 0.5;
128        let mut landweber = Landweber::new(omega);
129        let state = IterState::new().param(vec![2.0, 4.0]);
130        let (state, kv) = landweber
131            .next_iter(&mut Problem::new(TestProblem::new()), state)
132            .unwrap();
133        assert!(kv.is_none());
134        let new_param = state.get_param().unwrap();
135        assert_relative_eq!(new_param[0], 1.0, epsilon = f64::EPSILON);
136        assert_relative_eq!(new_param[1], 2.0, epsilon = f64::EPSILON);
137    }
138}