-- SPDX-License-Identifier: MPL-2.0

{- | Matrix constructors, products, spaces, and matrix functions. Spectral
functions are Double precision. The simple log/sqrt entries use symmetric
spectral algorithms; General variants accept nonsymmetric invertible real
matrices off the nonpositive-real spectral branch cut, with finite budgets.
Exponential uses scaling and squaring.
-}
module HMath.Matrix.Advanced (
    generateMatrix,
    triangularMatrix,
    symmetricMatrix,
    hermitianMatrix,
    permutationMatrix,
    blockMatrix,
    bandedMatrix,
    conjugateTranspose,
    zipMatrixWith,
    batchedMultiply,
    kronecker,
    khatriRao,
    matrixPower,
    trace,
    matrixNorm1,
    matrixNormInfinity,
    frobeniusNorm,
    spectralNorm,
    pseudoinverse,
    complexPseudoinverse,
    matrixRank,
    conditionNumber,
    nullSpace,
    columnSpace,
    rowSpace,
    changeCoordinates,
    matrixExp,
    matrixLog,
    matrixSqrt,
    matrixLogGeneral,
    matrixSqrtGeneral,
    polar,
) where

import Control.Monad (foldM)
import Data.Complex (Complex ((:+)), conjugate, imagPart, realPart)
import Data.List (sort)
import HMath.Error
import HMath.Matrix
import HMath.Matrix.Internal
import HMath.Matrix.LAPACK
import HMath.Matrix.LU (inverse, solve)
import Prelude hiding (log, sqrt)
import qualified Prelude as P

{- | Given row count, column count, then zero-based (row,column) generator, construct a
row-major matrix. Negative dimensions or Int-size element product overflow return Left
MathError.
-}
generateMatrix :: Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
generateMatrix :: forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
generateMatrix = Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build

{- | Given lower flag then matrix, preserve entries on/below the diagonal if True, on/above
if False, and zero the rest. Rectangular matrices are allowed; return checked
construction result.
-}
triangularMatrix :: (Num a) => Bool -> Matrix a -> Either MathError (Matrix a)
triangularMatrix :: forall a. Num a => Bool -> Matrix a -> Either MathError (Matrix a)
triangularMatrix Bool
lower Matrix a
a = Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
r Int
c (\Int
i Int
j -> if (Bool
lower Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
i) Bool -> Bool -> Bool
|| (Bool -> Bool
not Bool
lower Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
i) then Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j else a
0)
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Return (A+transpose(A))/2 for a square matrix. Non-square input returns
NonSquareMatrix; arithmetic follows Fractional.
-}
symmetricMatrix :: (Fractional a) => Matrix a -> Either MathError (Matrix a)
symmetricMatrix :: forall a. Fractional a => Matrix a -> Either MathError (Matrix a)
symmetricMatrix Matrix a
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | Bool
otherwise = Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
r Int
c (\Int
i Int
j -> (Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
+ Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
j Int
i) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2)
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Return (A+conjugateTranspose(A))/2 for a square complex matrix. Non-square input fails;
this symmetrizes rather than merely asserting the original is Hermitian.
-}
hermitianMatrix :: (RealFloat a) => Matrix (Complex a) -> Either MathError (Matrix (Complex a))
hermitianMatrix :: forall a.
RealFloat a =>
Matrix (Complex a) -> Either MathError (Matrix (Complex a))
hermitianMatrix Matrix (Complex a)
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError (Matrix (Complex a))
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | Bool
otherwise = Int
-> Int
-> (Int -> Int -> Complex a)
-> Either MathError (Matrix (Complex a))
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
r Int
c (\Int
i Int
j -> (Matrix (Complex a) -> Int -> Int -> Complex a
forall a. Matrix a -> Int -> Int -> a
get Matrix (Complex a)
a Int
i Int
j Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
+ Complex a -> Complex a
forall a. Num a => Complex a -> Complex a
conjugate (Matrix (Complex a) -> Int -> Int -> Complex a
forall a. Matrix a -> Int -> Int -> a
get Matrix (Complex a)
a Int
j Int
i)) Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/ Complex a
2)
  where
    (Int
r, Int
c) = Matrix (Complex a) -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix (Complex a)
a

{- | Given a permutation p of [0..n-1], construct n-by-n P with P[i,p[i]]=1. Duplicate,
missing, or out-of-range indices return Left InvalidArgument.
-}
permutationMatrix :: (Num a) => [Int] -> Either MathError (Matrix a)
permutationMatrix :: forall a. Num a => [Int] -> Either MathError (Matrix a)
permutationMatrix [Int]
p
    | [Int] -> [Int]
