linearAlgebra.jl

function typeofMatrix(
array::Union{AnyMatrix, AnyMatrixVector, AnyMatrixVector₂})

alias: typeofMat

Return the type of a matrix, either Hermitian, Diagonal, LowerTriangular, or Matrix. Argument array may be a matrix of one of these types, but also one of the following:

ℍVector, ℍVector₂, 𝔻Vector, 𝔻Vector₂, 𝕃Vector, 𝕃Vector₂, 𝕄Vector, 𝕄Vector₂.

Those are Array of Matrices types. See also aliases for the symbols ℍ, 𝔻, 𝕃 and 𝕄.

Note that this function is different from Julia function typeof, which returns the concrete type (see example below), thus cannot be used for typecasting matrices.

Examples

using LinearAlgebra, PosDefManifold
P=randP(3) # generate a 3x3 Hermitian matrix
typeofMatrix(P) # returns `Hermitian`
typeof(P) # returns `Hermitian{Float64,Array{Float64,2}}`
# typecast P as a `Matrix` M
M=Matrix(P)
# typecast M as a matrix of the same type as P and write the result in A
A=typeofMatrix(P)(M)

Pset=randP(3, 4) # generate a set of 4 3x3 Hermitian matrix
# Pset is an ℍVector type
typeofMatrix(Pset) # again returns `Hermitian`

function typeofVector(
array::Union{AnyMatrix, AnyMatrixVector, AnyMatrixVector₂})

alias: typeofVec

Return the type of a Vector, either HermitianVector, DiagonalVector, LowerTriangularVector, or MatrixVector. The aliases of those are, respectvely, ℍVector, 𝔻Vector, 𝕃Vector and 𝕄Vector. Argument array may be a vector of one of these types, but also one of the following:

ℍ, 𝔻, 𝕃 and 𝕄, ℍVector₂, 𝔻Vector₂, 𝕃Vector₂, 𝕄Vector₂.

See aliases for the symbols ℍ, 𝔻, 𝕃 and 𝕄. The last four are Array of Matrices types.

Note that this function is different from Julia function typeof only in that it returns the vector type also if array is not of the ℍVector, 𝔻Vector, 𝕃Vector or 𝕄Vector type.

Examples

using LinearAlgebra, PosDefManifold
P=randP(3, 4) # generate 4 3x3 Hermitian matrix
typeofMatrix(P) # returns `Array{Hermitian,1}`
typeof(P) # also returns `Array{Hermitian,1}`

typeofMatrix(P[1]) # returns `Array{Hermitian,1}`
typeof(P[1]) # returns `Hermitian{Float64,Array{Float64,2}}`

(1) function dim(X::AnyMatrix, [d])
(2) function dim(vector::AnyMatrixVector, [d])
(3) function dim(vector₂::AnyMatrixVector₂, [d])

(1) $X$ is a real or complex Matrix, Diagonal, LowerTriangular or Hermitian matrix. Return a 2-tuple containing the dimensions of $X$, which is two times the same dimension for all possible types of $X$ with the exception of the Matrix type, which can be rectangular. Optionally you can specify a dimension (1 or 2) to get just the length of that dimension.

(2) vector is an 𝕄Vector, 𝔻Vector, 𝕃Vector or ℍVector type (see AnyMatrixVector type). Return a 3-tuple containing the number of matrices it holds (dimension 1) and their dimensions (dimension 2 and 3). Optionally you can specify a dimension (1, 2, or 3) to get just the length of that dimension.

(3) vector₂ is an 𝕄Vector₂, 𝔻Vector₂, 𝕃Vector₂ or ℍVector₂ type (see AnyMatrixVector type). Return a 4-tuple containing

• the number of vectors of matrices it holds (dimension 1),
• a vector holding the number of matrices in each vector of matrices (dimensions 2),
• the two dimensions of the matrices (dimension 3 and 4).

Optionally you can specify a dimension (1, 2, 3 or 4) to get just the length of that dimension.

vector and vector₂ are Array of Matrices types. See also aliases for the symbols ℍ, 𝔻, 𝕃 and 𝕄.

This function replaces Julia size function, which cannot be used to retrive dimension for matrix vectors. It is not possible to overload the size function for matrix vectors since this causes problems to other Julia functions.

Examples

using LinearAlgebra, PosDefManifold
# (1)
M=randn(3, 4) # generate a 3x4 `Matrix`
dim(M) # returns (3, 4)
dim(M, 1) # returns 3
dim(M, 2) # returns 4
dim(M, 3) # out of range: returns 0

# (2)
Pset=randP(3, 4) # generate an ℍVector holding 4 3x3 Hermitian matrices
dim(Pset) # returns (4, 3, 3)
dim(Pset, 1) # returns 4
dim(Pset, 2) # returns 3
dim(Pset, 3) # returns 3

# (3)
# Generate a set of 4 random 3x3 SPD matrices
Pset=randP(3, 4)
# Generate a set of 40 random 4x4 SPD matrices
Qset=randP(3, 40)
A=ℍVector₂([Pset, Qset])
dim(A) # return (2, [4, 40], 3, 3)
dim(A, 1) # return 2
dim(A, 2) # return [4, 40]
dim(A, 2)[1] # return 4
dim(A, 3) # return 3
dim(A, 4) # return 3
dim(A, 5) # out of range: return 0

# note: to create an ℍVector₂ object holding k ℍVector objects use
sets=ℍVector₂(undef, k) # and then fill them

function remove(X::Union{Vector, Matrix}, what::Union{Int, Vector{Int}};
				dims=1)

Remove one or more elements from a vector or one or more columns or rows from a matrix.

If X is a Matrix, dims=1 (default) remove rows, dims=2 remove columns.

If X is a Vector, dims has no effect.

The second argument is either an integer or a vector of integers.

Examples

a=randn(5)
b=remove(a, 2)
b=remove(a, collect(1:3)) # remove rows 1 to 3
A=randn(3, 3)
B=remove(A, 2)
B=remove(A, 2; dims=2)
A=randn(5, 5)
B=remove(A, collect(1:2:5)) # remove rows 1, 3 and 5
C=remove(A, [1, 4])
A=randn(10, 10)
A=remove(A, [collect(2:3); collect(8:10)]; dims=2)

