comparison: src/linops.rs
src/linops.rs
- branch
- dev
- changeset 110
- a1278320be26
- parent 104
- e7f1cb4bec78
- child 124
- 6aa955ad8122
equal
deleted
inserted
replaced
| 109:943c6b3b9414 | 110:a1278320be26 |
|---|---|
| 1 /*! | 1 /*! |
| 2 Abstract linear operators. | 2 Abstract linear operators. |
| 3 */ | 3 */ |
| 4 | 4 |
| 5 use crate::direct_product::Pair; | 5 use crate::direct_product::Pair; |
| 6 use crate::error::DynResult; | |
| 6 use crate::instance::Instance; | 7 use crate::instance::Instance; |
| 7 pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space}; | 8 pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space}; |
| 8 use crate::norms::{Linfinity, Norm, NormExponent, PairNorm, L1, L2}; | 9 use crate::norms::{Linfinity, Norm, NormExponent, PairNorm, L1, L2}; |
| 9 use crate::types::*; | 10 use crate::types::*; |
| 10 use numeric_literals::replace_float_literals; | 11 use numeric_literals::replace_float_literals; |
| 82 { | 83 { |
| 83 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. | 84 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. |
| 84 /// This is not expected to be the norm, just any bound on it that can be | 85 /// This is not expected to be the norm, just any bound on it that can be |
| 85 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm | 86 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm |
| 86 /// in `X`, and `codexp` in the codomain. | 87 /// in `X`, and `codexp` in the codomain. |
| 87 fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> F; | 88 /// |
| 89 /// This may fail with an error if the bound is for some reason incalculable. | |
| 90 fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> DynResult<F>; | |
| 88 } | 91 } |
| 89 | 92 |
| 90 // Linear operator application into mutable target. The [`AsRef`] bound | 93 // Linear operator application into mutable target. The [`AsRef`] bound |
| 91 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; | 94 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; |
| 92 // the former is assumed to be e.g. a view into the latter. | 95 // the former is assumed to be e.g. a view into the latter. |
| 181 where | 184 where |
| 182 X: Space + Clone + Norm<F, E>, | 185 X: Space + Clone + Norm<F, E>, |
| 183 F: Num, | 186 F: Num, |
| 184 E: NormExponent, | 187 E: NormExponent, |
| 185 { | 188 { |
| 186 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> F { | 189 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> { |
| 187 F::ONE | 190 Ok(F::ONE) |
| 188 } | 191 } |
| 189 } | 192 } |
| 190 | 193 |
| 191 impl<X: Clone + Space> Adjointable<X, X> for IdOp<X> { | 194 impl<X: Clone + Space> Adjointable<X, X> for IdOp<X> { |
| 192 type AdjointCodomain = X; | 195 type AdjointCodomain = X; |
| 265 Y: AXPY<F> + Clone + Norm<F, E2>, | 268 Y: AXPY<F> + Clone + Norm<F, E2>, |
| 266 F: Num, | 269 F: Num, |
| 267 E1: NormExponent, | 270 E1: NormExponent, |
| 268 E2: NormExponent, | 271 E2: NormExponent, |
| 269 { | 272 { |
| 270 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> F { | 273 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> { |
| 271 F::ZERO | 274 Ok(F::ZERO) |
| 272 } | 275 } |
| 273 } | 276 } |
| 274 | 277 |
| 275 impl<'a, F: Num, X, XD, Y, Yprime: Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F> | 278 impl<'a, F: Num, X, XD, Y, Yprime: Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F> |
| 276 where | 279 where |
| 349 Yexp: NormExponent, | 352 Yexp: NormExponent, |
| 350 Zexp: NormExponent, | 353 Zexp: NormExponent, |
| 351 T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>, | 354 T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>, |
| 352 S: BoundedLinear<Z, Zexp, Yexp, F>, | 355 S: BoundedLinear<Z, Zexp, Yexp, F>, |
| 353 { | 356 { |
| 354 fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> F { | 357 fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> DynResult<F> { |
| 355 let zexp = self.intermediate_norm_exponent; | 358 let zexp = self.intermediate_norm_exponent; |
| 356 self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp) | 359 Ok(self.outer.opnorm_bound(zexp, yexp)? * self.inner.opnorm_bound(xexp, zexp)?) |
| 357 } | 360 } |
| 358 } | 361 } |
| 359 | 362 |
| 360 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. | 363 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. |
| 361 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)] | 364 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)] |
| 689 { | 692 { |
| 690 fn opnorm_bound( | 693 fn opnorm_bound( |
| 691 &self, | 694 &self, |
| 692 PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>, | 695 PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>, |
| 693 expr: ExpR, | 696 expr: ExpR, |
| 694 ) -> F { | 697 ) -> DynResult<F> { |
| 695 // An application of the triangle inequality bounds the norm by the maximum | 698 // An application of the triangle inequality bounds the norm by the maximum |
| 696 // of the individual norms. A simple observation shows this to be exact. | 699 // of the individual norms. A simple observation shows this to be exact. |
| 697 let na = self.0.opnorm_bound(expa, expr); | 700 let na = self.0.opnorm_bound(expa, expr)?; |
| 698 let nb = self.1.opnorm_bound(expb, expr); | 701 let nb = self.1.opnorm_bound(expb, expr)?; |
| 699 na.max(nb) | 702 Ok(na.max(nb)) |
| 700 } | 703 } |
| 701 } | 704 } |
| 702 | 705 |
| 703 impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> | 706 impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> |
| 704 for ColOp<S, T> | 707 for ColOp<S, T> |
| 713 { | 716 { |
| 714 fn opnorm_bound( | 717 fn opnorm_bound( |
| 715 &self, | 718 &self, |
| 716 expa: ExpA, | 719 expa: ExpA, |
| 717 PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>, | 720 PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>, |
| 718 ) -> F { | 721 ) -> DynResult<F> { |
| 719 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, | 722 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, |
| 720 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} | 723 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} |
| 721 let ns = self.0.opnorm_bound(expa, exps); | 724 let ns = self.0.opnorm_bound(expa, exps)?; |
| 722 let nt = self.1.opnorm_bound(expa, expt); | 725 let nt = self.1.opnorm_bound(expa, expt)?; |
| 723 ns.max(nt) | 726 Ok(ns.max(nt)) |
| 724 } | 727 } |
| 725 } | 728 } |
| 726 }; | 729 }; |
| 727 } | 730 } |
| 728 | 731 |
