module Math.LinearMap.Category (
LinearFunction (..), type (-+>)(), Bilinear
, lfun
, LinearMap (..), type (+>)()
, (⊕), (>+<)
, adjoint
, (<.>^), (-+|>)
, Tensor (..), type (⊗)(), (⊗)
, SymmetricTensor(..), squareV, squareVs
, type (⊗〃+>)(), currySymBilin
, Norm(..), Seminorm
, spanNorm
, euclideanNorm
, (|$|)
, normSq
, (<$|)
, scaleNorm
, normSpanningSystem
, normSpanningSystem'
, Variance, spanVariance, (|&>), varianceSpanningSystem
, dualNorm, dualNorm', dependence
, densifyNorm, wellDefinedNorm
, (\$), pseudoInverse, roughDet
, linearRegressionW, linearRegressionWVar
, eigen
, constructEigenSystem
, roughEigenSystem
, finishEigenSystem
, Eigenvector(..)
, LSpace
, TensorSpace (..)
, LinearSpace (..)
, SemiInner (..), cartesianDualBasisCandidates, embedFreeSubspace
, FiniteDimensional (..)
, addV, scale, inner, flipBilin, bilinearFunction
, (.⊗)
, DualSpace, riesz, coRiesz, showsPrecAsRiesz, (.<)
, HilbertSpace, SimpleSpace
, Num'(..)
, Fractional'
, RealFrac', RealFloat', LinearShowable
, ClosedScalarWitness(..), ScalarSpaceWitness(..), DualSpaceWitness(..)
, LinearManifoldWitness(..)
, relaxNorm, transformNorm, transformVariance
, findNormalLength, normalLength
, summandSpaceNorms, sumSubspaceNorms
, sharedNormSpanningSystem, sharedSeminormSpanningSystem
, sharedSeminormSpanningSystem'
, convexPolytopeHull
, convexPolytopeRepresentatives
) where
import Math.LinearMap.Category.Class
import Math.LinearMap.Category.Instances
import Math.LinearMap.Asserted
import Math.VectorSpace.Docile
import Data.Tree (Tree(..), Forest)
import Data.List (sortBy, foldl')
import qualified Data.Set as Set
import Data.Set (Set)
import Data.Ord (comparing)
import Data.List (maximumBy)
import Data.Maybe (catMaybes)
import Data.Foldable (toList)
import Data.Semigroup
import Data.VectorSpace
import Data.Basis
import Prelude ()
import qualified Prelude as Hask
import Control.Category.Constrained.Prelude hiding ((^))
import Control.Arrow.Constrained
import Linear ( V0(V0), V1(V1), V2(V2), V3(V3), V4(V4)
, _x, _y, _z, _w )
import Data.VectorSpace.Free
import Math.VectorSpace.ZeroDimensional
import qualified Linear.Matrix as Mat
import qualified Linear.Vector as Mat
import Control.Lens ((^.))
import qualified Data.Vector.Unboxed as UArr
import Numeric.IEEE
import qualified GHC.Exts as GHC
import qualified Data.Type.Coercion as GHC
infixr 7 -+|>
(-+|>) :: ( EnhancedCat f (LinearFunction s)
, LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s
, Object f u, Object f v )
=> DualVector u -> v -> f u v
du-+|>v = arr . LinearFunction $ (v^*) . (du<.>^)
spanNorm :: ∀ v . LSpace v => [DualVector v] -> Seminorm v
spanNorm = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness
-> \dvs -> Norm . LinearFunction $ \v -> sumV [dv ^* (dv<.>^v) | dv <- dvs]
spanVariance :: ∀ v . LSpace v => [v] -> Variance v
spanVariance = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness -> spanNorm
relaxNorm :: ∀ v . SimpleSpace v => Norm v -> [v] -> Norm v
relaxNorm = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness
-> \me vs -> let vs' = normSpanningSystem' me
in dualNorm . spanVariance $ vs' ++ vs
scaleNorm :: ∀ v . LSpace v => Scalar v -> Norm v -> Norm v
scaleNorm = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness -> \μ (Norm n) -> Norm $ μ^2 *^ n
newtype Norm v = Norm {
applyNorm :: v -+> DualVector v
}
type Seminorm v = Norm v
instance LSpace v => Semigroup (Norm v) where
Norm m <> Norm n = Norm $ m^+^n
instance LSpace v => Monoid (Seminorm v) where
mempty = Norm zeroV
mappend = (<>)
type Variance v = Norm (DualVector v)
euclideanNorm :: HilbertSpace v => Norm v
euclideanNorm = Norm id
adhocNorm :: FiniteDimensional v => Norm v
adhocNorm = Norm uncanonicallyToDual
dualNorm :: SimpleSpace v => Norm v -> Variance v
dualNorm = spanVariance . normSpanningSystem'
dualNorm' :: ∀ v . SimpleSpace v => Variance v -> Norm v
dualNorm' = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness -> spanNorm . normSpanningSystem'
transformNorm :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w)
=> (v+>w) -> Norm w -> Norm v
transformNorm = case ( dualSpaceWitness :: DualSpaceWitness v
, dualSpaceWitness :: DualSpaceWitness w ) of
(DualSpaceWitness, DualSpaceWitness)
-> \f (Norm m) -> Norm . arr $ (adjoint $ f) . (fmap m $ f)
transformVariance :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w)
=> (v+>w) -> Variance v -> Variance w
transformVariance = case ( dualSpaceWitness :: DualSpaceWitness v
, dualSpaceWitness :: DualSpaceWitness w ) of
(DualSpaceWitness, DualSpaceWitness)
-> \f (Norm m) -> Norm . arr $ f . (fmap m $ adjoint $ f)
infixl 6 ^%
(^%) :: (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
v ^% Norm m = v ^/ sqrt ((m-+$>v)<.>^v)
findNormalLength :: ∀ s . RealFrac' s => Norm s -> Maybe s
findNormalLength (Norm m) = case ( closedScalarWitness :: ClosedScalarWitness s
, m-+$>1 ) of
(ClosedScalarWitness, o) | o > 0 -> Just . sqrt $ recip o
| otherwise -> Nothing
normalLength :: ∀ s . RealFrac' s => Norm s -> s
normalLength (Norm m) = case ( closedScalarWitness :: ClosedScalarWitness s
, m-+$>1 ) of
(ClosedScalarWitness, o) | o >= 0 -> sqrt $ recip o
| o < 0 -> error "Norm fails to be positive semidefinite."
| otherwise -> error "Norm yields NaN."
infixr 0 <$|, |$|
(<$|) :: LSpace v => Norm v -> v -> DualVector v
Norm m <$| v = m-+$>v
normSq :: LSpace v => Seminorm v -> v -> Scalar v
normSq (Norm m) v = (m-+$>v)<.>^v
(|$|) :: (LSpace v, Floating (Scalar v)) => Seminorm v -> v -> Scalar v
(|$|) m = sqrt . normSq m
infixl 1 |&>
(|&>) :: LSpace v => DualVector v -> Variance v -> v
dv |&> Norm m = GHC.sym coerceDoubleDual $ m-+$>dv
densifyNorm :: ∀ v . LSpace v => Norm v -> Norm v
densifyNorm = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness
-> \(Norm m) -> Norm . arr $ sampleLinearFunction $ m
wellDefinedNorm :: ∀ v . LinearSpace v => Norm v -> Maybe (Norm v)
wellDefinedNorm = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness
-> \(Norm m) -> Norm <$> wellDefinedVector m
data OrthonormalSystem v = OrthonormalSystem {
orthonormalityNorm :: Norm v
, orthonormalVectors :: [v]
}
orthonormaliseFussily :: ∀ v . (LSpace v, RealFloat (Scalar v))
=> Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily = onf dualSpaceWitness
where onf :: DualSpaceWitness v -> Scalar v -> Norm v -> [v] -> [v]
onf DualSpaceWitness fuss me = go []
where go _ [] = []
go ws (v₀:vs)
| mvd > fuss = let μ = 1/sqrt mvd
v = vd^*μ
in v : go ((v,dvd^*μ):ws) vs
| otherwise = go ws vs
where vd = orthogonalComplementProj' ws $ v₀
dvd = applyNorm me $ vd
mvd = dvd<.>^vd
orthogonalComplementProj' :: LSpace v => [(v, DualVector v)] -> (v-+>v)
orthogonalComplementProj' ws = LinearFunction $ \v₀
-> foldl' (\v (w,dw) -> v ^-^ w^*(dw<.>^v)) v₀ ws
orthogonalComplementProj :: LSpace v => Norm v -> [v] -> (v-+>v)
orthogonalComplementProj (Norm m)
= orthogonalComplementProj' . map (id &&& (m-+$>))
data Eigenvector v = Eigenvector {
ev_Eigenvalue :: Scalar v
, ev_Eigenvector :: v
, ev_FunctionApplied :: v
, ev_Deviation :: v
, ev_Badness :: Scalar v
}
deriving instance (Show v, Show (Scalar v)) => Show (Eigenvector v)
constructEigenSystem :: (LSpace v, RealFloat (Scalar v))
=> Norm v
-> Scalar v
-> (v-+>v)
-> [v]
-> [[Eigenvector v]]
constructEigenSystem me ε₀ f = iterate (
sortBy (comparing $
negate . abs . ev_Eigenvalue)
. map asEV
. orthonormaliseFussily (1/4) me
. newSys)
. map (asEV . (^%me))
where newSys [] = []
newSys (Eigenvector λ v fv dv ε : evs)
| ε>ε₀ = case newSys evs of
[] -> [fv^/λ, dv^/sqrt(ε+ε₀)]
vn:vns -> fv^/λ : vn : dv^/sqrt(ε+ε₀) : vns
| ε>=0 = v : newSys evs
| otherwise = newSys evs
asEV v = Eigenvector λ v fv dv ε
where λ² = fv'<.>^fv
λ = fv'<.>^v
ε = normSq me dv
fv = f $ v
fv' = me<$|fv
dv | λ²>0 = v ^-^ fv^*(λ/λ²)
| otherwise = zeroV
finishEigenSystem :: ∀ v . (LSpace v, RealFloat (Scalar v))
=> Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem me = go . sortBy (comparing $ negate . ev_Eigenvalue)
where go [] = []
go [v] = [v]
go vs@[Eigenvector λ₀ v₀ fv₀ _dv₀ _ε₀, Eigenvector λ₁ v₁ fv₁ _dv₁ _ε₁]
| λ₀>λ₁ = [ asEV v₀' fv₀', asEV v₁' fv₁' ]
| otherwise = vs
where
v₀' = v₀^*μ₀₀ ^+^ v₁^*μ₀₁
fv₀' = fv₀^*μ₀₀ ^+^ fv₁^*μ₀₁
v₁' = v₀^*μ₁₀ ^+^ v₁^*μ₁₁
fv₁' = fv₀^*μ₁₀ ^+^ fv₁^*μ₁₁
fShift₁v₀ = fv₀ ^-^ λ₁*^v₀
(μ₀₀,μ₀₁) = normalised ( λ₀ λ₁
, (me <$| fShift₁v₀)<.>^v₁ )
(μ₁₀,μ₁₁) = (μ₀₁, μ₀₀)
go evs = lo'' ++ upper'
where l = length evs
lChunk = l`quot`3
(loEvs, (midEvs, hiEvs)) = second (splitAt $ l 2*lChunk)
$ splitAt lChunk evs
(lo',hi') = splitAt lChunk . go $ loEvs++hiEvs
(lo'',mid') = splitAt lChunk . go $ lo'++midEvs
upper' = go $ mid'++hi'
asEV v fv = Eigenvector λ v fv dv ε
where λ = (me<$|v)<.>^fv
dv = v^*λ ^-^ fv
ε = normSq me dv / λ^2
normalised (x,y) = (x/r, y/r)
where r = sqrt $ x^2 + y^2
roughEigenSystem :: (FiniteDimensional v, IEEE (Scalar v))
=> Norm v
-> (v+>v)
-> [Eigenvector v]
roughEigenSystem me f = go fBas 0 [[]]
where go [] _ (_:evs:_) = evs
go [] _ (evs:_) = evs
go (v:vs) oldDim (evs:evss)
| normSq me vPerp > fpε = case evss of
evs':_ | length evs' > oldDim
-> go (v:vs) (length evs) evss
_ -> let evss' = tail . constructEigenSystem me fpε (arr f)
$ map ev_Eigenvector (head $ evss++[evs]) ++ [vPerp]
in go vs (length evs) evss'
| otherwise = go vs oldDim (evs:evss)
where vPerp = orthogonalComplementProj me (ev_Eigenvector<$>evs) $ v
fBas = (^%me) <$> snd (decomposeLinMap id) []
fpε = epsilon * 8
eigen :: (FiniteDimensional v, HilbertSpace v, IEEE (Scalar v))
=> (v+>v) -> [(Scalar v, v)]
eigen f = map (ev_Eigenvalue &&& ev_Eigenvector)
$ iterate (finishEigenSystem euclideanNorm) (roughEigenSystem euclideanNorm f) !! 2
roughDet :: (FiniteDimensional v, IEEE (Scalar v)) => (v+>v) -> Scalar v
roughDet = roughEigenSystem adhocNorm >>> map ev_Eigenvalue >>> product
orthonormalityError :: LSpace v => Norm v -> [v] -> Scalar v
orthonormalityError me vs = normSq me $ orthogonalComplementProj me vs $ sumV vs
normSpanningSystem :: SimpleSpace v
=> Seminorm v -> [DualVector v]
normSpanningSystem me@(Norm m)
= catMaybes . map snd . orthonormaliseDuals 0
. map (id&&&(m-+$>)) $ normSpanningSystem' me
normSpanningSystem' :: (FiniteDimensional v, IEEE (Scalar v))
=> Seminorm v -> [v]
normSpanningSystem' me = orthonormaliseFussily 0 me $ enumerateSubBasis entireBasis
varianceSpanningSystem :: ∀ v . SimpleSpace v => Variance v -> [v]
varianceSpanningSystem = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness -> normSpanningSystem
sharedNormSpanningSystem :: SimpleSpace v
=> Norm v -> Seminorm v -> [(DualVector v, Scalar v)]
sharedNormSpanningSystem nn@(Norm n) nm
= first (n-+$>) <$> sharedNormSpanningSystem' 0 (nn, dualNorm nn) nm
sharedNormSpanningSystem' :: ∀ v . SimpleSpace v
=> Int -> (Norm v, Variance v) -> Seminorm v -> [(v, Scalar v)]
sharedNormSpanningSystem' = snss dualSpaceWitness
where snss :: DualSpaceWitness v -> Int -> (Norm v, Variance v)
-> Seminorm v -> [(v, Scalar v)]
snss DualSpaceWitness nRefine (nn@(Norm n), Norm n') (Norm m)
= sep =<< iterate (finishEigenSystem nn)
(roughEigenSystem nn $ arr n' . arr m) !! nRefine
sep (Eigenvector λ v λv _ _)
| λ>=0 = [(v, sqrt λ)]
| otherwise = []
sharedSeminormSpanningSystem :: ∀ v . SimpleSpace v
=> Seminorm v -> Seminorm v -> [(DualVector v, Maybe (Scalar v))]
sharedSeminormSpanningSystem nn nm
= finalise dualSpaceWitness
<$> sharedNormSpanningSystem' 1 (combined, dualNorm combined) nn
where combined = densifyNorm $ nn<>nm
finalise :: DualSpaceWitness v -> (v, Scalar v) -> (DualVector v, Maybe (Scalar v))
finalise DualSpaceWitness (v, μn)
| μn^2 > epsilon = (v'^*μn, Just $ sqrt (max 0 $ 1 μn^2)/μn)
| otherwise = (v', Nothing)
where v' = combined<$|v
sharedSeminormSpanningSystem' :: ∀ v . SimpleSpace v
=> Seminorm v -> Seminorm v -> [v]
sharedSeminormSpanningSystem' nn nm
= fst <$> sharedNormSpanningSystem' 1 (combined, dualNorm combined) nn
where combined = densifyNorm $ nn<>nm
dependence :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u~Scalar v)
=> Variance (u,v) -> (u+>v)
dependence = case ( dualSpaceWitness :: DualSpaceWitness u
, dualSpaceWitness :: DualSpaceWitness v ) of
(DualSpaceWitness,DualSpaceWitness)
-> \(Norm m) -> fmap ( snd . m . (id&&&zeroV) )
$ pseudoInverse (arr $ fst . m . (id&&&zeroV))
summandSpaceNorms :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v)
=> Norm (u,v) -> (Norm u, Norm v)
summandSpaceNorms = case ( dualSpaceWitness :: DualSpaceWitness u
, dualSpaceWitness :: DualSpaceWitness v ) of
(DualSpaceWitness,DualSpaceWitness)
-> \nuv -> let spanSys = normSpanningSystem nuv
in ( densifyNorm $ spanNorm (fst<$>spanSys)
, densifyNorm $ spanNorm (snd<$>spanSys) )
sumSubspaceNorms :: ∀ u v . (LSpace u, LSpace v, Scalar u~Scalar v)
=> Norm u -> Norm v -> Norm (u,v)
sumSubspaceNorms = case ( dualSpaceWitness :: DualSpaceWitness u
, dualSpaceWitness :: DualSpaceWitness v ) of
(DualSpaceWitness,DualSpaceWitness)
-> \(Norm nu) (Norm nv) -> Norm $ nu *** nv
instance (SimpleSpace v, Show (DualVector v)) => Show (Norm v) where
showsPrec p n = showParen (p>9) $ ("spanNorm "++) . shows (normSpanningSystem n)
type LinearShowable v = (Show v, RieszDecomposable v)
convexPolytopeHull :: ∀ v . SimpleSpace v => [v] -> [DualVector v]
convexPolytopeHull vs = case dualSpaceWitness :: DualSpaceWitness v of
DualSpaceWitness
-> [dv^/η | (dv,η) <- candidates, all ((<=η) . (dv<.>^)) vs]
where vrv = spanVariance vs
nmv = dualNorm' vrv
candidates :: [(DualVector v, Scalar v)]
candidates = [ (dv, dv<.>^v) | v <- vs
, let dv = nmv<$|v ]
convexPolytopeRepresentatives :: ∀ v . SimpleSpace v => [DualVector v] -> [v]
convexPolytopeRepresentatives dvs
= [v^/η | ((v,η),dv) <- zip candidates dvs
, all (\(w,ψ) -> dv<.>^w <= ψ) candidates]
where nmv :: Norm v
nmv = spanNorm dvs
vrv = dualNorm nmv
candidates :: [(v, Scalar v)]
candidates = [ (v, dv<.>^v) | dv <- dvs
, let v = dv|&>vrv ]
linearRegressionW :: ∀ s x m y
. ( LinearSpace x, FiniteDimensional y, SimpleSpace m
, Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s )
=> Norm y -> (x -> (m +> y)) -> [(x,y)] -> m
linearRegressionW σy modelMap = fst . linearRegressionWVar modelMap . map (second (,σy))
linearRegressionWVar :: ∀ s x m y
. ( LinearSpace x, FiniteDimensional y, SimpleSpace m
, Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s )
=> (x -> (m +> y)) -> [(x, (y, Norm y))] -> (m, [DualVector m])
linearRegressionWVar = lrw (dualSpaceWitness, dualSpaceWitness)
where lrw :: (DualSpaceWitness y, DualSpaceWitness m)
-> (x -> (m +> y)) -> [(x, (y, Norm y))] -> (m, [DualVector m])
lrw (DualSpaceWitness, DualSpaceWitness) modelMap dataxy
= ( leastSquareSol, deviations )
where leastSquareSol = (lfun $ forward' . zipWith ((<$|) . snd . snd) dataxy
. forward)
\$ forward' [σy<$|y | (_,(y,σy)) <- dataxy]
forward :: m -> [y]
forward m = [modelMap x $ m | (x,_)<-dataxy]
forward' :: [DualVector y] -> DualVector m
forward' = sumV . zipWith ($) modelGens
modelGens :: [DualVector y +> DualVector m]
modelGens = ((adjoint$) . modelMap . fst)<$>dataxy
deviations = [ m $ dy ^/ ψ
| (m,(dy,ψ)) <- zip modelGens ddys
, ψ > 0
]
ddys = [ (dy, ψ) | (x,(yd,σy)) <- dataxy
, let ym = modelMap x $ leastSquareSol
δy = yd ^-^ ym
dy = σy<$|δy
ψ = dy<.>^δy
]