Logo der Physikalisch-Technischen Bundesanstalt

Numerische Methoden

Arbeitsgruppe 8.43

Software Toolbox für rechenintensive Probleme

PyThia ist ein open-source Software Paket, dass entwickelt wurde um Polynomchaos-Surrogate für hochdimensionale Funktionen auf eine nicht-intrusive Weise zu erzeugen. Die Software wurde im Rahmen des ZIM Projektes ZF4014017RR7 zur Bestimmung von Unsicherheiten von Formparametern für Nanostrukturen, wie sie etwa bei Photolithographiemasken auftreten, entwickelt. Allerdings erlaubt die Allgemeinheit der zugrundeliegenden mathematischen Konzepte PyThia auch ohne Anpassungen auf andere Problemstellungen anzuwenden. Weitere Informationen zu dem PyThia Software Paket stehen auf der offiziellen Homepage zur Verfügung.

Zugang und Installation

PyThia ist öffentlich zugänglich und kann frei verwendet werden. Das Git Repository beschreibt außerdem wie man zu der Weiterentwicklung der Software beitragen kann. Dies ist beispielsweise dann praktisch, wenn neue Funktionen für die eigene Anwendung benötigt werden.

Mit einem Internetzugang kann PyThia einfach via pip installiert werden:

pip install pythia-uq 

Als Alternative, falls kein Internetzugang vorhanden ist oder neue, experimentelle Funktionen benötigt werden, kann PyThia auch direkt aus der Paketquelle installiert werden. Dies geht ganz einfach mit den folgenden drei Kommandozeilen:

git clone gitlab1.ptb.de/pythia/pythia.git # Herunterladen des Software Pakets 
cd pythia # Wechsel in den neuen Ordner 
pip install . # Installation von PyThia via pip

 

 

PyThias Konzepte

Allgemeine Polynomchaos Expansion

Die Polynomchaos Methode verwendet eine gewichtete Summe von orthonormalen Polynomen um eine Zielgröße zu approximieren, die von mehreren stochastischen Parametern abhängt. Angenommen wir haben einen Satz von $M$ Parametern $y=(y_1,\dots,y_M)\in\Gamma\subseteq\mathbb{R}^M$ für den jeder Parameter $y_m$ nach einem Maß $\pi_m$ verteilt ist. Dann ist es möglich jede quadrat-integrierbare Funktion $f\in L^2(\Gamma,\pi;\mathbb{R}^J)$ in eine Reihe

$$ f(y) = \sum_{\mu \in \mathbb{N}_0^M}f_\mu(y)P_\mu \qquad \text{wobei} \qquad f_\mu = \int_\Gamma f(y) P_\mu(y)\, \mathrm{d}\pi(y)$$

zu entwickeln. Hierbei nehmen wir an, dass die multivariate Parameterverteilung durch $\pi=\prod_{m=1}^M\pi_m$ gegeben ist und dass ${P_\mu}_{\mu\in\mathbb{N}_0^M}$ eine orthogonale und normierte Polynombasis in $L^2(\Gamma,\pi;\mathbb{R}^J)$ ist. Solche Reihenentwicklungen werden (allgemeine) polynomchaos Expansionen (PCE) genannt. Eine Approximation (Surrogat) von $f$ kann auf diese Weise einfach erhalten werden indem man eine Menge $\Lambda\subseteq\mathbb{N}_0^M$ wählt und die PCE entsprechend abschneidet

$$ f(y) \approx f^{\mathrm{PC}}(y) = \sum_{\mu\in\Lambda}f_\mu P_\mu (y). $$

Dieses Surrogat hat den Vorteil, dass Funktionsauswertungen $f^{\mathrm{PC}}(y^{(i)})$ in beliebigen Punkten $y^{(i)}\in\Gamma$ nur Auswertungen der Basispolynome $P_\mu$ benötigen und daher möglicherweise zeitaufwändige Auswertungen der eigentlichen Funktion $f$ umgehen. Außerdem ermöglicht die polynomielle Struktur des Surrogates $f^{\mathrm{PC}}$ Ableitungen und Integrationen bezüglich der Parameter analytisch zu berechnen.

 

 

