Is there an "Algebra of basis-functions"

In my research, I often work with the following kind of problem that essentially reduces to approximating a conditional expectation, say \mathbb{E}[g(X)|X=x], by a linear combination of basis-functions, \mathbb{E}[g(X)|X=x] \approx \sum_{i} a_i \psi_i(x), where the coefficients a_i are e.g. estimated using ordinary least squares. To apply least-squares, we want to compute for a sample x^1,...,x^n the matrix \mathbf{X} = (\psi_i(x^n))_{i,n}, while for pseudo regression, we may wish to compute the matrix \Psi = (\mathbb{E}[\psi_i(X)\psi_j(X)])_{i,j}, etc.

So that raises the question: is there a package that allow you to generate some Basis object based on simple rules (such as: take polynomials in X_1 of degrees 1:d, take Laguerre polynomials in X_2 of degrees 1:d, and take all products such that the degree of the resulting polynomials are at most of degree d) and then constructs and computes the functions \psi_i (and possibly the matrix \Psi, if possible)? (This can then be seen as an algebra of basis-functions in the sense that we join together sets of basis-functions, take their products, etc.)

For my own research I’ve built quite a bit of code to do just this and although I have spent some time optimizing for speed, I’m not using any GPU things and the like, while I think some kernel abstractions are possible here. I’ve also written API for various uses, e.g.

  • to automatically generate a Basis that consists of the orthonormal/orthogonal basis functions for a (product)distribution that you specify
  • to check that a certain subset of the Basis satisfies a certain condition on their moments
  • to split the basis to compute only the basis-functions corresponding to certain variables

If there already is a package that does this (ideally with high performance), I’d like to check if I can adapt my API for that and use it in my projects (and perhaps contribute my code there). If there is not yet a package for this, I’d like to consider making my code into a stand-alone package once my PhD-schedule quiets down a bit.

1 Like

The problem is with the sampling. Normal least-square systems were decided with a finite data point. Maybe the function properties offer a nice solution. Should you sample points or do some fancier calculus tricks to get a more accurate least-square with less cost? That depends on a lot of factors like the function chosen. Numerical integration on its own, the art of sampling points and perhaps doing other fancier tricks to get a descriptive picture of a fixed function, let alone fitting the functions with a least-square. I think that there would probably be some nicely constructed functions where you could perform a nice integration and then perform a minimization on the integral itself instead of having to sample. I’m not that good at calculus though.

It depends on whether or not you want a quick β€œgood enough” solution or you want the best and so on.

2 Likes

I see your point. My question is actually concerned less with the regression part (which serves more as the motivation for wanting this) and more with the algebra of basis-functions part of my post: just like we can make fancy models by composing parts in MTK and nice plots using AlgebraOfGraphics, I was wondering if there is a way to make complicated sets of basis functions from easier ones (and then, if you have a unified API for this, you can think about fitting OLS (simply sampling points and regressing), pseudo-regression (computing the necessary expectations, if possible), or other things).

I think StatsModels.jl comes with this out of the box?

julia> using StatsModels, DataFrames

julia> using StableRNGs; rng = StableRNG(1);

julia> f = @formula(y ~ 1 + a + b + c + b&c)
FormulaTerm
Response:
  y(unknown)
Predictors:
  1
  a(unknown)
  b(unknown)
  c(unknown)
  b(unknown) & c(unknown)

julia> df = DataFrame(y = rand(rng, 9), a = 1:9, b = rand(rng, 9), c = repeat(["d","e","f"], 3))
9Γ—4 DataFrame
 Row β”‚ y          a      b         c
     β”‚ Float64    Int64  Float64   String
─────┼────────────────────────────────────
   1 β”‚ 0.585195       1  0.236782  d
   2 β”‚ 0.0773379      2  0.943741  e
   3 β”‚ 0.716628       3  0.445671  f
   4 β”‚ 0.320357       4  0.763679  d
   5 β”‚ 0.653093       5  0.145071  e
   6 β”‚ 0.236639       6  0.021124  f
   7 β”‚ 0.709684       7  0.152545  d
   8 β”‚ 0.557787       8  0.617492  e
   9 β”‚ 0.05079        9  0.481531  f

julia> f = apply_schema(f, schema(f, df))
FormulaTerm
Response:
  y(continuous)
Predictors:
  1
  a(continuous)
  b(continuous)
  c(DummyCoding:3β†’2)
  b(continuous) & c(DummyCoding:3β†’2)

julia> resp, pred = modelcols(f, df);

julia> pred
9Γ—7 Matrix{Float64}:
 1.0  1.0  0.236782  0.0  0.0  0.0       0.0
 1.0  2.0  0.943741  1.0  0.0  0.943741  0.0
 1.0  3.0  0.445671  0.0  1.0  0.0       0.445671
 1.0  4.0  0.763679  0.0  0.0  0.0       0.0
 1.0  5.0  0.145071  1.0  0.0  0.145071  0.0
 1.0  6.0  0.021124  0.0  1.0  0.0       0.021124
 1.0  7.0  0.152545  0.0  0.0  0.0       0.0
 1.0  8.0  0.617492  1.0  0.0  0.617492  0.0
 1.0  9.0  0.481531  0.0  1.0  0.0       0.481531

julia> coefnames(f)
("y", ["(Intercept)", "a", "b", "c: e", "c: f", "b & c: e", "b & c: f"])
1 Like

Thanks, I didn’t think about that yet; so you could do all the logic for generating a set of functions (adding sets of functions together, adding cross-products, etc.) on the FormulaTerms, and then parse things along the FormulaTerm to do fancy things (e.g. compute expectations when the data follows a product-distribution, etc.).

I’ll keep this in mind and revisit this once I have time to check if this would support all the things that I would want to do and how I can write API for those things!

1 Like

I use Chebyshev polynomials (with transformations) most of the time. I have coded them with derivatives, Smolyak combinations, etc in

You might want to take a look at ApproxFun.jl if you haven’t already. Not sure if it supports everything you need, but it’s the first-stop shop for working with function spaces via orthogonal bases (similar to chebfun.org for Matlab, whose thorough and pedagogical documentation you might find helpful even if you’re using Julia and ApproxFun).

2 Likes

Thank you, that is a good pointer. The way I understand it is that you should provide it with a function already, which is kind of the problem in my case: I’m interested in estimating a conditional expectation by sampled points, so I want to be able to do some form of regression, such as least-squares, and to be able to go up to higher dimensions. For my typical use-case I’m looking at roughly 6 dimensions and in the order of 100-1000 basis-functions. (I know this is not ideal, but the algorithms NEED basis functions with certain properties for this, so fancier approximation methods are off the table for me, unfortunately)

If you’re doing least-squares, ApproxFun.jl can help build the Vandermonde matrix for your orthogonal function basis. See Frequently Asked Questions Β· ApproxFun.jl (the example starts with interpolation, but progresses to regression by the third code block and multivariate regression by the fourth).