AbstractFloat support (#84)

* Float32 support for linear shrinkages

* Float32 support for nonlinear shrinkage

* fix unit tests

* cast matrix types to float as suggested by @mateuszbaran

* bump package version

Author: biona001

Project.toml

@@ -1,7 +1,7 @@
 name = "CovarianceEstimation"
 uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
 authors = ["Mateusz Baran <[email protected]>", "Thibaut Lienart"]
-version = "0.2.7"
+version = "0.2.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/linearshrinkage.jl

@@ -255,8 +255,9 @@ function linear_shrinkage(::DiagonalUnitVariance, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F   = I
+    T   = float(eltype(S))
     κ   = n - Int(corrected)
-    γ   = κ/n
+    γ   = T(κ/n)
     Xc² = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
@@ -265,15 +266,15 @@ function linear_shrinkage(::DiagonalUnitVariance, Xc::AbstractMatrix,
         λ  /= κ * (ΣS² - 2tr(S) + p)
     elseif λ == :ss
         # use the standardised data matrix
-        d   = 1.0 ./ vec(sum(Xc², dims=1))
+        d   = one(T) ./ vec(sum(Xc², dims=1))
         S̄   = rescale(S, sqrt.(d)) # this has diagonal 1/κ
         ΣS̄² = sumij2(S̄, with_diag=true)
         λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
-        λ  /= κ * ΣS̄² - p / κ
+        λ  /= T(κ * ΣS̄² - p / κ)
     else
         throw(ArgumentError("Unsupported shrinkage method for target DiagonalUnitVariance: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink(F, S, λ)
 end
 
@@ -298,8 +299,9 @@ function linear_shrinkage(::DiagonalCommonVariance, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F   = target_B(S, p)
+    T   = float(eltype(F))
     κ   = n - Int(corrected)
-    γ   = κ/n
+    γ   = T(κ/n)
     Xc² = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
@@ -309,28 +311,28 @@ function linear_shrinkage(::DiagonalCommonVariance, Xc::AbstractMatrix,
         λ  /= κ * (ΣS² - p*v^2)
     elseif λ == :ss
         # use the standardised data matrix
-        d   = 1.0 ./ vec(sum(Xc², dims=1))
+        d   = one(T) ./ vec(sum(Xc², dims=1))
         S̄   = rescale(S, sqrt.(d)) # this has diagonal 1/κ
         v̄   = κ # tr(S̄)/p
         ΣS̄² = sumij2(S̄, with_diag=true)
         λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
-        λ  /= κ * ΣS̄² - p/κ
+        λ  /= T(κ * ΣS̄² - p/κ)
     elseif λ == :rblw
         # https://arxiv.org/pdf/0907.4698.pdf equations 17, 19
         trS² = sum(abs2, S)
         tr²S = tr(S)^2
         # note: using corrected or uncorrected S does not change λ
-        λ = ((n-2)/n * trS² + tr²S) / ((n+2) * (trS² - tr²S/p))
+        λ = T(((n-2)/n * trS² + tr²S) / ((n+2) * (trS² - tr²S/p)))
     elseif λ == :oas
         # https://arxiv.org/pdf/0907.4698.pdf equation 23
         trS² = sum(abs2, S)
         tr²S = tr(S)^2
         # note: using corrected or uncorrected S does not change λ
-        λ = ((1.0-2.0/p) * trS² + tr²S) / ((n+1.0-2.0/p) * (trS² - tr²S/p))
+        λ = ((one(T)-T(2.0)/p) * trS² + tr²S) / ((n+one(T)-T(2.0)/p) * (trS² - tr²S/p))
     else
         throw(ArgumentError("Unsupported shrinkage method for target DiagonalCommonVariance: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink(F, S, λ)
 end
 
@@ -355,8 +357,9 @@ function linear_shrinkage(::DiagonalUnequalVariance, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F   = target_D(S)
+    T   = float(eltype(F))
     κ   = n - Int(corrected)
-    γ   = κ/n
+    γ   = T(κ / n)
     Xc² = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
@@ -365,14 +368,14 @@ function linear_shrinkage(::DiagonalUnequalVariance, Xc::AbstractMatrix,
         λ  /= κ * ΣS²
     elseif λ == :ss
         # use the standardised data matrix
-        d   = 1.0 ./ vec(sum(Xc², dims=1))
+        d   = one(T) ./ vec(sum(Xc², dims=1))
         ΣS̄² = sumij2(rescale(S, sqrt.(d)))
         λ   = sumij(rescale!(uccov(Xc²), d)) / γ^2 - ΣS̄²
         λ  /= κ * ΣS̄²
     else
         throw(ArgumentError("Unsupported shrinkage method for target DiagonalUnequalVariance: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink(F, S, λ)
 end
 
@@ -405,24 +408,25 @@ function linear_shrinkage(::CommonCovariance, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F, v, c = target_C(S, p)
+    T   = float(eltype(F))
     κ   = n - Int(corrected)
-    γ   = κ/n
+    γ   = T(κ/n)
     Xc² = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
         ΣS² = sumij2(S, with_diag=true)
         λ   = sumij(uccov(Xc²), with_diag=true) / γ^2 - ΣS²
         λ  /= κ * (ΣS² - p*(p-1)*c^2 - p*v^2)
     elseif λ == :ss
-        d   = 1.0 ./ vec(sum(Xc², dims=1))
+        d   = one(T) ./ vec(sum(Xc², dims=1))
         S̄   = rescale(S, sqrt.(d))
         ΣS̄² = sumij2(S̄, with_diag=true)
         λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
         λ  /= κ * ΣS̄² - p/κ - κ * sumij(S̄; with_diag=false)^2 / (p * (p - 1))
     else
         throw(ArgumentError("Unsupported shrinkage method for target CommonCovariance: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink!(F, S, λ)
 end
 
@@ -451,8 +455,9 @@ function linear_shrinkage(::PerfectPositiveCorrelation, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F   = target_E(S)
+    T   = float(eltype(F))
     κ   = n - Int(corrected)
-    γ   = κ/n
+    γ   = T(κ/n)
     Xc² = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
@@ -461,7 +466,7 @@ function linear_shrinkage(::PerfectPositiveCorrelation, Xc::AbstractMatrix,
         λ  -= sum_fij(Xc, S, n, κ)
         λ  /= sumij2(S - F)
     elseif λ == :ss
-        d   = 1.0 ./ vec(sum(Xc², dims=1))
+        d   = one(T) ./ vec(sum(Xc², dims=1))
         s   = sqrt.(d)
         S̄   = rescale(S, s)
         ΣS̄² = sumij2(S̄)
@@ -472,7 +477,7 @@ function linear_shrinkage(::PerfectPositiveCorrelation, Xc::AbstractMatrix,
     else
         throw(ArgumentError("Unsupported shrinkage method for target PerfectPositiveCorrelation: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink!(F, S, λ)
 end
 
@@ -505,8 +510,9 @@ function linear_shrinkage(::ConstantCorrelation, Xc::AbstractMatrix,
                           corrected::Bool)
 
     F, r̄ = target_F(S, p)
+    T    = float(eltype(F))
     κ    = n - Int(corrected)
-    γ    = κ/n
+    γ    = T(κ/n)
     Xc²  = Xc.^2
     # computing the shrinkage
     if λ ∈ [:auto, :lw]
@@ -515,7 +521,7 @@ function linear_shrinkage(::ConstantCorrelation, Xc::AbstractMatrix,
         λ  -= r̄ * sum_fij(Xc, S, n, κ)
         λ  /= sumij2(S - F)
     elseif λ == :ss
-        d    = 1.0 ./ vec(sum(Xc², dims=1))
+        d    = one(T) ./ vec(sum(Xc², dims=1))
         s    = sqrt.(d)
         S̄    = rescale(S, s)
         F̄, r̄ = target_F(S̄, p)
@@ -526,6 +532,6 @@ function linear_shrinkage(::ConstantCorrelation, Xc::AbstractMatrix,
     else
         throw(ArgumentError("Unsupported shrinkage method for target ConstantCorrelation: $λ."))
     end
-    λ = clamp(λ, 0.0, 1.0)
+    λ = clamp(λ, zero(T), one(T))
     return linshrink!(F, S, λ)
 end

src/nonlinearshrinkage.jl

@@ -27,24 +27,24 @@ const EPAN_3 = 0.3 * INVPI
 
 Return the Epanechnikov kernel evaluated at `x`.
 """
-epanechnikov(x::Real) = EPAN_1 * max(0.0, 1.0 - x^2/5.0)
+epanechnikov(x::T) where T<:Real = float(T)(EPAN_1 * max(0.0, 1.0 - x^2/5.0))
 
 """
     epnanechnikov_HT(x)
 
 Return the Hilbert Transform of the Epanechnikov kernel evaluated at `x`
 if `|x|≂̸√5`.
 """
-function epanechnikov_HT1(x::Real)
-    -EPAN_3 * x + EPAN_2 * (1.0 - x^2/5.0) * log(abs((SQRT5 - x)/(SQRT5 + x)))
+function epanechnikov_HT1(x::T) where T <: Real
+    float(T)(-EPAN_3 * x + EPAN_2 * (1.0 - x^2/5.0) * log(abs((SQRT5 - x)/(SQRT5 + x))))
 end
 
 """
     epnanechnikov_HT2(x)
 Return the Hilbert Transform of the Epanechnikov kernel evaluated at `x`
 if `|x|=√5`.
 """
-epanechnikov_HT2(x::Real) = -EPAN_3*x
+epanechnikov_HT2(x::T) where T <: Real = float(T)(-EPAN_3*x)
 
 """
     analytical_nonlinear_shrinkage(S, n, p; decomp)
@@ -66,19 +66,20 @@ function analytical_nonlinear_shrinkage(S::AbstractMatrix{<:Real},
     sample_perm = @view perm[max(1, (p - η) + 1):p]
     λ    = @view F.values[sample_perm]
     U    = F.vectors[:, perm]
+    T    = float(eltype(F))
 
     # dominant cost forming of S or eigen(S) --> O(max{np^2, p^3})
 
     # compute analytical nonlinear shrinkage kernel formula
     L = repeat(λ, outer=(1, min(p, η)))
 
     # Equation (4.9)
-    h = η^(-1//3)
+    h = T(η^(-1//3))
     H = h * L'
     x = (L .- L') ./ H
 
     # additional useful definitions
-    γ  = p/η
+    γ  = T(p/η)
     πλ = π * λ
 
     # Equation (4.7) in http://www.econ.uzh.ch/static/wp/econwp264.pdf
@@ -96,17 +97,17 @@ function analytical_nonlinear_shrinkage(S::AbstractMatrix{<:Real},
     if p <= η
         # Equation (4.3)
         πγλ = γ * πλ
-        denom = @. (πγλ * f̃)^2 + (1.0 - γ - πγλ * Hf̃)^2
+        denom = @. (πγλ * f̃)^2 + (one(T) - γ - πγλ * Hf̃)^2
         d̃ = λ ./ denom
     else
         # Equation (C.8)
-        hs5  = h * SQRT5
-        Hf̃0  = (0.3/h^2 + 0.75/hs5 * (1.0 - 0.2/h^2) * log((1+hs5)/(1-hs5)))
-        Hf̃0 *= mean(1.0 ./ πλ)
+        hs5  = T(h * SQRT5)
+        Hf̃0  = T((0.3/h^2 + 0.75/hs5 * (1.0 - 0.2/h^2) * log((1+hs5)/(1-hs5))))
+        Hf̃0 *= mean(one(T) ./ πλ)
         # Equation (C.5)
-        d̃0 = INVPI / ((γ - 1.0) * Hf̃0)
+        d̃0 = T(INVPI / ((γ - one(T)) * Hf̃0))
         # Eq. (C.4)
-        d̃1 = @. 1.0 / (π * πλ * (f̃^2 + Hf̃^2))
+        d̃1 = @. one(T) / (π * πλ * (f̃^2 + Hf̃^2))
         d̃  = [fill(d̃0, (p - η, 1)); d̃1]
     end
 

src/utils.jl

@@ -1,12 +1,12 @@
 function linshrink(F::AbstractMatrix, S::AbstractMatrix, λ::Real)
-    return Symmetric((1.0 .- λ).*S .+ λ.*F)
+    return Symmetric((one(λ) .- λ).*S .+ λ.*F)
 end
 
 function linshrink(F::UniformScaling, S::AbstractMatrix, λ::Real)
-    return Symmetric((1.0 .- λ).*S + λ.*F)
+    return Symmetric((one(λ) .- λ).*S + λ.*F)
 end
 
 function linshrink!(F::AbstractMatrix, S::AbstractMatrix, λ::Real)
-    F .= (1.0 .- λ).*S .+ λ.*F
+    F .= (one(λ) .- λ).*S .+ λ.*F
     return Symmetric(F)
 end

test/runtests.jl

@@ -75,4 +75,5 @@ end
     include("test_biweight.jl")
     include("test_linearshrinkage.jl")
     include("test_nonlinearshrinkage.jl")
+    include("test_float32.jl")
 end

test/test_float32.jl

@@ -0,0 +1,26 @@
+@testset "Float32 matrices" begin
+    # linear shrinkages
+    for X in test_matrices
+        x = convert(Matrix{Float32}, X)
+        for target in [
+            DiagonalUnitVariance(),
+            DiagonalCommonVariance(),
+            DiagonalUnequalVariance(),
+            CommonCovariance(),
+            PerfectPositiveCorrelation(),
+            ConstantCorrelation()
+            ]
+            for shrinkage in [:lw, :ss]
+                @test eltype(cov(LinearShrinkage(target, shrinkage), x)) == Float32
+            end
+        end
+    end
+
+    # nonlinear shrinkages
+    ANS = AnalyticalNonlinearShrinkage()
+    for X in test_matrices
+        size(X, 1) < 12 && continue
+        x = convert(Matrix{Float32}, X)
+        @test eltype(cov(ANS, x)) == Float32
+    end
+end