Nicht-Intrusive Approximation mittels Linearer Regression

Lineare Regression stellt eine sinnvolle Methode dar um den Satz an PC Koeffizienten ${f_\mu}_{\mu\in\Lambda}$ gleichzeitig zu berechnen. Wir setzen

$$ \mathcal{V}_\Lambda = { v_N = \sum_{\mu\in\Lambda} v_\mu P_\mu \;\;  \vert\, v\mu\in\mathbb{R}^J \text{ für alle }\mu\in\Lambda} $$

als Approximationsraum und definieren das Residuum $\mathcal{R}(v_N) = f - v_N$ für alle $v_N\in\mathcal{V}_\Lambda$. Angenommen $0\in\Lambda$, dann impliziert $f_0 = \mathbb{E}_\pi[f]$, dass $\mathcal{R}(v_N)$ eine zentrierte Zufallsvariable ist. Wenn wir nun versuchen die Varianz von $\mathcal{R}$ zu minimieren erhalten wir das Least-Squares Minimierungsproblem

$$ u_N = \operatorname*{arg\, min}_{v_N\in\mathcal{V}\Lambda} \Vert \mathcal{R}(v_N)\Vert{L^2(\Gamma,\pi;\mathbb{R}^J)}^2. $$

Die Norm $\Vert \mathcal{R}(v_N)\Vert_{L^2(\Gamma,\pi;\mathbb{R}^J)}$ kann im Allgemeinen nicht analytisch berechnet werden, da keine Stammfunktion von $f$ bekannt ist. Deterministische Quadraturformeln um das integral zu approximieren sind für niedrige $M$ möglich, werden allerdings für moderat große Dimensionen ($M>3$) bereits untragbar, da die Anzahl der Quadraturpunkte exponentiell von $M$ abhängt. Um dies zu umgehen ist es möglich die Norm empirisch im Monte Carlo Sinn zu approximieren, d.h. wir verwenden die Seminorm

$$ \Vert \mathcal{R}(v_N)\Vert_{L^2(\Gamma,\pi;\mathbb{R}^J)}^2 \approx \Vert \mathcal{R}(v_N)\Vert_n^2 := \sum_{i=1}^n w^{(i)} \Vert \mathcal{R}(v_N(y^{(i)})) \Vert_2^2, $$

wobei $y^{(1)},\dots,y^{(n)}$ zufällige Realisierungen der Parameter $y\sim\pi$ und $w^{(i)}>0$ entsprechende Gewichte für die empirische Quadratur sind. Das Optimalitätskriterium erster Ordnung für das empirische Minimierungsproblem

$$ u_N = \operatorname*{arg\, min}_{v_N\in\mathcal{V}\Lambda} \Vert \mathcal{R}(v_N)\Vert{n}^2 $$

führt dann zu einem Gleichungssystem $\mathbf{G}_N \mathbf{u}_N = \mathbf{F}_N$, dass wir nach den PC Koeffizienten $\mathbf{u}_N$ von $u_N$ lösen können. Die Matrix $\mathbf{G}_N\in\mathbb{R}^{\vert\Lambda\vert\times\vert\Lambda\vert}$ ist die empirische Approximation der Gram-Matrix und $\mathbf{F}_N\in\mathbb{R}^{\vert\Lambda\vert\times J}$ enthält die Auswertungen der Zielfunktion $f$

$$ \mathbf{G}N[k,\ell] = \sum{i=1}^n w^{(i)} P_{\iota(k)}(y^{(i)}) P_{\iota(\ell)}(y^{(i)})$$  und

$$\qquad \mathbf{F}N[\ell,j] = \sum{i=1}^n w^{(i)} f_j(y^{(i)}) P_{\iota(\ell)}(y^{(i)}), $$

