Least-Squares von 3D-Daten

Numerische Berechnungen von Least-Squares Approximationen 
mittels DesignMatrix und SingularValueDecomposition

Prof. Dr. Bernhard Zgraggen, HSR

Import Data and Plots 3D → {ListPointPlot3D, ListPlot}

Import of "Table" data

data=Import[ToFileName[NotebookDirectory[],"pointsall.txt"],"Table"];
Length[data]

Scatter plot 3D:

ListPointPlot3D[data,PlotStyle->{Directive[Red,PointSize[Tiny]]}]

The 3D data points are functional (no different z-values for the same pair (x,y)). This is checked by the following code.

Select[Tally[data,({#1[[1]],#1[[2]]}=={#2[[1]],#2[[2]]})&],#[[2]]>1&]

The projection of the 3D-points onto the x-y-plane yields a irregular grid. The points here consist only of the first two coordinates of the 3D-points.

ListPlot[({#[[1]],#[[2]]}&)/@data,PlotStyle->{Darker[Red],PointSize[Small]}]

The ListCurvePathPlot command tries to reconstruct curves from the x-y-coordinates.

ListCurvePathPlot[({#[[1]],#[[2]]}&)/@data]

Some computations of ranges of data values (min - max):

{Min[#[[1]]&/@data],Max[#[[1]]&/@data]}
{Min[#[[2]]&/@data],Max[#[[2]]&/@data]}
{Min[#[[3]]&/@data],Max[#[[3]]&/@data]}

Least-Squares Approximations → {DesignMatrix, Fit, SingularValueDecomposition}

Firstly, a set of basis functions of a certain degree is generated (complete "ansatz"):

Clear[x,y]
degree=16;
basis=Flatten@Join[Table[List@@Expand[(x+y)^i],{i,0,degree}]];
Length[basis]

basis1=Flatten@Join[Table[Table[x^k y^(i-k),{k,0,i}],{i,0,degree}]];

basis1=basis1/.{x->(x-Mean[#[[1]]&/@data])/StandardDeviation[#[[1]]&/@data], y->(y-Mean[#[[2]]&/@data])/StandardDeviation[#[[2]]&/@data]};

Mathematica can do least-squares approximations internally:

f=Fit[data,basis,{x,y}]//Simplify;

Design matrix (rather large here):

dm=DesignMatrix[data,basis1,{x,y}, IncludeConstantBasis->False];
Dimensions[dm]

Design matrix (rather large here):

{u,d,v}=SingularValueDecomposition[dm,Length[basis1],Tolerance ->0];

Checking the diagonal part of matrix d (as many rows as there are basis functions) and the dimensions of u, d, v:

{Dimensions[u],Dimensions[d],Dimensions[v]}

Take[d,Length[basis]]//MatrixForm
Take[v,Length[basis]]//MatrixForm

Combinations of data plots with plots of least-squares approximation:

z=(#[[3]])&/@data;  (* 3rd coordinate of data = height measurement *)

as=v.((Transpose[u].z)/Diagonal[d])

fs=as.basis1//Expand;

Show[ListPointPlot3D[data,PlotStyle->{Directive[Red,PointSize[Small]]}],Plot3D[fs,{x,0,202.335},{y,0,148.453},PlotRange->{-10-31.10029`,10+0.198112`},PlotStyle->{Opacity[0.33]}]]

Show[ListPointPlot3D[data,PlotStyle->{Directive[Red,PointSize[Small]]}],Plot3D[f,{x,0,202.335},{y,0,148.453},PlotRange->{-10-31.10029`,10+0.198112`},PlotStyle->{Opacity[0.33]}]]

Computing the sum of squared residuals ( = least-squares error) and maximal error of approximation:

s=Norm[Transpose[u].z]  (* sum of squared residuals = squared length of transformed z-vector *)

Max[Abs[#[[3]]-fs/.{x->#[[1]],y->#[[2]]}]&/@data]

Scatter plot for the least-squares approximation and combination with scatter plot of data:

ListPointPlot3D[{#[[1]],#[[2]],fs/.{x->#[[1]],y->#[[2]]}}&/@data]

Show[ListPointPlot3D[data,PlotStyle->{Directive[Red,PointSize[Small]]}],ListPointPlot3D[{#[[1]],#[[2]],fs/.{x->#[[1]],y->#[[2]]}}&/@data,PlotRange->{-31.10029`,0.198112`}]]