{-# LANGUAGE BangPatterns #-}

{- | Symmetric eigenproblems in pure Haskell: cyclic Jacobi for full
decomposition (PCA covariance, @m×m@ kernels) and power iteration for the
dominant eigenpair (FISTA step sizes). Deterministic, sign-canonicalised output.
-}
module DataFrame.LinearAlgebra.Eigen (
    jacobiEigenSym,
    powerIterTop,
) where

import Control.Monad (forM_, when)
import Control.Monad.ST (runST)
import Data.List (sortBy)
import Data.Ord (Down (..), comparing)
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU
import qualified Data.Vector.Unboxed.Mutable as VUM
import DataFrame.LinearAlgebra (Matrix, dot, matVec, scaleV)

{- | Cyclic Jacobi eigendecomposition of a symmetric matrix. Eigenvalues are
returned in descending order paired with eigenvectors as rows, each
sign-canonicalised (largest-magnitude component positive) for unique output.
-}
jacobiEigenSym :: Matrix -> (VU.Vector Double, Matrix)
jacobiEigenSym :: Matrix -> (Vector Double, Matrix)
jacobiEigenSym Matrix
a0
    | Matrix -> Bool
forall a. Vector a -> Bool
V.null Matrix
a0 = (Vector Double
forall a. Unbox a => Vector a
VU.empty, Matrix
forall a. Vector a
V.empty)
    | Bool
otherwise = (forall s. ST s (Vector Double, Matrix)) -> (Vector Double, Matrix)
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Double, Matrix))
 -> (Vector Double, Matrix))
-> (forall s. ST s (Vector Double, Matrix))
-> (Vector Double, Matrix)
forall a b. (a -> b) -> a -> b
$ do
        MVector s Double
a <- Int -> ST s (MVector (PrimState (ST s)) Double)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
VUM.new (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d)
        [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i ->
            [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j ->
                MVector (PrimState (ST s)) Double -> Int -> Double -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector s Double
MVector (PrimState (ST s)) Double
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) ((Matrix
a0 Matrix -> Int -> Vector Double
forall a. Vector a -> Int -> a
V.! Int
i) Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
VU.! Int
j)
        MVector s Double
v <- Int -> Double -> ST s (MVector (PrimState (ST s)) Double)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
VUM.replicate (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d) Double
0
        [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> MVector (PrimState (ST s)) Double -> Int -> Double -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector s Double
MVector (PrimState (ST s)) Double
v (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) Double
1
        MVector (PrimState (ST s)) Double
-> MVector (PrimState (ST s)) Double -> Integer -> ST s ()
forall {f :: * -> *}.
PrimMonad f =>
MVector (PrimState f) Double
-> MVector (PrimState f) Double -> Integer -> f ()
sweep MVector s Double
MVector (PrimState (ST s)) Double
a MVector s Double
MVector (PrimState (ST s)) Double
v Integer
0
        Vector Double
afrozen <- MVector (PrimState (ST s)) Double -> ST s (Vector Double)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
VU.freeze MVector s Double
MVector (PrimState (ST s)) Double
a
        Vector Double
vmat <- MVector (PrimState (ST s)) Double -> ST s (Vector Double)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
MVector (PrimState m) a -> m (Vector a)
VU.freeze MVector s Double
MVector (PrimState (ST s)) Double
v
        let diag :: Vector Double
diag = Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
VU.generate Int
d (\Int
i -> Vector Double
afrozen Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
VU.! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i))
            vecs :: Matrix
vecs =
                Int -> (Int -> Vector Double) -> Matrix
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
d ((Int -> Vector Double) -> Matrix)
-> (Int -> Vector Double) -> Matrix
forall a b. (a -> b) -> a -> b
$ \Int
col ->
                    Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
VU.generate Int
d ((Int -> Double) -> Vector Double)
-> (Int -> Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ \Int
row -> Vector Double
vmat Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
VU.! (Int
row Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
col)
            paired :: [(Double, Vector Double)]
paired =
                ((Double, Vector Double) -> (Double, Vector Double) -> Ordering)
-> [(Double, Vector Double)] -> [(Double, Vector Double)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy
                    (((Double, Vector Double) -> Down Double)
-> (Double, Vector Double) -> (Double, Vector Double) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Double -> Down Double
forall a. a -> Down a
Down (Double -> Down Double)
-> ((Double, Vector Double) -> Double)
-> (Double, Vector Double)
-> Down Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Vector Double) -> Double
forall a b. (a, b) -> a
fst))
                    ([Double] -> [Vector Double] -> [(Double, Vector Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Vector Double -> [Double]
forall a. Unbox a => Vector a -> [a]
VU.toList Vector Double
diag) (Matrix -> [Vector Double]
forall a. Vector a -> [a]
V.toList Matrix
vecs))
        (Vector Double, Matrix) -> ST s (Vector Double, Matrix)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure
            ( [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
VU.fromList (((Double, Vector Double) -> Double)
-> [(Double, Vector Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Vector Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Vector Double)]
paired)
            , [Vector Double] -> Matrix
forall a. [a] -> Vector a
V.fromList (((Double, Vector Double) -> Vector Double)
-> [(Double, Vector Double)] -> [Vector Double]
forall a b. (a -> b) -> [a] -> [b]
map (Vector Double -> Vector Double
canonicalSign (Vector Double -> Vector Double)
-> ((Double, Vector Double) -> Vector Double)
-> (Double, Vector Double)
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Vector Double) -> Vector Double
forall a b. (a, b) -> b
snd) [(Double, Vector Double)]
paired)
            )
  where
    d :: Int
d = Matrix -> Int
forall a. Vector a -> Int
V.length Matrix
a0
    maxSweeps :: Integer
maxSweeps = Integer
100
    tol :: Double
tol = Double
1e-12
    sweep :: MVector (PrimState f) Double
-> MVector (PrimState f) Double -> Integer -> f ()
sweep MVector (PrimState f) Double
a MVector (PrimState f) Double
v Integer
s
        | Integer
s Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
maxSweeps = () -> f ()
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
        | Bool
otherwise = do
            Double
off <- MVector (PrimState f) Double -> f Double
forall {f :: * -> *} {t}.
(PrimMonad f, Unbox t, Num t) =>
MVector (PrimState f) t -> f t
offNorm MVector (PrimState f) Double
a
            Bool -> f () -> f ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
off Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
tol) (f () -> f ()) -> f () -> f ()
forall a b. (a -> b) -> a -> b
$ do
                [Int] -> (Int -> f ()) -> f ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2] ((Int -> f ()) -> f ()) -> (Int -> f ()) -> f ()