wobei $\iota\colon\mathbb{N}^{\vert\Lambda\vert}\to\Lambda$ eine beliebige Nummerierung der Menge $\Lambda$ ist. Als Resultat benötigen wir lediglich Trainingsdaten Paare ${y^{(i)}, f(y^{(i)})}{i=1}^n$ und entsprechende Gewichte ${w^{(i)}}$ um die Koeffizienten ${u\mu}$ des PC Surrogates zu berechnen.

Das folgende Animation zeigt beispielhaft wie die Approximation einer Funktion $f=(f_1,f_2)\colon[-7,7]\to\mathbb{R}^2$ für verschiedene Polynomgrade durch die beschriebene lineare Regression in PyThia umgesetzt wird.

Die Approximation kann mit den folgenden Zeilen Code generiert werden, wobei wir annehmen, dass die Trainingsdaten y_trainw_train und f_train bereits erzeugt wurden.

import pythia as pt # Definition des stochastischen Parameters 
 
param = pt.parameter.RandomParameter(name="y", domain=[-7, 7], ...
... distribution="uniform") # Definition der Multiindexmenge 
 
Lambda Lambda = pt.index.PCMultiIndex(1) 
Lambda.mdx = mdx.mdxFromShape([deg]) 
Lambda.refresh() # Erzeugen des PC Surrogates 
 
surrogate = pt.chaos.PolynomialChaos([param], Lambda, y_train, w_train, ...
... f_train) # Auswertung des Surrogates in 100 zufälligen Punkten  
 
y = np.random.uniform(-7, 7, (100, 1)) val = surrogate.approximation(y)

 

 

Erzeugen von Trainingsdaten

n vielen Anwendungen haben wir Zugang zu Punktauswertungen der Zielfunktion, allerdings kann es sehr zeitintensiv sein $f$ zu evaluieren. Dies ist beispielsweise der Fall, wenn das Auswerten von $f$ das Lösen von Differentialgleichungen verlangt. In diesen Fällen ist es vorteilhaft, wenn die Anzahl der Trainingsdaten so gering wie möglich ist um eine lange offline Phase zur Erzeugung der Daten zu vermeiden. Auf der anderen Seite müssen wir sicherstellen, dass die Datenpunkte $y^{(i)}$ dafür sorgen, dass die empirische Gram-Matrix $\mathbf{G}_N$ invertierbar und gut konditioniert ist um Probleme mit der numerischen Stabilität zu vermeiden. Da wir orthonormale Polynome als Basis verwenden konvergiert die Gram-Matrix gegen die Identität für $n\to\infty$. Wir müssen daher also sicherstellen, dass die Samples ${y^{(i)}}_{i=1}^n$ clever gewählt sind und dafür sorgen, dass $\Vert \mathbf{G}_N-I\Vert_2$ klein ist.

Natürlich gibt es viele verschiedene zufällige oder pseudo-zufällige Sampling Methoden, die je nach Zielfunktion $f$ unterschiedlich gute Ergebnisse liefern. Da wir allerdings keine zusätzlichen Informationen über $f$ annehmen wollen, verlassen wir uns ausschließlich auf Samples, die wir mit Hilfe der Parameterverteilung generieren können. Der offensichtliche Ansatz die Samplepunkte $y^{(i)}$ direkt nach der Verteilung $\pi$ zu ziehen und uniforme Gewichte $w^{(i)}=1/n$ zu wählen, hat den Nachteil, dass möglicherweise sehr viele Samples benötigt werden. In der Tat gibt es Abschätzungen die zeigen, dass die Anzahl der Samples $n$ wie $\vert\Lambda\vert^2\log(\vert\Lambda\vert)$ skalieren muss um zu garantieren, dass $\Vert \mathbf{G}_N-I\Vert_2$ klein ist. Das Hauptproblem stammt dabei daher, dass $\pi$ nicht notwendigerweise dafür ausgelegt ist, die inneren Produkte über die Polynombasen effizient empirisch zu approximieren.

Als eine Alternative schlägt Cohen & Migliorati 2017 vor, die Samples aus einer umgewichteten Dichte zu ziehen. Diese abgewandelte Dichte beinhaltet Informationen über die gewählte Polynombasis:

