Thread Subject: 2D Orthogonal Polynomial Surface Fitting

Hello everyone,

I would really appreciate all of your help. I've scoured the file exchange, looking for a way of fitting a 2D surface with polynomials. After using two (both of which are great! Links at bottom of post) functions on the file exchange, I discovered that orthogonal polynomials are needed, such as Chebyshev or Legendre.

However, I can't for the life of me find any code which does this already. Has anyone else approached this problem, or can advise me on how to proceed?

Thanks a million,

Fran


http://www.mathworks.com/matlabcentral/fileexchange/10065-polyfitn
http://www.mathworks.com/matlabcentral/fileexchange/13719-2d-weighted-polynomial-fitting-and-evaluation

"Fran " <fp04@ic.ac.uk> wrote in message <ho0f6o$10l$1@fred.mathworks.com>...
> Hello everyone,
>
> I would really appreciate all of your help. I've scoured the file exchange, looking for a way of fitting a 2D surface with polynomials. After using two (both of which are great! Links at bottom of post) functions on the file exchange, I discovered that orthogonal polynomials are needed, such as Chebyshev or Legendre.
>
> However, I can't for the life of me find any code which does this already. Has anyone else approached this problem, or can advise me on how to proceed?
====

Well, as I'm sure you noticed, there are lots of file exchange submissions that do this for 1D, e.g.

http://www.mathworks.com/matlabcentral/fileexchange/20932-legendre-polynomial-fitting

What exactly are your needs in terms of speed? The extension of these tools to 2D is to assume that your 2D orthogonal polynomials are tensor products of 1D polynomials.
So, if the matrix A are the samples of your 2D surface, a simple way to generalize these tools to 2D would be to use a for-loop to do a 1D fit to all the columns of A, followed by the same thing on all the columns.

Unfortunately I need the fitting to be relatively quick - I will be utilising the code in an optimisation loop with (probably) hundreds of iterations.

I would find it very surprising if this kind of thing wasn't out there, somewhere, already. I think it wold be such a common problem!

Thanks for your help, using a for loop will be the first place I start if nothing else becomes available from this thread! :)

Fran

"Fran " <fp04@ic.ac.uk> wrote in message <ho0f6o$10l$1@fred.mathworks.com>...
> Hello everyone,
>
> I would really appreciate all of your help. I've scoured the file exchange, looking for a way of fitting a 2D surface with polynomials. After using two (both of which are great! Links at bottom of post) functions on the file exchange, I discovered that orthogonal polynomials are needed, such as Chebyshev or Legendre.
>
> However, I can't for the life of me find any code which does this already. Has anyone else approached this problem, or can advise me on how to proceed?
>

Why are they "needed"?

If you need to fit such a high order polynomial model
that orthogonal polynomials are necessary, then I'll
contend that you are trying to fit too high an order
model anyway. What will you get out of all of those
dozens of coefficients?

Just use a tool that can fit your surface from data. So
use splines, gridfit (on the FEX, and really a form of
spline model itself), radial basis functions, neural nets,
etc.

http://www.mathworks.com/matlabcentral/fileexchange/8998

http://www.mathworks.com/products/splines/

http://www.mathworks.com/products/neuralnet/

http://www.mathworks.com/products/curvefitting/

John

"Fran " <fp04@ic.ac.uk> wrote in message <ho0f6o$10l$1@fred.mathworks.com>...
> Hello everyone,
>
> I would really appreciate all of your help. I've scoured the file exchange, looking for a way of fitting a 2D surface with polynomials. After using two (both of which are great! Links at bottom of post) functions on the file exchange, I discovered that orthogonal polynomials are needed, such as Chebyshev or Legendre.

It seems to me those are separate issue: fitting and orthogonal functions. Fitting is projection of a function on a (polynomial) subspace; orthogonal polynomials are a requirement on a selected basis of a subspace. You can fit with whatever subspace, and represent this subspace with orthogonal functions or not.

The first problem call for fitting routine, the second requires change of the basis. They are not directly related.

You can start out with *any* basis, and do some sort of Gram-schmidt process to transform into orthogonal basis. Legendre and Chebychev are known as orthogonal because their analytic expression can be conveniently expressed, that's all.

Bruno

"Fran " <fp04@ic.ac.uk> wrote in message <ho13l2$m42$1@fred.mathworks.com>...
> Unfortunately I need the fitting to be relatively quick - I will be utilising the code in an optimisation loop with (probably) hundreds of iterations.
>
> I would find it very surprising if this kind of thing wasn't out there, somewhere, already. I think it wold be such a common problem!
>
> Thanks for your help, using a for loop will be the first place I start if nothing else becomes available from this thread! :)
=================

If you're still convinced you need this, in spite of the comments above, one way to go would be to convert the 1D fitting operator to matrix form, which we know MATLAB can manipulate very fast.

Namely, if your surface samples form an NxN matrix A, you would fit each column of
eye(N) resulting in corresponding columns of a matrix K. The 2D fit would then be given by

