[Enhancement]: Generalize approximate_kernel_expectation to non-square matrices

Feature/Enhancement Description

Considering the current implementation (below), we note that gbar is assumed to be square of dims equal to distribution, but this is not always the case.

function approximate_kernel_expectation(method::AbstractApproximationMethod, g::Function, m::AbstractVector{T}, P::AbstractMatrix{T}) where {T <: Real}
    ndims = length(m)

    weights = getweights(method, m, P)
    points  = getpoints(method, m, P)

    gbar = zeros(ndims, ndims)
    foreach(zip(weights, points)) do (weight, point)
        axpy!(weight, g(point), gbar) # gbar = gbar + weight * g(point)
    end

    return gbar
end

Motivation / Use Case

RxGP.jl forms kernel function expectations that are non-square.

Proposed Solution

Per @HoangMHNguyen's work, the below seems to work well.

function approximate_kernel_expectation(method::AbstractApproximationMethod, g::Function, m::AbstractVector{T}, P::AbstractMatrix{T}) where {T <: Real}
    weights = getweights(method, m, P)
    points  = getpoints(method, m, P)

    gbar = g(m) .* 0.0
    foreach(zip(weights, points)) do (weight, point)
        axpy!(weight, g(point), gbar) # gbar = gbar + weight * g(point)
    end
    return gbar
end

Alternatives Considered

No response

Example Use Cases

No response

Priority (from your perspective)

Critical for my use case

Related Issues / Discussions

No response

Additional Context

No response

Start at approximate_kernel_expectation and inspect how getweights, getpoints, and g(point) are used. Check the current tests or call sites for square-matrix assumptions, then verify that a non-square result shaped like g(m) is accumulated correctly without changing existing square-matrix behavior.