Best fit of a surface

Dear all, i have a matrix A(188x215) but is has many pics as you can see in the image

i would like to visualize a plane that best fits this matrix.

So far i made it with some for cycles but i do think that there can be a better way to solve that.

Any help would be appreciated

Answer 1

[arich82]

If you're just looking for a best fit plane (for some definition of 'best'), you can use the backslash operator, just like finding the the least-squares fit for a line.

For a line, you have data x and y, but unknown coeffcients m abnd b:

y1 = m*x1 + b;
y2 = m*x2 + b;
y3 = m*x3 + b;
...

This represents an overdetermined system X*M = Y, where we can solve M = [b; m]

    X = [ ...
          1, x1; ...
          1, x2; ...
          1, x3; ...
        ];

    Y = [y1; y2; y3;];

    M = X\Y;

For a plane, we have the general equation

0 = a*(x - x0) + b*(y - y0) + c*(z - z0)
z = (1/c)*(a*x0 + b*y0 + c*z0) - ((a/c)*x + *b/c)*y);

Let C = (1/c)*(a*x0 + b*y0 + c*z0), A = -a/c, B = -b/c, then for M = [C, A, B], you again have the overdetermined equation Z = X*M, which can be solved via backslash:

X = [ ...
      1, x1, y1; ...
      1, x2, y2; ...
      1, x3, y3; ...
    ];

Z = [z1; z2; z3;];

M = X\Z;

In general, if you have vectors of data x(:), y(:), z(:) of length n, then you can create the data variables as X = [ones(n, 1), x(:), y(:)]. Below is a full example:

m = 188; 
n = 215; 
[y, x] = meshgrid(1:n, 1:m);
z = sind(10*x + 5*y) + log(x.*y); % filler data; this would be your A matrix
X = [ones(m*n, 1), x(:), y(:)];
M = X\z(:);

figure;
surf(x, y, z, 'EdgeColor', 'none'); % original wavy data
hold all;
surf(x, y, reshape(X*M, m, n), 'EdgeColor', 'none', 'FaceColor', 'b'); % blue fitted-plane

Alternatively, if you just want to get rid of the peaks, you might consider filtering your data (e.g. filter2).

Hope this helps.

--Andy

[Matt J]

The numerical analysts on here may contradict me, but I have seen it said that solving by QR decomp is more stable for polynomial fitting than direct use of backslash. This is what is done inside polyfit() at any rate.

[Michele]

Thank you! Also for explaining it so clearly and making also undertand why it works like this!

[arich82]

Matt J--

Do you have a reference thread for the QR stability comment? Or a pathological case where mldivide would fail in fitting a simple plane? I'm genuinely interested.

Once upon a time, the Matlab help documentation described mldivide as using QR for non-square matrices, though perhaps there was some option for a more stable decomp that wasn't used by default. (It now appears there's just a brief comment under "Systems of Linear Equations" about mldivide and qr; it would be nice MW didn't gut their help docs...) The mldivide doc no longer offers any real info, and the new example in qr shows a "Q-less" solution to a system of linear equations, but doesn't offer any insights on what might be different with the backslash operator.

[Matt J]

Once upon a time, the Matlab help documentation described mldivide as using QR for non-square matrices,

It still does,

http://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#brs10rz-1

so I guess it doesn't matter if all you want is the fit. If you wanted an error analysis of the fit, like polyfit provides for example, it helps to do the QR decomp explicitly and reuse it for the error analysis...

Answer 2

[Matt J]

See this recent, similar thread on multivariate polynomial fitting

http://www.mathworks.com/matlabcentral/answers/113898-polynomial-fitting-with-multiple-independent-variables

[Michele]

Thank you for the link! I had a look and there are many interersting infos!