function isSquare(X::Matrix)=size(X, 1)==size(X, 2)

Return true if matrix X is square, false otherwise.

function det1(X::AnyMatrix; <tol::Real=0.>)

Return the argument matrix $X$ normalized so as to have unit determinant. For square positive definite matrices this is the best approximant from the set of matrices in the special linear group - see Bhatia and Jain (2014)🎓.

$X$ can be a real or complex Diagonal, LowerTriangular, Matrix, or Hermitian matrix. (see AnyMatrix type)

If the determinant is not greater than tol (which defalts to zero) a warning is printed and $X$ is returned.

See Julia det function.

See also: tr1.

Examples

using LinearAlgebra, PosDefManifold
P=randP(5) # generate a random real positive definite matrix 5x5
Q=det1(P)
det(Q) # must be 1
# using a tolerance
Q=det1(P; tol=1e-12)

tr1(X::AnyMatrix; tol::Real=0.)

Return the argument matrix $X$ normalized so as to have unit trace.

$X$ can be a real or complex Diagonal, LowerTriangular, Matrix or Hermitian matrix (see AnyMatrix type). Its trace must be real. If the absolute value of its imaginary part is greater than tol (which defalts to zero) a warning is printed and $X$ is returned. Also, if the trace is not greater than tol a warning is printed and $X$ is returned.

See: Julia trace function.

See also: tr, det1.

Examples

using LinearAlgebra, PosDefManifold

P=randP(5) # generate a random real positive definite matrix 5x5
Q=tr1(P)
tr(Q)  # must be 1
# using a tolerance
Q=tr1(P; tol=1e-12)

Pc=randP(ComplexF64, 5) # generate a random real positive definite matrix 5x5
Qc=tr1(Pc)
tr(Qc)  # must be 1

nearestPosDef(X::Union{𝔻, 𝕄}; tol::Real=0.)

Return the nearest symmetric/Hermitian positive semi-definite matrix of a diagonal or of an arbitary square matrix X according to the Frobenius norm. If the eigenvalues of the symmetric part of X are all non-negative, the result is positive definite and will be flagged as Hermitian, otherwise it is positive semi-definite and will not be flagged. The nearest matrix is given by

$(Y+H)/2$

where

$Y=(X+X^H)/2$

is the symmetric part of $X$, and $H$ is the symmetric polar factor of $Y$. See Higham(1988)🎓 for details and for the way it is computed.

See also: det1, procrustes.

Examples

using LinearAlgebra, PosDefManifold
X=randn(5, 5) # generate an arbitrary 5x5 matrix
S=nearestPosDef(X)

P=randP(5) # generate a random real positive definite 5x5 matrix
S=nearestPosDef(Matrix(P)) # typecasting an Hermitian matrix as a `Matrix`
# Since P is a positive definite matrix S must be equal to P
S ≈ P ? println(" ⭐ ") : println(" ⛔ ")

nearestOrthogonal(X::AnyMatrix)

alias: nearestOrth

Return the nearest orthogonal matrix of a square Hermitian, LowerTriangular, Diagonal or generic Matrix X (see AnyMatrix type). This is given by

$UV^H$,

where

$\textrm(SVD)=UΛV^H$.

If X is Diagonal, return X.

See also: nearestPosDef, procrustes.

Examples

using PosDefManifold
U=nearestOrth(randn(5, 5))

(1) normalizeCol!(X::𝕄{T}, j::Int)
(2) normalizeCol!(X::𝕄{T}, j::Int, by::Number)
(3) normalizeCol!(X::𝕄{T}, range::UnitRange)
(4) normalizeCol!(X::𝕄{T}, range::UnitRange, by::Number)
for all above: where T<:RealOrComplex

Given a Matrix type $X$ comprised of real or complex elements,

• (1) normalize the $j^{th}$ column to unit norm
• (2) divide the elements of the $j^{th}$ column by number $by$
• (3) normalize the columns in $range$ to unit norm
• (4) divide the elements of columns in $range$ by number $by$.

$by$ is a number of abstract supertype Number. It should be an integer, real or complex number. For efficiency, it should be of the same type as the elements of $X$.

$range$ is a UnitRange type.

Methods (1) and (3) call the BLAS.nrm2 routine for computing the norm of concerned columns. See Threads.

See norm and also randn for the example below.

See also: colNorm, colProd.

Examples

using PosDefManifold
X=randn(10, 20)
normalizeCol!(X, 2)                  # (1) normalize columns 2
normalizeCol!(X, 2, 10.0)            # (2) divide columns 2 by 10.0
normalizeCol!(X, 2:4)                # (3) normalize columns 2 to 4
X=randn(ComplexF64, 10, 20)
normalizeCol!(X, 3)                  # (1) normalize columns 3
normalizeCol!(X, 3:6, (2.0 + 0.5im)) # (4) divide columns 3 to 5 by (2.0 + 0.5im)

    (1) ispos(λ::Vector{T};
	<
	tol::Real=0,
	rev=true,
	🔔=true,
	msg="">)

    (2) ispos(Λ::𝔻{T};
	< same optional keyword arguments as in (1) > )

	for all above: where T<:Real

Return $true$ if all numbers in (1) real vector $λ$ or in (2) real Diagonal matrix $Λ$ are not inferior to $tol$, otherwise return $false$. This is used, for example, in spectral functions to check that all eigenvalues are positive.

The following are <optional keyword arguments>:

• If $rev=true$ the (1) elements in $λ$ or (2) the diagonal elements

in $Λ$ will be chacked in reverse order. This is done for allowing a very fast check when the elements are sorted and it is known from where is best to start checking.

If the result is $false$:

• if $🔔=true$ a bell character will be printed. In most systems this will ring a bell on the computer.
• if string $msg$ is provided, a warning will print $msg$ followed by:

"at position pos", where pos is the position where the first non-positive element has been found.

Examples

using PosDefManifold
a=[1, 0, 2, 8]
ispos(a, msg="non-positive element found")

