Fix cov2cor and cor2cov with Hermitian and Symmetric matrices (#953)

* Fix `cov2cor` and `cor2cov` with `Hermitian` and `Symmetric` matrices

* Fix types and write only to active triangle

* Fix tests and docs

* Revert workaround for Julia 1.0

Author: devmotion

Project.toml

@@ -1,7 +1,7 @@
 name = "StatsBase"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 authors = ["JuliaStats"]
-version = "0.34.4"
+version = "0.34.5"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"

docs/src/cov.md

@@ -13,7 +13,9 @@ std(::CovarianceEstimator, ::AbstractVector)
 cor
 mean_and_cov
 cov2cor
+StatsBase.cov2cor!
 cor2cov
+StatsBase.cor2cov!
 CovarianceEstimator
 SimpleCovariance
 ```

src/StatsBase.jl

@@ -18,7 +18,7 @@ using SparseArrays
 import Random: rand, rand!
 import LinearAlgebra: BlasReal, BlasFloat
 import Statistics: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
-                   cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
+                   cor, corm, unscaled_covzm, quantile, sqrt!,
                    median, middle
 using StatsAPI: StatisticalModel, RegressionModel
 import StatsAPI: pairwise, pairwise!, params, params!,

src/cov.jl

@@ -128,41 +128,124 @@ end
     cov2cor(C::AbstractMatrix, [s::AbstractArray])
 
 Compute the correlation matrix from the covariance matrix `C` and, optionally, a vector
-of standard deviations `s`. Use `StatsBase.cov2cor!` for an in-place version.
+of standard deviations `s`. Use [`StatsBase.cov2cor!`](@ref) for an in-place version.
 """
-cov2cor(C::AbstractMatrix, s::AbstractArray=sqrt.(view(C, diagind(C)))) = cov2cor!(copy(C), s)
+function cov2cor(C::AbstractMatrix, s::AbstractArray = map(sqrt, view(C, diagind(C))))
+    zs = zero(eltype(s))
+    T = typeof(zero(eltype(C)) / (zs * zs))
+    return cov2cor!(copyto!(similar(C, T), C), s)
+end
+
+# Original implementation: https://github.com/JuliaStats/Statistics.jl/blob/22dee82f9824d6045e87aa4b97e1d64fe6f01d8d/src/Statistics.jl#L633-L657
+"""
+    cov2cor!(C::AbstractMatrix, [s::AbstractArray])
+
+Convert the covariance matrix `C` to a correlation matrix in-place, optionally using a vector of
+standard deviations `s`.
+"""
+function cov2cor!(C::AbstractMatrix, s::AbstractArray = map(sqrt, view(C, diagind(C))))
+    Base.require_one_based_indexing(C, s)
+    n = length(s)
+    size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
+    for j = 1:n
+        sj = s[j]
+        for i = 1:(j-1)
+            C[i,j] = adjoint(C[j,i])
+        end
+        C[j,j] = oneunit(C[j,j])
+        for i = (j+1):n
+            C[i,j] = _clampcor(C[i,j] / (s[i] * sj))
+        end
+    end
+    return C
+end
+_clampcor(x::Real) = clamp(x, -1, 1)
+_clampcor(x) = x
+
+# Preserve structure of Symmetric and Hermitian covariance matrices
+function cov2cor!(C::Union{Symmetric{<:Real},Hermitian}, s::AbstractArray)
+    n = length(s)
+    size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
+    A = parent(C)
+    if C.uplo === 'U'
+        for j = 1:n
+            sj = s[j]
+            for i = 1:(j-1)
+                A[i,j] = _clampcor(A[i,j] / (s[i] * sj))
+            end
+            A[j,j] = oneunit(A[j,j])
+        end
+    else
+        for j = 1:n
+            sj = s[j]
+            A[j,j] = oneunit(A[j,j])
+            for i = (j+1):n
+                A[i,j] = _clampcor(A[i,j] / (s[i] * sj))
+            end
+        end
+    end
+    return C
+end
 
 """
     cor2cov(C, s)
 
 Compute the covariance matrix from the correlation matrix `C` and a vector of standard
-deviations `s`. Use `StatsBase.cor2cov!` for an in-place version.
+deviations `s`. Use [`StatsBase.cor2cov!`](@ref) for an in-place version.
 """
-cor2cov(C::AbstractMatrix, s::AbstractArray) = cor2cov!(copy(C), s)
+function cor2cov(C::AbstractMatrix, s::AbstractArray)
+    zs = zero(eltype(s))
+    T = typeof(zero(eltype(C)) * (zs * zs))
+    return cor2cov!(copyto!(similar(C, T), C), s)
+end
 
 """
     cor2cov!(C, s)
 
-Converts the correlation matrix `C` to a covariance matrix in-place using a vector of
+Convert the correlation matrix `C` to a covariance matrix in-place using a vector of
 standard deviations `s`.
 """
 function cor2cov!(C::AbstractMatrix, s::AbstractArray)
+    Base.require_one_based_indexing(C, s)
     n = length(s)
     size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
-    for i in CartesianIndices(size(C))
-        @inbounds C[i] *= s[i[1]] * s[i[2]]
+    for j in 1:n
+        sj = s[j]
+        for i in 1:(j-1)
+            C[i,j] = adjoint(C[j,i])
+        end
+        C[j,j] = sj^2
+        for i in (j+1):n
+            C[i,j] *= s[i] * sj
+        end
     end
     return C
 end
 
-"""
-    _cov2cor!(C)
-
-Compute the correlation matrix from the covariance matrix `C`, in-place.
-
-The leading diagonal is used to determine the standard deviations by which to normalise.
-"""
-_cov2cor!(C::AbstractMatrix) = cov2cor!(C, sqrt.(diag(C)))
+# Preserve structure of Symmetric and Hermitian correlation matrices
+function cor2cov!(C::Union{Symmetric{<:Real},Hermitian}, s::AbstractArray)
+    n = length(s)
+    size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
+    A = parent(C)
+    if C.uplo === 'U'
+        for j in 1:n
+            sj = s[j]
+            for i in 1:(j-1)
+                A[i,j] *= s[i] * sj
+            end
+            A[j,j] = sj^2
+        end
+    else
+        for j in 1:n
+            sj = s[j]
+            A[j,j] = sj^2
+            for i in (j+1):n
+                A[i,j] *= s[i] * sj
+            end
+        end
+    end
+    return C
+end
 
 """
     CovarianceEstimator
@@ -255,10 +338,10 @@ The keyword argument `mean` can be:
     of size `(N,1)`.
 """
 function cor(ce::CovarianceEstimator, X::AbstractMatrix; kwargs...)
-    return _cov2cor!(cov(ce, X; kwargs...))
+    return cov2cor!(cov(ce, X; kwargs...))
 end
 function cor(ce::CovarianceEstimator, X::AbstractMatrix, w::AbstractWeights; kwargs...)
-    return _cov2cor!(cov(ce, X, w; kwargs...))
+    return cov2cor!(cov(ce, X, w; kwargs...))
 end
 
 """

test/cov.jl

@@ -6,6 +6,27 @@ struct EmptyCovarianceEstimator <: CovarianceEstimator end
 @testset "StatsBase.Covariance" begin
 weight_funcs = (weights, aweights, fweights, pweights)
 
+function test_isapprox_preserves_symherm_structure(f::F, x::AbstractMatrix, y::AbstractMatrix, args...) where F
+    for wrapper in (identity, x -> Symmetric(x, :U), x -> Symmetric(x, :L), x -> Hermitian(x, :U), x -> Hermitian(x, :L))
+        A = wrapper(copy(x))
+        fA = @inferred(f(A, args...))
+        @test fA ≈ y
+        if f === StatsBase.cov2cor! || f === StatsBase.cor2cov!
+            @test fA === A
+            if A isa Union{Symmetric,Hermitian}
+                @test parent(fA) != fA # only active triangle is written to
+            end
+        else
+            @test fA !== A
+            if A isa Union{Symmetric,Hermitian}
+                @test fA isa (A isa Symmetric ? Symmetric : Hermitian)
+                @test fA.uplo == A.uplo
+                @test parent(fA) != fA # only active triangle is written to
+            end
+        end
+    end
+end
+
 @testset "$f" for f in weight_funcs
     X = randn(3, 8)
 
@@ -120,18 +141,32 @@ weight_funcs = (weights, aweights, fweights, pweights)
             cor2 = cor(X, wv2, 2)
 
             @testset "cov2cor" begin
-                @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
-                @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
-                @test cov2cor(cov1) ≈ cor1
-                @test cov2cor(cov2) ≈ cor2
-                @test cov2cor(cov1, std1) ≈ cor1
-                @test cov2cor(cov2, std2) ≈ cor2
+                test_isapprox_preserves_symherm_structure(cov2cor, cov(X, dims = 1), cor(X, dims = 1), std(X, dims = 1))
+                test_isapprox_preserves_symherm_structure(cov2cor, cov(X, dims = 2), cor(X, dims = 2), std(X, dims = 2))
+                test_isapprox_preserves_symherm_structure(cov2cor, cov1, cor1)
+                test_isapprox_preserves_symherm_structure(cov2cor, cov2, cor2)
+                test_isapprox_preserves_symherm_structure(cov2cor, cov1, cor1, std1)
+                test_isapprox_preserves_symherm_structure(cov2cor, cov2, cor2, std2)
+            end
+            @testset "StatsBase.cov2cor!" begin
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov(X, dims = 1), cor(X, dims = 1), std(X, dims = 1))
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov(X, dims = 2), cor(X, dims = 2), std(X, dims = 2))
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov1, cor1)
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov2, cor2)
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov1, cor1, std1)
+                test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov2, cor2, std2)
             end
             @testset "cor2cov" begin
-                @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
-                @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
-                @test cor2cov(cor1, std1) ≈ cov1
-                @test cor2cov(cor2, std2) ≈ cov2
+                test_isapprox_preserves_symherm_structure(cor2cov, cor(X, dims = 1), cov(X, dims = 1), std(X, dims = 1))
+                test_isapprox_preserves_symherm_structure(cor2cov, cor(X, dims = 2), cov(X, dims = 2), std(X, dims = 2))
+                test_isapprox_preserves_symherm_structure(cor2cov, cor1, cov1, std1)
+                test_isapprox_preserves_symherm_structure(cor2cov, cor2, cov2, std2)
+            end
+            @testset "StatsBase.cor2cov!" begin
+                test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor(X, dims = 1), cov(X, dims = 1), std(X, dims = 1))
+                test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor(X, dims = 2), cov(X, dims = 2), std(X, dims = 2))
+                test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor1, cov1, std1)
+                test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor2, cov2, std2)
             end
         end
     end
@@ -198,41 +233,24 @@ weight_funcs = (weights, aweights, fweights, pweights)
                 cor2 = cor(X, wv2, 2)
 
                 @testset "cov2cor" begin
-                    @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
-                    @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
-                    @test cov2cor(cov1, std1) ≈ cor1
-                    @test cov2cor(cov2, std2) ≈ cor2
+                    test_isapprox_preserves_symherm_structure(cov2cor, cov1, cor1)
+                    test_isapprox_preserves_symherm_structure(cov2cor, cov2, cor2)
+                    test_isapprox_preserves_symherm_structure(cov2cor, cov1, cor1, std1)
+                    test_isapprox_preserves_symherm_structure(cov2cor, cov2, cor2, std2)
                 end
-
-                @testset "cov2cor!" begin
-                    tmp_cov1 = copy(cov1)
-                    @test !(tmp_cov1 ≈ cor1)
-                    StatsBase.cov2cor!(tmp_cov1, std1)
-                    @test tmp_cov1 ≈ cor1
-
-                    tmp_cov2 = copy(cov2)
-                    @test !(tmp_cov2 ≈ cor2)
-                    StatsBase.cov2cor!(tmp_cov2, std2)
-                    @test tmp_cov2 ≈ cor2
+                @testset "StatsBase.cov2cor!" begin
+                    test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov1, cor1)
+                    test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov2, cor2)
+                    test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov1, cor1, std1)
+                    test_isapprox_preserves_symherm_structure(StatsBase.cov2cor!, cov2, cor2, std2)
                 end
-
                 @testset "cor2cov" begin
-                    @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
-                    @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
-                    @test cor2cov(cor1, std1) ≈ cov1
-                    @test cor2cov(cor2, std2) ≈ cov2
+                    test_isapprox_preserves_symherm_structure(cor2cov, cor1, cov1, std1)
+                    test_isapprox_preserves_symherm_structure(cor2cov, cor2, cov2, std2)
                 end
-
-                @testset "cor2cov!" begin
-                    tmp_cor1 = copy(cor1)
-                    @test !(tmp_cor1 ≈ cov1)
-                    StatsBase.cor2cov!(tmp_cor1, std1)
-                    @test tmp_cor1 ≈ cov1
-
-                    tmp_cor2 = copy(cor2)
-                    @test !(tmp_cor2 ≈ cov2)
-                    StatsBase.cor2cov!(tmp_cor2, std2)
-                    @test tmp_cor2 ≈ cov2
+                @testset "StatsBase.cor2cov!" begin
+                    test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor1, cov1, std1)
+                    test_isapprox_preserves_symherm_structure(StatsBase.cor2cov!, cor2, cov2, std2)
                 end
             end
         end