forall a. Ord a => [a] -> [a]
sort [Int]
p [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
/= [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"permutationMatrix: indices must be a permutation")
    | Bool
otherwise = Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> if [Int]
p [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
1 else a
0)
  where
    n :: Int
n = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
p

{- | Assemble a rectangular grid of matrix blocks. Require equal heights within each block
row and equal widths within each block column. Ragged grids or total-size overflow
fail; empty grids yield a 0-by-0 matrix.
-}
blockMatrix :: [[Matrix a]] -> Either MathError (Matrix a)
blockMatrix :: forall a. [[Matrix a]] -> Either MathError (Matrix a)
blockMatrix [] = Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
0 Int
0 []
blockMatrix [[Matrix a]]
blocks
    | ([Matrix a] -> Bool) -> [[Matrix a]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Matrix a] -> Int) -> [Matrix a] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Matrix a]]
blocks = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"blockMatrix: ragged block grid")
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
0 Int
0 []
    | ([Matrix a] -> Bool) -> [[Matrix a]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\[Matrix a]
row -> (Matrix a -> Bool) -> [Matrix a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape ([Matrix a]
row [Matrix a] -> Int -> Matrix a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0))) (Int -> Bool) -> (Matrix a -> Int) -> Matrix a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Int) -> Int
forall a b. (a, b) -> a
fst ((Int, Int) -> Int) -> (Matrix a -> (Int, Int)) -> Matrix a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape) [Matrix a]
row) [[Matrix a]]
blocks =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"blockMatrix: inconsistent block-row heights")
    | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Int
j -> ([Matrix a] -> Bool) -> [[Matrix a]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape (([[Matrix a]]
blocks [[Matrix a]] -> Int -> [Matrix a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Matrix a] -> Int -> Matrix a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))) (Int -> Bool) -> ([Matrix a] -> Int) -> [Matrix a] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Int) -> Int
forall a b. (a, b) -> b
snd ((Int, Int) -> Int)
-> ([Matrix a] -> (Int, Int)) -> [Matrix a] -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape (Matrix a -> (Int, Int))
-> ([Matrix a] -> Matrix a) -> [Matrix a] -> (Int, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Matrix a] -> Int -> Matrix a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) [[Matrix a]]
blocks) [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"blockMatrix: inconsistent block-column widths")
    | Integer
totalR Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) Bool -> Bool -> Bool
|| Integer
totalC Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"blockMatrix: shape overflow")
    | Bool
otherwise =
        Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked
            (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
totalR)
            (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
totalC)
            ( [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
                [ [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
block [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Matrix a
block <- [Matrix a]
row] | Int
i <- [Int
0 .. (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape ([Matrix a]
row [Matrix a] -> Int -> Matrix a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                | [Matrix a]
row <- [[Matrix a]]
blocks
                ]
            )
  where
    n :: Int
n = [Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Matrix a]]
blocks [[Matrix a]] -> Int -> [Matrix a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
    totalR :: Integer
totalR = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> Integer
forall a. Integral a => a -> Integer
toInteger ((Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape ([Matrix a]
row [Matrix a] -> Int -> Matrix a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0))) | [Matrix a]
row <- [[Matrix a]]
blocks]
    totalC :: Integer
totalC = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> Integer
forall a. Integral a => a -> Integer
toInteger ((Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
b)) | Matrix a
b <- [[Matrix a]]
blocks [[Matrix a]] -> Int -> [Matrix a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0]

{- | Given rows, columns, then (diagonal offset,values) pairs, construct a matrix with zeros
elsewhere. Offset k means column-row=k. Each list must match that diagonal's length;
duplicates add. Invalid dimensions/counts fail.
-}
bandedMatrix :: (Num a) => Int -> Int -> [(Int, [a])] -> Either MathError (Matrix a)
bandedMatrix :: forall a.
Num a =>
Int -> Int -> [(Int, [a])] -> Either MathError (Matrix a)
bandedMatrix Int
r Int
c [(Int, [a])]
diagonals
    | Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NegativeDimensions (Int
r, Int
c))
    | Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
r Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
c Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
ElementCountOverflow (Int
r, Int
c))
    | ((Int, [a]) -> Bool) -> [(Int, [a])] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\(Int
k, [a]
xs) -> Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Int -> Integer
forall a. Integral a => a -> Integer
diagonalLength Int
k) [(Int, [a])]
diagonals =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"bandedMatrix: diagonal length mismatch")
    | Bool
otherwise =
        Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build
            Int
r
            Int
c
            ( \Int
i Int
j -> [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
0 (-Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
k)) | (Int
k, [a]
xs) <- [(Int, [a])]
diagonals, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k]
            )
  where
    diagonalLength :: a -> Integer
diagonalLength a
k = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
0 (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
r) (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- a -> Integer
forall a. Integral a => a -> Integer
toInteger a
k) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
0 (-a -> Integer
forall a. Integral a => a -> Integer
toInteger a
k))

