From 473582d73b5296e0f66fd81a290dcb66645e4ef0 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 17 Oct 2023 00:05:48 +0200 Subject: [PATCH] Ensure that the output of `X_A_Xt` etc. is symmetric --- src/chol.jl | 16 ++++++++-------- src/pdiagmat.jl | 16 ++++++++-------- src/pdmat.jl | 16 ++++++++-------- src/pdsparsemat.jl | 33 ++++++++++++++------------------- src/scalmat.jl | 32 ++++++++++++++++---------------- test/specialarrays.jl | 8 ++++---- test/testutils.jl | 20 ++++++++++++-------- 7 files changed, 70 insertions(+), 71 deletions(-) diff --git a/src/chol.jl b/src/chol.jl index 9b97602..ec3ec0c 100644 --- 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 index c705b09..a14cf3f 100644 --- 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 index ef8eeaa..bc09e3e 100644 --- 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 index c7c4266..032e0ed 100644 --- 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 index 713cac7..feabe8a 100644 --- 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 index ac1b20f..e896d72 100644 --- 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 index 16a8860..96843c9 100644 --- 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