Theoretische und methodische Grundlagen des Visual Computing – Approximation

Lineare Least-Squares-Approximation

Least Squares in 1D: Gegeben sind \(x_1, \dotsc , x_N \in \real \), \(f_1, \dotsc , f_N \in \real \) und \(m \in \natural \).
Gesucht ist eine Funktion \(g \in \PP _m(\real )\) mit \(F := \sum _{i=1}^N (f_i - g(x_i))^2 \to \min \).
Sei zunächst \(\{p_k\}_{k=0}^m\) eine Basis von \(\PP _m(\real )\) (z. B. Monome). Dann erhält man mit
\(g(x) = \sum _{k=0}^m c_k p_k(x)\), dass \(0 \overset {!}{=} \partial _{c_j} F = -\sum _{i=1}^N 2(f_i - \sum _{k=0}^m c_k p_k(x_i)) p_j(x_i)\) für \(j = 0, \dotsc , m\) erfüllt sein muss. Die hinreichende Bedingung für ein Minimum, dass die Hesse-Matrix
\(F’’ = (\partial _{c_i} \partial _{c_j} F)_{i,j} = (\sum _{k=1}^N p_i(x_k) p_j(x_k))_{i,j}\) positiv definit ist, ist erfüllt, weil die \(p_i\) l. u. sind.
Damit bekommt man \(A (c_0, \dotsc , c_m)^\tp = ((f, p_0), \dotsc , (f, p_m))^\tp \) mit \(A := ((p_i, p_j))_{i,j=0}^m\) und \((p_i, p_j) := \sum _{k=1}^N p_i(x_k) p_j(x_k)\).

lineare Least-Squares-Approximation:
Gegeben seien \(C \in \real ^{N \times k}\) mit \(N > k\) und vollem Spaltenrang \(k\) sowie \(\vec {d} \in \real ^N\).
Gesucht ist \(\vec {y} \in \real ^k\) mit \(|\vec {r}|^2 \to \min \) und Residuum \(\vec {r} := C\vec {y} - \vec {d}\) (da \(C\vec {y} = \vec {d}\) überbest.).
Es gilt \(|\vec {r}|^2 = \vec {r}^\tp \vec {r} = \vec {y}^\tp A \vec {y} - 2\vec {b}^\tp \vec {y} + \vec {d}^\tp \vec {d}\) mit \(A := C^\tp C \in \real ^{k \times k}\) und \(\vec {b} = C^\tp \vec {d} \in \real ^k\).
\(A\) ist positiv definit, weil \(C\) vollen Spaltenrang besitzt. Durch \(\forall _{i=1,\dotsc ,k}\; \partial _{y_i} \vec {r}^\tp \vec {r} \overset {!}{=} 0\) erhält man die Normalengleichungen \(A\vec {y} = \vec {b}\). Damit erhält man \(\vec {y} = A^{-1} \vec {b} = (C^\tp C)^{-1} C^\tp \vec {d}\) (Hesse-Matrix ist \(H = 2A\), d. h. positiv definit).
Ein Nachteil der Normalengleichungen ist, dass die Konditionszahl groß sein kann.

Zusammenhang mit Regression: Seien nun \(\vecs {x}{1}, \dotsc , \vecs {x}{N} \in \real ^d\) und \(f_1, \dotsc , f_N \in \real \) gegeben. Gesucht ist \(g \in \PP _m(\real ^d)\) wie vorhin. Ist \(\{p_k\}_{k=1}^M\) eine Basis von \(\PP _m(\real ^d)\) (mit \(M := \binom {d+m}{m}\)), so erhält man das für \(N > M\) überbestimmte LGS \(\forall _{i=1,\dotsc ,N}\; g(\vecs {x}{i}) = \vec {p}(\vecs {x}{i})^\tp \vec {c} = f_i\) mit \(\vec {p}(\vec {x}) := (p_1(\vec {x}), \dotsc , p_M(\vec {x}))^\tp \). Damit erhält man die Normalengleichungen \(A\vec {c} = \vec {b}\) mit
\(A := \sum _{i=1}^N \vec {p}(\vecs {x}{i}) \vec {p}(\vecs {x}{i})^\tp \) und \(\vec {b} := \sum _{i=1}^N f_i \vec {p}(\vecs {x}{i})\) (\(A = C^\tp C\) mit \(C_{i,j} = p_j(\vecs {x}{i})\)).

