Ausgleichsrechnung (Least-Squares)

Computer-Algebra und Visualisierungen der Ausgleichsrechnung

Prof. Dr. Bernhard Zgraggen, HSR

Einführung in dieses Notebook

Dieses Dokument ist ein "selbsterklärendes" Mathematica-Notebook, das interaktiv bearbeitet werden kann ("Hands-On"). Input-Zellen können mittels der Tastenkombination ⇧⌤ (in Worten: Shift-Enter) ausgeführt werden. Falls Textzellen geschlossen sind, so können sie durch Klick auf das Icon ▼ oder auf die blaue Zellklammer rechterhand geöffnet werden.

Wenn Sie in Mathematica auf Befehle oder Begriffe stossen, die Sie nicht verstehen, so steht Ihnen die mächtige Hilfe/Doc Center zur Verfügung; diese ist interaktiv.

Allfällige Aufgaben/Übungen befinden sich am Ende des Dokuments.

Linear Least-Squares: Diskrete Daten → {Fit, LeastSquares, LinearModelFit}

Symbolische Daten, Designmatrix und Normalgleichungen → {DesignMatrix}

Definition of symbolic data in a list of variable length.

datasymb[n_]:=Table[{x[i],y[i]},{i,0,n}]

datasymb[4]//TableForm

datasymb[4]

Definition of Newton polynomials (as one possible set of basis functions).

newtonpoly[k_,var_]:=Product[var-x[i],{i,0,k-1}]

Table[newtonpoly[k,x],{k,0,4}]//TableForm

The design matrix has entries (Subscript[g, j])(x[i]) in row i and column j. It depends of the data and the basis functions (Subscript[g, j]).

The next commands generate an abstract design matrix and one for the basis of Newton polynomials. The latter has a triangular structure.

n=4; (* number of basis functions minus 1 *)
nn=8; (* number of data minus 1: the symbol N is occupied by Mathematica *)
DesignMatrix[datasymb[nn],Table[Subscript["g",k-1][t],{k,1,n+1}],t,IncludeConstantBasis->False];
%//MatrixForm

n=2; (* number of basis functions minus 1 *)
nn=5; (* number of data minus 1: the symbol N is occupied by Mathematica *)
DesignMatrix[datasymb[nn],Table[newtonpoly[k,t],{k,0,n}],t];
%//MatrixForm

The normal equations

   G denotes the design matrix.
    <(Subscript[g, j]),(Subscript[g, k])>  denotes the canonical inner product (Underoverscript[∑ , (i=0), N])(Subscript[g, j])(x[i])(Subscript[g, k])(x[i]) of basis functions.
    
    
                                                                              (G^tr)·G a|  ==   (G^tr)y|
                                                                              
    
    (({{<(Subscript[g, 0]),(Subscript[g, 0])>, <(Subscript[g, 1]),(Subscript[g, 0])>, ⋯, <(Subscript[g, n]),(Subscript[g, 0])>}, {<(Subscript[g, 0]),(Subscript[g, 1])>, <(Subscript[g, 1]),(Subscript[g, 1])>, ⋯, <(Subscript[g, n]),(Subscript[g, 1])>}, {⋮, ⋮, ⋱, ⋮}, {<(Subscript[g, 0]),(Subscript[g, n])>, <(Subscript[g, 1]),(Subscript[g, n])>, …, <(Subscript[g, n]),(Subscript[g, n])>}}) ) (({{(Subscript[a, 0])}, {(Subscript[a, 1])}, {⋮}, {(Subscript[a, n])}}))== (({{<y,(Subscript[g, 0])>}, {<y,(Subscript[g, 1])>}, {⋮}, {<y,(Subscript[g, n])>}}) )
    
    
    In many cases the square matrix (G^tr)·G is ill-conditioned. Using matrix decompositions like CholeskyDecomposition or QRDecomposition helps here. 
    Another powerful approach  is the use of orthogonal basis functions, such that  <(Subscript[g, j]), (Subscript[g, k])> == 0 whenever j ≠ k. This means that the matrix (G^tr)·G is diagonal.

Diskrete Orthogonalisierung → {Orthogonalize}

The Mathematica command Orthogonalize allows orthogonalization incl. normalization (by default with the Gram-Schmidt procedure) with respect to an inner product. Here it is used for the "discrete" inner product defined in the section above. The output formulas normally are very large (due to normalization) and symbolic Simplify consumes much time and memory.