{- | Transpose a complex matrix and conjugate every entry, exchanging its row/column
dimensions.
-}
conjugateTranspose :: (RealFloat a) => Matrix (Complex a) -> Matrix (Complex a)
conjugateTranspose :: forall a. RealFloat a => Matrix (Complex a) -> Matrix (Complex a)
conjugateTranspose = (Complex a -> Complex a)
-> Matrix (Complex a) -> Matrix (Complex a)
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries Complex a -> Complex a
forall a. Num a => Complex a -> Complex a
conjugate (Matrix (Complex a) -> Matrix (Complex a))
-> (Matrix (Complex a) -> Matrix (Complex a))
-> Matrix (Complex a)
-> Matrix (Complex a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Complex a) -> Matrix (Complex a)
forall a. Matrix a -> Matrix a
transposeMatrix

{- | Given elementwise function and two identically shaped matrices, apply it to
corresponding entries and return a new matrix. Mismatched shapes return ShapeMismatch.
-}
zipMatrixWith :: (a -> b -> c) -> Matrix a -> Matrix b -> Either MathError (Matrix c)
zipMatrixWith :: forall a b c.
(a -> b -> c)
-> Matrix a -> Matrix b -> Either MathError (Matrix c)
zipMatrixWith a -> b -> c
f Matrix a
a Matrix b
b
    | Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
/= Matrix b -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix b
b = MathError -> Either MathError (Matrix c)
forall a b. a -> Either a b
Left ((Int, Int) -> (Int, Int) -> MathError
ShapeMismatch (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a) (Matrix b -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix b
b))
    | Bool
otherwise = Int -> Int -> [c] -> Either MathError (Matrix c)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
r Int
c ((a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> b -> c
f (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a) (Matrix b -> [b]
forall a. Matrix a -> [a]
toList Matrix b
b))
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Multiply pairs of matrices from equally long lists, preserving batch order. Different
batch lengths or any checked multiplication failure return Left MathError.
-}
batchedMultiply :: (Num a) => [Matrix a] -> [Matrix a] -> Either MathError [Matrix a]
batchedMultiply :: forall a.
Num a =>
[Matrix a] -> [Matrix a] -> Either MathError [Matrix a]
batchedMultiply [Matrix a]
as [Matrix a]
bs
    | [Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Matrix a]
as Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Matrix a]
bs = MathError -> Either MathError [Matrix a]
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch ([Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Matrix a]
as) ([Matrix a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Matrix a]
bs))
    | Bool
otherwise = [Either MathError (Matrix a)] -> Either MathError [Matrix a]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence ((Matrix a -> Matrix a -> Either MathError (Matrix a))
-> [Matrix a] -> [Matrix a] -> [Either MathError (Matrix a)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked [Matrix a]
as [Matrix a]
bs)

{- | Return the block Kronecker product of matrices A and B, of shape
(Arows*Brows,Acols*Bcols). Reject dimension/element-count overflow.
-}
kronecker :: (Num a) => Matrix a -> Matrix a -> Either MathError (Matrix a)
kronecker :: forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
kronecker Matrix a
a Matrix a
b
    | Integer
rows Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) Bool -> Bool -> Bool
|| Integer
cols Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"kronecker: dimension overflow")
    | Bool
otherwise =
        Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build
            (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
rows)
            (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
cols)
            (\Int
i Int
j -> Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a (Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
br) (Int
j Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
bc) a -> a -> a
forall a. Num a => a -> a -> a
* Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
b (Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
br) (Int
j Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
bc))
  where
    (Int
ar, Int
ac) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    (Int
br, Int
bc) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
b
    rows :: Integer
rows = Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
ar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
br
    cols :: Integer
cols = Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
ac Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
bc

{- | Return columnwise Kronecker product, requiring equal column counts. Shape is
(Arows*Brows,columns); mismatched columns or dimension overflow fail.
-}
khatriRao :: (Num a) => Matrix a -> Matrix a -> Either MathError (Matrix a)
khatriRao :: forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
khatriRao Matrix a
a Matrix a
b
    | Int
ac Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
bc = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
InnerDimensionMismatch Int
ac Int
bc)
    | Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
ar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
br Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) =
        MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"khatriRao: dimension overflow")
    | Bool
otherwise = Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build (Int
ar Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
br) Int
ac (\Int
i Int
j -> Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a (Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
br) Int
j a -> a -> a
forall a. Num a => a -> a -> a
* Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
b (Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
br) Int
j)
  where
    (Int
ar, Int
ac) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a; (Int
br, Int
bc) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
b

{- | Given Integer exponent then square real matrix, compute power by squaring. Zero yields
identity; negative exponent inverts first. Non-square/singular input or multiplication
errors propagate.
-}
matrixPower :: (RealFloat a) => Integer -> Matrix a -> Either MathError (Matrix a)
matrixPower :: forall a.
RealFloat a =>
Integer -> Matrix a -> Either MathError (Matrix a)
matrixPower Integer
power Matrix a
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | Integer
power Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = Matrix a -> Either MathError (Matrix a)
forall a. RealFloat a => Matrix a -> Either MathError (Matrix a)
inverse Matrix a
a Either MathError (Matrix a)
-> (Matrix a -> Either MathError (Matrix a))
-> Either MathError (Matrix a)
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Integer -> Matrix a -> Either MathError (Matrix a)
forall a.
RealFloat a =>
Integer -> Matrix a -> Either MathError (Matrix a)
matrixPower (-Integer
power)
    | Bool