$$ \mathrm{d}\rho\Lambda = w\Lambda^{-1} \mathrm{d}\pi $$  für inverse Gewichte

$$w\Lambda^{-1}(y) = \frac{1}{\vert\Lambda\vert}\sum{k=0}^{\vert\Lambda\vert-1} \vert P_{\iota(k)}(y)\vert^2. $$

Verwenden wir Samples $y^{(i)}\sim\rho\Lambda$ und die entsprechenden Gewichte $w^{(i)}=w\Lambda(y^{(i)})$ um das Surrogat zu trainieren, reicht es die Anzahl der Trainingsdaten $n$ mit $\Lambda\vert\log(\vert\Lambda\vert)$ zu skalieren um sicherzustellen, dass $\mathbf{G}_N$ mit hoher Wahrscheinlichkeit dicht an der Identität ist. Da in jedem Fall wenigstens $n=\vert\Lambda\vert$ Samples erforderlich sind um zu garantieren, dass $\mathbf{G}_N$ (mit hoher Wahrscheinlichkeit) regulär ist, lässt sich sogar zeigen, dass diese obere Schranke an $n$ nicht verbessert werden kann.

Das folgende Bild zeigt die optimal gewichteten Dichten $\rho_\Lambda$ für verschiedene Standardverteilungen und Mengen $\Lambda$.

 

 

 

Um Samples $y^{(i)}\sim\rho_\Lambda$ mit PyThia zu erzeugen braucht es nur einige wenige Zeilen Code. Die folgenden Zeilen zeigen beispielhaft wie die optimalen Samples für $\pi=\mathcal{U}([0,1]\times[-1,5])$ generiert werden können:

import pythia as pt # Definition der stochastischen Parameter 
 
param1 = pt.parameter.RandomParameter(name="y1", ... 
... domain=[0, 1], distribution="uniform") 
 
param2 = pt.parameter.RandomParameter(name="y2", ... 
... domain=[-1, 5], distribution="uniform") 
 
params = [param1, param2] # Definition der optimalen ...
... Dichte für Lambda = {0,...,3} x {0,...,5} 
 
 
deg1 = 3 # maximaler Polynomgrad für den ersten Parameter 
 
 
deg2 = 5 # maximaler Polynomgrad für den zweiten Parameter 
 
sampler = pt.sampler.WLSTensorSampler ... 
... (params, [deg1, deg2]) # Anwenden des Samplers 
 
 
y_train = sampler.sample(100_000) # Samples erzeugen 
 
 
w_train = sampler.weight(y_train) # entsprechende Gewichte berechnen 
 
 
val = sampler.pdf(y_train) # PDF auswerten

 

 

Globale Sensitivitätsanalyse

Während der Modellierung eines physikalischen Prozesses ist es oft nicht klar wie sensitiv das Modell auf Änderungen in den Eingangsparametern ist. Daher kann es nötig sein eine Sensitivitätsanalyse durchzuführen um zu sehen wie sich Änderungen von Parametern auf die Zielfunktion auswirken. Unglücklicherweise ist die Berechnung von lokalen oder globalen Parametersensitivitäten oft sehr aufwändig da entweder partielle Ableitungen von $f$ approximiert oder hochdimensionale Integrale ausgewertet werden müssen. Wenn wir allerdings Zugang zu einem PC Surrogat $f^{\mathrm{PC}}$ von $f$ haben ist es möglich die Resultate von Sudret 2008 zu verwenden um alle Sobol Indizes ohne zusätzlichen Rechenaufwand anzunähern.

Die Sobol Indizes $\operatorname{Sob}_{j_1,\dots,j_s}$ beschreiben den Anteil, den die Änderung einer Kombination der Parametern $j_1,\dots,j_s\in{1,\dots,M}$ auf die totale Varianz von $f$ hat. Mit den Teilmengen von $\mathbb{N}_0^M$, bei denen genau die Komponenten $j_1,\dots,j_s$ von Null verschieden sind, d.h.

