diff --git a/Math/LinearMap/Category.hs b/Math/LinearMap/Category.hs --- a/Math/LinearMap/Category.hs +++ b/Math/LinearMap/Category.hs @@ -196,27 +196,6 @@ --- | For real matrices, this boils down to 'transpose'. --- For free complex spaces it also incurs complex conjugation. --- --- The signature can also be understood as --- --- @ --- adjoint :: (v +> w) -> (DualVector w +> DualVector v) --- @ --- --- Or --- --- @ --- adjoint :: (DualVector v +> DualVector w) -> (w +> v) --- @ --- --- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is --- equivalent, via 'riesz' isomorphism). -adjoint :: (LSpace v, LSpace w, Scalar v ~ Scalar w) - => (v +> DualVector w) -+> (w +> DualVector v) -adjoint = arr fromTensor . transposeTensor . arr asTensor - diff --git a/Math/VectorSpace/Docile.hs b/Math/VectorSpace/Docile.hs --- a/Math/VectorSpace/Docile.hs +++ b/Math/VectorSpace/Docile.hs @@ -34,6 +34,7 @@ import Data.Set (Set) import Data.Ord (comparing) import Data.List (maximumBy, unfoldr) +import qualified Data.Vector as Arr import Data.Foldable (toList) import Data.Semigroup @@ -91,15 +92,15 @@ -- the functional list. -> ([(Int,v)] -> Forest (Int, DualVector v)) -- ^ Suitable definition of 'dualBasisCandidates'. -cartesianDualBasisCandidates dvs abss vcas = go 0 sorted +cartesianDualBasisCandidates dvs abss vcas = go 0 0 sorted where sorted = sortBy (comparing $ negate . snd . snd) [ (i, (av, maximum av)) | (i,v)<-vcas, let av = abss v ] - go k ((i,(av,_)):scs) - | k [a] -> [a] + bringToFront i l = case splitAt i l of + (_,[]) -> [] + (f,s:l') -> s : f++l' instance (Fractional'' s, SemiInner s) => SemiInner (ZeroDim s) where dualBasisCandidates _ = [] @@ -117,26 +123,36 @@ (<.>^) :: LSpace v => DualVector v -> v -> Scalar v f<.>^v = (applyDualVector$f)$v -orthonormaliseDuals :: (SemiInner v, LSpace v, Fractional'' (Scalar v)) - => [(v, DualVector v)] -> [(v,DualVector v)] -orthonormaliseDuals [] = [] -orthonormaliseDuals ((v,v'₀):ws) - = (v,v') : [(w, w' ^-^ (w'<.>^v)*^v') | (w,w')<-wssys] - where wssys = orthonormaliseDuals ws - v'₁ = foldl' (\v'i (w,w') -> v'i ^-^ (v'i<.>^w)*^w') v'₀ wssys - v' = v'₁ ^/ (v'₁<.>^v) +orthonormaliseDuals :: (SemiInner v, LSpace v, RealFrac' (Scalar v)) + => Scalar v -> [(v, DualVector v)] -> [(v,DualVector v)] +orthonormaliseDuals _ [] = [] +orthonormaliseDuals ε ((v,v'₀):ws) + | abs ovl > ε = (v,v') : [(w, w' ^-^ (w'<.>^v)*^v') | (w,w')<-wssys] + | otherwise = (v,zeroV) : wssys + where wssys = orthonormaliseDuals ε ws + v'₁ = foldl' (\v'i (w,w') -> v'i ^-^ (v'i<.>^w)*^w') (v'₀ ^/ (v'₀<.>^v)) wssys + v' = v'₁ ^/ ovl + ovl = v'₁<.>^v -dualBasis :: (SemiInner v, LSpace v, Fractional'' (Scalar v)) => [v] -> [DualVector v] -dualBasis vs = snd <$> orthonormaliseDuals (zip' vsIxed candidates) +dualBasis :: (SemiInner v, LSpace v, RealFrac' (Scalar v)) => [v] -> [DualVector v] +dualBasis vs = snd <$> orthonormaliseDuals epsilon (zip' vsIxed candidates) where zip' ((i,v):vs) ((j,v'):ds) | i^(lookupArr Arr.! i) /= 0 + , Just best' <- findBest (n-1) bv' + = Just $ (i,v') : best' + | otherwise = findBest n alts vsIxed = zip [0..] vs + lookupArr = Arr.fromList vs + n = Arr.length lookupArr instance SemiInner ℝ where dualBasisCandidates = fmap ((`Node`[]) . second recip) @@ -181,16 +197,13 @@ combineBaseis _ forbidden ([], bv) = combineBaseis True forbidden ([],bv) -instance ∀ s u v . ( LSpace u, FiniteDimensional (DualVector u), SemiInner (DualVector u) - , SemiInner v, FiniteDimensional v - , Scalar u ~ s, Scalar v ~ s, RealFrac' s ) +instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s ) => SemiInner (Tensor s u v) where dualBasisCandidates = map (fmap (second $ arr transposeTensor . arr asTensor)) . dualBasisCandidates . map (second $ arr asLinearMap) -instance ∀ s u v . ( SemiInner u, FiniteDimensional u, Scalar u ~ s - , SemiInner v, FiniteDimensional v, Scalar v ~ s, RealFrac' s ) +instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s ) => SemiInner (LinearMap s u v) where dualBasisCandidates = sequenceForest . map (second pseudoInverse) -- this is not efficient @@ -254,7 +267,6 @@ -- library). uncanonicallyFromDual :: DualVector v -+> v uncanonicallyToDual :: v -+> DualVector v - instance (Num''' s) => FiniteDimensional (ZeroDim s) where @@ -307,7 +319,7 @@ #define FreeFiniteDimensional(V, VB, dimens, take, give) \ instance (Num''' s, LSpace s) \ => FiniteDimensional (V s) where { \ - data SubBasis (V s) = VB; \ + data SubBasis (V s) = VB deriving (Show); \ entireBasis = VB; \ enumerateSubBasis VB = toList $ Mat.identity; \ subbasisDimension VB = dimens; \ @@ -392,7 +404,7 @@ = [ u⊗v | u <- enumerateSubBasis bu, v <- enumerateSubBasis bv ] subbasisDimension (TensorBasis bu bv) = subbasisDimension bu * subbasisDimension bv decomposeLinMap muvw = case decomposeLinMap $ curryLinearMap $ muvw of - (bu, mvwsg) -> first (TensorBasis bu) . go id $ mvwsg [] + (bu, mvwsg) -> first (TensorBasis bu) . go $ mvwsg [] where (go, _) = tensorLinmapDecompositionhelpers decomposeLinMapWithin (TensorBasis bu bv) muvw = case decomposeLinMapWithin bu $ curryLinearMap $ muvw of @@ -418,13 +430,13 @@ tensorLinmapDecompositionhelpers :: ( FiniteDimensional v, LSpace w , Scalar v~s, Scalar w~s ) - => ( DList w -> [v+>w] -> (SubBasis v, DList w) + => ( [v+>w] -> (SubBasis v, DList w) , SubBasis v -> DList w -> [v+>w] -> DList (v+>w) -> (Bool, (SubBasis v, DList w)) ) tensorLinmapDecompositionhelpers = (go, goWith) - where go _ [] = decomposeLinMap zeroV - go prevdc (mvw:mvws) = case decomposeLinMap mvw of - (bv, cfs) -> snd (goWith bv prevdc mvws (mvw:)) + where go [] = decomposeLinMap zeroV + go (mvw:mvws) = case decomposeLinMap mvw of + (bv, cfs) -> snd (goWith bv cfs mvws (mvw:)) goWith bv prevdc [] prevs = (False, (bv, prevdc)) goWith bv prevdc (mvw:mvws) prevs = case decomposeLinMapWithin bv mvw of Right cfs -> goWith bv (prevdc . cfs) mvws (prevs . (mvw:)) @@ -439,6 +451,9 @@ \\nif it cannot decompose in the given basis, do so in a proper\ \\nsuperbasis of the given one (so that any vector that could be\ \\ndecomposed in the old basis can also be decomposed in the new one)." + +deriving instance (Show (SubBasis u), Show (SubBasis v)) + => Show (SubBasis (Tensor s u v)) instance ∀ s u v . @@ -475,7 +490,10 @@ uncanonicallyFromDual = fmap uncanonicallyFromDual >>> arr asTensor >>> transposeTensor >>> arr fromTensor >>> fmap uncanonicallyFromDual +deriving instance (Show (SubBasis (DualVector u)), Show (SubBasis v)) + => Show (SubBasis (LinearMap s u v)) + infixr 0 \$ -- | Inverse function application, aka solving of a linear system: @@ -486,8 +504,10 @@ -- f '$' f '\$' u ≡ u -- @ -- --- If @f@ does not have full rank, the behaviour is undefined (but we expect --- it to be reasonably well-behaved or even give a least-squares solution). +-- If @f@ does not have full rank, the behaviour is undefined. However, it +-- does not need to be a proper isomorphism: the +-- first of the above equations is still fulfilled if only @f@ is /injective/ +-- (overdetermined system) and the second if it is /surjective/. -- -- If you want to solve for multiple RHS vectors, be sure to partially -- apply this operator to the linear map, like @@ -502,18 +522,53 @@ -- @ -- [f '\$' v₁, f '\$' v₂, ...] -- @ -(\$) :: ( FiniteDimensional u, FiniteDimensional v, SemiInner v - , Scalar u ~ Scalar v, Fractional' (Scalar v) ) +(\$) :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v -> u -(\$) m = fst . \v -> recomposeSB mbas [v'<.>^v | v' <- v's] - where v's = dualBasis $ mdecomp [] - (mbas, mdecomp) = decomposeLinMap m +(\$) m + | du > dv = (unsafeRightInverse m $) + | du < dv = (unsafeLeftInverse m $) + | otherwise = let v's = dualBasis $ mdecomp [] + (mbas, mdecomp) = decomposeLinMap m + in fst . \v -> recomposeSB mbas [v'<.>^v | v' <- v's] + where du = subbasisDimension (entireBasis :: SubBasis u) + dv = subbasisDimension (entireBasis :: SubBasis v) -pseudoInverse :: ( FiniteDimensional u, FiniteDimensional v, SemiInner v - , Scalar u ~ Scalar v, Fractional' (Scalar v) ) +pseudoInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v+>u -pseudoInverse m = recomposeContraLinMap (fst . recomposeSB mbas) v's +pseudoInverse m + | du > dv = unsafeRightInverse m + | du < dv = unsafeLeftInverse m + | otherwise = unsafeInverse m + where du = subbasisDimension (entireBasis :: SubBasis u) + dv = subbasisDimension (entireBasis :: SubBasis v) + +-- | If @f@ is injective, then +-- +-- @ +-- unsafeLeftInverse f . f ≡ id +-- @ +unsafeLeftInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) + => (u+>v) -> v+>u +unsafeLeftInverse m = unsafeInverse (m' . (fmap uncanonicallyToDual $ m)) + . m' . arr uncanonicallyToDual + where m' = adjoint $ m :: DualVector v +> DualVector u + +-- | If @f@ is surjective, then +-- +-- @ +-- f . unsafeRightInverse f ≡ id +-- @ +unsafeRightInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) + => (u+>v) -> v+>u +unsafeRightInverse m = (fmap uncanonicallyToDual $ m') + . unsafeInverse (m . (fmap uncanonicallyToDual $ m')) + where m' = adjoint $ m :: DualVector v +> DualVector u + +-- | Invert an isomorphism. For other linear maps, the result is undefined. +unsafeInverse :: ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) + => (u+>v) -> v+>u +unsafeInverse m = recomposeContraLinMap (fst . recomposeSB mbas) v's where v's = dualBasis $ mdecomp [] (mbas, mdecomp) = decomposeLinMap m @@ -525,7 +580,7 @@ let (bas, compos) = decomposeLinMap $ sampleLinearFunction $ applyDualVector $ dv in fst . recomposeSB bas $ compos [] -sRiesz :: (FiniteDimensional v, InnerSpace v) => DualSpace v -+> v +sRiesz :: FiniteDimensional v => DualSpace v -+> v sRiesz = LinearFunction $ \dv -> let (bas, compos) = decomposeLinMap $ dv in fst . recomposeSB bas $ compos [] @@ -543,6 +598,7 @@ . showsPrec 7 (sRiesz$dv) instance Show (LinearMap ℝ (V0 ℝ) ℝ) where showsPrec = showsPrecAsRiesz +instance Show (LinearMap ℝ ℝ ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V1 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V2 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V3 ℝ) ℝ) where showsPrec = showsPrecAsRiesz @@ -561,22 +617,22 @@ instance Show (LinearMap s v (V0 s)) where show _ = "zeroV" -instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v) +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V1 ℝ)) where showsPrec p m = showParen (p>6) $ ("ex .< "++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) -instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v) +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V2 ℝ)) where showsPrec p m = showParen (p>6) $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m) -instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v) +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V3 ℝ)) where showsPrec p m = showParen (p>6) $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m) . (" ^+^ ez.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._z)) $ m) -instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v) +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V4 ℝ)) where showsPrec p m = showParen (p>6) $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) @@ -635,3 +691,24 @@ unsafeFromFullUnboxVect arrv = arr (unsafeFromFullUnboxVect arrv :: LinearMap s u v) + +-- | For real matrices, this boils down to 'transpose'. +-- For free complex spaces it also incurs complex conjugation. +-- +-- The signature can also be understood as +-- +-- @ +-- adjoint :: (v +> w) -> (DualVector w +> DualVector v) +-- @ +-- +-- Or +-- +-- @ +-- adjoint :: (DualVector v +> DualVector w) -> (w +> v) +-- @ +-- +-- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is +-- equivalent, via 'riesz' isomorphism). +adjoint :: (LSpace v, LSpace w, Scalar v ~ Scalar w) + => (v +> DualVector w) -+> (w +> DualVector v) +adjoint = arr fromTensor . transposeTensor . arr asTensor diff --git a/linearmap-category.cabal b/linearmap-category.cabal --- a/linearmap-category.cabal +++ b/linearmap-category.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: linearmap-category -version: 0.1.0.0 +version: 0.1.0.1 synopsis: Native, complete, matrix-free linear algebra. description: The term /numerical linear algebra/ is often used almost synonymous with /matrix modifications/. However, what's interesting