# it will print:
# ┌ Warning: non-positive element found at position 2
# └ @ [here julie will point to the line of code issuing the warning]

(1) colProd(X::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int)
(2) colProd(X::Union{𝕄{T}, ℍ{T}}, Y::Union{𝕄{T}, ℍ{T}}, j::Int, l::Int)
for all above: where T<:RealOrComplex

(1) Given a real or complex Matrix or Hermitian matrix $X$, return the dot product of the $j^{th}$ and $l^{th}$ columns, defined as,

$\sum_{i=1}^{r} \big(x_{ij}^*x_{il}\big),$

where $r$ is the number of rows of $X$ and $^*$ denotes complex conjugate (nothing if the matrix is real).

(2) Given real or complex Matrix or Hermitian matrices $X$ and $Y$, return the dot product of the $j^{th}$ column of $X$ and the $l^{th}$ column of $Y$, defined as,

$\sum_{i=1}^{r} \big(x_{ij}^*y_{il}\big),$

where $r$ is the number of rows of $X$ and of $Y$ and $^*$ is as above.

Arguments $j$ and $l$ must be positive integers in range

• (1) j,l in 1:size(X, 2),
• (2) j in 1:size(X, 2), l in 1:size(Y, 2).

See also: normalizeCol!, colNorm.

Examples

using PosDefManifold
X=randn(10, 20)
p=colProd(X, 1, 3)
Y=randn(10, 30)
q=colProd(X, Y, 2, 25)

(1) sumOfSqr(A::Array)
(2) sumOfSqr(H::ℍ{T})
(3) sumOfSqr(L::𝕃{T})
(4) sumOfSqr(D::𝔻{T})
(5) sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, j::Int)
(6) sumOfSqr(X::Union{𝕄{T}, ℍ{T}}, range::UnitRange)
for (1)-(6) above: where T<:RealOrComplex

alias: ss

Return

• (1) the sum of squares of the elements in an array $A$ of any dimensions.
• (2) as in (1), but for an Hermitian matrix $H$, using only the lower triangular part.
• (3) as in (1), but for a LowerTriangular matrix $L$.
• (4) as in (1), but for a Diagonal matrix $D$ (sum of squares of diagonal elements).
• (5) the sum of square of the $j^{th}$ column of a Matrix or Hermitian $X$.
• (6) the sum of square of the columns of a Matrix or Hermitian $X$ in a given range.

All methods support real and complex matrices.

Only method (1) works for arrays of any dimensions.

Methods (1)-(4) return the square of the Frobenius norm.

For method (5), $j$ is a positive integer in range 1:size(X, 1).

For method (6), $range$ is a UnitRange type.

See also: colNorm, sumOfSqrDiag, sumOfSqrTril.

Examples

using PosDefManifold
X=randn(10, 20)
sum2=sumOfSqr(X)        # (1) sum of squares of all elements
sum2=sumOfSqr(X, 1)     # (2) sum of squares of elements in column 1
sum2=sumOfSqr(X, 2:4)   # (3) sum of squares of elements in column 2 to 4

sumOfSqrDiag(X::AnyMatrix)

alias: ssd

Sum of squares of the diagonal elements in real or complex Matrix, Diagonal, Hermitian or LowerTriangular matrix $X$. If $X$ is rectangular (which can be only if it is of the Matrix type), the main diagonal is considered.

See AnyMatrix type

See also: sumOfSqr, sumOfSqrTril.

Examples

using LinearAlgebra, PosDefManifold
X=randn(10, 20)
sumDiag2=sumOfSqrDiag(X) # (1)
sumDiag2=sumOfSqrDiag(𝔻(X)) # (2)
# 𝔻=LinearAlgebra.Diagonal is declated in the main module

colNorm(X::Union{𝕄{T}, ℍ{T}}, j::Int) where T<:RealOrComplex

Given a real or complex Matrix or Hermitian matrix $X$, return the Euclidean norm of its $j^{th}$ column.

This function calls the BLAS.nrm2 routine. See Threads.

See also: normalizeCol!, colProd, sumOfSqr.

Examples

using PosDefManifold
X=randn(10, 20)
normOfSecondColumn=colNorm(X, 2)

sumOfSqrTril(X::AnyMatrix, k::Int=0)

alias: sst

Given a real or complex Matrix, Diagonal, Hermitian or LowerTriangular matrix $X$ (see AnyMatrix type), return the sum of squares of the elements in its lower triangle up to the $k^{th}$ underdiagonal.

Matrix $X$ may be rectangular.

$k$ must be in range

• 1-size(X, 1):c-1 for $X$ Matrix, Diagonal or Hermitian,
• 1-size(X, 1):0 for $X$ LowerTriangular.

For $X$ Diagonal the result is

• $0$ if $k<0$,
• the sum of the squares of the diagonal elements otherwise.

See julia tril(M, k::Integer) function for numbering of diagonals.

See also: sumOfSqr, sumOfSqrDiag.

Examples

using PosDefManifold
A=[4. 3.; 2. 5.; 1. 2.]
#3×2 Array{Float64,2}:
# 4.0  3.0
# 2.0  5.0
# 1.0  2.0

s=sumOfSqrTril(A, -1)
# 9.0 = 1²+2²+2²

s=sumOfSqrTril(A, 0)
# 50.0 = 1²+2²+2²+4²+5²

(1) tr(P::ℍ{T}, Q::ℍ{T})
(2) tr(P::ℍ{T}, M::𝕄{T})
(3) tr(D::𝔻{T}, H::Union{ℍ{T}, 𝕄{T}})
(4) tr(H::Union{ℍ{T}, 𝕄{T}}, D::𝔻{T})
for all above: where T<:RealOrComplex

Given (1) two Hermitian positive definite matrix $P$ and $Q$, return the trace of the product $PQ$. This is real even if $P$ and $Q$ are complex.

$P$ must always be flagged as Hermitian. See typecasting matrices.

In (2) $Q$ is a Matrix object, in which case return

• a real trace if the product $PQ$ is real or if it has all positive real eigenvalues.
• a complex trace if the product $PQ$ is not real and has complex eigenvalues.

