Lower time complexity of 10 for loop code

I am attempting to approximate 52 data points to the following equation:

F(x,y,z)= C_1*x^2 + C_2*x + C_3*y^2+C_4*y+C_5*z^2+C_6*z+ C_7*x*y+C_8*x*z+C_9*y*z+C_10

Since I know the values for x,y,z, and f(x,y,z) I decided to create 10 nested for loops. Each for loop iterating variable represents an index for a point which I use to set up a linear system. (I.e. A(1,:) = [x(i)^2, x(i), y(i)^2, y(i), z(i)^2,...]) I do this for all ten rows of A. The ten for loops will allow me to keep points constant and eventually get every combination. I then solve C = A\f. And then determine the average error for this solution by plugging in the value for C and the x,y,z into the equation and determining the difference between it and the actual value for f. Eventually, with every combination of points being tested, I hope to find a minimal average error. Unfortunately, while I expect this loop to work eventually it is very time inefficient. Does anyone have any advice to minimize the amount of time it takes? Note that I have code to prevent repeated indices being used for one system of equations. Also note that the order of the indices is irrelevant. For example indices[1 2 3 4 5 6 7 8 9 10] is exactly the same as [10 9 8 7 6 5 1 2 3 4].

Answer 1

[Walter Roberson]

X = x(:); Y = y(:); Z = z(:);
A = [X.^2, X, Y.^2, Y, Z.^2, Z, X.*Y, X.*Z, Y.*Z, ones(length(X),1)];
F = f(:);
C = A\F;

No loop needed.

[Thompson Bliss]

Unfortunately this code does not work. C must be a vector of length 10. This code produces a 52 x 520 matrix for C.

[Walter Roberson]

No it doesn't.

>> A = rand(52,10);
>> F = rand(52,1);
>> C = A\F;
>> size(C)

ans =

      10     1

Look at the documentation:

http://www.mathworks.com/help/matlab/ref/mldivide.html#output_argument_x

"Solution, returned as a vector, full matrix, or sparse matrix. If A is an m-by-n matrix and B is an m-by-p matrix, then x is an n-by-p matrix, including the case when p==1."

I used F = f(:); so F is going to be a something-by-1 matrix, so p must be 1 in the above notation. According to the documentation, the output must then be somevalue-by-1 . Your claimed output of 52 x 520 is not possible with the code I show, no matter how many elements were in x.

In order to get a 52 x 520 output, your A would have had to be something-by-52 and your B would have had to have been something-by-520. Your B is derived from your f values, so you would have had to have had something-by-520 f values. That disagrees with your commentary that you have 52 points to approximate.

Here is a more complete example of the code for interpolating 52 points:

x = rand(1,52); y = rand(1,52); z = rand(1,52); f = rand(1,52);
X = x(:); Y = y(:); Z = z(:);
A = [X.^2, X, Y.^2, Y, Z.^2, Z, X.*Y, X.*Z, Y.*Z, ones(length(X),1)];
F = f(:);
C = A\F;

>> size(C)

ans =

      10     1

Notice the complete lack of any loop to build anything, let alone a matrix 520 wide.