Tools

function eig(A)

function eig(A, B)

Call Julia function eigen and return its output sorted by descending order of eigenvalues.

function nonDiagonality(C::Union{Matrix, Diagonal, SorH})

Measure of deviancy from diagonality of $n⋅n$ square matrix C, defined as (Congedo et al., 2008)🎓.

$\frac{\sum_{i≠j}|c_{ij}|^2}{(n-1)\sum_{i}|c_{ii}|^2}$

It is equal to $0$ if $C$ is diagonal, equal to $1$ if $C$ is perfectly uniform.

Examples:

using Diagonalizations
C=ones(10, 10)                   # uniform matrix
nd=nonDiagonality(C)             # must be 1
D=Diagonal(abs.(randn(10, 10)))  # diagonal matrix
nd=nonDiagonality(D)             # must be 0

function spForm(P::Union{Mat, Real, Complex})

Measure of deviancy from scaled permutation form of $n⋅n$ square matrix P, defined as

$\frac{1}{2(n-1)}\bigg(\sum_{row}1-\beta(row)+\sum_{col}1-\beta(col)\bigg)$,

where for each row and column of P, β is the maximum of the absolute values divided by the sum of the absolute values.

This index is equal to $0$ if in each row and column $P$ has only one non-zero element, that is, if $P$ is a scaled permutation matrix. The larger the index, the farther away $P$ is from this form.

This measure and several existing variants are well-known in the blind source separation / independent component analysis community, where it is used to compare approximate joint diagonalization algorithms on simulated data. In fact, if $A$ is the inverse of the approximate joint diagonalizer that is used to generate the data and $B$ the approximate joint diagonalizer estimated by an algorithm, $P=BA$ must be as close as possible to a scaled permutation matrix (see scale and permutation).

Return 0.0 (zero) if P is a real of complex number.

Examples:

using Diagonalizations, PosDefManifold
# create 20 random commuting matrices
# they all have the same eigenvectors
Cset=randP(3, 20; eigvalsSNR=Inf, commuting=true)
# estimate the approximate joint diagonalizer (ajd)
a=ajd(Cset)
# the ajd must be equivalent to the eigenvector matrix of
# any of the matrices in Cset
spForm(a.F'*eigvecs(Cset[1]))+1.0≈1.0 ? println(" ⭐ ") : println(" ⛔ ")

(1)
function genDataMatrix(t::Int, n::Int, A=nothing)

(2)
function genDataMatrix(::Type{Complex{T}},
                       t::Int, n::Int, A=○) where {T<:AbstractFloat}

(1) Generate a $t⋅n$ random data matrix as $XA$, where $X$ is a $t⋅n$ matrix with entries randomly drawn from a Gaussian distribution and $A$ a $n⋅n$ symmetric matrix, which, if not provided as argument A, will be generated with entries randomly drawn from a uniform distribution ∈[-1, 1].

(2) as (1), but $X$ is generated randomly from a complex Gaussian distribution and $A$is Hermitian (complex) which, if not provided as argument A, will be generated with entries randomly drawn from a uniform distribution ∈[-1-i1, 1+i1].

Examples:

A=genDataMatrix(100, 20) # (1), real
A=genDataMatrix(ComplexF64, 100, 20) # (2), complex