n=2; (* number of basis functions minus 1 *)
nn=4; (* number of data minus 1: the symbol N is occupied by Mathematica *)
orthonewton1=Orthogonalize[Table[newtonpoly[k,t],{k,0,n}],Sum[(#1/.{t->datasymb[nn][[i]][[1]]})* ( #2/.{t->datasymb[nn][[i]][[1]]}),{i,1,nn+1}]&];

(* orthonewton1 = Simplify[orthonewton1]; *)

The matrix (G^tr)·G now becomes diagonal (G is the DesignMatrix with respect to the orthogonalized basis).

To avoid normalization a version of the Gram-Schmidt procedure only using Projections is defined:

myGramSchmidt[vecs_,ip___]:=Module[{ovecs=vecs},
Do[ovecs[[i]]-=Projection[ovecs[[i]],ovecs[[j]],ip],{i,2,Length[vecs]},{j,1,i-1}];
ovecs]

This modified orthogonalization is used as before but Simplify does not consume too much time etc:

n=2; (* number of basis functions minus 1 *)
nn=4; (* number of data minus 1: the symbol N is occupied by Mathematica *)
orthonewton2=myGramSchmidt[Table[newtonpoly[k,t],{k,0,n}],Sum[(#1/.{t->datasymb[nn][[i]][[1]]})* ( #2/.{t->datasymb[nn][[i]][[1]]}),{i,1,nn+1}]&];

orthonewton1/.{x[k_]:>k}
orthonewton2/.{x[k_]:>k}

The matrix (G^tr)·G now becomes diagonal (G is the DesignMatrix with respect to the orthogonalized basis).

Table[If[j≠k, Simplify@Sum[(orthonewton2[[j]]/.{t->datasymb[nn][[i]][[1]]})* ( orthonewton2[[k]]/.{t->datasymb[nn][[i]][[1]]}),{i,1,nn+1}],
"*"],{j,1,n+1},{k,1,n+1}
]//MatrixForm

The orthogonalization functions may be used for continous data, too. Discrete summation is replaced by integration. Here are two important examples (Legendre and Chebyshev polynomials up to normalization).

n=5;
Orthogonalize[Table[x^k,{k,0,n}],Integrate[#1*  #2,{x,-1,1}]&]
myGramSchmidt[Table[x^k,{k,0,n}],Integrate[#1*  #2,{x,-1,1}]&]

n=5;
Orthogonalize[Table[x^k,{k,0,n}],Integrate[((#1*  #2)/(Sqrt[1-(x^2)])),{x,-1,1}]&]
myGramSchmidt[Table[x^k,{k,0,n}],Integrate[((#1*  #2)/(Sqrt[1-(x^2)])),{x,-1,1}]&]

Least-Squares: Kontinuierliche Daten → {Orthogonalize}

Nichtlinearer Least-Squares → {FindFit, NonlinearModelFit}

Aufgabe(n)

(1)
Für ganzzahlige x-Daten {0, 1, 2, ..., N} lassen sich orthogonale Polynome durch die Formel

(Subscript[P, m,N])(t)=(∑_(i=0)^m)(((-1))^i)(({{m}, {i}}))(({{m+i}, {i}}))((t^((i)))/(N^((i))))       wobei    (t^((i)))= t(t-1)(t-2)⋯(t-i+1)

explizit darstellen.

Überprüfen Sie, dass die Polynome (Subscript[P, m,N])(t)für verschiedene m orthogonal sind.

Hinweis: Die Polynome können wie folgt in Mathematica codiert werden:

p[m_,nn_,t_]:=(∑_(i=0)^m)(((-1))^i)Binomial[m,i]Binomial[m+i,i]((FactorialPower[t,i])/(FactorialPower[nn,i]))

Lösung : Die Normalmatrix besteht aus allen Skalarprodukten von Spalten der Designmatrix

n=4;
normalmatrix=Table[(∑_(t=0)^n)p[m1,n,t]*p[m2,n,t],{m1,0,n},{m2,0,n}];

normalmatrix//MatrixForm

Experimental (Example 1.4 Script)

{p[0,4,t],p[1,4,t],p[2,4,t]}//FunctionExpand

{p[0,4,x-3],p[1,4,x-3],p[2,4,x-3]}//FunctionExpand

Table[p[i,9,t]//FunctionExpand,{i,0,9}]//Column

ip[g1_,g2_]:=Sum[(g1/.{t->i})*(g2/.{t->i}),{i,0,4}]
ip1[g1_,y_]:=Sum[(g1/.{t->i})*(y[[i+1]]),{i,0,4}]

MatrixForm@Table[ip[p[i,4,x-3],p[j,4,x-3]],{i,0,2},{j,0,2}]

y={1.7,2,2.26,2.42,2.7}

Table[ip1[p[i,4,t],y],{i,0,2}]

Table[ip1[p[i,4,t],y]/ip[p[i,4,t],p[i,4,t]],{i,0,2}]

%.{p[0,4,x-3],p[1,4,x-3],p[2,4,x-3]}//FunctionExpand//ExpandAll