src/linops.rs

branch
dev
changeset 151
402d717bb5c0
parent 150
c4e394a9c84c
child 152
dab30b331f49
equal deleted inserted replaced
150:c4e394a9c84c 151:402d717bb5c0
34 + Add<Self, Output = Self::Owned> 34 + Add<Self, Output = Self::Owned>
35 + Add<Self::Owned, Output = Self::Owned> 35 + Add<Self::Owned, Output = Self::Owned>
36 + Sub<Self, Output = Self::Owned> 36 + Sub<Self, Output = Self::Owned>
37 + Sub<Self::Owned, Output = Self::Owned> 37 + Sub<Self::Owned, Output = Self::Owned>
38 + Neg 38 + Neg
39 + for<'b> Add<&'b Self, Output = <Self as VectorSpace>::Owned>
40 + for<'b> Sub<&'b Self, Output = <Self as VectorSpace>::Owned>
39 { 41 {
40 type Field: Num; 42 type Field: Num;
41 type Owned: ClosedSpace 43 type Owned: ClosedSpace
42 + AXPY< 44 + AXPY<
43 Self, 45 Self,
67 + DivAssign<Self::Field> 69 + DivAssign<Self::Field>
68 + AddAssign<Self> 70 + AddAssign<Self>
69 + AddAssign<Self::Owned> 71 + AddAssign<Self::Owned>
70 + SubAssign<Self> 72 + SubAssign<Self>
71 + SubAssign<Self::Owned> 73 + SubAssign<Self::Owned>
74 + for<'b> AddAssign<&'b Self>
75 + for<'b> SubAssign<&'b Self>
72 where 76 where
73 X: Space, 77 X: Space,
74 { 78 {
75 /// Computes `y = βy + αx`, where `y` is `Self`. 79 /// Computes `y = βy + αx`, where `y` is `Self`.
76 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field); 80 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field);
86 } 90 }
87 91
88 /// Set self to zero. 92 /// Set self to zero.
89 fn set_zero(&mut self); 93 fn set_zero(&mut self);
90 } 94 }
95
96 pub trait ClosedVectorSpace: Instance<Self> + VectorSpace<Owned = Self> {}
97 impl<X: Instance<X> + VectorSpace<Owned = Self>> ClosedVectorSpace for X {}
91 98
92 /// Efficient in-place application for [`Linear`] operators. 99 /// Efficient in-place application for [`Linear`] operators.
93 #[replace_float_literals(F::cast_from(literal))] 100 #[replace_float_literals(F::cast_from(literal))]
94 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> { 101 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> {
95 /// Computes `y = αAx + βy`, where `A` is `Self`. 102 /// Computes `y = αAx + βy`, where `A` is `Self`.
139 pub trait Adjointable<X, Yʹ>: Linear<X> 146 pub trait Adjointable<X, Yʹ>: Linear<X>
140 where 147 where
141 X: Space, 148 X: Space,
142 Yʹ: Space, 149 Yʹ: Space,
143 { 150 {
144 type AdjointCodomain: Space; 151 type AdjointCodomain: ClosedSpace;
145 type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain> 152 type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
146 where 153 where
147 Self: 'a; 154 Self: 'a;
148 155
149 /// Form the adjoint operator of `self`. 156 /// Form the adjoint operator of `self`.
161 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint 168 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
162 /// have the full space as the codomain, etc. 169 /// have the full space as the codomain, etc.
163 pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>: 170 pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
164 Linear<X> 171 Linear<X>
165 { 172 {
166 type PreadjointCodomain: Space; 173 type PreadjointCodomain: ClosedSpace;
167 type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain> 174 type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain>
168 where 175 where
169 Self: 'a; 176 Self: 'a;
170 177
171 /// Form the adjoint operator of `self`. 178 /// Form the adjoint operator of `self`.
188 IdOp(PhantomData) 195 IdOp(PhantomData)
189 } 196 }
190 } 197 }
191 198
192 impl<X: Space> Mapping<X> for IdOp<X> { 199 impl<X: Space> Mapping<X> for IdOp<X> {
193 type Codomain = X::OwnedVariant; 200 type Codomain = X::OwnedSpace;
194 201
195 fn apply<I: Instance<X>>(&self, x: I) -> X { 202 fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
196 x.own() 203 x.own().into_owned()
197 } 204 }
198 } 205 }
199 206
200 impl<X: Space> Linear<X> for IdOp<X> {} 207 impl<X: Space> Linear<X> for IdOp<X> {}
201 208
224 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> { 231 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> {
225 Ok(F::ONE) 232 Ok(F::ONE)
226 } 233 }
227 } 234 }
228 235
229 impl<X: Clone + Space> Adjointable<X, X> for IdOp<X> { 236 impl<X: Clone + Space> Adjointable<X, X::OwnedSpace> for IdOp<X> {
230 type AdjointCodomain = X; 237 type AdjointCodomain = X::OwnedSpace;
231 type Adjoint<'a> 238 type Adjoint<'a>
232 = IdOp<X> 239 = IdOp<X::OwnedSpace>
233 where 240 where
234 X: 'a; 241 X: 'a;
235 242
236 fn adjoint(&self) -> Self::Adjoint<'_> { 243 fn adjoint(&self) -> Self::Adjoint<'_> {
237 IdOp::new() 244 IdOp::new()
238 } 245 }
239 } 246 }
240 247
241 impl<X: Clone + Space> Preadjointable<X, X> for IdOp<X> { 248 impl<X: Clone + Space> Preadjointable<X, X::OwnedSpace> for IdOp<X> {
242 type PreadjointCodomain = X; 249 type PreadjointCodomain = X::OwnedSpace;
243 type Preadjoint<'a> 250 type Preadjoint<'a>
244 = IdOp<X> 251 = IdOp<X::OwnedSpace>
245 where 252 where
246 X: 'a; 253 X: 'a;
247 254
248 fn preadjoint(&self) -> Self::Preadjoint<'_> { 255 fn preadjoint(&self) -> Self::Preadjoint<'_> {
249 IdOp::new() 256 IdOp::new()
295 302
296 impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp 303 impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
297 where 304 where
298 F: Num, 305 F: Num,
299 X: VectorSpace<Field = F> + HasDual<F>, 306 X: VectorSpace<Field = F> + HasDual<F>,
300 X::DualSpace: VectorSpace<Owned = X::DualSpace>, 307 X::DualSpace: ClosedVectorSpace,
301 { 308 {
302 type AdjointCodomain = X::DualSpace; 309 type AdjointCodomain = X::DualSpace;
303 type Adjoint<'b> 310 type Adjoint<'b>
304 = SimpleZeroOp 311 = SimpleZeroOp
305 where 312 where
458 for ZeroOp<X, Y, OY, OXprime, F> 465 for ZeroOp<X, Y, OY, OXprime, F>
459 where 466 where
460 X: HasDual<F, DualSpace = Xprime>, 467 X: HasDual<F, DualSpace = Xprime>,
461 Y: HasDual<F, DualSpace = Yprime>, 468 Y: HasDual<F, DualSpace = Yprime>,
462 F: Float, 469 F: Float,
463 Xprime: VectorSpace<Field = F, Owned = Xprime>, 470 Xprime: ClosedVectorSpace<Field = F>,
464 Xprime::Owned: AXPY<Field = F>, 471 //Xprime::Owned: AXPY<Field = F>,
465 Yprime: Space + Instance<Yprime>, 472 Yprime: ClosedSpace,
466 OY: OriginGenerator<Y>, 473 OY: OriginGenerator<Y>,
467 OXprime: OriginGenerator<X::DualSpace>, 474 OXprime: OriginGenerator<X::DualSpace>,
468 { 475 {
469 type AdjointCodomain = Xprime; 476 type AdjointCodomain = Xprime;
470 type Adjoint<'c> 477 type Adjoint<'c>
538 A: Space, 545 A: Space,
539 B: Space, 546 B: Space,
540 S: Mapping<A>, 547 S: Mapping<A>,
541 T: Mapping<B>, 548 T: Mapping<B>,
542 S::Codomain: Add<T::Codomain>, 549 S::Codomain: Add<T::Codomain>,
543 <S::Codomain as Add<T::Codomain>>::Output: Space, 550 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
544 { 551 {
545 type Codomain = <S::Codomain as Add<T::Codomain>>::Output; 552 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
546 553
547 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain { 554 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
548 x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b)) 555 x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b))
554 A: Space, 561 A: Space,
555 B: Space, 562 B: Space,
556 S: Linear<A>, 563 S: Linear<A>,
557 T: Linear<B>, 564 T: Linear<B>,
558 S::Codomain: Add<T::Codomain>, 565 S::Codomain: Add<T::Codomain>,
559 <S::Codomain as Add<T::Codomain>>::Output: Space, 566 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
560 { 567 {
561 } 568 }
562 569
563 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> 570 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
564 where 571 where
691 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T> 698 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
692 where 699 where
693 A: Space, 700 A: Space,
694 Xʹ: Space, 701 Xʹ: Space,
695 Yʹ: Space, 702 Yʹ: Space,
696 R: Space + ClosedAdd, 703 R: ClosedSpace + ClosedAdd,
697 S: Adjointable<A, Xʹ, AdjointCodomain = R>, 704 S: Adjointable<A, Xʹ, AdjointCodomain = R>,
698 T: Adjointable<A, Yʹ, AdjointCodomain = R>, 705 T: Adjointable<A, Yʹ, AdjointCodomain = R>,
699 Self: Linear<A>, 706 Self: Linear<A>,
700 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 707 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
701 // Pair<Xʹ,Yʹ>, 708 // Pair<Xʹ,Yʹ>,
716 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T> 723 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
717 where 724 where
718 A: Space, 725 A: Space,
719 Xʹ: Space, 726 Xʹ: Space,
720 Yʹ: Space, 727 Yʹ: Space,
721 R: Space + ClosedAdd, 728 R: ClosedSpace + ClosedAdd,
722 S: Preadjointable<A, Xʹ, PreadjointCodomain = R>, 729 S: Preadjointable<A, Xʹ, PreadjointCodomain = R>,
723 T: Preadjointable<A, Yʹ, PreadjointCodomain = R>, 730 T: Preadjointable<A, Yʹ, PreadjointCodomain = R>,
724 Self: Linear<A>, 731 Self: Linear<A>,
725 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>, 732 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
726 { 733 {
800 where 807 where
801 A: Space, 808 A: Space,
802 B: Space, 809 B: Space,
803 Xʹ: Space, 810 Xʹ: Space,
804 Yʹ: Space, 811 Yʹ: Space,
805 R: Space, 812 R: ClosedSpace,
806 S: Adjointable<A, Xʹ>, 813 S: Adjointable<A, Xʹ>,
807 T: Adjointable<B, Yʹ>, 814 T: Adjointable<B, Yʹ>,
808 Self: Linear<Pair<A, B>>, 815 Self: Linear<Pair<A, B>>,
809 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>, 816 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
810 { 817 {
823 where 830 where
824 A: Space, 831 A: Space,
825 B: Space, 832 B: Space,
826 Xʹ: Space, 833 Xʹ: Space,
827 Yʹ: Space, 834 Yʹ: Space,
828 R: Space, 835 R: ClosedSpace,
829 S: Preadjointable<A, Xʹ>, 836 S: Preadjointable<A, Xʹ>,
830 T: Preadjointable<B, Yʹ>, 837 T: Preadjointable<B, Yʹ>,
831 Self: Linear<Pair<A, B>>, 838 Self: Linear<Pair<A, B>>,
832 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>, 839 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
833 { 840 {
854 A: Space + Norm<ExpA, F>, 861 A: Space + Norm<ExpA, F>,
855 B: Space + Norm<ExpB, F>, 862 B: Space + Norm<ExpB, F>,
856 S: BoundedLinear<A, ExpA, ExpR, F>, 863 S: BoundedLinear<A, ExpA, ExpR, F>,
857 T: BoundedLinear<B, ExpB, ExpR, F>, 864 T: BoundedLinear<B, ExpB, ExpR, F>,
858 S::Codomain: Add<T::Codomain>, 865 S::Codomain: Add<T::Codomain>,
859 <S::Codomain as Add<T::Codomain>>::Output: Space, 866 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
860 ExpA: NormExponent, 867 ExpA: NormExponent,
861 ExpB: NormExponent, 868 ExpB: NormExponent,
862 ExpR: NormExponent, 869 ExpR: NormExponent,
863 { 870 {
864 fn opnorm_bound( 871 fn opnorm_bound(
910 pub struct Scaled<F: Float>(pub F); 917 pub struct Scaled<F: Float>(pub F);
911 918
912 impl<Domain, F> Mapping<Domain> for Scaled<F> 919 impl<Domain, F> Mapping<Domain> for Scaled<F>
913 where 920 where
914 F: Float, 921 F: Float,
915 Domain: Space + ClosedMul<F>, 922 Domain: Space + Mul<F>,
923 <Domain as Mul<F>>::Output: ClosedSpace,
916 { 924 {
917 type Codomain = <Domain as Mul<F>>::Output; 925 type Codomain = <Domain as Mul<F>>::Output;
918 926
919 /// Compute the value of `self` at `x`. 927 /// Compute the value of `self` at `x`.
920 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain { 928 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
923 } 931 }
924 932
925 impl<Domain, F> Linear<Domain> for Scaled<F> 933 impl<Domain, F> Linear<Domain> for Scaled<F>
926 where 934 where
927 F: Float, 935 F: Float,
928 Domain: Space + ClosedMul<F>, 936 Domain: Space + Mul<F>,
929 { 937 <Domain as Mul<F>>::Output: ClosedSpace,
930 } 938 {
939 }

mercurial