otherwise = do i <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
r; go power a i
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    go :: t -> Matrix a -> Matrix a -> Either MathError (Matrix a)
go t
0 Matrix a
_ Matrix a
result = Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right Matrix a
result
    go t
k Matrix a
base Matrix a
result = do
        next <- if t -> Bool
forall a. Integral a => a -> Bool
odd t
k then Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked Matrix a
result Matrix a
base else Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right Matrix a
result
        if k == 1 then Right next else do squared <- multiplyChecked base base; go (k `div` 2) squared next

{- | Return sum of main diagonal of a square matrix; non-square input returns
NonSquareMatrix. A 0-by-0 matrix has trace zero.
-}
trace :: (Num a) => Matrix a -> Either MathError a
trace :: forall a. Num a => Matrix a -> Either MathError a
trace Matrix a
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError a
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | Bool
otherwise = a -> Either MathError a
forall a b. b -> Either a b
Right ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
i | Int
i <- [Int
0 .. Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Return maximum absolute column sum (induced one norm). Empty dimensions give zero; no
finiteness check is added.
-}
matrixNorm1 :: (Ord a, Num a) => Matrix a -> a
matrixNorm1 :: forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix a
a = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [[a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a -> a
forall a. Num a => a -> a
abs (Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j) | Int
i <- [Int
0 .. Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
j <- [Int
0 .. Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Return maximum absolute row sum (induced infinity norm). A matrix with no elements
gives zero.
-}
matrixNormInfinity :: (Ord a, Num a) => Matrix a -> a
matrixNormInfinity :: forall a. (Ord a, Num a) => Matrix a -> a
matrixNormInfinity Matrix a
a
    | [a] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a) = a
0
    | Bool
otherwise = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> ([a] -> [a]) -> [a] -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs) (Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a))

{- | Return scaled Euclidean norm of all entries, reducing intermediate overflow; empty
matrix norm is zero. No Either result is added.
-}
frobeniusNorm :: (RealFloat a) => Matrix a -> a
frobeniusNorm :: forall a. RealFloat a => Matrix a -> a
frobeniusNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm ([a] -> a) -> (Matrix a -> [a]) -> Matrix a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> [a]
forall a. Matrix a -> [a]
toList

{- | Return largest singular value via SVD, or zero for an empty singular spectrum. Invalid
input/backend decomposition failures return Left MathError.
-}
spectralNorm :: Matrix Double -> Either MathError Double
spectralNorm :: Matrix Double -> Either MathError Double
spectralNorm Matrix Double
a = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Double] -> Double) -> (SVD -> [Double]) -> SVD -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
:) ([Double] -> [Double]) -> (SVD -> [Double]) -> SVD -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SVD -> [Double]
singularValues (SVD -> Double) -> Either MathError SVD -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Matrix Double -> Either MathError SVD
svd Matrix Double
a

validTol :: Double -> Either MathError ()
validTol :: Double -> Either MathError ()
validTol Double
t
    | Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
t Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
t = MathError -> Either MathError ()
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()

{- | Given finite nonnegative relative tolerance then matrix, count positive singular values
whose ratio to the largest exceeds tolerance. Zero tolerance counts all positive
singular values. Invalid tolerance or SVD failure propagates.
-}
matrixRank :: Double -> Matrix Double -> Either MathError Int
matrixRank :: Double -> Matrix Double -> Either MathError Int
matrixRank Double
tol Matrix Double
a = do
    Double -> Either MathError ()
validTol Double
tol
    f <- Matrix Double -> Either MathError SVD
svd Matrix Double
a
    let s = SVD -> [Double]
singularValues SVD
f; biggest = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
s)
    Right (length [x | x <- s, x > 0, tol == 0 || x / biggest > tol])

{- | Given relative singular-value tolerance then m-by-n matrix, return n-by-m Moore-Penrose
pseudoinverse via thresholded SVD. Rectangular/rank-deficient inputs are supported;
invalid tolerance or numeric/decomposition failures propagate.
-}
pseudoinverse :: Double -> Matrix Double -> Either MathError (Matrix Double)
pseudoinverse :: Double -> Matrix Double -> Either MathError (Matrix Double)
pseudoinverse Double
tol Matrix Double
a = do
    Double -> Either MathError ()
validTol Double
tol
    f <- Matrix Double -> Either MathError SVD
svd Matrix Double
a
    let s = SVD -> [Double]