Weighted Least Squares (WLS)

Weighted Least Squares (WLS): Gegeben seien \(\vecs {x}{1}, \dotsc , \vecs {x}{N} \in \real ^d\) und \(f_1, \dotsc , f_N \in \real \).
Die WLS-Methode sucht nun für \(\vec {\xi } \in \real ^d\) fest ein Polynom \(g \in \PP _m(\real ^d)\) mit
\(F := \sum _{i=1}^N \theta (|\vec {\xi } - \vecs {x}{i}|) (f_i - g(x_i))^2 \to \min \) mit einer Gewichtsfunktion \(\theta \colon [0, \infty ) \to \real \).
Mit dem Ansatz \(g(\vec {x}) = \vec {p}(\vec {x})^\tp \vec {c}\) und \(\vec {\nabla } F(\vec {c}) = \vec {0}\) erhält man das LGS
\(\sum _{i=1}^N \theta (r_i) \vec {p}(\vecs {x}{i}) \vec {p}(\vecs {x}{i})^\tp \cdot \vec {c} = \sum _{i=1}^N \theta (r_i) f_i \vec {p}(\vecs {x}{i})\) mit \(r_i := |\vec {\xi } - \vecs {x}{i}|\).
Gebräuchliche Gewichtsfunktionen sind \(\theta (r) := e^{-r^2/h^2}\) für \(h > 0\) (Gauß),
\(\theta (r) := (1 - \frac {r}{h})^4 (\frac {4r}{h} + 1)\) für \(r \in [0, h]\) und \(h > 0\) (Wendland) und \(\theta (r) := \frac {1}{r^2 + \varepsilon ^2}\) für \(\varepsilon > 0\).
Wegen numerischen Instabilitäten kann es hilfreich sein, \(\vec {\xi }\) in den Ursprung zu verschieben.

WLS in 1D: Seien \(x_1, \dotsc , x_N \in \real \), \(f_1, \dotsc , f_N \in \real \) und \(\xi \in \real \).

linearer Ansatz: \(g(x) := \vec {p}(x)^\tp \vec {c}\) mit \(\vec {p}(x) := (1, x)^\tp \) führt zu \(\vec {c} = A^{-1} \vec {b}\) mit
\(A := \sum _{i=1}^N \theta (|\xi - x_i|) \smallpmatrix {1&x_i\\x_i&x_i^2}\) und \(\vec {b} := \sum _{i=1}^N \theta (|\xi - x_i|) f_i \smallpmatrix {1\\x_i}\)
quadratischer Ansatz: \(g(x) := \vec {p}(x)^\tp \vec {c}\) mit \(\vec {p}(x) := (1, x, x^2)^\tp \) führt zu \(\vec {c} = A^{-1} \vec {b}\) mit
\(A := \sum _{i=1}^N \theta (|\xi - x_i|) \smallpmatrix {1&x_i&x_i^2\\x_i&x_i^2&x_i^3\\x_i^2&x_i^3&x_i^4}\) und \(\vec {b} := \sum _{i=1}^N \theta (|\xi - x_i|) f_i \smallpmatrix {1\\x_i\\x_i^2}\)

Moving Least Squares (MLS)

