1
$\begingroup$

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.

$\endgroup$
6
  • 3
    $\begingroup$ 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. $\endgroup$ Aug 13, 2015 at 11:23
  • 2
    $\begingroup$ Is $A(2,1)$ wrong? It should be $\sum{x_i}$, instead of $\sum{x^n_i}$, right ? $\endgroup$ Aug 13, 2015 at 11:56
  • $\begingroup$ Right, my bad. Corrected. $\endgroup$ Aug 13, 2015 at 15:27
  • $\begingroup$ @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. $\endgroup$ Aug 13, 2015 at 15:32
  • 2
    $\begingroup$ 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 $\endgroup$ Aug 14, 2015 at 1:55

2 Answers 2

4
$\begingroup$

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)$.

$\endgroup$
8
$\begingroup$

$\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.

$\endgroup$
9
  • 2
    $\begingroup$ 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. $\endgroup$
    – Bill Barth
    Aug 13, 2015 at 12:46
  • $\begingroup$ 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. $\endgroup$ Aug 13, 2015 at 15:36
  • 2
    $\begingroup$ 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? $\endgroup$
    – Kirill
    Aug 13, 2015 at 17:05
  • 3
    $\begingroup$ 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. $\endgroup$ Aug 13, 2015 at 18:22
  • 3
    $\begingroup$ 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. $\endgroup$ Aug 15, 2015 at 11:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.