singularValues SVD
f; biggest = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
s); k = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s
    d <- diagonalChecked [if x > 0 && (tol == 0 || x / biggest > tol) then 1 / x else 0 | x <- s]
    utd <- multiplyChecked d (transposeMatrix (svdU f))
    result <- multiplyChecked (transposeMatrix (svdVT f)) utd
    if k == 0
        then zeroMatrixChecked (snd (shape a)) (fst (shape a))
        else finiteMatrix result >> Right result

{- | Return largest/smallest singular-value ratio; singular matrices yield positive
infinity, empty spectrum yields 1. Input/SVD failures return Left MathError.
-}
conditionNumber :: Matrix Double -> Either MathError Double
conditionNumber :: Matrix Double -> Either MathError Double
conditionNumber Matrix Double
a = do
    f <- Matrix Double -> Either MathError SVD
svd Matrix Double
a
    case singularValues f of
        [] -> Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
1
        [Double]
xs -> Double -> Either MathError Double
forall a b. b -> Either a b
Right (if [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 else [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs)

{- | Given relative SVD rank tolerance then matrix, return columns of left singular vectors
spanning its numerical column space. Invalid tolerance/decomposition errors propagate.
-}
columnSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
columnSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
columnSpace Double
tol Matrix Double
a = do rank <- Double -> Matrix Double -> Either MathError Int
matrixRank Double
tol Matrix Double
a; f <- svd a; sliceMatrixChecked (0, 0) (fst (shape a), rank) (svdU f)

{- | Given relative rank tolerance then matrix, return a column basis for its row space by
applying columnSpace to its transpose.
-}
rowSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
rowSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
rowSpace Double
tol Matrix Double
a = Double -> Matrix Double -> Either MathError (Matrix Double)
columnSpace Double
tol (Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a
transposeMatrix Matrix Double
a)

{- | Given relative singular-value tolerance then m-by-n matrix, return an n-by-(n-rank)
basis for its numerical right null space using full SVD. Invalid
tolerance/decomposition failures propagate.
-}
nullSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
nullSpace :: Double -> Matrix Double -> Either MathError (Matrix Double)
nullSpace Double
tol Matrix Double
a = do
    Double -> Either MathError ()
validTol Double
tol
    f <- Matrix Double -> Either MathError SVD
fullSVD Matrix Double
a
    rank <- matrixRank tol a
    let n = (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a)
    sliceMatrixChecked (0, rank) (n, n - rank) (transposeMatrix (svdVT f))

{- | Given source basis matrix, target square basis matrix, then source coordinate list, map
to physical coordinates and solve in target basis. Shape mismatch, singular target, and
checked solve errors propagate.
-}
changeCoordinates :: (RealFloat a) => Matrix a -> Matrix a -> [a] -> Either MathError [a]
changeCoordinates :: forall a.
RealFloat a =>
Matrix a -> Matrix a -> [a] -> Either MathError [a]
changeCoordinates Matrix a
source Matrix a
target [a]
coordinates = do
    rhs <- Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
coordinates) Int
1 [a]
coordinates
    physical <- multiplyChecked source rhs
    solve target (toList physical)

{- | Given positive finite tolerance, positive iteration budget, then square finite matrix,
approximate exp(A) by scaling, Taylor series and squaring. Return NonConvergence on
exhausted series budget; input/shape/numerical failures propagate.
-}
matrixExp :: (RealFloat a) => a -> Int -> Matrix a -> Either MathError (Matrix a)
matrixExp :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> Either MathError (Matrix a)
matrixExp a
tol Int
budget Matrix a
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"matrixExp: positive iteration budget required")
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
        let norm :: a
norm = Matrix a -> a
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix a
a
        if a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
norm Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
norm
            then MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
            else do
                let scaling :: Int
scaling = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (if a
norm a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0.5 then Int
0 else a -> Int
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (a -> a -> a
forall a. Floating a => a -> a -> a
P.logBase a
2 a
norm) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 :: Int)
                    small :: Matrix a
small = a -> Matrix a -> Matrix a
forall a. Num a => a -> Matrix a -> Matrix a
scaleMatrix (Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
1 (-Int
scaling)) Matrix a
a
                i <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
r
                result <- series small 1 i i
                foldM (\Matrix a
x Int
_ -> Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked Matrix a
x Matrix a
x Either MathError (Matrix a)
-> (Matrix a -> Either MathError (Matrix a))
-> Either MathError (Matrix a)
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Matrix a -> Either MathError (Matrix a)
forall a. RealFloat a => Matrix a -> Either MathError (Matrix a)
checked) result [1 .. scaling]
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    checked :: Matrix a -> Either MathError (Matrix a)
checked Matrix a
x = Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
x Either MathError ()
-> Either MathError (Matrix a) -> Either MathError (Matrix a)
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right Matrix a
x
    series :: Matrix a
