[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. ```julia 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. ```julia 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_