{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.Function (
SqRt, sqrt, sqrtSchur, sqrtDenmanBeavers,
Exp, exp, expRealHermitian,
Log, log, logUnipotentUpper,
LiftReal, liftReal,
) where
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Matrix.Lazy.UpperTriangular as LazyUpper
import qualified Numeric.LAPACK.Matrix.Array.Mosaic as ArrMosaic
import qualified Numeric.LAPACK.Matrix.Mosaic.Basic as Mosaic
import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Quadratic as Quad
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Diagonal as Diagonal
import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Scalar as Scalar
import qualified Numeric.LAPACK.Shape.Private as ShapePriv
import qualified Numeric.LAPACK.Shape as ExtShape
import Numeric.LAPACK.Matrix.Array.Mosaic
(UnitUpperP, LowerP, UpperP, FlexUpperP,
HermitianP, HermitianPosSemidefP)
import Numeric.LAPACK.Matrix.Array.Private (Quadratic, ArrayMatrix, packTag)
import Numeric.LAPACK.Matrix.Square (Square)
import Numeric.LAPACK.Matrix
((#!), (#+#), (#-#), (#\##), (#*\), (#*##), (##*#))
import Numeric.LAPACK.Vector ((|+|), (|-|))
import qualified Numeric.Netlib.Class as Class
import qualified Type.Data.Bool as Bool
import Type.Data.Bool (False, True)
import Foreign.Storable (Storable)
import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape
import qualified Data.NonEmpty as NonEmpty
import qualified Data.Complex as Complex
import qualified Data.List.HT as ListHT
import qualified Data.Stream as Stream
import Data.Complex (Complex)
import Data.Semigroup ((<>))
import Data.Function.HT (nest)
import Data.Tuple.HT (mapFst, swap)
import Data.Stream (Stream)
import qualified Prelude as P
import Prelude hiding (sqrt, exp, log)
class (MatrixShape.Property property) => SqRt property where
sqrt ::
(Layout.Packing pack,
MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
Shape.C sh, Class.Real a) =>
Quadratic pack property lower upper sh a ->
Quadratic pack property lower upper sh a
instance SqRt Omni.Unit where
sqrt a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerIdentity -> a
Omni.PowerUpperTriangular -> sqrtUpper (const Scalar.one) a
Omni.PowerLowerTriangular ->
Triangular.transpose .
sqrtUpper (const Scalar.one) .
Triangular.transpose $ a
Omni.PowerDiagonal -> error "non-unit diagonal impossible"
Omni.PowerFull -> error "unit full matrix impossible"
instance SqRt Omni.Arbitrary where
sqrt m =
case Omni.powerSingleton $ ArrMatrix.shape m of
Omni.PowerDiagonal -> liftDiagonal P.sqrt m
Omni.PowerUpperTriangular -> sqrtUpper P.sqrt m
Omni.PowerLowerTriangular ->
Triangular.transpose . sqrtUpper P.sqrt . Triangular.transpose $ m
Omni.PowerFull ->
if Shape.size (Quad.size m) <= 2
then
ArrMatrix.Array $ Array.fromList (ArrMatrix.shape m) $
case Array.toList $ ArrMatrix.unwrap m of
[qA,qB,qC,qD] ->
let ((a,b),(c,d)) = sqrt2 ((qA,qB),(qC,qD))
in [a,b,c,d]
[qA] -> [P.sqrt qA]
_ -> []
else
case Scalar.precisionOfFunctor m of
Scalar.Float -> sqrtDenmanBeavers m
Scalar.Double -> sqrtDenmanBeavers m
instance SqRt Omni.Symmetric where
sqrt a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerSymmetric ->
Symmetric.fromHermitian . sqrtHermitian . Hermitian.fromSymmetric
$ a
_ -> error "Symmetric.sqrt: impossible shape"
instance
(neg ~ False, Bool.C zero, Bool.C pos) =>
SqRt (Omni.Hermitian neg zero pos) where
sqrt a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerDiagonal -> liftHermitianDiagonal P.sqrt a
Omni.PowerHermitian -> sqrtHermitian a
_ -> error "Hermitian.sqrt: impossible shape"
sqrtHermitian ::
(Layout.Packing pack,
Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
Shape.C sh, Class.Real a) =>
Quadratic pack herm Layout.Filled Layout.Filled sh a ->
Quadratic pack herm Layout.Filled Layout.Filled sh a
sqrtHermitian a =
case Scalar.precisionOfFunctor a of
Scalar.Float -> liftHermitian P.sqrt a
Scalar.Double -> liftHermitian P.sqrt a
sqrtUpper ::
(Layout.Packing pack, Omni.TriDiag diag, Shape.C sh, Class.Floating a) =>
(a -> a) -> FlexUpperP pack diag sh a -> FlexUpperP pack diag sh a
sqrtUpper sqrtF a =
Triangular.adaptOrder a $ ArrMatrix.lift1 Mosaic.reunpack $
LazyUpper.toStorable $ LazyUpper.sqrt sqrtF $ LazyUpper.fromStorable a
sqrt2 :: (Floating a) => ((a,a),(a,a)) -> ((a,a),(a,a))
sqrt2 ((a,b),(c,d)) =
let x = P.sqrt $ a+d + 2 * P.sqrt (a*d - b*c)
y = (d-a)/x
in (((x-y)/2, b/x), (c/x, (x+y)/2))
sqrtSchur :: (Shape.C sh, Class.Real a) => Square sh a -> Square sh a
sqrtSchur = applyUnchecked $ applyPermutable $ \a ->
let (q,u) = Square.schur $ checkZeroOffDiagonal a
in q <> sqrtUpper P.sqrt (Triangular.takeUpper u) #*## Square.adjoint q
checkZeroOffDiagonal ::
(Shape.Indexed sh, Class.Floating a) =>
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
checkZeroOffDiagonal a =
if all (\ij -> Scalar.isZero (a#!ij)) $
ListHT.mapAdjacent (flip (,)) $ Shape.indices $ Quad.size a
then a
else error "sqrtSchur: non-real eigenvalues"
sqrtDenmanBeavers ::
(ArrMatrix.Homogeneous prop, ArrMatrix.Additive prop) =>
(Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
(Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
sqrtDenmanBeavers = applyUnchecked $ \a ->
limit (Scalar.selectReal 1e-6 1e-14) $ fmap fst $
Stream.iterate
(\(b,c) ->
(Matrix.scaleReal 0.5 (b #+# Matrix.inverse c),
Matrix.scaleReal 0.5 (Matrix.inverse b #+# c)))
(a, Matrix.identityFrom a)
limit ::
(Extent.Measure meas, Extent.C vert, Extent.C horiz,
Shape.C height, Eq height, Shape.C width, Eq width,
Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar,
ArrMatrix.ArrayMatrix
pack prop lower upper meas vert horiz height width ~ m) =>
ar -> Stream (m a) -> m a
limit eps =
ArrMatrix.Array . snd . Stream.head .
Stream.dropWhile
(\(b0,b1) ->
Vector.normInf (b0|-|b1) > 0.5*eps * Vector.normInf (b0|+|b1)) .
streamMapAdjacent (,) . fmap ArrMatrix.unwrap
streamMapAdjacent :: (a -> a -> b) -> Stream a -> Stream b
streamMapAdjacent f xs = Stream.zipWith f xs (Stream.tail xs)
class (MatrixShape.Property property) => Exp property where
exp ::
(Layout.Packing pack,
MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
Shape.C sh, Class.Floating a) =>
Quadratic pack property lower upper sh a ->
Quadratic pack property lower upper sh a
instance Exp Omni.Arbitrary where
exp a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerDiagonal ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> liftDiagonal P.exp a
Scalar.Complex -> liftDiagonal P.exp a
Omni.PowerFull ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real ->
ArrMatrix.liftUnpacked1 id $
applyUnchecked (expScaledPade (#\##)) $ Matrix.toFull a
Scalar.Complex ->
ArrMatrix.liftUnpacked1 id $
applyUnchecked (expScaledPade (#\##)) $ Matrix.toFull a
Omni.PowerUpperTriangular ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> applyUnchecked (expScaledPade solveUpper) a
Scalar.Complex -> applyUnchecked (expScaledPade solveUpper) a
Omni.PowerLowerTriangular ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> applyUnchecked (expScaledPade solveLower) a
Scalar.Complex -> applyUnchecked (expScaledPade solveLower) a
solveUpper ::
(Layout.Packing pack, Shape.C sh, Eq sh, Class.Floating a) =>
UpperP pack sh a -> UpperP pack sh a -> UpperP pack sh a
solveUpper d n =
ArrMosaic.assureMirrored $ d #\## Triangular.toSquare n
solveLower ::
(Layout.Packing pack, Shape.C sh, Eq sh, Class.Floating a) =>
LowerP pack sh a -> LowerP pack sh a -> LowerP pack sh a
solveLower d n =
ArrMosaic.assureMirrored $ d #\## Triangular.toSquare n
instance
(neg ~ True, zero ~ True, pos ~ True) =>
Exp (Omni.Hermitian neg zero pos) where
exp a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerDiagonal ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> liftHermitianDiagonal P.exp a
Scalar.Complex -> liftHermitianDiagonal P.exp a
Omni.PowerHermitian ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> liftHermitian P.exp a
Scalar.Complex -> liftHermitian P.exp a
_ -> error "Hermitian.exp: impossible shape"
instance Exp Omni.Symmetric where
exp a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerSymmetric ->
case Scalar.complexSingletonOfFunctor a of
Scalar.Real ->
Symmetric.fromHermitian . liftHermitian P.exp .
Hermitian.fromSymmetric $ a
Scalar.Complex ->
applyUnchecked
(expScaledPade
(\d n ->
Symmetric.assureSymmetry $
d #\## Symmetric.toSquare n))
a
_ -> error "Symmetric.exp: impossible shape"
expScaledPade ::
(ArrMatrix.Homogeneous prop, ArrMatrix.Subtractive prop) =>
(Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
(Shape.C sh, Eq sh,
Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
(Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a) ->
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
expScaledPade solve a0 =
let s = max 0 $ (1+) $ ceiling $ logBase 2 $ norm1 a0
a = Matrix.scaleReal ((1/2)^s) a0
(odds,evens) = deinterleave $ NonEmpty.tail expPadeCoefficients
Stream.Cons eye as = Matrix.powers a
(oddPowers,evenPowers) = deinterleave $ Stream.toList as
v = foldl (#+#) eye $ zipWith Matrix.scaleReal evens evenPowers
u =
fmap (NonEmpty.foldl1 (#+#)) $ NonEmpty.fetch $
zipWith Matrix.scaleReal odds oddPowers
op vm f um = maybe vm (f vm) um
in nest s Matrix.square $ solve (op v (#-#) u) (op v (#+#) u)
deinterleave :: [a] -> ([a],[a])
deinterleave =
let go [] = ([],[])
go (x:xs) = mapFst (x:) $ swap $ go xs
in go
expPadeCoefficients :: (Class.Real a) => NonEmpty.T [] a
expPadeCoefficients =
let eps = Scalar.selectReal 1e-8 1e-16
q = expPadeOrder eps
coeff k = fromIntegral (q-k+1) / fromIntegral (k*(2*q-k+1))
in NonEmpty.scanl (*) (1 `asTypeOf` eps) $ map coeff [1..q]
expPadeOrder :: (Class.Real a) => a -> Int
expPadeOrder eps =
let factorials = scanl (*) 1 $ iterate (1+) 1
in subtract 1 $ length $ takeWhile id $
zipWith3 (\num den twoPower -> num > den*twoPower*eps)
factorials (ListHT.sieve 2 factorials) (iterate (2*) 1)
expRealHermitian ::
(Layout.Packing pack, Shape.C sh, Class.Real a) =>
HermitianP pack sh a -> HermitianPosSemidefP pack sh a
expRealHermitian =
Hermitian.relaxSemidefinite .
Hermitian.assurePositiveDefiniteness .
exp
class (MatrixShape.Property property) => Log property where
log ::
(Layout.Packing pack,
MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
Shape.C sh, Class.Real a) =>
Quadratic pack property lower upper sh a ->
Quadratic pack property lower upper sh a
instance Log Omni.Arbitrary where
log = liftReal P.log
instance
(neg ~ True, zero ~ True, pos ~ True) =>
Log (Omni.Hermitian neg zero pos) where
log = liftReal P.log
instance Log Omni.Symmetric where
log = liftReal P.log
logUnipotentUpper ::
(Layout.Packing pack, Shape.C sh, Class.Floating a) =>
UnitUpperP pack sh a -> UpperP pack sh a
logUnipotentUpper = logUnipotent . Triangular.relaxUnitDiagonal
logUnipotent ::
(Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
(ArrMatrix.Scale prop, ArrMatrix.Subtractive prop) =>
(Shape.C sh, Class.Floating a) =>
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
logUnipotent a =
case Scalar.complexSingletonOfFunctor a of
Scalar.Real -> logUnipotentAux a
Scalar.Complex -> logUnipotentAux a
logUnipotentAux ::
(Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
(ArrMatrix.Scale prop, ArrMatrix.Subtractive prop) =>
(Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
logUnipotentAux = applyUnchecked $ \a ->
let b = a #-# Matrix.identityFrom a
in foldl (#+#) (ArrMatrix.zero $ ArrMatrix.shape b) $
zipWith ArrMatrix.scale
(zipWith (/) (cycle [1,-1]) (iterate (1+) 1))
(Stream.takeWhile ((0<) . normInf) $ Matrix.powers1 b)
class (MatrixShape.Property property) => LiftReal property where
liftReal ::
(Layout.Packing pack,
MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
Shape.C sh, Class.Real a) =>
(a -> a) ->
Quadratic pack property lower upper sh a ->
Quadratic pack property lower upper sh a
instance LiftReal Omni.Arbitrary where
liftReal f a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerDiagonal -> liftDiagonal f a
Omni.PowerUpperTriangular -> flip applyUnchecked a $ \b ->
let (vr,d,vlAdj) = Triangular.eigensystem b
scal = Triangular.takeDiagonal $ vlAdj <> vr
in ArrMosaic.assureMirrored $
Triangular.toSquare vr
#*\ Vector.divide (Array.map f d) scal
##*# vlAdj
Omni.PowerLowerTriangular -> flip applyUnchecked a $ \b ->
let (vr,d,vlAdj) = Triangular.eigensystem b
scal = Triangular.takeDiagonal $ vlAdj <> vr
in ArrMosaic.assureMirrored $
Triangular.toSquare vr
#*\ Vector.divide (Array.map f d) scal
##*# vlAdj
Omni.PowerFull ->
case Scalar.precisionOfFunctor a of
Scalar.Float -> liftRealFull f a
Scalar.Double -> liftRealFull f a
liftRealFull ::
(MatrixShape.Property prop,
MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
Shape.C sh, Class.Real a, Scalar.RealOf a ~ a) =>
(a -> a) ->
Quadratic Layout.Unpacked prop lower upper sh a ->
Quadratic Layout.Unpacked prop lower upper sh a
liftRealFull f = applyPermutable $ applyUnchecked $ \a ->
let (vr,d,vlAdj) = Square.eigensystem $ Matrix.toFull a
vrR = matrixRealPart vr
vlAdjR = matrixRealPart vlAdj
dR = Array.map Complex.realPart d
scal = Square.takeDiagonal $ vlAdjR <> vrR
in if Scalar.isZero $ Vector.normInf $ Array.map Complex.imagPart d
then ArrMatrix.liftUnpacked1 id $
vrR #*\ Vector.divide (Array.map f dR) scal ##*# vlAdjR
else error "liftReal: non-real eigenvalues"
matrixRealPart ::
(ArrayMatrix pack property lower upper meas vert horiz height width ~ matrix,
Extent.Measure meas, Extent.C vert, Extent.C horiz,
Shape.C height, Shape.C width, Class.Real a) =>
matrix (Complex a) -> matrix a
matrixRealPart =
ArrMatrix.Array . Array.map Complex.realPart . ArrMatrix.unwrap
instance
(neg ~ True, zero ~ True, pos ~ True) =>
LiftReal (Omni.Hermitian neg zero pos) where
liftReal f a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerDiagonal -> liftHermitianDiagonal f a
Omni.PowerHermitian -> liftHermitianReal f a
_ -> error "Hermitian.liftReal: impossible shape"
instance LiftReal Omni.Symmetric where
liftReal f a =
case Omni.powerSingleton $ ArrMatrix.shape a of
Omni.PowerSymmetric ->
Symmetric.fromHermitian .
liftHermitianReal f . Hermitian.fromSymmetric $ a
_ -> error "Symmetric.liftReal: impossible shape"
liftHermitianReal ::
(Layout.Packing pack, Shape.C sh, Class.Real a) =>
(a -> a) -> HermitianP pack sh a -> HermitianP pack sh a
liftHermitianReal f a =
case Scalar.precisionOfFunctor a of
Scalar.Float -> liftHermitian f a
Scalar.Double -> liftHermitian f a
norm1 ::
(Extent.Measure meas, Extent.C vert, Extent.C horiz,
Shape.C height, Shape.C width, Class.Floating a) =>
ArrMatrix.ArrayMatrix pack prop lower upper meas vert horiz height width a ->
Scalar.RealOf a
norm1 = Vector.norm1 . ArrMatrix.unwrap
normInf ::
(Extent.Measure meas, Extent.C vert, Extent.C horiz,
Shape.C height, Shape.C width, Class.Floating a) =>
ArrMatrix.ArrayMatrix pack prop lower upper meas vert horiz height width a ->
Scalar.RealOf a
normInf = Vector.normInf . ArrMatrix.unwrap
liftDiagonal ::
(Layout.Packing pack, Shape.C sh, Class.Floating a, Class.Floating b) =>
(a -> b) ->
Quadratic pack Omni.Arbitrary Layout.Empty Layout.Empty sh a ->
Quadratic pack Omni.Arbitrary Layout.Empty Layout.Empty sh b
liftDiagonal f = Diagonal.lift $ Array.map f
liftHermitianDiagonal ::
(Layout.Packing pack,
Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
Shape.C sh, Class.Floating a, Class.Floating b) =>
(a -> b) ->
Quadratic pack herm Layout.Empty Layout.Empty sh a ->
Quadratic pack herm Layout.Empty Layout.Empty sh b
liftHermitianDiagonal f a =
case packTag a of
Layout.Packed -> ArrMatrix.lift1 (Array.map f) a
Layout.Unpacked ->
let b = Square.diagonal $ Array.map f $ Quad.takeDiagonal a
in ArrMatrix.liftUnpacked1 id $
if ArrMatrix.order a == ArrMatrix.order b
then b
else Square.transpose b
liftHermitian ::
(Layout.Packing pack,
Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Storable ar) =>
(ar -> ar) ->
Quadratic pack herm Layout.Filled Layout.Filled sh a ->
Quadratic pack herm Layout.Filled Layout.Filled sh a
liftHermitian f = applyPermutable $ applyUnchecked $ \a ->
let (q,d) = Hermitian.eigensystem a
in ArrMatrix.lift1 id $
Hermitian.congruenceDiagonalAdjoint (Square.toFull q) (Array.map f d)
applyPermutable ::
(Shape.C sh, Class.Floating a) =>
(Quadratic pack prop lower upper (ExtShape.IntIndexed sh) a ->
Quadratic pack prop lower upper (ExtShape.IntIndexed sh) a) ->
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
applyPermutable f =
Quad.mapSize ExtShape.deconsIntIndexed .
f .
Quad.mapSize ExtShape.IntIndexed
applyUnchecked ::
(Shape.C sh, Class.Floating a) =>
(Quadratic pack prop lower upper (ShapePriv.Unchecked sh) a ->
Quadratic pack prop lower upper (ShapePriv.Unchecked sh) a) ->
Quadratic pack prop lower upper sh a ->
Quadratic pack prop lower upper sh a
applyUnchecked f =
Quad.mapSize ShapePriv.deconsUnchecked .
f .
Quad.mapSize ShapePriv.Unchecked