From 18e7c953f4e1efbf008a40495eee454c650713c5 Mon Sep 17 00:00:00 2001
From: Marco-Congedo <marco.congedo@gmail.com>
Date: Sat, 28 Dec 2019 09:53:05 +0100
Subject: [PATCH] implemented `frf` and `invfrf` (whitening) functions

---
 docs/src/linearAlgebra.md |  4 +++
 src/PosDefManifold.jl     |  2 ++
 src/linearAlgebra.jl      | 73 +++++++++++++++++++++++++++++++++++++++
 src/test.jl               | 12 +++++++
 4 files changed, 91 insertions(+)

diff --git a/docs/src/linearAlgebra.md b/docs/src/linearAlgebra.md
index 1c6aafc..25bc865 100644
--- a/docs/src/linearAlgebra.md
+++ b/docs/src/linearAlgebra.md
@@ -141,6 +141,8 @@ congruence
 | Function   | Description |
 |:---------- |:----------- |
 | [`evd`](@ref) | Eigenvalue-Eigenvector decomposition of a matrix in ``UΛU'=P`` form|
+| [`frf`](@ref) | Full-rank factorization of an Hermitian matrix |
+| [`invfrf`](@ref) | Inverse of the full-rank factorization of an Hermitian matrix (whitening) |
 | [`spectralFunctions`](@ref) | Mother function for creating spectral functions of eigenvalues|
 | [`pow`](@ref)| Power of a positive matrix for any number of exponents in one pass|
 | [`invsqrt`](@ref)| Principal square root inverse (whitening) of a positive matrix|
@@ -151,6 +153,8 @@ congruence
 
 ```@docs
 evd
+frf
+invfrf
 spectralFunctions
 pow
 invsqrt
diff --git a/src/PosDefManifold.jl b/src/PosDefManifold.jl
index 93ca0f7..b21558a 100644
--- a/src/PosDefManifold.jl
+++ b/src/PosDefManifold.jl
@@ -137,6 +137,8 @@ export
     fVec,
     congruence, cong,
     evd,
+    frf,
+    invfrf,
     spectralFunctions,
     ispos,
     pow,
diff --git a/src/linearAlgebra.jl b/src/linearAlgebra.jl
index 73459a3..e21c42b 100644
--- a/src/linearAlgebra.jl
+++ b/src/linearAlgebra.jl
@@ -1451,6 +1451,79 @@ function evd(S::Union{𝕄{T}, ℍ{T}}) where T<:RealOrComplex # return tuple (
 end
 
 
+
+"""
+```
+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`](@ref).
+
+## Examples
+```
+using LinearAlgebra, PosDefManifold
+P=randP(3)
+F = frf(P)
+F*F'≈P ? println(" ⭐ ") : println(" ⛔ ")
+```
+"""
+function frf(P::ℍ{T}) where T<:RealOrComplex
+   #size(P, 1)≠size(A, 2) && throw(ArgumentError(📌*", frf function: input matrix must be square"))
+   λ, U=eigen(P)
+   return U*Diagonal(sqrt.(λ))
+end
+
+
+"""
+```
+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`](@ref).
+
+## Examples
+```
+using LinearAlgebra, PosDefManifold
+P=randP(3)
+F = invfrf(P)
+F*P*F'≈I ? println(" ⭐ ") : println(" ⛔ ")
+```
+"""
+function invfrf(P::ℍ{T}) where T<:RealOrComplex
+   #size(P, 1)≠size(A, 2) && throw(ArgumentError(📌*", frf function: input matrix must be square"))
+   Λ, U=evd(P)
+   return invsqrt(Λ)*U'
+end
+
+
 """
     (1) spectralFunctions(P::ℍ{T}, func) where T<:RealOrComplex
     (2) spectralFunctions(D::𝔻{S}, func) where S<:Real
diff --git a/src/test.jl b/src/test.jl
index f7ae893..a8c6892 100644
--- a/src/test.jl
+++ b/src/test.jl
@@ -280,6 +280,18 @@ function tests();
     U*Λ*U'≈PC ? OK() : OH(name*" Complex Input")
 
 
+    name="frf"; newTest(name)
+    F = frf(P)
+    F*F'≈P ? OK() : OH(name*" Real Input")
+    F = frf(PC)
+    F*F'≈PC ? OK() : OH(name*" Complex Input")
+
+    name="invfrf"; newTest(name)
+    F = invfrf(P)
+    F*P*F'≈I ? OK() : OH(name*" Real Input")
+    F = invfrf(PC)
+    F*PC*F'≈I ? OK() : OH(name*" Complex Input")
+
     name="spectralFunctions"; newTest(name);
     spectralFunctions(P, x->x+1); RUN()
     spectralFunctions(PC, abs2); RUN()