forall a b. (a -> b) -> a -> b
$ \Int
p ->
                    [Int] -> (Int -> f ()) -> f ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> f ()) -> f ()) -> (Int -> f ()) -> f ()
forall a b. (a -> b) -> a -> b
$ \Int
q -> MVector (PrimState f) Double
-> MVector (PrimState f) Double -> Int -> Int -> f ()
forall {m :: * -> *} {a}.
(PrimMonad m, Unbox a, Ord a, Floating a) =>
MVector (PrimState m) a
-> MVector (PrimState m) a -> Int -> Int -> m ()
rotate MVector (PrimState f) Double
a MVector (PrimState f) Double
v Int
p Int
q
                MVector (PrimState f) Double
-> MVector (PrimState f) Double -> Integer -> f ()
sweep MVector (PrimState f) Double
a MVector (PrimState f) Double
v (Integer
s Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)
    offNorm :: MVector (PrimState f) t -> f t
offNorm MVector (PrimState f) t
a = Int -> t -> f t
forall {f :: * -> *}.
(PrimState f ~ PrimState f, PrimMonad f) =>
Int -> t -> f t
go Int
0 t
0
      where
        go :: Int -> t -> f t
go Int
i !t
acc
            | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
d = t -> f t
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure t
acc
            | Bool
otherwise = do
                t
r <- Int -> Int -> t -> f t
forall {f :: * -> *}.
(PrimState f ~ PrimState f, PrimMonad f) =>
Int -> Int -> t -> f t
goRow Int
i (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
acc
                Int -> t -> f t
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
r
        goRow :: Int -> Int -> t -> f t
goRow Int
i Int
j !t
acc
            | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
d = t -> f t
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure t
acc
            | Bool
otherwise = do
                t
x <- MVector (PrimState f) t -> Int -> f t
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState f) t
MVector (PrimState f) t
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j)
                Int -> Int -> t -> f t
goRow Int
i (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (t
acc t -> t -> t
forall a. Num a => a -> a -> a
+ t
x t -> t -> t
forall a. Num a => a -> a -> a
* t
x)
    rotate :: MVector (PrimState m) a
-> MVector (PrimState m) a -> Int -> Int -> m ()
rotate MVector (PrimState m) a
a MVector (PrimState m) a
v Int
p Int
q = do
        a
apq <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q)
        Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (a -> a
forall a. Num a => a -> a
abs a
apq a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
1e-300) (m () -> m ()) -> m () -> m ()
forall a b. (a -> b) -> a -> b
$ do
            a
app <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p)
            a
