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

isproper(::NaturalParametersSpace, ::Type{MvNormalMeanCovariance}, η, conditioner) is not correct #145

Open
Nimrais opened this issue Nov 15, 2023 · 3 comments

Comments

@Nimrais
Copy link
Member

Nimrais commented Nov 15, 2023

It's not possible to pass the isproper check during the creation phase for valid natural multivariate normal parameters.

begin
	using ExponentialFamily
	using LinearAlgebra
	
	weighted_mean = [0.5929739312294535, -0.5765312236076422, -0.27502558591534837, 0.7921902812250886]
	negative_half_precision = [-1.4988377222596585 0.24017412940042754 -0.25091748712776246 0.021297606627674978; 0.24017412940042765 -1.3055257012352184 0.12334134123122117 -0.048756680716141254; -0.25091748712776263 0.1233413412312212 -1.5108823225676407 0.06676100253121688; 0.021297606627674905 -0.04875668071614124 0.06676100253121682 -1.3471915876466825]
	packed_parameters = [0.5929739312294535, -0.5765312236076422, -0.27502558591534837, 0.7921902812250886, -1.4988377222596585, 0.24017412940042765, -0.25091748712776263, 0.021297606627674905, 0.24017412940042754, -1.3055257012352184, 0.1233413412312212, -0.04875668071614124, -0.25091748712776246, 0.12334134123122117, -1.5108823225676407, 0.06676100253121682, 0.021297606627674978, -0.048756680716141254, 0.06676100253121688, -1.3471915876466825]

	(η₁, η₂) = ExponentialFamily.unpack_parameters(
		MvNormalMeanCovariance,
		packed_parameters
	)

	k = div(-1 + isqrt(1 + 4 * length(packed_parameters)), 2)

	# correct sizes pass
	@show length(packed_parameters) < 2 || (length(packed_parameters) !== (k + k^2))
	@show length(η₁) === size(η₂, 1)
	@show size(η₂, 1) === size(η₂, 2)
	
	#matrix is not posdef
	@show isposdef(η₂)

	#but it's acutllay posdef
	@show LinearAlgebra.eigen(-η₂).values

	# this fails
	ExponentialFamily.ExponentialFamilyDistribution(MvNormalMeanCovariance, packed_parameters)
end

outputs

length(packed_parameters) < 2 || length(packed_parameters) !== k + k ^ 2 = false
length(η₁) === size(η₂, 1) = true
size(η₂, 1) === size(η₂, 2) = true
isposdef(η₂) = false
(LinearAlgebra.eigen(-η₂)).values = [1.1309033948541316, 1.296832849037286, 1.3503981343797664, 1.884302955438015]

Error

Parameter vector [0.5929739312294535, -0.5765312236076422, -0.27502558591534837, 0.7921902812250886, -1.4988377222596585, 0.24017412940042765, -0.25091748712776263, 0.021297606627674905, 0.24017412940042754, -1.3055257012352184, 0.1233413412312212, -0.04875668071614124, -0.25091748712776246, 0.12334134123122117, -1.5108823225676407, 0.06676100253121682, 0.021297606627674978, -0.048756680716141254, 0.06676100253121688, -1.3471915876466825] is not a valid natural parameter for distribution ExponentialFamily.MvNormalMeanCovariance.

error(::String)@error.jl:35
@Nimrais Nimrais changed the title isproper(::NaturalParametersSpace, ::Type{MvNormalMeanCovariance}, η, conditioner) not correct isproper(::NaturalParametersSpace, ::Type{MvNormalMeanCovariance}, η, conditioner) is not correct Nov 15, 2023
@bvdmitri
Copy link
Member

This is a known problem (limitation) of the isposdef function from LinearAlgebra. The check for positive definitness is not stable and is affected y small perturbations in the off-diagonal elements. I'm not sure we can do anything about it and IMO it is not really an issue with the ExponentialFamily.jl.

As a workaround, convert your matrix to Hermitian.

@bvdmitri
Copy link
Member

bvdmitri commented Jun 4, 2024

This might be a potential issue for the manifold projection. Recording it here so we don't forget, but those small differences in the off-diagonal entries do indeed lead to errors during the optimization :(

@bvdmitri
Copy link
Member

bvdmitri commented Jun 4, 2024

Perhaps, we can modify our check to something like

return isnothing(conditioner) &&
               length(η₁) === size(η₂, 1) &&
               (size(η₂, 1) === size(η₂, 2)) &&
               isapprox(norm(η₂ - transpose(η₂)), 0.0; atol=1e-10) &&
               isposdef(Hermitian(-η₂))

the extra norm check if the matrix if symmetric or not, and then Hermitian uses only one part of the matrix. WDYT @Nimrais . The check is already expensive and this makes it even more expensive. Another solution would be to bypass this check during the optimization entirely. I'll look into it.

EDIT: I bypassed the check in the optimization as for now

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants