Best multivariate polynomial regression

I have a dataset (x,y) where x is a n-dimensional vector and y is an m-dimensional vector. (m=3, n>2) My goal is to find the best polynomial in x fitting the (x,y) dataset.

The dimension of x is pretty big (right now it is 25), and I don't want to enter manually all the possibilities (ie x1*x3*x5, x1*x4*x6, ...). I can use Matlab, Mathematica and R. How can I do this?

Also, I would be interested in hearing your suggestions about the following problem: how can I choose from the result the most relevant coefficients? (maybe x1*x2 is more relevant than x2*x3)

Thank you

• I don't think you meant to refer to mutinomials did you? It sounds like you are interested in fitting a polynomial in up to 25 variables to fit y. But you say y is 3 dimensio0nal. So does that mean that you are trying to separately fit polynomials to the three components of Y? There needs to be some clarification. Also you haven't mentioned how much data you have. Overfitting can be a major concern if you have 25 variables and are looking at powers and cross products for each of these terms. Jun 19, 2012 at 22:12
• Sorry I meant multivariate. I generate the data myself so I have virtually as much data as I want. Let's say that the size of my dataset is about 200 and that I could decrease the number of variables (hence using m and n). I am not sure if I should fit the three components separately or together, but I think we can start with the simplest way (ie separately). I also suspect that the degree of the polynom is quite low (2 or 3 should suffice). Jun 19, 2012 at 22:38
• I am now getting more puzzled. Why would you be fitting a model to data that you generate yourself? If you generated it yourself wouldn't you alreayd know the model? Jun 19, 2012 at 23:09
• I am analyzing the moments of 2D light distributions that I generate from a path tracer. Jun 20, 2012 at 16:05

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.

• I had the first part (beta = X\Y) but the second part is very useful and exactly what I was looking for, thanks! Jun 22, 2012 at 13:44
• @Flav FWIW, the Eureqa, A.K.A. nutonian or formulize, application is pretty good for this kind of thing: formulize.nutonian.com . For your problem, it will not deal with vector Y very well, so you will have to run it thrice. Jun 22, 2012 at 16:20