I want to fit a polynomial surface to some 3d points. First in a loop, matrix R is calculated as follows:
for j=1:size(indices,1)
A(j,1:3)=Temp.Location(indices(j),:);
X = Temp.Location(indices(j),1);
Y = Temp.Location(indices(j),2);
R(j, 1) = X*X;
R(j, 2) = Y*Y;
R(j, 3) = Y;
R(j, 4) = X;
R(j, 5) = X*Y;
R(j, 6) = 1;
Z(j, 1) = Temp.Location(indices(j),3);
end
And then coeficients are calculated as follows:
Cof=inv(R'*R)*R'*Z;
The results (z and predicted z) are as follows:
enter image description here
Is the method wrong?
How least square is used as surface fitiing?
Thanks
points: [ x ,y ,z ,predicted z]
[
]4
NOT AN ANSWER.
This isn't an answer to your question, but it's a cleanup of your code, because ... because it was needed.
You wrote:
In doing this, you're re-allocating the matrix $R$ each time through the loop. No big deal for 40 points, but a total mess when there are 100000 points. So let's just fix that a bit by allocating $R$ in advance:
Now let's get more idomatic: when you want to extract items 1, 2, 3 of an array, you can write
or you can just write
U([1 2 3]) = X([1 2 3])
which is the more "matlab-like" way to do things. Since you're taking each of the items in
Temp.Location, in the order specified byindices, you can rewrite most of your code like this:That's looking better, but it's still only partway there. I'm going to make a vector the same length as
Xcontaining only the number $1$, and call this vectorU, for "unit":And now I'm going to stick $X, Y,$ and $U$ all into an array $M$, as its columns. I've gotten rid of $A$ because you never used it.
And finally, I'm going to look at
M'M. The first entry of this will be the row-vector $X'$ multiplied by the column vector $X$, i.e., it'll be $$ x_1^2 + x_2 ^2 + \ldots + x_n^2. $$ Now that might seem odd, because your code just computes $x_1^2$ and puts it into $R(j, 1)$, and similarly for the others. But if you're doing least-squares fitting right, you need this sum, not the individual values. Your code, later, perhaps sumsR(j, ...)over allj--- it almost certainly should do so. So now my (conjectured) version of the right code looks like this:At this point
Ris a symmetric matrix containing all the values you need to complete the least-squares process. Of course, I may have mangled your intentions somewhat --- just looking at your code fragment, it's not clear where you're going with the matrixRthat you built --- but I've written far more idiomatic and readable matlab code, which will run hugely faster if you ever get around to processing large amounts of data.Finally, now that we compute
Ronly once, that first line, where we allocateRas an array full of zeroes ... that can be removed!N.B.:
In your coefficient computation, you're computing
inv(R'*R) * R' * Z. That's not going to get you the coefficient you need, because you're squaring a bunch of things that you've already squared. It's also computationally a terrible idea. We seldom solve $Ax = b$ by writingInstead, you get a huge efficiency improvement if you write
You might want to read the documentation for the backslash operator; it could substantially simplify your program.
More important though: I want to suggest that you go back and re-read whatever presentation of least-squares fitting you've seen, until you understand it deeply. Then, and only then, start writing code. Randomly searching through the space of all possible programs for the one that does what you want is almost never a good strategy, even if the random search is clustered "near" some correct code: the search space is simply too large to ever result in a correct answer.
If you think you've already done that deep understanding perhaps you can write out here the mathematics that you're trying to emulate with your program.