Methods (3) and (4) return the trace of the product $DH$ or $HD$, where $D$ is a Diagonal matrix and $H$ an $Hermitian$ or $Matrix$ object. The result is of the same type as the input matrices.

For all methods all arguments must be of the same type.

Math

Let $P$ and $Q$ be Hermitian matrices, using the properties of the trace (e.g., the cyclic property and the similarity invariance) you can use this function to fast compute the trace of several expressions. For example:

$\textrm{tr}(PQ)=\textrm{tr}(P^{1/2}QP^{1/2})$

and

$\textrm{tr}(PQP)=\textrm{tr}(P^{2}Q)$ (see example below).

See: trace.

See also: DiagOfProd, tr1.

Examples

using PosDefManifold
P=randP(ComplexF64, 5) # generate a random complex positive definite matrix 5x5
Q=randP(ComplexF64, 5) # generate a random complex positive definite matrix 5x5
tr(P, Q) ≈ tr(P*Q) ? println(" ⭐ ") : println(" ⛔ ")
tr(P, Q) ≈ tr(sqrt(P)*Q*sqrt(P)) ? println(" ⭐ ") : println(" ⛔ ")
tr(sqr(P), Q) ≈ tr(P*Q*P) ? println(" ⭐ ") : println(" ⛔ ")

(1) quadraticForm(v::Vector{T}, P::ℍ{T}) where T<:Real
(2) quadraticForm(v::Vector{T}, L::𝕃{T}) where T<:Real
(3) quadraticForm(v::Vector{T}, X::𝕄{T}, forceLower::Bool=false) where T<:Real
(4) quadraticForm(v::Vector{S}, X::Union{𝕄{S}, ℍ{S}, 𝕃{S}}) where S<:Complex

alias: qf

(1) Given a real vector $v$ and a real Hermitian matrix $P$, compute the quadratic form

$v^TPv$,

where the superscript T denotes transpose. It uses only the lower triangular part of $P$.

(2) As in (1), given a real vector $v$ and a LowerTriangular matrix $L$.

(3) As in (1), given a real vector $v$ and a real generic Matrix $M$, if forceLower=true. If forceLower=false, the product $v^TMv$ is evaluated instead using the whole matrix $M$.

(4) Quadratic form $v^HPv$, where superscript H denotes complex conjugate and transpose, for a complex vector v and a complex Matrix, LowerTrianglar or Hermitian matrix. The whole matrix is used.

Math

For $v$ and $X$ real and $X$ symmetric, the quadratic form is

$\sum_i(v_i^2x_{ii})+\sum_{i>j}(2v_iv_jx_{ij})$.

For $L$ lower triangular is

$\sum_i(v_i^2x_{ii})+\sum_{i>j}(v_iv_jx_{ij})$.

These formula are used in methods (1), (2) and (3).

Examples

using PosDefManifold
P=randP(5) # generate a random real positive definite matrix 5x5
v=randn(5)
q1=quadraticForm(v, P) # or q1=qf(v, P)
q2=v'*P*v
q1 ≈ q2 ? println(" ⭐ ") : println(" ⛔ ")

