Ensure that the output of X_A_Xt etc. is symmetric (#191)

JuliaStats · Oct 17, 2023 · 283d865 · 283d865
1 parent 187f741
commit 283d865
Show file tree

Hide file tree

Showing 7 changed files with 70 additions and 71 deletions.
diff --git a/src/chol.jl b/src/chol.jl
@@ -93,26 +93,26 @@ end
 
 # tri products
 
-function X_A_Xt(A::Cholesky, X::AbstractMatrix)
+function X_A_Xt(A::Cholesky, X::AbstractMatrix{<:Real})
     @check_argdims size(A, 1) == size(X, 2)
     Z = X * chol_lower(A)
-    return Z * transpose(Z)
+    return Symmetric(Z * transpose(Z))
 end
 
-function Xt_A_X(A::Cholesky, X::AbstractMatrix)
+function Xt_A_X(A::Cholesky, X::AbstractMatrix{<:Real})
     @check_argdims size(A, 1) == size(X, 1)
     Z = chol_upper(A) * X
-    return transpose(Z) * Z
+    return Symmetric(transpose(Z) * Z)
 end
 
-function X_invA_Xt(A::Cholesky, X::AbstractMatrix)
+function X_invA_Xt(A::Cholesky, X::AbstractMatrix{<:Real})
     @check_argdims size(A, 1) == size(X, 2)
     Z = X / chol_upper(A)
-    return Z * transpose(Z)
+    return Symmetric(Z * transpose(Z))
 end
 
-function Xt_invA_X(A::Cholesky, X::AbstractMatrix)
+function Xt_invA_X(A::Cholesky, X::AbstractMatrix{<:Real})
     @check_argdims size(A, 1) == size(X, 1)
     Z = chol_lower(A) \ X
-    return transpose(Z) * Z
+    return Symmetric(transpose(Z) * Z)
 end
diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl
@@ -166,28 +166,28 @@ end
 
 ### tri products
 
-function X_A_Xt(a::PDiagMat, x::AbstractMatrix)
+function X_A_Xt(a::PDiagMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
     z = a.diag .* transpose(x)
-    return x * z
+    return Symmetric(x * z)
 end
 
-function Xt_A_X(a::PDiagMat, x::AbstractMatrix)
+function Xt_A_X(a::PDiagMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
     z = a.diag .* x
-    return transpose(x) * z
+    return Symmetric(transpose(x) * z)
 end
 
-function X_invA_Xt(a::PDiagMat, x::AbstractMatrix)
+function X_invA_Xt(a::PDiagMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
     z = transpose(x) ./ a.diag
-    return x * z
+    return Symmetric(x * z)
 end
 
-function Xt_invA_X(a::PDiagMat, x::AbstractMatrix)
+function Xt_invA_X(a::PDiagMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
     z = x ./ a.diag
-    return transpose(x) * z
+    return Symmetric(transpose(x) * z)
 end
 
 ### Specializations for `Array` arguments with reduced allocations

diff --git a/src/pdmat.jl b/src/pdmat.jl
@@ -160,28 +160,28 @@ end
 
 ### tri products
 
-function X_A_Xt(a::PDMat, x::AbstractMatrix)
+function X_A_Xt(a::PDMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
     z = x * chol_lower(a.chol)
-    return z * transpose(z)
+    return Symmetric(z * transpose(z))
 end
 
-function Xt_A_X(a::PDMat, x::AbstractMatrix)
+function Xt_A_X(a::PDMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
     z = chol_upper(a.chol) * x
-    return transpose(z) * z
+    return Symmetric(transpose(z) * z)
 end
 
-function X_invA_Xt(a::PDMat, x::AbstractMatrix)
+function X_invA_Xt(a::PDMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
     z = x / chol_upper(a.chol)
-    return z * transpose(z)
+    return Symmetric(z * transpose(z))
 end
 
-function Xt_invA_X(a::PDMat, x::AbstractMatrix)
+function Xt_invA_X(a::PDMat, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
     z = chol_lower(a.chol) \ x
-    return transpose(z) * z
+    return Symmetric(transpose(z) * z)
 end
 
 ### Specializations for `Array` arguments with reduced allocations

diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl
@@ -156,33 +156,28 @@ end
 
 ### tri products
 
-function X_A_Xt(a::PDSparseMat, x::AbstractMatrix)
-    # `*` is not defined for `PtL` factor components,
-    # so we can't use `x * chol_lower(a.chol)`
-    C = a.chol
-    PtL = sparse(C.L)[C.p, :]
-    z = x * PtL
-    z * transpose(z)
+function X_A_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    @check_argdims a.dim == size(x, 2)
+    z = a.mat * transpose(x)
+    return Symmetric(x * z)
 end
 
 
-function Xt_A_X(a::PDSparseMat, x::AbstractMatrix)
-    # `*` is not defined for `UP` factor components,
-    # so we can't use `chol_upper(a.chol) * x`
-    # Moreover, `sparse` is only defined for `L` factor components
-    C = a.chol
-    UP = transpose(sparse(C.L))[:, C.p]
-    z = UP * x
-    transpose(z) * z
+function Xt_A_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    @check_argdims a.dim == size(x, 1)
+    z = a.mat * x
+    return Symmetric(transpose(x) * z)
 end
 
 
-function X_invA_Xt(a::PDSparseMat, x::AbstractMatrix)
+function X_invA_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    @check_argdims a.dim == size(x, 2)
     z = a.chol \ collect(transpose(x))
-    x * z
+    return Symmetric(x * z)
 end
 
-function Xt_invA_X(a::PDSparseMat, x::AbstractMatrix)
+function Xt_invA_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    @check_argdims a.dim == size(x, 1)
     z = a.chol \ x
-    transpose(x) * z
+    return Symmetric(transpose(x) * z)
 end
diff --git a/src/scalmat.jl b/src/scalmat.jl
@@ -147,43 +147,43 @@ end
 
 ### tri products
 
-function X_A_Xt(a::ScalMat, x::AbstractMatrix)
+function X_A_Xt(a::ScalMat, x::AbstractMatrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 2)
-    a.value * (x * transpose(x))
+    return Symmetric(a.value * (x * transpose(x)))
 end
 
-function Xt_A_X(a::ScalMat, x::AbstractMatrix)
+function Xt_A_X(a::ScalMat, x::AbstractMatrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 1)
-    a.value * (transpose(x) * x)
+    return Symmetric(a.value * (transpose(x) * x))
 end
 
-function X_invA_Xt(a::ScalMat, x::AbstractMatrix)
+function X_invA_Xt(a::ScalMat, x::AbstractMatrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 2)
-    (x * transpose(x)) / a.value
+    return Symmetric((x * transpose(x)) / a.value)
 end
 
-function Xt_invA_X(a::ScalMat, x::AbstractMatrix)
+function Xt_invA_X(a::ScalMat, x::AbstractMatrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 1)
-    (transpose(x) * x) / a.value
+    return Symmetric((transpose(x) * x) / a.value)
 end
 
 # Specializations for `x::Matrix` with reduced allocations
-function X_A_Xt(a::ScalMat, x::Matrix)
+function X_A_Xt(a::ScalMat, x::Matrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 2)
-    lmul!(a.value, x * transpose(x))
+    return Symmetric(lmul!(a.value, x * transpose(x)))
 end
 
-function Xt_A_X(a::ScalMat, x::Matrix)
+function Xt_A_X(a::ScalMat, x::Matrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 1)
-    lmul!(a.value, transpose(x) * x)
+    return Symmetric(lmul!(a.value, transpose(x) * x))
 end
 
-function X_invA_Xt(a::ScalMat, x::Matrix)
+function X_invA_Xt(a::ScalMat, x::Matrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 2)
-    _rdiv!(x * transpose(x), a.value)
+    return Symmetric(_rdiv!(x * transpose(x), a.value))
 end
 
-function Xt_invA_X(a::ScalMat, x::Matrix)
+function Xt_invA_X(a::ScalMat, x::Matrix{<:Real})
     @check_argdims LinearAlgebra.checksquare(a) == size(x, 1)
-    _rdiv!(transpose(x) * x, a.value)
+    return Symmetric(_rdiv!(transpose(x) * x, a.value))
 end
diff --git a/test/specialarrays.jl b/test/specialarrays.jl
@@ -70,16 +70,16 @@ using StaticArrays
             @test invquad(A, Y) isa SVector{10, Float64}
             @test invquad(A, Y) ≈ diag(Matrix(Y)' * (Matrix(A) \ Matrix(Y)))
 
-            @test X_A_Xt(A, X) isa SMatrix{10, 10, Float64}
+            @test X_A_Xt(A, X) isa Symmetric{Float64,<:SMatrix{10, 10, Float64}}
             @test X_A_Xt(A, X) ≈ Matrix(X) * Matrix(A) *  Matrix(X)'
 
-            @test X_invA_Xt(A, X) isa SMatrix{10, 10, Float64}
+            @test X_invA_Xt(A, X) isa Symmetric{Float64,<:SMatrix{10, 10, Float64}}
             @test X_invA_Xt(A, X) ≈ Matrix(X) * (Matrix(A) \ Matrix(X)')
 
-            @test Xt_A_X(A, Y) isa SMatrix{10, 10, Float64}
+            @test Xt_A_X(A, Y) isa Symmetric{Float64,<:SMatrix{10, 10, Float64}}
             @test Xt_A_X(A, Y) ≈ Matrix(Y)' * Matrix(A) * Matrix(Y)
 
-            @test Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64}
+            @test Xt_invA_X(A, Y) isa Symmetric{Float64,<:SMatrix{10, 10, Float64}}
             @test Xt_invA_X(A, Y) ≈ Matrix(Y)' * (Matrix(A) \ Matrix(Y))
         end
     end

diff --git a/test/testutils.jl b/test/testutils.jl
@@ -297,23 +297,27 @@ function pdtest_triprod(C, Cmat::Matrix, Imat::Matrix, X::Matrix, verbose::Int)
     Xt = copy(transpose(X))
 
     _pdt(verbose, "X_A_Xt")
-    # default tolerance in isapprox is different on 0.4. rtol argument can be deleted
-    # ≈ form used when 0.4 is no longer supported
-    lhs, rhs = X_A_Xt(C, Xt), Xt * Cmat * X
-    @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))
+    M = X_A_Xt(C, Xt)
+    @test M ≈ Xt * Cmat * X
+    @test issymmetric(M)
     @test_throws DimensionMismatch X_A_Xt(C, rand(n, d + 1))
 
     _pdt(verbose, "Xt_A_X")
-    lhs, rhs = Xt_A_X(C, X), Xt * Cmat * X
-    @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))
+    M = Xt_A_X(C, X)
+    @test M ≈ Xt * Cmat * X
+    @test issymmetric(M)
     @test_throws DimensionMismatch Xt_A_X(C, rand(d + 1, n))
 
     _pdt(verbose, "X_invA_Xt")
-    @test X_invA_Xt(C, Xt) ≈ Xt * Imat * X
+    M = X_invA_Xt(C, Xt)
+    @test M ≈ Xt * Imat * X
+    @test issymmetric(M)
     @test_throws DimensionMismatch X_invA_Xt(C, rand(n, d + 1))
 
     _pdt(verbose, "Xt_invA_X")
-    @test Xt_invA_X(C, X) ≈ Xt * Imat * X
+    M = Xt_invA_X(C, X)
+    @test M ≈ Xt * Imat * X
+    @test issymmetric(M)
     @test_throws DimensionMismatch Xt_invA_X(C, rand(d + 1, n))
 end