Moving Least Squares (MLS): Gegeben seien \(\vecs {x}{1}, \dotsc , \vecs {x}{N} \in \real ^d\) und \(f_1, \dotsc , f_N \in \real \).
Die MLS-Methode sucht nun für jedes \(\vec {x} \in \real ^d\) ein Polynom \(g_{\vec {x}} \in \PP _m(\real ^d)\) mit
\(F_{\vec {x}} := \sum _{i=1}^N \theta (|\vec {x} - \vecs {x}{i}|) (f_i - g(x_i))^2 \to \min \) mit einer Gewichtsfunktion \(\theta \colon [0, \infty ) \to \real \).
Mit dem Ansatz \(g_{\vec {x}}(\vec {y}) = \vec {p}(\vec {y})^\tp \vecs {c}{\vec {x}}\) und \(\vec {\nabla } F_{\vec {x}}(\vecs {c}{\vec {x}}) = \vec {0}\) erhält man das LGS
\(\sum _{i=1}^N \theta (r_i) \vec {p}(\vecs {x}{i}) \vec {p}(\vecs {x}{i})^\tp \cdot \vecs {c}{\vec {x}} = \sum _{i=1}^N \theta (r_i) f_i \vec {p}(\vecs {x}{i})\) mit \(r_i := |\vec {x} - \vecs {x}{i}|\).
Die Approximation ist dann gegeben durch \(h\colon \real ^d \to \real \), \(h(\vec {x}) := g_{\vec {x}}(\vec {x})\), d. h. für jeden Auswertungspunkt \(\vec {x}\) muss ein WLS-Problem gelöst werden.

Bilddeformation mit MLS:
Seien \(\Omega := [0, 1]^2\) und \(\vecs {p}{i}, \vecs {q}{i} \in \Omega \) für \(i = 1, \dotsc , N\) Kontrollpunkte.
Gesucht ist nun eine Deformationsfunktion \(\vec {f}\) mit \(\vec {f}(\vecs {p}{i}) = \vecs {q}{i}\) für \(i = 1, \dotsc , N\) (Interpolation), \(\vec {f}\) glatt und \([\forall _{i=1,\dotsc ,N}\; \vecs {p}{i} = \vecs {q}{i}] \implies \vec {f} = \id _\Omega \).
Dazu benutzt man MLS mit affinen Transformationen, d. h. man setzt \(\vec {f}(\vec {x}) = \vecs {g}{\vec {x}}(\vec {x})\) an, wobei \(\vecs {g}{\vec {x}}(\vec {y}) := M_{\vec {x}} \vec {y} + \vecs {t}{\vec {x}}\), \(F_{\vec {x}} := \sum _{i=1}^N w_i |\vecs {g}{\vec {x}}(\vecs {p}{i}) - \vecs {q}{i}|^2 \to \min \) und \(w_i := \frac {1}{|\vecs {p}{i} - \vec {x}|^{2\alpha }}\).
Durch \(\vecs {\nabla }{\vec {t}} F \overset {!}{=} \vec {0}\) erhält man \(\vec {t} = \vecs {q}{\ast } - M\vecs {p}{\ast }\) mit \(\vecs {p}{\ast } := \frac {\sum _{i=1}^N w_i\vecs {p}{i}}{\sum _{i=1}^N w_i}\) und \(\vecs {q}{\ast } := \frac {\sum _{i=1}^N w_i\vecs {q}{i}}{\sum _{i=1}^N w_i}\).
Ist \(M = (m_{i,j})_{i,j=1}^2\), so bekommt man mit \(\partial _{m_{i,j}} F \overset {!}{=} 0\), \(\vecs {\varrho }{i} = \vecs {p}{i} - \vecs {p}{\ast }\) und \(\vecs {\sigma }{i} = \vecs {q}{i} - \vecs {q}{\ast }\) die Bedingung \(MR = S\) mit \(R := \sum _{i=1}^N w_i \vecs {\varrho }{i} \vecs {\varrho }{i}^\tp \) und \(S := \sum _{i=1}^N w_i \vecs {\varrho }{i} \vecs {\sigma }{i}^\tp \), d. h. \(M = SR^{-1}\).
Somit bekommt man \(\vecs {g}{\vec {x}}(\vec {y}) = M(\vec {y} - \vecs {p}{\ast }) + \vecs {q}{\ast }\).