Review

[olivierverdier]

JuliaCon/proceedings-review#179

The idea of this package is clear: use Manifolds.jl to represent statistical manifolds, in this case exponential families.

There is one fundamental flaw with this package: a manifold is independent of the representation. The choice of coordinates does not matter.

To be perfectly clear, the "natural parameters" (relative to a choice of sufficient statistics) are simply coordinates of a manifold, and should be implemented as such.

For instance, there is just one half line (possibly with different metrics, possibly with different coordinates). Therefore, the implementations of "shifted negative numbers", or "shifted positive numbers" make no sense. It leads to code duplication and wheel reinvention. You should only ever need to use PositiveNumbers from Manifolds.jl (adding different coordinate systems, or custom metrics).

Another (less important) issue that I have is the simplicity of the implemented families: there are no multivariate normal distributions, for instance. Interesting cases could be normal distributions with isotropic variance. Again, that should not be difficult to do using the existing bricks provided in Manifolds.jl.

[Nimrais]

@olivierverdier Thank you for the review. Before I address other parts of your review, there is some misunderstanding about what the package already implements. The paper directly provides examples of two different multivariate normal families. Also, you can find that they are indeed implemented in the code https://github.com/ReactiveBayes/ExponentialFamilyManifolds.jl/blob/main/src/natural_manifolds/normal.jl#L27-L32 (full covariance) and https://github.com/ReactiveBayes/ExponentialFamilyManifolds.jl/blob/main/src/natural_manifolds/normal.jl#L51-L56 (a form of an isotropic gaussian).

[olivierverdier]

Fair enough, sorry I overlooked those!

[Nimrais]

@olivierverdier Thank you for this insightful review about manifold representation independence. However, I agree with your main point: a manifold is independent of the representation. However, from my point of view, the real manifold is a set of exponential distributions specified with base measure, sufficient statistics, and domain. So, however, as you correctly noted, shifted negative numbers makes no sense; representing them only through manifold from Manifolds.jl also makes little sense. It is a coincidence that some coordinates are the same, and I can reuse them.

There are two reasons for that: (1) I want to avoid the transformation of natural parameters into some manifolds-friendly form (for example, for multivariate normal, it becomes too slow I need to invert the matrix on each iteration, and (2) Fisher information matrix in natural parameters is just a hessian of the log partition, so to evaluate I do not need to integrate.

The key insight here is that while geometrically, these might be the same manifolds, the statistical structure (Fisher metric, exponential family structure) provides additional context that makes specific coordinate systems computationally privileged. Natural parameters aren't just another coordinate system - they're the coordinate system where the Fisher information takes its most straightforward form from the computational point of view.

So, to accommodate these two goals (don't change coordinates too much) and be able to compute fisher information only through derivation, my refactoring proposal for this issue is the following one:

# The manifold is the space of distributions
struct ExponentialFamilyManifold{𝔽,T,D,C} <: AbstractManifold{𝔽}
    dims::D
    conditioner::C
end

# Natural parameters form an atlas
struct NaturalParametersAtlas{T} <: AbstractAtlas{ℝ} end

# Fisher Information as the metric
struct FisherInformationMetric <: RiemannianMetric end

# Decorate with the Fisher metric
const StatisticalManifold{𝔽,T,D,C} = MetricManifold{𝔽,ExponentialFamilyManifold{𝔽,T,D,C},FisherInformationMetric}

# Implementation of chart operations
function Manifolds.get_parameters(M::ExponentialFamilyManifold, A::NaturalParametersAtlas, i, p)
    # Convert distribution p to its natural parameters
    return natural_parameters(p)
end

function Manifolds.get_point(M::ExponentialFamilyManifold, A::NaturalParametersAtlas, i, η)
    # Convert natural parameters η back to the distribution
    return distribution_from_natural(M, η)
end

# Implement the Fisher metric - it's straightforward in natural parameters
function Manifolds.inner(
    M::StatisticalManifold,
    p,
    X,
    Y,
    A::NaturalParametersAtlas = get_default_atlas(M)
)
    i = get_chart_index(M, A, p)
    η = get_parameters(M, A, i, p)
    # In natural parameters, Fisher metric is just the Hessian of log partition function
    return dot(X, hessian_log_partition(M, η) * Y)
end

This approach allows the implementation of another set of coordinates if needed and, as an extra bonus, allows the implementation of an approximation to an exponential map on the StatisticalManifold:

struct ChartFirstOrderRetraction{TA<:AbstractAtlas} <: AbstractRetractionMethod
    A::TA
end

# Second-order retraction in charts
struct ChartSecondOrderRetraction{TA<:AbstractAtlas} <: AbstractRetractionMethod
    A::TA
end

# Implementation of statistical manifold
function retract_chart_first_order!(M::StatisticalManifold, q, p, X, m::ChartFirstOrderRetraction)
    i = get_chart_index(M, m.A, p)
    p_chart = get_parameters(M, m.A, i, p)
    X_chart = get_coordinates(M, p, X, induced_basis(M, m.A, i))
    # Simple addition in chart coordinates
    q_chart = p_chart + X_chart
    return get_point!(M, q, m.A, i, q_chart)
end

For a rough draft of this implementation, you can check #26 for one specific distribution, Normal-Gamma.

So, would such an implementation be able to resolve your concern?

[olivierverdier]

Not really. You have to change coordinates back and forth.
(it's not much, right? either $x \mapsto x - b$ (ShiftedPositiveNumber), or $x \mapsto -x -b$ (ShiftedNegativeNumber), or $x \mapsto -x$ (SymmetricNegativeDefinite). These operations are pretty much for free.).

As long as these three classes are removed I have no further objections.

---

You can work with natural coordinates all you want, but the actual point in the manifold may be encoded with a (very simple) change of variable.

The implementation of the exponential distribution could be for instance:

"""
    get_natural_manifold_base(::Type{Exponential}, ::Tuple{}, conditioner=nothing)

Get the natural manifold base for the `Exponential` distribution.
"""
function get_natural_manifold_base(::Type{Exponential}, ::Tuple{}, conditioner=nothing)
    return PositiveNumbers()
end

"""
    partition_point(::Type{Exponential}, ::Tuple{}, p, conditioner=nothing)

Converts the `point` to a compatible representation for the natural manifold of type `Exponential`.
"""
function partition_point(::Type{Exponential}, ::Tuple{}, p, conditioner=nothing)
    return -p
end

If you really need some other manifold for some reason, sure, go ahead, but the three ones that you created should be removed.

[Nimrais]

@olivierverdier Hi, sorry for the month delay, but this package is a foundation package for many packages in the ecosystem, so changes here broke many other packages. I needed to prepare changes, or at least the possibility of their implementation, beforehand. Your main concern is addressed in this PR. If you agree, you can approve the PR and close the issue once the PR is merged.

But anyway, thanks for such a careful review! Your suggestion turned out to be a really good one!

[Nimrais]

@olivierverdier We need to move forward with updates in other repos. Also, because this release is breaking, we would need to adapt all future PRs. So I suggest that once you are ready, close the issue, not the PR.