-> Int -> Matrix a -> Matrix a -> Either MathError (Matrix a)
series Matrix a
small Int
k Matrix a
term Matrix a
total
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
budget = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"matrixExp" Int
budget)
        | Bool
otherwise = do
            productTerm <- Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked Matrix a
term Matrix a
small
            let next = a -> Matrix a -> Matrix a
forall a. Num a => a -> Matrix a -> Matrix a
scaleMatrix (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k) Matrix a
productTerm
            totalNext <- addMatrixChecked total next >>= checked
            if matrixNorm1 next <= tol * max 1 (matrixNorm1 totalNext)
                then Right totalNext
                else series small (k + 1) next totalNext

symmetricSpectral ::
    (Double -> Either MathError Double) -> Matrix Double -> Either MathError (Matrix Double)
symmetricSpectral :: (Double -> Either MathError Double)
-> Matrix Double -> Either MathError (Matrix Double)
symmetricSpectral Double -> Either MathError Double
f Matrix Double
a
    | Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
c = MathError -> Either MathError (Matrix Double)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
r, Int
c))
    | ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any
        (\(Double
x, Double
y) -> Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-12 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
x) (Double -> Double
forall a. Num a => a -> a
abs Double
y)))
        ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
a) (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList (Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a
transposeMatrix Matrix Double
a))) =
        MathError -> Either MathError (Matrix Double)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"spectral function: symmetric matrix required")
    | Bool
otherwise = do
        eig <- Matrix Double -> Either MathError Eigen
symmetricEigen Matrix Double
a
        values <- traverse (f . realPart) (eigenvalues eig)
        if any ((> 1e-10) . abs . imagPart) (eigenvalues eig)
            then Left (DomainError "spectral function: nonreal spectrum")
            else do
                let q = (Complex Double -> Double)
-> Matrix (Complex Double) -> Matrix Double
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries Complex Double -> Double
forall a. Complex a -> a
realPart (Eigen -> Matrix (Complex Double)
eigenvectors Eigen
eig)
                d <- diagonalChecked values
                qd <- multiplyChecked q d
                result <- multiplyChecked qd (transposeMatrix q)
                finiteMatrix result >> Right result
  where
    (Int
r, Int
c) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a

{- | Compute principal real matrix logarithm of a symmetric positive-definite Double matrix
using symmetric eigendecomposition. Reject non-square/asymmetric matrices, nonpositive
spectrum, non-real decomposition, or numeric failure.
-}
matrixLog :: Matrix Double -> Either MathError (Matrix Double)
matrixLog :: Matrix Double -> Either MathError (Matrix Double)
matrixLog =
    (Double -> Either MathError Double)
-> Matrix Double -> Either MathError (Matrix Double)
symmetricSpectral
        ( \Double
x ->
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
                then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"matrixLog: positive definite matrix required")
                else Double -> Either MathError Double
forall a b. b -> Either a b
Right (Double -> Double
forall a. Floating a => a -> a
P.log Double
x)
        )

{- | Compute symmetric positive-semidefinite matrix square root by spectral decomposition;
tolerate eigenvalues down to -1e-12 and clamp them to zero. Invalid
shape/symmetry/spectrum or decomposition failure returns Left MathError.
-}
matrixSqrt :: Matrix Double -> Either MathError (Matrix Double)
matrixSqrt :: Matrix Double -> Either MathError (Matrix Double)
matrixSqrt =
    (Double -> Either MathError Double)
-> Matrix Double -> Either MathError (Matrix Double)
symmetricSpectral
        ( \Double
x ->
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
1e-12
                then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"matrixSqrt: positive semidefinite matrix required")
                else Double -> Either MathError Double
forall a b. b -> Either a b
Right (Double -> Double
forall a. Floating a => a -> a
P.sqrt (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
x))
        )

{- | Return (U,P) from SVD so A=U*P, with P positive semidefinite; rectangular input uses
the thin polar factors. Validation/decomposition/multiplication errors propagate.
-}
polar :: Matrix Double -> Either MathError (Matrix Double, Matrix Double)
polar :: Matrix Double -> Either MathError (Matrix Double, Matrix Double)
polar Matrix Double
a = do
    f <- Matrix Double -> Either MathError SVD
svd Matrix Double
a
    u <- multiplyChecked (svdU f) (svdVT f)
    d <- diagonalChecked (singularValues f)
    vd <- multiplyChecked (transposeMatrix (svdVT f)) d
    p <- multiplyChecked vd (svdVT f)
    Right (u, p)

{- | Complex128 Moore-Penrose pseudoinverse, thresholded relative to the largest
singular value. Both rectangular and rank-deficient inputs are supported.
-}
complexPseudoinverse ::
    Double -> Matrix (Complex Double) -> Either MathError (Matrix (Complex Double))
