# Polynomial approximation

Is there any universal method to fill this matrix for any $n$ value:

$\textbf{A} = \left[ \matrix{n & \sum x_i & \sum x_i^2 & \cdots & \sum x_i^n \cr \sum x_i & \sum x_i^2 & \sum x_i^3 & \cdots & \sum x_i^{n+1} \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr \sum x_i^n & \sum x_i^{n+1} & \sum x_i^{n+2} &\cdots & \sum x_i^{n+n}} \right]$

Or maybe, it is easier to carry out with Vandermonde matrix to calculate the polynomial coefficients?

I would be grateful for any suggestions and Fortran code.

• Are you hoping for something other than for loops? Use a loop to compute the first row then use another loop to compute the rest of the rows by copying values from the row above and then compute the last value in the row. – Doug Lipinski Aug 13 '15 at 11:23
• Is $A(2,1)$ wrong? It should be $\sum{x_i}$, instead of $\sum{x^n_i}$, right ? – Tolga Birdal Aug 13 '15 at 11:56
• Right, my bad. Corrected. – Beginner in fort Aug 13 '15 at 15:27
• @DougLipinski My intention is to find out how can I fill this matrix for any n degree, because I can not imagine any universal algorythm for that - only for particular degree, not any. I have calculated polynomial coefficients using the Vandermonde matrix and appended it with y data values and then I have used Gauss elimination. But I would also like to know how can I do it using another method or to know, whether it is not efficient. – Beginner in fort Aug 13 '15 at 15:32
• I don't understand what you mean. Why will loops not work for "any n"? Just use n to define your loop limits. Here's a fully functional solution for MATLAB (just define n and x first): A = zeros(n+1); for j=0:n, A(1,j+1)=sum(x.^j); end, for i=1:n, A(i+1,1:n)=A(i,2:n+1); A(i+1,n+1)=sum(x.^(n+i)); end – Doug Lipinski Aug 14 '15 at 1:55

You can easily do this with for loops. Just use n to define your loop limits. Here's a fully functional solution for MATLAB (just define n and x first):

%pre-allocate A
A = zeros(n+1);

%first row:
for j=0:n
A(1,j+1)=sum(x.^j);
end

%rows 2 through n
for i=1:n
A(i+1,1:n)=A(i,2:n+1); %copy from previous row
A(i+1,n+1)=sum(x.^(n+i)); %compute last entry in row
end


This method requires (2*n+1)*(2*numel(x)) operations so if x has about n values the method is $\mathcal{O}(n^2)$.

$\mathbf{A}$ is an $(n+1) \times (n+1)$ matrix. It can be obtained as follows:

$\textbf{A} = \left[ \matrix{1 & 1 & 1 & \cdots & 1 \cr x_0 & x_1 & x_2 & \cdots & x_{n} \cr \vdots & \vdots & \vdots & \cdots & \vdots \cr x_0^n & x_1^n & x_2^n & \cdots & x_n^{n}}\right] \left[ \matrix{1 & x_0 & x_0^2 & \cdots & x_0^n \cr 1 & x_1 & x_1^2 & \cdots & x_1^{n} \cr \vdots & \vdots & \vdots & \cdots & \vdots \cr 1 & x_n & x_n^2 & \cdots & x_i^{n}}\right]$

Note that the decomposition doesn't involve any summation and can be easily computed using simple indexing and power operations. The matrix multiplication is an effective operator to compute the final matrix $\mathbf{A}$. This is both efficient and easy to implement.

The matrix on the right is simply a Vandermonde matrix with an added column (right-most). Here is a MATLAB 2-liner to implement this:

% x is the input column vector, n is the number of elements e.g. n=length(x)
W = [ fliplr(vander(x)) (x.^n)];
A = W'*W;


This should work faster than for-looping in MATLAB since vander function doesn't contain any loops. It is also definitely shorter and cleaner.

I believe that making a quick remark is important here. Unlike common belief, when properly implemented matrix multiplication is not $O(N^3)$, but $<O(N^{2.4})$. Many available packages are already very optimized. Using vector instruction sets, along with parallelization brings further speed-up. However, the theoretical complexity never comes to $O(N^2)$ which is achieved via for-looping. This approach is only neater and faster in MATLAB for small sized matrices (up to ~ N<1500).

Also note that, even if the matrix notation is used here, the multiplication could still benefit from the special structure of these matrices. If one is to utilize this property, similar complexity can be obtained.

• Depending on what OP needs to do with this matrix, they may be better off just leaving it in this form. I.e. if the need to form matrix-vector products, they're better off not forming $A$ at all. – Bill Barth Aug 13 '15 at 12:46
• Yes, I understand that I can do it using Vandermonde matrix and then using Gauss elimination, or even simplier inverse matrices method, but my intention is to know how to fill this specified array I have written in the post or to know whether it is not efficient and the Vandermonde one is better to use. – Beginner in fort Aug 13 '15 at 15:36
• This matrix multiplication takes $O(n^3)$ time, but the matrix as originally written takes $O(n^2)$ to evaluate (there are $2n$ distinct components, each taking $O(n)$ time). Isn't this form more expensive? – Kirill Aug 13 '15 at 17:05
• It might be more efficient for certain software (e.g. MATLAB), which requires for-less implementations. Both matrices can be formed in a vectorized fashion and multiplication is done through BLAS, which is already optimized. – Tolga Birdal Aug 13 '15 at 18:22
• Okay, I have run the tests on 2 other computers, including a Macbook Air. It seems that before n~1500, vectorized code outperforms. They seem to equate around that point. Then for-loops start to beat. My $x$ is a uniformly distributed random vector. – Tolga Birdal Aug 15 '15 at 11:25