$$ \mathcal{F}_{j_1,\dots,j_s} = { \mu\in\mathbb{N}_0^M \,\vert\, \mu_j \neq 0 \Leftrightarrow j=j_1,\dots,j_s}, $$

können wir die Sobol Indizes mittels der Koeffizienten ${f_\mu}$ der PC Expansion beschreiben:

$$ \operatorname{Sob}{j_1,\dots,j_s} = \frac{\sum{\mu\in\mathcal{F}_{j_1,\dots,j_s}} f\mu^2}{\sum{\mu\in\mathbb{N}_{0}^M\setminus{0}} f\mu^2} \approx$$

$$ \frac{\sum{\mu\in\mathcal{F}{j_1,\dots,j_s}\cap\Lambda} f\mu^2}{\sum{\mu\in\Lambda\setminus{0}} f\mu^2} =: \operatorname{Sob}_{j_1,\dots,j_s}^\Lambda. $$

Folglich ist es sehr einfach die Sobol Indizes durch quadrieren und anschließendes aufsummieren der entsprechenden PC Koeffizienten zu approximieren, was ohne signifikanten Aufwand geschehen kann. Es ist sogar möglich den Approximationsfehler der Sobol Indizes durch das Residuum $\mathcal{R}(u_N)$ zu beschränken:

$$ \vert\operatorname{Sob}_{j_1,\dots,j_s} - \operatorname{Sob}_{j_1,\dots,j_s}^\Lambda\vert \leq \frac{\Vert \mathcal{R}(u_N)\Vert_2^2}{\operatorname{Var}[u_N]}. $$

Da PyThia das PC Surrogat berechnet indem $\Vert \mathcal{R}(u_N)\Vert_2^2$ (empirisch) minimiert wird, führt dies automatisch zu einer guten Approximation der Sobol Indizes.

PyThia berechnet die Sobol Indizes automatisch während der Initialisierung eines Polynomchaos Objektes. Der Zugriff auf die Sobol Indizes ist z.B. in den folgenden Codezeilen beschrieben:

import pythia as pt # Definition der stochastischen Parameter 
 
param1 = pt.parameter.RandomParameter(name="y1", ... 
... domain=[0, 1], distribution="uniform")
 
param2 = pt.parameter.RandomParameter(name="y2", ... 
... domain=[-1, 5], distribution="uniform") 
 
param3 = pt.parameter.RandomParameter (name="y2", ...
...domain=[-1, 5], distribution="uniform")
 
 
params = [param1, param2, param3] # Definition der Multiindexmenge Lambda 
 
Lambda = pt.index.PCMultiIndex(len(params)) 
Lambda.mdx = mdx.mdxFromShape([5, 5, 5]) # Lambda = {0, ..., 5}^3
 
 Lambda.refresh() # Erzeugen des PC Surrogates
 
 surrogate = pt.chaos.PolynomialChaos([param], Lambda, ...
... y_train, w_train, f_train) 
# Zugriff auf die Sobol Indizes 
# Bemerkung: Die Sobol Indizes sind entsprechend der 
Subskript-Tupel in Lambda.sdx geordnet surrogate.sobolCoeff[0]
# Sob_{1} surrogate.sobolCoeff[2] 
# Sob_{3} surrogate.sobolCoeff[4] 
# Sob_{1, 3} surrogate.sobolCoeff[6]
# Sob_{1, 2, 3} 

Nach oben

Anwendungen

 

Als ein simples Benchmark Beispiel um die Performance von PyThia zu visualisieren versuchen wir die Indikatorfunktion, die durch das PyThia Logo gegeben ist, zu approximieren. Aufgrund der vielen Sprünge in dieser Funktion ist dies ein sehr schweres Problem für eine globale Approximation durch Polynome und verlangt daher hohe Polynomgrade für eine adäquate Annäherung. Indem wir die optimale Sampling Dichte verwenden um die Trainingsdaten zu erzeugen, ist es jedoch möglich trotz des hohen benötigten Polynomgrades in wenigen Minuten ein Surrogat für die Logo-Funktion mit PyThia zu generieren.

