Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove Symmetric(::Diagonal) type piracy #115

Merged
merged 1 commit into from
Mar 29, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "AbstractGPs"
uuid = "99985d1d-32ba-4be9-9821-2ec096f28918"
authors = ["JuliaGaussianProcesses Team"]
version = "0.2.19"
version = "0.2.20"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
13 changes: 3 additions & 10 deletions src/abstract_gp/finite_gp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ true
"""
function Random.rand(rng::AbstractRNG, f::FiniteGP, N::Int)
m, C_mat = mean_and_cov(f)
C = cholesky(Symmetric(C_mat))
C = cholesky(_symmetric(C_mat))
return m .+ C.U' * randn(rng, promote_type(eltype(m), eltype(C)), length(m), N)
end
Random.rand(f::FiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
Expand Down Expand Up @@ -223,7 +223,7 @@ Distributions.loglikelihood(f::FiniteGP, Y::AbstractMatrix{<:Real}) = sum(logpdf

function Distributions.logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})
m, C_mat = mean_and_cov(f)
C = cholesky(Symmetric(C_mat))
C = cholesky(_symmetric(C_mat))
T = promote_type(eltype(m), eltype(C), eltype(Y))
return -((size(Y, 1) * T(log(2π)) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- m)) ./ 2
end
Expand Down Expand Up @@ -285,19 +285,12 @@ function dtc(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
return first(_compute_intermediates(f, y, u))
end

# Small bit of indirection to work around a cholesky-related bug whereby the interaction
# between `FillArrays` and `Diagonal` and `Cholesky` causes problems.
_cholesky(X) = cholesky(X)
function _cholesky(X::Diagonal{<:Real, <:FillArrays.AbstractFill})
return cholesky(Diagonal(collect(diag(X))))
end

# Factor out computations common to the `elbo` and `dtc`.
function _compute_intermediates(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
consistency_check(f, y, u)
chol_Σy = _cholesky(f.Σy)

A = cholesky(Symmetric(cov(u))).U' \ (chol_Σy.U' \ cov(f, u))'
A = cholesky(_symmetric(cov(u))).U' \ (chol_Σy.U' \ cov(f, u))'
Λ_ε = cholesky(Symmetric(A * A' + I))
δ = chol_Σy.U' \ (y - mean(f))

Expand Down
14 changes: 3 additions & 11 deletions src/posterior_gp/approx_posterior_gp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,6 @@ struct ApproxPosteriorGP{Tapprox, Tprior, Tdata} <: AbstractGP
data::Tdata
end

# If a matrix is `Diagonal`, we generally don't need to wrap it in a `Symmetric`, because
# it's already symmetric. This is used in a couple of places to avoid precisely this and
# having to add specialised methods of e.g. `_cholesky` for complicated wrapped types.
_symmetric(X) = Symmetric(X)
_symmetric(X::Diagonal) = X

"""
approx_posterior(::VFE, fx::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

Expand All @@ -25,14 +19,14 @@ Intelligence and Statistics. 2009.
"""
function approx_posterior(::VFE, fx::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
U_y = _cholesky(_symmetric(fx.Σy)).U
U = cholesky(Symmetric(cov(u))).U
U = cholesky(_symmetric(cov(u))).U

B_εf = U' \ (U_y' \ cov(fx, u))'

b_y = U_y' \ (y - mean(fx))

D = B_εf * B_εf' + I
Λ_ε = cholesky(Symmetric(D))
Λ_ε = cholesky(_symmetric(D))

m_ε = Λ_ε \ (B_εf * b_y)

Expand Down Expand Up @@ -117,7 +111,7 @@ function update_approx_posterior(
)
U11 = f_post_approx.data.U
C12 = cov(u.f, f_post_approx.data.z, u.x)
C22 = Symmetric(cov(u))
C22 = _symmetric(cov(u))
U = update_chol(Cholesky(U11,'U', 0), C12, C22).U
U22 = U[end-length(u)+1:end, end-length(u)+1:end]
U12 = U[1:length(f_post_approx.data.z), end-length(u)+1:end]
Expand Down Expand Up @@ -154,8 +148,6 @@ function update_approx_posterior(
return ApproxPosteriorGP(VFE(), f_post_approx.prior, cache)
end

LinearAlgebra.Symmetric(X::Diagonal) = X

# AbstractGP interface implementation.

function Statistics.mean(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
Expand Down
2 changes: 1 addition & 1 deletion src/posterior_gp/posterior_gp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ where `m` and `k` are the mean function and kernel of `fx.f` respectively.
"""
function posterior(fx::FiniteGP, y::AbstractVector{<:Real})
m, C_mat = mean_and_cov(fx)
C = cholesky(Symmetric(C_mat))
C = cholesky(_symmetric(C_mat))
δ = y - m
α = C \ δ
return PosteriorGP(fx.f, (α=α, C=C, x=fx.x, δ=δ))
Expand Down
12 changes: 12 additions & 0 deletions src/util/common_covmat_ops.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,15 @@
# If a matrix is `Diagonal`, we generally don't need to wrap it in a `Symmetric`, because
# it's already symmetric. This is used in a couple of places to avoid precisely this and
# having to add specialised methods of e.g. `_cholesky` for complicated wrapped types.
_symmetric(X) = Symmetric(X)
_symmetric(X::Diagonal) = X

# Small bit of indirection to work around a cholesky-related bug whereby the interaction
# between `FillArrays` and `Diagonal` and `Cholesky` causes problems.
_cholesky(X) = cholesky(X)
function _cholesky(X::Diagonal{<:Real, <:FillArrays.AbstractFill})
return cholesky(Diagonal(collect(diag(X))))
end

"""
update_chol(chol::Cholesky, C12::AbstractMatrix, C22::AbstractMatrix)
Expand Down
6 changes: 0 additions & 6 deletions test/abstract_gp/finite_gp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,12 +169,6 @@ end
# atol=1e-8, rtol=1e-8,
# )

# Supporting utility function.
@testset "_cholesky" begin
D = Diagonal(Fill(5.0, 10))
@test AbstractGPs._cholesky(D).U ≈ AbstractGPs._cholesky(collect(D)).U
end

# Ensure that the elbo is close to the logpdf when appropriate.
@test elbo(y, ŷ, fx) isa Real
@test elbo(y, ŷ, fx) ≈ logpdf(y, ŷ)
Expand Down
16 changes: 0 additions & 16 deletions test/posterior_gp/approx_posterior_gp.jl
Original file line number Diff line number Diff line change
@@ -1,25 +1,9 @@
@testset "approx_posterior_gp" begin

@testset "_symmetric" begin
@testset "Matrix" begin
X = randn(5, 5)
@test AbstractGPs._symmetric(X) isa Symmetric
@test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
end
@testset "Diagonal" begin
X = Diagonal(randn(5))
@test AbstractGPs._symmetric(X) isa Diagonal
@test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
end
end

rng = MersenneTwister(123456)
N_cond = 3
N_a = 5
N_b = 6

@test Symmetric(Diagonal(randn(rng, 5))) isa Diagonal

# Specify prior.
f = GP(sin, Matern32Kernel())

Expand Down
19 changes: 19 additions & 0 deletions test/util/common_covmat_ops.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,23 @@
@testset "common_covmat_ops" begin
# Supporting utility functions.
@testset "_cholesky" begin
D = Diagonal(Fill(5.0, 10))
@test AbstractGPs._cholesky(D).U ≈ AbstractGPs._cholesky(collect(D)).U
end

@testset "_symmetric" begin
@testset "Matrix" begin
X = randn(5, 5)
@test AbstractGPs._symmetric(X) isa Symmetric
@test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
end
@testset "Diagonal" begin
X = Diagonal(randn(5))
@test AbstractGPs._symmetric(X) isa Diagonal
@test collect(AbstractGPs._symmetric(X)) == collect(Symmetric(X))
end
end

@testset "update_chol" begin
X = rand(5)
y = rand(5)
Expand Down