complexPseudoinverse :: Double
-> Matrix (Complex Double)
-> Either MathError (Matrix (Complex Double))
complexPseudoinverse Double
tol Matrix (Complex Double)
a = do
    Double -> Either MathError ()
validTol Double
tol
    f <- Matrix (Complex Double) -> Either MathError ComplexSVD
complexSVD Matrix (Complex Double)
a
    let values = ComplexSVD -> [Double]
complexSingularValues ComplexSVD
f; biggest = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
values)
    d <-
        diagonalChecked [if x > 0 && (tol == 0 || x / biggest > tol) then (1 / x) :+ 0 else 0 | x <- values]
    du <- multiplyChecked d (conjugateTranspose (complexU f))
    result <- multiplyChecked (conjugateTranspose (complexVT f)) du
    finiteValues (concatMap (\Complex Double
z -> [Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z, Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z]) (toList result))
    Right result

{- | Principal square root for invertible real matrices whose spectra avoid
the nonpositive real axis. This nonsymmetric reference method uses scaled
Denman-Beavers iteration with private inverse solves. The budget bounds the
number of coupled iterations; failure to meet the reconstruction tolerance
returns NonConvergence. Ill-conditioned inverse steps may report a LU error.
The existing matrixSqrt additionally handles symmetric singular PSD inputs.
-}
matrixSqrtGeneral :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
matrixSqrtGeneral :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
matrixSqrtGeneral Double
tolerance Int
budget Matrix Double
a = do
    Double -> Int -> Either MathError ()
generalConfiguration Double
tolerance Int
budget
    Matrix Double -> Either MathError ()
principalDomain Matrix Double
a
    Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
squareRootCore Double
tolerance Int
budget Matrix Double
a

{- | Real principal matrix logarithm on the same spectral domain as
matrixSqrtGeneral. Inverse scaling by successive principal square roots makes
the Mercator series converge in matrix 1-norm. The budget bounds each phase:
root reductions, each inner Denman-Beavers solve, series terms, and the final
exponential reconstruction check. There is no unbounded convergence loop.
-}
matrixLogGeneral :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
matrixLogGeneral :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
matrixLogGeneral Double
tolerance Int
budget Matrix Double
a = do
    Double -> Int -> Either MathError ()
generalConfiguration Double
tolerance Int
budget
    Matrix Double -> Either MathError ()
principalDomain Matrix Double
a
    ident <- Int -> Either MathError (Matrix Double)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
n
    if n == 0
        then Right a
        else do
            (scaled, rootCount) <- reduce ident 0 a
            delta <- subtractMatrixChecked scaled ident
            logarithm <- series delta 1 delta (scaleMatrix 0 delta)
            let result = (Double -> Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries (Int -> Double -> Double
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
rootCount) Matrix Double
logarithm
            finiteMatrix result
            reconstructed <- matrixExp (min 1e-14 (tolerance / 16)) budget result
            residual <- relativeMatrixResidual reconstructed a
            if residual <= tolerance
                then Right result
                else Left (NonConvergence "matrixLogGeneral: exponential residual" budget)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a)
    innerTolerance :: Double
innerTolerance = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1e-14 (Double
tolerance Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
32)
    reduce :: Matrix Double
-> Int -> Matrix Double -> Either MathError (Matrix Double, Int)
reduce Matrix Double
ident Int
k Matrix Double
current = do
        delta <- Matrix Double -> Matrix Double -> Either MathError (Matrix Double)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
subtractMatrixChecked Matrix Double
current Matrix Double
ident
        let distanceFromIdentity = Matrix Double -> Double
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix Double
delta
        if isNaN distanceFromIdentity || isInfinite distanceFromIdentity
            then Left NonFiniteResult
            else
                if distanceFromIdentity <= 0.5
                    then Right (current, k)
                    else
                        if k >= budget
                            then Left (NonConvergence "matrixLogGeneral: inverse scaling" budget)
                            else squareRootCore innerTolerance budget current >>= reduce ident (k + 1)
    series :: Matrix Double
-> Int
-> Matrix Double
-> Matrix Double
-> Either MathError (Matrix Double)
series Matrix Double
delta Int
k Matrix Double
power Matrix Double
total
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
budget = MathError -> Either MathError (Matrix Double)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"matrixLogGeneral: Mercator series" Int
budget)
        | Bool
otherwise = do
            let coefficient :: Double