result=(K*(K*A).').'


and you could reuse K in each iteration of your optimization loop.

On 20 Mar, 09:23, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:

> orthogonal polynomials are a requirement on a selected basis of a subspace.

Wrong.

Orthogonal bases are simpler or more convenient to deal with
than non-orthogonal bases, but there are no requirements anywhere
that any one basis needs to be orthogonal. Once you have a
vector space, you can choose to deal with any bais you want,
orthgonal or not. The key property to look out for, is whether
or not the basis is complete.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <3bb7e6cd-8320-47a2-8acf-c4749762fa3c@g10g2000yqh.googlegroups.com>...

> Orthogonal bases are simpler or more convenient to deal with
> than non-orthogonal bases, but there are no requirements anywhere
> that any one basis needs to be orthogonal. Once you have a
> vector space, you can choose to deal with any bais you want,
> orthgonal or not. The key property to look out for, is whether
> or not the basis is complete.

Completeness isn't the only consideration. There are numerical considerations as well.
For example, the basis 1,x,x^2,x^3, is a complete basis for the space of cubic polynomials, but it is more numerically stable to fit using orthogonal polynomials.

Rune Allnor <allnor@tele.ntnu.no> wrote in message <3bb7e6cd-8320-47a2-8acf-c4749762fa3c@g10g2000yqh.googlegroups.com>...
> On 20 Mar, 09:23, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
>
> > orthogonal polynomials are a requirement on a selected basis of a subspace.
>
> Wrong.
>

I think you just miss-understood what I wrote.

Bruno

On 20 Mar, 19:01, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message <3bb7e6cd-8320-47a2-8acf-c4749762f...@g10g2000yqh.googlegroups.com>...
> > On 20 Mar, 09:23, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
>
> > > orthogonal polynomials are a requirement on a selected basis of a subspace.
>
> > Wrong.
>
> I think you just miss-understood what I wrote.

No, I don't. I already *know* what you mean, so I don't have
to *read* your posts. Sounds familiar?

Rune

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ho20lc$jdg$1@fred.mathworks.com>...

> You can start out with *any* basis, and do some sort of Gram-schmidt process to transform into orthogonal basis. Legendre and Chebychev are known as orthogonal because their analytic expression can be conveniently expressed, that's all.
>

Just an addition note that an even better way to carry out orthogonalization is (i) to build a Gram matrix - a matrix where each element is of the scalar products of two functions, then using (ii) SVD/eig to build the orthogonal basis.

Bruno

On 20 Mar, 19:12, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> "Bruno Luong" <b.lu...@fogale.findmycountry> wrote in message <ho20lc$jd...@fred.mathworks.com>...
> > You can start out with *any* basis, and do some sort of Gram-schmidt process to transform into orthogonal basis. Legendre and Chebychev are known as orthogonal because their analytic expression can be conveniently expressed, that's all.
>
> Just an addition note that an even better way to carry out orthogonalization is (i) to build a Gram matrix - a matrix where each element is of the scalar products of two functions, then using (ii) SVD/eig to build the orthogonal basis.

Why?

There is no need to compute SVDs or EVD if all you want is
an orthogonal basis.

Your approach (i) is little more than a QR decomposition.
Once you have that, it's no need to use an additional SVD.
In fact, the SVD and EIG algorithms is based on chains of
repeated QR decompositions.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <d68aa588-4272-4618-a5f0-acccd15ca142@b30g2000yqd.googlegroups.com>...
> On 20 Mar, 19:12, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> > "Bruno Luong" <b.lu...@fogale.findmycountry> wrote in message <ho20lc$jd...@fred.mathworks.com>...
> > > You can start out with *any* basis, and do some sort of Gram-schmidt process to transform into orthogonal basis. Legendre and Chebychev are known as orthogonal because their analytic expression can be conveniently expressed, that's all.
> >
> > Just an addition note that an even better way to carry out orthogonalization is (i) to build a Gram matrix - a matrix where each element is of the scalar products of two functions, then using (ii) SVD/eig to build the orthogonal basis.
>
> Why?
>
> There is no need to compute SVDs or EVD if all you want is
> an orthogonal basis.
>
> Your approach (i) is little more than a QR decomposition.
> Once you have that, it's no need to use an additional SVD.
> In fact, the SVD and EIG algorithms is based on chains of
> repeated QR decompositions.
>

Sorry Rune, I have desire in explaining to you. My post is for OP.

Bruno

Thanks for your help everyone. It is really appreciated!

John:

I do need a large number of polynomial coefficients. This is because I am trying to reconstruct a smooth surface from an undersampled dataset by exploiting the fact that it is sparse in the polynomial domain. I need to calculate a large number of polynomial coefficients in order to 'see' that the majority of coefficients are zero, and am therefore at the correct surface solution.

Bruno:

You are correct about the fitting and choice of basis being separate issues, as I can always project from the non-orthogonal to orthogonal basis. For convenience, I'll just start directly my fitting to Legendre or Chebyshev's.

Matt J:

Very helpful, this matrix approach is the obvious direction to take.

Thanks again!