aqq <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
q Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q)
            let theta :: a
theta = (a
aqq a -> a -> a
forall a. Num a => a -> a -> a
- a
app) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
apq)
                s' :: a
s' = if a
theta a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
1 else a -> a
forall a. Num a => a -> a
signum a
theta
                t :: a
t = a
s' a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a -> a
forall a. Num a => a -> a
abs a
theta a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
sqrt (a
theta a -> a -> a
forall a. Num a => a -> a -> a
* a
theta a -> a -> a
forall a. Num a => a -> a -> a
+ a
1))
                c :: a
c = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
t a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)
                sn :: a
sn = a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
c
            [Int] -> (Int -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> m ()) -> m ()) -> (Int -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
                a
aip <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p)
                a
aiq <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p) (a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
aip a -> a -> a
forall a. Num a => a -> a -> a
- a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
aiq)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
a (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q) (a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
aip a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
aiq)
            [Int] -> (Int -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> m ()) -> m ()) -> (Int -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
                a
apj <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j)
                a
aqj <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
a (Int
q Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
a (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) (a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
apj a -> a -> a
forall a. Num a => a -> a -> a
- a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
aqj)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
a (Int
q Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) (a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
apj a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
aqj)
            [Int] -> (Int -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] ((Int -> m ()) -> m ()) -> (Int -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
                a
vip <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
v (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p)
                a
viq <- MVector (PrimState m) a -> Int -> m a
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
VUM.read MVector (PrimState m) a
v (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
v (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p) (a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
vip a -> a -> a
forall a. Num a => a -> a -> a
- a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
viq)
                MVector (PrimState m) a -> Int -> a -> m ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
VUM.write MVector (PrimState m) a
v (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q) (a
sn a -> a -> a
forall a. Num a => a -> a -> a
* a
vip a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
viq)

canonicalSign :: VU.Vector Double -> VU.Vector Double
canonicalSign :: Vector Double -> Vector Double
canonicalSign Vector Double
vec =
    let idx :: Int
idx = Vector Double -> Int
forall a. (Unbox a, Ord a) => Vector a -> Int
VU.maxIndex ((Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map Double -> Double
forall a. Num a => a -> a
abs Vector Double
vec)
     in if Vector Double
vec Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
VU.! Int
idx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then (Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map Double -> Double
forall a. Num a => a -> a
negate Vector Double
vec else Vector Double
vec

{- | Dominant eigenvalue and eigenvector of a symmetric PSD matrix via power
iteration with a deterministic all-ones start.
-}
powerIterTop :: Int -> Matrix -> (Double, VU.Vector Double)
powerIterTop :: Int -> Matrix -> (Double, Vector Double)
powerIterTop Int
iters Matrix
a
    | Matrix -> Bool
forall a. Vector a -> Bool
V.null Matrix
a = (Double
0, Vector Double
forall a. Unbox a => Vector a
VU.empty)
    | Bool
otherwise = Int -> Vector Double -> (Double, Vector Double)
forall {t}.
(Eq t, Num t) =>
t -> Vector Double -> (Double, Vector Double)
go Int
iters (Vector Double -> Vector Double
normalize (Int -> Double -> Vector Double
forall a. Unbox a => Int -> a -> Vector a
VU.replicate Int
d Double
1))
  where
    d :: Int
d = Matrix -> Int
forall a. Vector a -> Int
V.length Matrix
a
    normalize :: Vector Double -> Vector Double
normalize Vector Double
v =
        let nrm :: Double
nrm = Double -> Double
forall a. Floating a => a -> a
sqrt (Vector Double -> Vector Double -> Double
dot Vector Double
v Vector Double
v) in if Double
nrm Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Vector Double
v else Double -> Vector Double -> Vector Double
scaleV (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nrm) Vector Double
v
    go :: t -> Vector Double -> (Double, Vector Double)
go t
0 Vector Double
v = (Vector Double -> Vector Double -> Double
dot Vector Double
v (Matrix -> Vector Double -> Vector Double
matVec Matrix
a Vector Double
v), Vector Double
v)
    go t
k Vector Double
v =
        let av :: Vector Double
av = Matrix -> Vector Double -> Vector Double
matVec Matrix
a Vector Double
v
            nrm :: Double
nrm = Double -> Double
forall a. Floating a => a -> a
sqrt (Vector Double -> Vector Double -> Double
dot Vector Double
av Vector Double
av)
         in if Double
nrm Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-300 then (Double
0, Vector Double
v) else t -> Vector Double -> (Double, Vector Double)
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) (Double -> Vector Double -> Vector Double
scaleV (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nrm) Vector Double
av)