coefficient = (if Int -> Bool
forall a. Integral a => a -> Bool
odd Int
k then Double
1 else -Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k
                term :: Matrix Double
term = Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> Matrix a -> Matrix a
scaleMatrix Double
coefficient Matrix Double
power
            next <- Matrix Double -> Matrix Double -> Either MathError (Matrix Double)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
addMatrixChecked Matrix Double
total Matrix Double
term
            finiteMatrix next
            let termNorm = Matrix Double -> Double
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix Double
term
            if isNaN termNorm || isInfinite termNorm
                then Left NonFiniteResult
                else
                    if termNorm <= innerTolerance * max 1 (matrixNorm1 next)
                        then Right next
                        else do
                            nextPower <- multiplyChecked power delta
                            finiteMatrix nextPower
                            series delta (k + 1) nextPower next

generalConfiguration :: Double -> Int -> Either MathError ()
generalConfiguration :: Double -> Int -> Either MathError ()
generalConfiguration Double
tolerance Int
budget
    | Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
tolerance Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
tolerance = MathError -> Either MathError ()
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"general matrix function: positive iteration budget required")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
principalDomain :: Matrix Double -> Either MathError ()
principalDomain :: Matrix Double -> Either MathError ()
principalDomain Matrix Double
a = do
    spectrum <- Matrix Double -> Either MathError Eigen
eigen Matrix Double
a
    if any (\Complex Double
value -> Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
value Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
&& Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
value Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0) (eigenvalues spectrum)
        then
            Left
                ( DomainError
                    "principal real matrix function: invertible matrix with no nonpositive-real-axis eigenvalues required"
                )
        else Right ()

squareRootCore :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
squareRootCore :: Double -> Int -> Matrix Double -> Either MathError (Matrix Double)
squareRootCore Double
tolerance Int
budget Matrix Double
a
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Matrix Double -> Either MathError (Matrix Double)
forall a b. b -> Either a b
Right Matrix Double
a
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
inputNorm Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
inputNorm = MathError -> Either MathError (Matrix Double)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Double
inputNorm Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = MathError -> Either MathError (Matrix Double)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"matrixSqrtGeneral: invertible matrix required")
    | Bool
otherwise = do
        let scaling :: Int
scaling = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Double -> Double
forall a. Floating a => a -> a -> a
P.logBase Double
4 Double
inputNorm) :: Int
            scaled :: Matrix Double
scaled = (Double -> Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries (Int -> Double -> Double
forall a. RealFloat a => Int -> a -> a
scaleFloat (-Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
scaling)) Matrix Double
a
        Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix Double
scaled
        ident <- Int -> Either MathError (Matrix Double)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
n
        root <- iterateRoot scaled 0 scaled ident
        let result = (Double -> Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries (Int -> Double -> Double
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
scaling) Matrix Double
root
        finiteMatrix result
        squaredResult <- multiplyChecked result result
        residual <- relativeMatrixResidual squaredResult a
        if residual <= max tolerance (32 * epsilon (0 :: Double))
            then Right result
            else Left (NonConvergence "matrixSqrtGeneral: reconstruction residual" budget)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a)
    inputNorm :: Double
inputNorm = Matrix Double -> Double
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix Double
a
    iterateRoot :: Matrix Double
-> Int
-> Matrix Double
-> Matrix Double
-> Either MathError (Matrix Double)
iterateRoot Matrix Double
target Int
k Matrix Double
y Matrix Double
z = do
        squared <- Matrix Double -> Matrix Double -> Either MathError (Matrix Double)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked Matrix Double
y Matrix Double
y
        residual <- relativeMatrixResidual squared target
        if residual <= tolerance
            then Right y
            else
                if k >= budget
                    then Left (NonConvergence "matrixSqrtGeneral: Denman-Beavers" budget)
                    else do
                        inverseZ <- inverse z
                        inverseY <- inverse y
                        nextY <- scaleMatrix 0.5 <$> addMatrixChecked y inverseZ
                        nextZ <- scaleMatrix 0.5 <$> addMatrixChecked z inverseY
                        finiteMatrix nextY
                        finiteMatrix nextZ
                        iterateRoot target (k + 1) nextY nextZ

-- Divide entries before subtraction to avoid overflowing a finite residual
-- when two very large finite matrices have opposite-signed coefficients.
relativeMatrixResidual :: Matrix Double -> Matrix Double -> Either MathError Double
relativeMatrixResidual :: Matrix Double -> Matrix Double -> Either MathError Double
relativeMatrixResidual Matrix Double
actual Matrix Double
expected = do
    Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix Double
actual
    Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix Double
expected
    let scale :: Double
scale = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
expected)))
    delta <- (Double -> Double -> Double)
-> Matrix Double
-> Matrix Double
-> Either MathError (Matrix Double)
forall a b c.
(a -> b -> c)
-> Matrix a -> Matrix b -> Either MathError (Matrix c)
zipMatrixWith (\Double
x Double
y -> Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Matrix Double
actual Matrix Double
expected
    let numerator = Matrix Double -> Double
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 Matrix Double
delta
        denominator = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) (Matrix Double -> Double
forall a. (Ord a, Num a) => Matrix a -> a
matrixNorm1 ((Double -> Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) Matrix Double
expected))
        residual = Double
numerator Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denominator
    if isNaN residual || isInfinite residual then Left NonFiniteResult else Right residual