--- a/src/sliding_fb.rs Thu Jan 23 23:39:40 2025 +0100
+++ b/src/sliding_fb.rs Fri Jan 31 18:00:17 2025 -0500
@@ -22,7 +22,7 @@
use crate::forward_model::{
ForwardModel,
AdjointProductBoundedBy,
- LipschitzValues,
+ BoundedCurvature,
};
//use crate::tolerance::Tolerance;
use crate::plot::{
@@ -99,6 +99,7 @@
/// Internal type of adaptive transport step length calculation
pub(crate) enum TransportStepLength<F : Float, G : Fn(F, F) -> F> {
/// Fixed, known step length
+ #[allow(dead_code)]
Fixed(F),
/// Adaptive step length, only wrt. maximum transport.
/// Content of `l` depends on use case, while `g` calculates the step length from `l`.
@@ -158,16 +159,16 @@
ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ);
}
},
- AdaptiveMax{ l : ℓ_v, ref mut max_transport, g : ref calculate_θ } => {
+ AdaptiveMax{ l : ℓ_F, ref mut max_transport, g : ref calculate_θ } => {
*max_transport = max_transport.max(γ1.norm(Radon));
- let θτ = τ * calculate_θ(ℓ_v, *max_transport);
+ let θτ = τ * calculate_θ(ℓ_F, *max_transport);
for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ);
}
},
- FullyAdaptive{ l : ref mut adaptive_ℓ_v, ref mut max_transport, g : ref calculate_θ } => {
+ FullyAdaptive{ l : ref mut adaptive_ℓ_F, ref mut max_transport, g : ref calculate_θ } => {
*max_transport = max_transport.max(γ1.norm(Radon));
- let mut θ = calculate_θ(*adaptive_ℓ_v, *max_transport);
+ let mut θ = calculate_θ(*adaptive_ℓ_F, *max_transport);
// Do two runs through the spikes to update θ, breaking if first run did not cause
// a change.
for _i in 0..=1 {
@@ -179,9 +180,9 @@
let n = g.norm2();
if n >= F::EPSILON {
// Estimate Lipschitz factor of ∇v
- let this_ℓ_v = (dv_x - v.differential(&ρ.x)).norm2() / n;
- *adaptive_ℓ_v = adaptive_ℓ_v.max(this_ℓ_v);
- θ = calculate_θ(*adaptive_ℓ_v, *max_transport);
+ let this_ℓ_F = (dv_x - v.differential(&ρ.x)).norm2() / n;
+ *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F);
+ θ = calculate_θ(*adaptive_ℓ_F, *max_transport);
changes = true
}
}
@@ -274,10 +275,9 @@
F : Float + ToNalgebraRealField,
I : AlgIteratorFactory<IterInfo<F, N>>,
A : ForwardModel<RNDM<F, N>, F>
- + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType=F>,
- //+ TransportLipschitz<L2Squared, FloatType=F>,
+ + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType=F>
+ + BoundedCurvature<FloatType=F>,
for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>,
- for<'b> A::Preadjoint<'b> : LipschitzValues<FloatType=F>,
A::PreadjointCodomain : DifferentiableRealMapping<F, N>,
RNDM<F, N> : SpikeMerging<F>,
Reg : SlidingRegTerm<F, N>,
@@ -301,12 +301,19 @@
//let ℓ = opA.transport.lipschitz_factor(L2Squared) * max_transport;
let ℓ = 0.0;
let τ = config.τ0 / opA.adjoint_product_bound(prox_penalty).unwrap();
- let calculate_θ = |ℓ_v, _| config.transport.θ0 / (τ*(ℓ + ℓ_v));
- let mut θ_or_adaptive = match opA.preadjoint().value_diff_unit_lipschitz_factor() {
- // We only estimate w (the uniform Lipschitz for of v), if we also estimate ℓ_v
- // (the uniform Lipschitz factor of ∇v).
- // We assume that the residual is decreasing.
- Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * b.norm2(), 0.0)),
+ let (maybe_ℓ_v0, maybe_transport_lip) = opA.curvature_bound_components();
+ let transport_lip = maybe_transport_lip.unwrap();
+ let calculate_θ = |ℓ_v, max_transport| {
+ let ℓ_F = ℓ_v + transport_lip * max_transport;
+ config.transport.θ0 / (τ*(ℓ + ℓ_F))
+ };
+ let mut θ_or_adaptive = match maybe_ℓ_v0 {
+ //Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * b.norm2(), 0.0)),
+ Some(ℓ_v0) => TransportStepLength::AdaptiveMax {
+ l: ℓ_v0 * b.norm2(), // TODO: could estimate computing the real reesidual
+ max_transport : 0.0,
+ g : calculate_θ
+ },
None => TransportStepLength::FullyAdaptive {
l : 10.0 * F::EPSILON, // Start with something very small to estimate differentials
max_transport : 0.0,