fidelity(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex

Given two positive definte Hermitian matrices $P$ and $Q$, return their fidelity:

$tr\big(P^{1/2}QP^{1/2}\big)^{1/2}.$

This is used in quantum physics and is related to the Wasserstein metric. See for example Bhatia, Jain and Lim (2019b)🎓.

Examples

using PosDefManifold
P=randP(5);
Q=randP(5);
f=fidelity(P, Q)

fDiag(func::Function, X::AnyMatrix, k::Int=0)

alias: 𝑓𝔻

Applies function func element-wise to the elements of the $k^{th}$ diagonal of real or complex Diagonal, LowerTriangular, Matrix or Hermitian matrix $X$ and return a diagonal matrix with these elements. $X$ must be square in all cases, but for the 𝕄=Matrix type argument, in which case it may be of dimension r⋅c, with r ≠ c.

See julia tril(M, k::Integer) function for numbering of diagonals.

Bt default the main diagonal is considered.

• If $X$ is Diagonal, $k$ is set automatically to zero (main diagonal).
• If $X$ is LowerTriangular, $k$ cannot be positive.

Note that if $X$ is rectangular the dimension of the result depends on the size of $X$ and on the chosen diagonal. For example,

• r ≠ c and $k$=0 (main diagonal), the result will be of dimension min(r,c)⋅min(r,c),
• $X$ 3⋅4 and $k=-1$, the result will be 2⋅2,
• $X$ 3⋅4 and $k=1$, the result will be 3⋅3, etc.

If the input matrix is complex, the function `func`
must be able to support complex arguments.

See also: DiagOfProd, tr.

Examples

using PosDefManifold
P=randP(5) # use P=randP(ComplexF64, 5) for generating an Hermitian matrix

# diagonal matrix with the inverse of the first sub-diagonal of P
D=fDiag(inv, P, -1)

(Λ, U) = evd(P) # Λ holds the eigenvalues of P, see evd

# diagonal matrix with the log of the eigenvalues
Δ=fDiag(log, Λ)

# using an anonymous function for the square of the eigenvalues
Δ=fDiag(x->x^2, Λ)

DiagOfProd(P::ℍ{T}, Q::ℍ{T}) where T<:RealOrComplex

alias: dop

Return the Diagonal matrix holding the diagonal of the product $PQ$ of two Hermitian matrices P and Q. Only the diagoanl part of the product is computed.

See also: tr, fDiag.

Examples

using PosDefManifold, LinearAlgebra
P, Q=randP(5), randP(5)
DiagOfProd(P, Q)≈Diagonal(P*Q) ? println("⭐ ") : println("⛔ ")

mgs(X::𝕄{T}, numCol::Int=0) where T<:RealOrComplex

Modified (stabilized) Gram-Schmidt orthogonalization of the columns of square or tall matrix $X$, which can be comprised of real or complex elements. The orthogonalized $X$ is returned by the function. $X$ is not changed.

All columns are orthogonalized by default. If instead argument numCol is provided, then only the first numCol columns of $X$ are orthogonalized. In this case only the firt numCol columns will be returned.

Examples

using LinearAlgebra, PosDefManifold
X=randn(10, 10);
U=mgs(X)        # result is 10⋅10
U=mgs(X, 3)     # result is 10⋅3
U'*U ≈ I ? println(" ⭐ ") : println(" ⛔ ")
# julia undertands also:
U'U ≈ I ? println(" ⭐ ") : println(" ⛔ ")

	(1) fVec(f::Function, 𝐏::AnyMatrixVector;
	<
	w::Vector=[],
	✓w=false,
	allocs=[])
	>

	(2) fVec(f::Function, g::Function, 𝐏::AnyMatrixVector;
	< same optional keyword arguments in (1) >)

Given a 1d array $𝐏={P_1,...,P_k}$ of $k$ matrices of the 𝕄Vector type, 𝔻Vector type, 𝕃Vector type or ℍVector type and an optional non-negative real weights vector $w={w_1,...,w_k}$, return expression

$(1)\hspace{6pt}f_{i=1}^{k}(w_iP_i)$,

or

$(2)\hspace{6pt}f_{i=1}^{k}(w_ig(P_i))$,

where $f$ is either the mean or the sum standard julia functions and $g$ is whatever matrix function applying to each matrix $P_k$, such as exp, log,sqrt`, etc, and anonymous functions.

This function is multi-threaded. It works by partitioning the $k$ operations required by the $f$ function in several groups, passing each group to a separate thread and combining the result of the intermediate operations. This function allows a gain in computational time only when the number of matrices (1) and/or their size (2) is high. Use mean and sum otherwise. The maximal gain is obtained when the number of matrices in 𝐏 is an exact multiple of the number of threads Julia is instructed to use. For this latter, see Threads.

!!! note "Nota Bene"

 Contrarily to Julia `mean` and `sum` function (v 1.1.0) the `fVec` function
 returns a matrix of the same type of the matrices in ``𝐏``.

<optional keword argument> allocs allows to pass pre-allocated memory for holding the intermediate result of each thread. Argument allocs must be a vector of as many matrices as threads and where the matrices have the same dimension as the the matrices in $𝐏$ (see the example here below). Using this option is worthwhile only if the size of the matrices is very high and/or when fVec is to be called repeatedly on many vector of matrices, where the matrices have always the same size, so that one allocation works for all calls.

If <optional keyword argument> ✓w=true is passed, the weights are normalized so as to sum up to 1, otherwise they are used as they are passed. This option is provided to allow calling this function repeatedly without normalizing the same weights vector each time. By default ✓w is false.

Examples

using LinearAlgebra, PosDefManifold
Pset=randP(4, 1000); # generate 1000 positive definite 4x4 matrices
mean(Pset) # arithmetic mean calling Julia function
Threads.nthreads() # check how many threads are available
fVec(mean, Pset) # multi-threaded arithmetic mean

inv(mean(inv, Pset)) # Harmonic mean calling Julia function
inv(fVec(mean, inv, Pset)) # multi-threaded Harmonic mean

exp(mean(log, Pset)) # log Euclidean mean calling Julia function
exp(fVec(mean, log, Pset)) # multi-threaded log Euclidean mean

# notice that Julia `exp` function has changed the type of the result
# to `Symmetric`. To obtain an `Hermitian` output use
ℍ(exp(fVec(mean, log, Pset)))

w=(randn(1000)).^2
w=w./sum(w)  		# generate normalized random weights

# weighted arithmetic mean calling Julia function
sum(Pset[i]*w[i] for i=1:length(w))
# multi-threaded weighted arithmetic mean
fVec(sum, Pset, w=w)

# weighted harmonic mean calling Julia function
inv(sum(inv(Pset[i])*w[i] for i=1:length(w)))
# multi-threaded weighted harmonic mean
inv(fVec(sum, inv, Pset, w=w))

# pre-allocating memory
Pset=randP(100, 1000); # generate 1000 positive definite 100x100 matrices
Qset=MatrixVector(repeat([similar(Pset[1])], Threads.nthreads()))
fVec(mean, log, Pset, allocs=Qset)

# How much computing time we save ?
# (example min time obtained with 4 threads & 4 BLAS threads)
using BenchmarkTools
# standard Julia function
@benchmark(mean(log, Pset)) 					# (5.271 s)
# fVec
@benchmark(fVec(mean, log, Pset))				# (1.540 s)

(1) congruence(B::AnyMatrix, P::AnyMatrix, matrixType)
(2) congruence(B::AnyMatrix, 𝐏::AnyMatrixVector, matrixVectorType)
(3) congruence(B::AnyMatrix, 𝑷::AnyMatrixVector₂, matrixVector₂Type)
(4) congruence(𝐁::AnyMatrixVector, 𝑷::AnyMatrixVector₂, matrixVector₂Type)

alias: cong

(1) Return the congruent transformation

$BPB^H$,

for $B$ and $P$ any combination of Hermitian, LowerTriangular, Diagonal or general Matrix type.

The result is of the matrixType argument, which must be provided and must be one of these four abstract type (not an instance of them). See aliases for shortening these type using symbols ℍ, 𝔻, 𝕃 and 𝕄.

(2) Return a vector of matrices holding the congruent transformations

$BP_kB^H$,

for all $k$ matrices in $𝐏={P_1,...,P_k}$, for $B$ and $𝐏$ any combination of matrix type Hermitian, LowerTriangular, Diagonal or Matrix ($B$) and vector of matrices type ℍVector, 𝔻Vector, 𝕃Vector and 𝕄Vector ($𝐏$). See Array of Matrices types.

The result is a vector of matrices of the matrixVectorType argument, which must be provided and must be one of the following abstract types: ℍVector, 𝔻Vector, 𝕃Vector or 𝕄Vector (and not an instance of these types).

(3) Return a vector of vector of matrices holding the congruent transformations

$BP_{mk}B^H$,

for all $m$ vectors of $k[m]$ vectors of matrices in $𝑷$, for $B$ and $𝑷$ any combination of matrix type Hermitian, LowerTriangular, Diagonal or Matrix ($B$) and vector of matrices type ℍVector₂, 𝔻Vector₂, 𝕃Vector₂ and 𝕄Vector₂ ($𝑷$). See Array of Matrices types.

The result is a vector of vector of matrices of the matrixVector₂Type argument, which must be provided and must be one of the following abstract types: ℍVector₂, 𝔻Vector₂, 𝕃Vector₂ or 𝕄Vector₂ (and not an instance of these types).

(4) Return a vector of vector of matrices holding the congruent transformations

$B_iP_{ij}B_j^H$, for $i,j∈[1,...,m]$.

for $𝐁$ holding $m$ matrices and $𝑷$ holding $m$ vectors holding $m$ matrices each. Note that, differently from method (3), here the vectors of $𝑷$ are all of the same length and this is eaxctly the length of $𝐁$. $𝐁$ and $𝑷$ may be any combination of matrix vector type ℍVector, 𝔻Vector, 𝕃Vector and 𝕄Vector ($𝐁$) and vector of matrices type ℍVector₂, 𝔻Vector₂, 𝕃Vector₂ and 𝕄Vector₂ ($𝑷$). See Array of Matrices types.

Note that this function computes the following algebraic expression:

$\begin{pmatrix} B_1 & \hspace{0.01cm} & 0 \\ \hspace{0.01cm} & \ddots & \hspace{0.01cm} \\ 0 & \hspace{0.01cm} & B_m \end{pmatrix} \begin{pmatrix} C_{11} & \cdots & C_{1m} \\ \vdots & \ddots & \vdots \\ C_{m1} & \cdots & C_{mm} \end{pmatrix} \begin{pmatrix}B_1^T & \hspace{0.01cm} & 0 \\ \hspace{0.01cm} & \ddots & \hspace{0.01cm} \\ 0 & \hspace{0.01cm} & B_m^T\end{pmatrix}$ .

The result is a vector of vector of matrices of the matrixVector₂Type argument, which must be provided and must be one of the following abstract types: ℍVector₂, 𝔻Vector₂, 𝕃Vector₂ or 𝕄Vector₂ (and not an instance of these types).

When you pass it to this function, make sure to typecast $𝐁$ as an ℍVector, 𝔻Vector, 𝕃Vector or 𝕄Vector type if it is not already created as one of these types. See the example here below and typecasting matrices.

Method (2), (3) and (4) are multi-threaded. See Threads.

`matrixVector₂Type` in (3)-(4).
For example, in (1) if ``B`` and ``P`` are `Hermitian`,
calling `cong(B, P, 𝔻)` will actually
return the diagonal part of ``B*P*B'`` and calling `cong(B, P, 𝕃)` will
actually return its lower triangular part. The full congruence can
be obtained as an `Hermitian` matrix by `cong(B, P, ℍ)` and as a generic
matrix object by `cong(B, P, 𝕄)`.

Examples

using LinearAlgebra, PosDefManifold

# (1)
P=randP(3) # generate a 3x3 positive matrix
M=randn(3, 3)
C=cong(M, P, ℍ) # equivalent to C=ℍ(M*P*M')

# (2)
Pset=randP(4, 100); # generate 100 positive definite 4x4 matrices
M=randn(4, 4)
Qset=cong(M, Pset, ℍVector) # = [M*Pset_1*M',...,M*Pset_k*M'] as an ℍVector type

# recenter the matrices in Pset to their Fisher mean:
Qset=cong(invsqrt(mean(Fisher, Pset)), Pset, ℍVector)

# as a check, the Fisher mean of Qset is now the identity
mean(Fisher, Qset)≈I ? println("⭐") : println("⛔")

# (3)
Pset1=randP(4, 10); # generate 10 positive definite 4x4 matrices
Pset2=randP(4, 8);
Pset=ℍVector₂([Pset1, Pset2]);
M=randn(4, 4)
Qset=cong(M, Pset, MatrixVector₂)
Qset[1][1]≈M*Pset[1][1]*M' ? println("⭐") : println("⛔")
Qset[1][5]≈M*Pset[1][5]*M' ? println("⭐") : println("⛔")
Qset[2][1]≈M*Pset[2][1]*M' ? println("⭐") : println("⛔")
Qset[2][4]≈M*Pset[2][4]*M' ? println("⭐") : println("⛔")

# (4)
Pset1=randP(4, 2); # generate 2 positive definite 4x4 matrices
Pset2=randP(4, 2);
Pset=ℍVector₂([Pset1, Pset2]);
U=𝕄Vector([randU(4), randU(4)])
Qset=cong(U, Pset, MatrixVector₂)
Qset[1][1]≈U[1]*Pset[1][1]*U[1]' ? println("⭐") : println("⛔")
Qset[1][2]≈U[1]*Pset[1][2]*U[2]' ? println("⭐") : println("⛔")
Qset[2][1]≈U[2]*Pset[2][1]*U[1]' ? println("⭐") : println("⛔")
Qset[2][2]≈U[2]*Pset[2][2]*U[2]' ? println("⭐") : println("⛔")

evd(S::Union{𝕄{T}, ℍ{T}}) where T<:RealOrComplex

Given a positive semi-definite matrix $S$, returns a 2-tuple $(Λ, U)$, where $U$ is the matrix holding in columns the eigenvectors and $Λ$ is the matrix holding the eigenvalues on the diagonal. This is the output of Julia eigen function in $UΛU'=S$ form.

As for the eigen function, the eigenvalues and associated eigenvectors are sorted by increasing values of eigenvalues.

$S$ may be real or complex and may be flagged by Julia as Hermitian (in this case PosDefManifold assumes it is positive definite).

See typecasting matrices.

See also: spectralFunctions.

Examples

using PosDefManifold
A=randn(3, 3);
S=A+A';
Λ, U=evd(S); # which is equivalent to (Λ, U)=evd(P)
(U*Λ*U') ≈ S ? println(" ⭐ ") : println(" ⛔ ")
# => UΛU'=S, UΛ=SU, ΛU'=U'S

frf(P::ℍ{T}) where T<:RealOrComplex

Full-rank factorization of Hermitian matrix P. It is given by

$F=UD^{1/2}$,

where

$\textrm{EVD}(P)=UDU^{H}$

is the eigenvalue-eigenvector decomposition of P. It verifies

$FF^H=P$,

thus $F^{-1}$ is a whitening matrix.

See also: invfrf.

Examples

using LinearAlgebra, PosDefManifold
P=randP(3)
F = frf(P)
F*F'≈P ? println(" ⭐ ") : println(" ⛔ ")

invfrf(P::ℍ{T}) where T<:RealOrComplex

Inverse of the full-rank factorization of Hermitian matrix P. It is given by

$F=D^{-1/2}U^H$,

where

$\textrm{EVD}(P)=UDU^{H}$

is the eigenvalue-eigenvector decomposition of P. It verifies

$FPF^H=I$,

thus $F$ is a whitening matrix.

See also: frf.

Examples

using LinearAlgebra, PosDefManifold
P=randP(3)
F = invfrf(P)
F*P*F'≈I ? println(" ⭐ ") : println(" ⛔ ")

(1) spectralFunctions(P::ℍ{T}, func) where T<:RealOrComplex
(2) spectralFunctions(D::𝔻{S}, func) where S<:Real

(1) This is the mother function for all spectral functions of eigenvalues implemented in this library, which are:

• pow (power),
• isqrt (inverse square root).

The function sqr (square) does not use it, as it can be obtained more efficiently by simple multiplication.

You can use this function if you need another spectral function of eigenvalues besides those and those already implemented in the standard package LinearAlgebra. In general, you won't call it directly.

func is the function that will be applied on the eigenvalues.

$P$ must be flagged as Hermitian. See typecasting matrices. It must be a positive definite or positive semi-definite matrix, depending on func.

A special method is provided for real Diagonal matrices (2).

Maths

The definition of spectral functions for a positive definite matrix $P$ is at it follows:

$f\big(P\big)=Uf\big(Λ\big)U^H,$

where $U$ is the matrix holding in columns the eigenvectors of $P$, $Λ$ is the matrix holding on diagonal its eigenvalues and $f$ is a function applying element-wise to the eigenvalues.

See also: evd.

Examples

using LinearAlgebra, PosDefManifold
n=5
P=randP(n) # P=randP(ComplexF64, 5) to generate an Hermitian complex matrix
noise=0.1;
Q=spectralFunctions(P, x->x+noise) # add white noise to the eigenvalues
tr(Q)-tr(P) ≈ noise*n ? println(" ⭐ ") : println(" ⛔ ")

(1) pow(P::ℍ{T}, args...) where T<:RealOrComplex
(2) pow(D::𝔻{S}, args...) where S<:Real

(1) Given a positive semi-definite Hermitian matrix $P$, return the power $P^{r_1}, P^{r_2},...$ for any number of exponents $r_1, r_2,...$. It returns a tuple comprising as many elements as arguments passed after $P$.

$P$ must be flagged as Hermitian. See typecasting matrices.

$arg1, arg2,...$ are real numbers.

A special method is provided for real Diagonal matrices (2).

See also: invsqrt.

Examples

using LinearAlgebra, PosDefManifold
P=randP(5);     # use P=randP(ComplexF64, 5) for generating an Hermitian matrix
Q=pow(P, 0.5);            # =>  QQ=P
Q, W=pow(P, 0.5, -0.5);
W*P*W ≈ I ? println(" ⭐ ") : println(" ⛔ ")
Q*Q ≈ P ? println(" ⭐ ") : println(" ⛔ ")
R, S=pow(P, 0.3, 0.7);
R*S ≈ P ? println(" ⭐ ") : println(" ⛔ ")

(1) invsqrt(P{T}::ℍ) where T<:RealOrComplex
(2) invsqrt(D{S}::𝔻) where S<:Real

Given a positive definite Hermitian matrix $P$, compute the inverse of the principal square root $P^{-1/2}$.

$P$ must be flagged as Hermitian. See typecasting matrices.

A special method is provided for real Diagonal matrices (2).

Maths

The principal square root of a positive definite matrix $P$ is the only symmetric (if $P$ is real) or Hermitian (if $P$ is complex) square root. Its inverse $P^{-1/2}$ is also named the whitening or sphering matrix since$P^{-1/2}PP^{-1/2}=I$.

See: typecasting matrices.

See also: pow.

Examples

using LinearAlgebra, PosDefManifold
P=randP(ComplexF64, 5);
Q=invsqrt(P);
Q*P*Q ≈ I ? println(" ⭐ ") : println(" ⛔ ")

(1) sqr(P::ℍ{T}) where T<:RealOrComplex
(2) sqr(X::Union{𝕄{T}, 𝕃{T}, 𝔻{S}}) where T<:RealOrComplex where S<:Real

(1) Given a positive semi-definite Hermitian matrix $P$, compute its square $P^{2}$.

$P$ must be flagged as Hermitian. See typecasting matrices.

A method is provided also for generic matrices of the Matrix type, LowerTriangular matrices and real Diagonal matrices (2). The output is of the same type as the input.

See also: pow.

Examples

using PosDefManifold
P=randP(5);
P²=sqr(P);  # =>  P²=PP
sqrt(P²)≈ P ? println(" ⭐ ") : println(" ⛔ ")

powerIterations(H::Union{ℍ{T}, 𝕄{T}}, q::Int;
<
evalues=false,
tol::Real=0,
maxiter::Int=300,
verbose=false>) where T<:RealOrComplex

powerIterations(L::𝕃{S}, q::Int;
< same optional keyword arguments in (1)>) where S<:Real

alias: powIter

(1) Compute the $q$ eigenvectors associated to the $q$ largest (real) eigenvalues of real or complex Hermitian or Matrix $H$ using the power iterations + Gram-Schmidt orthogonalization as suggested by Strang. The eigenvectors are returned with the same type as the elements of $H$.

$H$ must have real eigenvalues, that is, it must be a symmetric matrix if it is real or an Hermitian matrix if it is complex.

(2) as in (1), but using only the LowerTriangular view $L$ of a matrix. This option is available only for real matrices (see below).

The following are <optional keyword arguments>:

• `tol is the tolerance for the convergence of the power method (see below),
• `maxiter is the maximum number of iterations allowed for the power method,
• if `verbose=true, the convergence of all iterations will be printed,
• if evalues=true, return the 4-tuple(Λ, U, iterations, covergence)`,
• if evalues=false return the 3-tuple(U, iterations, covergence)`.

See also: mgs.

Examples

using LinearAlgebra, PosDefManifold
# Generate an Hermitian (complex) matrix
H=randP(ComplexF64, 10);
# 3 eigenvectors and eigenvalues
U, iterations, convergence=powIter(H, 3, verbose=true)
# all eigenvectors
Λ, U, iterations, convergence=powIter(H, size(H, 2), evalues=true, verbose=true);
U'*U ≈ I && U*Λ*U'≈H ? println(" ⭐ ") : println(" ⛔ ")

# passing a `Matrix` object
Λ, U, iterations, convergence=powIter(Matrix(H), 3, evalues=true)

# passing a `LowerTriangular` object (must be a real matrix in this case)
L=𝕃(randP(10))
Λ, U, iterations, convergence=powIter(L, 3, evalues=true)

(1) choL(P::ℍ{T}) where T<:RealOrComplex
(2) choL(D::𝔻{S}) where S<:Real

(1) Given a real or complex positive definite Hermitian matrix $P$, return the Cholesky lower triangular factor $L$ such that $LL^H=P$. To obtain $L^H$ or both $L$ and $L^H$, use instead julia function cholesky.

On output, $L$ is of type LowerTriangular.

(2) For a real Diagonal matrix $D$, return $D^{1/2}$.

See also: choInv.

Examples

using PosDefManifold
P=randP(5);
L=choL(P);
L*L'≈ P ? println(" ⭐ ") : println(" ⛔ ")

choInv(P::AbstractArray{T};
	kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex

For a real or complex positive definite matrix $P$, let $P=LL^H$ be its Cholesky decomposition and $P=L_1DL_1^H$ the related LDLt decomposition. In the above, $L$ is a lower triangular matrix, $D$ a positive-definite diagonal matrix and $L_1$ a unit lower triangular matrix. Return:

• if kindis :LLt (default), the 2-tuple $L$, $L^{-H}$
• if kindis :LDLt, the 3-tuple $L_1$, $D$, $L_1^{-H}$.

Those are obtained in one pass and for small matrices this is faster then calling Julia's chelosky function and inverting the lower factor unless you set

 BLAS.set_num_threads(1).

Input matrix P may be of type Matrix or Hermitian. Since only the lower triangle is used, P may also be a LowerTriangular view of a positive definite matrix. If P is real, it can also be of the Symmetric type.

The algorithm is a multiplicative Gaussian elimination. If run completely, in input matrix P there will be the Identity at the end.

Notes: Output $L^{-H}$ is an inverse square root (whitening matrix) of $P$, since $L^{-1}PL^{-H}=I$. It therefore yields the inversion of $P$ as $P^{-1}=L^{-H}L^{-1}$. It is the fastest whitening matrix to be computed, however it yields poor numerical precision, especially for large matrices.

The following relations holds:

• $L=PL^{-H}$
• $L^{H}=L^{-1}P$
• $L^{-H}=P^{-1}L$
• $L^{-1}=L^{H}P^{-1}$.

We also have

• $L^{H}L=L^{-1}P^{2}L^{-H}=UPU^H$, with $U$ orthogonal (see below) and
• $L^{-1}L^{-H}=L^{H}P^{-2}L=UP^{-1}U^H$.

$LL^{H}$ and $L^{H}L$ are unitarily similar, that is,

$ULL^{H}U^H=L^{H}L$,

where $U=L^{-1}P^{1/2}$, with $P^{1/2}=H$ the principal (unique symmetric) square root of $P$. This is seen writing $PP^{-1}=HHL^{-H}L^{-1}$; multiplying both sides on the left by $L^{-1}$ and on the right by $L$ we obtain

$L^{-1}PP^{-1}L=L^{-1}HHL^{-H}=I=(L^{-1}H)(L^{-1}H)^H$

and since $L^{-1}H$ is square it must be unitary.

From these expressions we have

• $H=LU=U^HL^H$
• $L=HU^H$
• $H^{-1}=U^HL^{-1}$
• $L^{-1}=UH^{-1}$.

$U$ is the polar factor of $L^{H}$, i.e., $L^{H}=UH$, since $LL^{H}=HU^HUH^H=H^2=P$.

From $L^{H}L=UCU^H$ we have $L^{H}LU=UC=ULL^{H}$ and from $U=L^{-1}H$ we have $L=HU^H$.

See also: choInv!, choL.

Examples

using PosDefManifold
n, t = 800, 6000
etol = 1e-9
Z=randn(t, n)
Y=Z'*Z
Yi=inv(Y)

A, B=choInv!(copy(Y))
norm(A*A'-Y)/√n < etol ? println(" ⭐ ") : println(" ⛔ ")
norm(B*B'-Yi)/√n < etol ? println(" ⭐ ") : println(" ⛔ ")

A, D, B=choInv!(copy(Y); kind=:LDLt)
norm(Y-A*D*A')/√n < etol ? println(" ⭐ ") : println(" ⛔ ")
norm(Yi-B*inv(D)*B')/√n < etol ? println(" ⭐ ") : println(" ⛔ ")

# repeat the test for complex matrices
Z=randn(ComplexF64, t, n)
Y=Z'*Z
Yi=inv(Y)

A, B=choInv!(copy(Y))
norm(A*A'-Y)/√n < etol ? println(" ⭐ ") : println(" ⛔ ")
norm(B*B'-Yi)/√n < etol ? println(" ⭐ ") : println(" ⛔ ")

A, D, B=choInv!(copy(Y); kind=:LDLt)
norm(Y-A*D*A')/√n < etol ? println(" ⭐ ") : println(" ⛔ ")
norm(Yi-B*inv(D)*B')/√n < etol ? println(" ⭐ ") : println(" ⛔ ")

choInv!(P::AbstractArray{T};
	kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex

The same thing as choInv, but destroys the input matrix. This function does nt require copying the input matrix, thus it is slightly faster.