Solved – Multivariate orthogonal polynomial regression

regression

As a means of motivating the question, consider a regresison problem where we seek to estimate $Y$ using observed variables $\{ a, b \}$

When doing multivariate polynomial regresison, I try to find the optimal paramitization of the function

$$f(y)=c_{1}a+c_{2}b+c_{3}a^{2}+c_{4}ab+c_{5}b^{2}+\cdots$$

which best fit the data in a least squared sense.

The problem with this, however, is that the parameters $c_i$ are not independent. Is there a way to do the regression on a different set of "basis" vectors which are orthogonal? Doing this has many obvious advantages

1) the coefficients are no longer correlated.
2) the values of the $c_i$'s themselves no longer depend on the degree of the coefficients.
3) This also has the computational advantage of being able to drop the higher order terms for a coarser but still accurate approximation to the data.

This is easily achieved in the single variable case using orthogonal polynomials, using a well studied set such as the Chebyshev Polynomials. It's not obvious however (to me anyway) how to generalize this! It occured to me that I could mutiply chebyshev polynomials pairwise, but I'm not sure if that is the mathematically correct thing to do.

Your help is appreciated

Best Answer

For completion's sake (and to help improve the stats of this site, ha) I have to wonder if this paper wouldn't also answer your question?

ABSTRACT: We discuss the choice of polynomial basis for approximation of uncertainty propagation through complex simulation models with capability to output derivative information. Our work is part of a larger research effort in uncertainty quantification using sampling methods augmented with derivative information. The approach has new challenges compared with standard polynomial regression. In particular, we show that a tensor product multivariate orthogonal polynomial basis of an arbitrary degree may no longer be constructed. We provide sufficient conditions for an orthonormal set of this type to exist, a basis for the space it spans. We demonstrate the benefits of the basis in the propagation of material uncertainties through a simplified model of heat transport in a nuclear reactor core. Compared with the tensor product Hermite polynomial basis, the orthogonal basis results in a better numerical conditioning of the regression procedure, a modest improvement in approximation error when basis polynomials are chosen a priori, and a significant improvement when basis polynomials are chosen adaptively, using a stepwise fitting procedure.

Otherwise, the tensor-product basis of one-dimensional polynomials is not only the appropriate technique, but also the only one I can find for this.

Related Solutions

Solved – Explicit formulation of predict() output of an orthogonal polynomial regression

This does not handle orthogonal polynomials and implements ordinary (not orthogonal which IMHO are not worth the effort with modern stat computing algorithms) polynomials and restricted cubic splines (natural splines), but the R rms package is useful in a similar context.

require(rms)
f <- ols(y ~ x1 + rcs(x2,4) + pol(x3,3)) # spline for x2 with 4 knots; cubic for x3
Function(f)   # derive R function containing algebraic expression for predictions
latex(f)      # compose LaTeX markup for pretty algebraic form of fitted model

Solved – Best multivariate polynomial regression

In Matlab, let Y be the $k\times 3$ matrix whose rows are observations of $y$, and similarly let X be the $k\times n$ matrix of corresponding $X$ observations. To regress Y against X, try

beta = X \ Y;

The variable beta will now be $n\times 3$ and have the regression coefficients.

This is for a model that is 'linear' in your input $X$ variables. To populate a matrix with the $n(n+1)/2$ pairs of products of the $X$ variables is a bit more tricky without using a bunch of for loops. One could do something like:

Z = bsxfun(@times,X,permute(X,[1,3,2]));  % this is k x n x n
Z = reshape(Z,size(Z,1),[]);              % this is k x (n^2)
o2n = 1:size(X,2);                        % elements 1 to n
keepidx = bsxfun(@ge,o2n,permute(o2n,[1,3,2]));  % keep only if i >= j
keepidx = keepidx(:)';
Z = Z(:,keepidx);

At the end of this, modulo my programming errors, Z should be $k \times n(n+1)/2$ and have the order 2 powers of elements of X. Then you can try

beta = [X,Z] \ Y;

Variable selection is another matter entirely. As a first pass, BIC seems to work OK for multivariate multiple regression.

Best Answer

Related Solutions

Solved – Explicit formulation of predict() output of an orthogonal polynomial regression

Solved – Best multivariate polynomial regression

Related Question