Die folgende Animation zeigt die $n$ verwendeten Samplepunkte, wobei wir $n=2\vert\Lambda\vert\log(\vert\Lambda\vert)$ wählen, sowie das entsprechende resultierende Surrogat für verschiedene totale Polynomgrade.

Natürlich wurde PyThia neben den hier beschriebenen Modellbeispielen auch verwendet um reale Anwendungen aus der Metrologie to modellieren. Einige dieser Anwendungen beschreiben wir im Folgenden.

Rekonstruktion von Nanostrukturen mittels Scatterometrie

Scatterometrie ist eine indirekte und nicht-destruktive optische Messtechnik, die verwendet wird um die Form von Nanostrukturen zu bestimmen. Es ist dabei sogar möglich die Struktur von Objekten zu bestimmen, die kleiner als die verwendete Wellenlänge sind. In dieser Anwendung wurden die relativen Intensitäten der reflektierten Welle für verschiedene Einfallswinkel, Azimuthe sowie s- und p-Polarisation der ausgehenden Welle modelliert. Um die Struktur des Querschnittes der Gitterlinien zu bestimmen wurde die Oberfläche durch 6 verschiedene Parameter charakterisiert. Ein Schema für den experimentellen Aufbau lässt sich in etwa so darstellen:

 

Da das Vorwärtsmodell das Lösen der zeit-harmonischen Maxwell Gleichungen mittels Finiter Elemente Methoden verlangt und daher sehr zeitintensiv ist, haben wir PyThia verwendet um ein PC Surrogat zu generieren. Dieses Surrogat wurde für eine Sensitivitätsanalyse verwendet, die die bereits erwartete Abhängigkeit des Vorwärtsmodells von den einzelnen Parametern bestätigt hat. Anschließend wurde das inverse Problem zur Bestimmung der wahrscheinlichsten Parameterwerte und deren assoziierte Unsicherheiten gelöst indem die Bayessche Posterior Verteilung durch einen Markov-Chain Monte Carlo Sampler, der auf dem Surrogat basiert, approximiert wurde. Das folgende Bild zeigt die Approximation der Messwerte im Vergleich mit dem in den entsprechenden Mittelwerten ausgewerteten PC Surrogat.

Rekonstruktion von Nanostrukturen mittels Röntgenfluoreszenz

Röntgenfluoreszenz unter streifendem Einfallswinkel (GIXRF) ist eine Messtechnik die verwendet wird um die atomare Zusammensetzung von Materialien zu bestimmen. Dabei werden Elektronen von inneren Schalen der Atome durch Röntgenstrahlen angeregt, was dafür sorgt, dass eine elementspezifische Fluoreszenzstrahlung freigesetzt wird, die detektiert werden kann. Durch multiple Messungen unter unterschiedlichen Einfallswinkeln nahe dem Winkel der totalen Reflektion ist es sogar möglich die Position der entsprechenden Atome zu bestimmen. Dies kann für die Bestimmung der allgemeinen Form der Nanostrukturen verwendet werden.

PyThia wurde verwendet um ein Surrogat abhängig von 7 Formparametern für das Vorwärtsmodell zu erzeugen, welches das Lösen der zeit-harmonischen Maxwell Gleichungen mittels Finiter Elemente Methoden benötigt und daher sehr zeitaufwändig ist. Das Surrogat wurde für eine globale Sensitivitätsanalyse verwendet. Weiterhin hat dieses Surrogat ermöglicht, die Verteilung der Parameter in einem Bayesschen inversen Problem mittels Markov-Chain Monte Carlo Sampling zu approximieren. Die Marginale dieser Verteilung sind im folgenden Bild dargestellt.

Publikationen

N. Farchmin, M. Hammerschmidt, P. I. Schneider, M. Wurm, B. Bodermann, M. Bär;S. Heidenreich
Journal of Micro/Nanolithography, MEMS, and MOEMS, 2(19),
1--11,
2020.
Export als:
BibTeX, XML