Funktionsfabriken in R

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Funktionen in R
      • Diagnose-Funktionen für Funktionen in R
      • Spezielle selbstdefinierte Funktionen in R
      • Funktionsfabriken in R
  • Graphik
    • Erzeugen von dreidimensionalen Graphiken mit scatterplot3d

In Diagnose-Funktionen für Funktionen in R wurde schon ein erster Einblick in die Bedeutung der Umgebung einer Funktion gegeben und wie man auf die Umgebung zugreift. Dies wird her benötigt und fortgeführt.

Und in Spezielle selbstdefinierte Funktionen in R wurden erste Konzepte der funktionalen Programmierung eingeführt; diese werden hier benötigt und erweitert.

Zu den hier vorgestellten Beispielen mit quadratischen Formen werden dreidimensionale Diagramme erstellt, die die Bibliothek scatterplot3d einsetzen. Kenntnisse über diese Bibliothek sind hilfreich, aber nicht nötig um dem eigentlichen Thema – Funktionsfabriken – zu folgen.

Einführung

In Spezielle selbstdefinierte Funktionen in R wurden bereits Funktionen höherer Ordnung vorgestellt und (teilweise) behandelt; dieser Begriff umfasst mehrere Arten von Funktionen:

1. Funktionen, die als Eingabewert wiederum eine Funktion besitzen; der Eingabewert wird dann auch als funktionaler Parameter bezeichnet. Ein Spezialfall davon liegt bei einem Funktional vor: jetzt ist der Rückgabewert ein Vektor.
2. Funktionen, die als Rückgabewert eine Funktion besitzen.
3. Funktionen die beides kombinieren.

Bekannte Beispiele für Erstere sind Funktionen wie outer() oder die Familie der apply()-Funktionen. Die zweite Gruppe wurden in Spezielle selbstdefinierte Funktionen in R noch nicht behandelt. Da dies sehr umfangreich ist und insbesondere ein tieferes Verständnis der Umgebung einer Funktion erfordert, ist dieser Gruppe ein eigener Artikel gewidmet.

Die dritte Gruppe muss dann nicht eigens behandelt werden, da man leicht die beiden Möglichkeiten kombinieren kann.

Hier werden:

• Beispiele für Funktionen mit als Rückgabewert gezeigt. Nämlich speziell – und das ist die häufigste Anwendung – Funktionsfabriken, die Funktionen herstellen, die untereinander eine ähnliche Struktur besitzen (hier werden quadratische Formen als Beispiel gewählt).
• Der Begriff der Umgebung wurde in Diagnose-Funktionen für Funktionen in R schon eingeführt und wird aus der Perspektive erläutert, wie er zum Einsatz von Funktionsfabriken nötig ist.
• Da die Implementierung einer Funktionsfabrik niemals zwingend nötig ist, sondern sich durch "herkömmliche" Funktionen ersetzen lässt, wird diskutiert, wann Funktionsfabriken eingesetzt werden sollen.

Beispiel zur Einführung: Die Implementierung einer quadratischen Form

Was ist eine quadratische Form?

Es gibt zahlreiche Anwendungen, in denen Terme der Art

½·a·x²

vorkommen; typische Beispiele sind der quadratische Terme in der Taylor-Entwicklung einer unendlich oft differenzierbaren Funktion

f(x) = f(x₀) + f'(x₀) · (x - x₀) + ½ · (x - x₀)² · f(x₀) + ...

oder die kinetische Energie T:

T = ½·m·v².

Mit einer quadratischen Form Q(x) werden Terme wie ½·a·x² auf mehrere Dimensionen verallgemeinert:

• Aus der Zahl a wird eine quadratische Matrix A und
• aus den beiden Faktoren x werden Vektoren.
• An die Stelle der Multiplikation von Zahlen tritt die Matrizenmultiplikation beziehungsweise das Skalarprodukt.

Damit diese wohldefiniert sind, muss man in der richtigen Reihenfolge multiplizieren:

Q(x) = x · (A · x).

Zuerst wird die Matrizenmultiplikation A·x ausgeführt, anschließend wird das Skalarprodukt mit x berechnet.

(Ob man den Faktor ½ in die Matrix A integriert oder explizit angibt, ist Geschmackssache – in der Literatur sollten Sie immer darauf achten, wie die quadratische Form definiert ist.)

Damit man die Multiplikationen oben ausführen kann, muss der zweite Faktor x ein Spaltenvektor sein; die Multiplikation A · x liefert dann ebenfalls einen Spaltenvektor. Um dann die erste Multiplikation wie ein Skalarprodukt zu berechnen, muss der erste Faktor ein Zeilenvektor sein. Man schreibt daher das Produkt x·A·x häufig wie ein Skalarprodukt ⟨ . | . ⟩:

Q(x) = ⟨x|A·x⟩

In dieser Schreibweise sollte dann auch klar sein, dass die Dimensionen übereinstimmen müssen:

• Der Vektor x hat n Komponenten.
• Die Matrix A ist eine nxn-Matrix.

In der Form Q(x) = ⟨x|A·x⟩ gibt es zahlreiche Anwendungen von mehrdimensionalen quadratischen Formen, um zwei Beispiele zu nennen:

1. Die Formel T = ½·m·v² von oben beschreibt die kinetische Energie einer Punktmasse m, die sich mit Geschwindigkeit v bewegt. Die Rotationsenergie eines ausgedehnten Körpers kann man sich als zusammengesetzt aus den kinetischen Energien von Massenelementen Δm vorstellen. Durch Einführung des sogenannten Trägheitstensors J kann man herleiten, dass die Rotationsenergie durch T_rot = ½· ⟨ω|J·ω⟩ berechnet wird, wobei ω der Vektor der Winkelgeschwindigkeit ist.
2. Die oben angegebene Formel für die Taylor-Entwicklung kann man auch so lesen, dass f(x) eine reelle Funktion ist, die für n-dimensionale Vektoren x definiert ist. Im linearen Term der Taylor-Entwicklung steht dann das Skalarprodukt aus dem Gradienten von f (an der Stelle x₀) und dem Vektor (x - x₀). Im quadratischen Term steht die quadratische Form, die mit Hilfe der Hesse-Matrix H gebildet wird:

½· ⟨(x - x₀)|H·(x - x₀)⟩.

Die mathematisch korrekte Definition einer quadratischen Form ist in Abbildung 1 gezeigt.

Abbildung 1: Definition einer quadratischen Form als einer reellen Abbildung, die durch eine quadratische Matrix induziert wird. Der Definitionsbereich sind n-dimensionale Vektoren. Unten wird für eine symmetrische Matrix die quadratische Form explizit angegeben. Das Beispiel wird unten fortgeführt.

Zur Veranschaulichung zeigt Abbildung 2 die quadratische Form aus dem Beispiel aus Abbildung 1 als "Gebirge" über der x₁-x₂-Ebene.

Abbildung 2: Darstellung einer zweidimensionalen quadratischen Form als Gebirge über der x₁-x₂-Ebene. Verwendet wird das Beispiel aus Abbildung 1.

Für Abbildung 3 werden die nur die Punkte auf dem Einheitskreis ausgewertet:

Q(x₁, x₂) = (x₁ - x₂)² für x₁² + x₂² = 1.

Genauer es werden zu den Winkeln

φ = 0°, 1°, ..., 359°

die Koordinaten der Punkte (x₁ = cos φ, x₂ = sin φ) in die quadratische Form eingesetzt und als Histogramm dargestellt. Man beachte, dass Abbildung 2 und 3 nur schwer miteinander uu vergleichen sind, da die z-Achse anders skaliert ist.

Zusätzlich werden in Abbildung 3 die beiden Eigenvektoren der Matrix in die x₁-x₂-Ebene eingetragen. Dieses Problem wird später verwendet, um die Umgebung (environment) einer Funktion besser zu verstehen.

Abbildung 3: Für die quadratische Form, die schon in Abbildung 2 dargestellt wurde, sind hier nur Funktionswerte für Punkte (x₁, x₂) auf dem Einheitskreis als Histogramm gezeichnet. In der x₁-x₂-Ebene sind zusätzlich die Eigenvektoren der zugehörigen Matrix A gezeigt. Unten wird dann – als Simulation einer Auswertung einer quadratischen Form – eine Funktion implementiert, die diesen Plot erstellt.

Implementierung der quadratischen Form als herkömmliche Funktion

Um besser zu verstehen, wie eine Funktionsfabrik arbeitet und wann sie eingesetzt werden soll, wird zunächst die quadratische Form als "herkömmliche" Funktion implementiert. Sie besitzt als Eingabewerte eine Matrix A und einen Vektor v, der Rückgabewert ist die reelle Zahl, die bei der Auswertung der quadratischen Form entsteht.

Damit die Eingabewerte eine wohldefinierte quadratische Form erzeugen, werden sie zunächst überprüft – durch den recycling-Mechanismus von R sind auch Kombinationen denkbar, bei denen mit ein Rückgabewert berechnet wird, obwohl die quadratische Form im Sinne der Mathematik nicht definiert ist. Die Prüfung der Eingabewerte umfasst:

1. Der Eingabewert A muss eine Matrix, v ein Vektor sein (Zeile 2).
2. Die Matrix A muss quadratisch sein, der Vektor v muss so viele Komponenten besitzen wie die Matrix A Zeilen (oder Spalten) besitzt (Zeile 3 und 4).

Das folgende Skript zeigt eine mögliche Implementierung:

quadrForm <- function(A, v){
  stopifnot(is.matrix(A), is.vector(v))
  N <- nrow(A)
  stopifnot( identical(N, ncol(A)), identical(N, length(v)) )
  return( sum( v * (A %*% v) ) )
}

Zeile 5: Auf den ersten Blick mag es befremdlich erscheinen, warum die Auswertung der quadratischen Form nicht ausschließlich mit dem Skalarprodukt erfolgt, etwa mit v %*% (A %*% v) anstelle von sum( v * (A %*% v) ) . Wird auch die Multiplikation von v (von links) als Skalarprodukt realisiert, wird insgesamt ein Feld mit den Dimensionen 1x1 erzeugt, was hier nicht sinnvoll ist, da der Rückgabewert der quadratischen Form wie eine Zahlen eingesetzt wird.

Einen Aufruf der Funktion quadrForm() zeigt das folgende Skript, man kann die Berechnung leicht mit dem Beispiel aus Abbildung 1 nachvollziehen:

A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE)
A
#      [,1] [,2]
# [1,]    1   -1
# [2,]   -1    1

quadrForm(A = A, v = c(1, -1))    # 4

Auswertung und Plot einer quadratischen Form

Als Beispiel für den Einsatz einer quadratischen Form soll Abbildung 3 dienen: Es soll eine Funktion plot_QF_Hist() implementiert werden, die folgende Anforderungen erfüllt (QF steht natürlich für quadratische Form und Hist für Histogramm):

1. Eingabewert ist eine quadratische Form qf, das heißt es soll ein Funktional im Sinne der funktionalen Programmierung entwickelt werden.
2. Die Koordinaten (x₁, x₂) der 360 Punkte auf dem Einheitskreis werden berechnet.
3. Für diese Punkte wird die quadratische Form qf ausgewertet.
4. In einem Dataframe werden die x₁-, x₂-Werte sowie die Werte qf(x₁, x₂) abgespeichert.
5. Dieses Dataframe ist zugleich der Rückgabewert der Funktion plot_QF_Hist(), damit man mit den Werten eventuell weitere Berechnungen anstellen kann.
6. Das Histogramm wie in Abbildung 3 wird erzeugt. Dazu wird die Bibliothek scatterplot3d verwendet.
7. Die Eigenvektoren werden noch nicht berechnet und in das Diagramm eingetragen – dies geschieht später.

library("scatterplot3d")

plot_QF_Hist <- function(qf){
  phis = (0:359) * pi / 180
  xs <- cos(phis)
  ys = sin(phis)
  N <- length(phis)
  clrs <- rainbow(n = N)
  
  vs <- mapply(FUN = c, xs, ys, SIMPLIFY = FALSE)
  zs <- sapply(X = vs, FUN = qf)
  
  df <- data.frame(x = xs, y = ys, z = zs)
  z.lim <- c( min(df$z, 0), max(df$z, 0) )
  
  sc <- scatterplot3d(x = 0, y = 0, z = 0,
                      color = "black",
                      type = "p",
                      xlim = c(-2, 2), ylim = c(-2, 2),
                      zlim = z.lim,
                      xlab = "x1", ylab = "x2", zlab = "z",
                      main = "z = Q(x)",
                      lwd = 2, angle = 72, box = TRUE)
  # Histogramm:
  # Alle Punkte aus df projizieren:
  df.proj <- sc$xyz.convert(x = df)
  # Alle Punkte auf dem Einheitskreis projizieren:
  df.0 <- cbind(df$x, df$y, rep(x = 0, times = N))
  df.0.proj <- sc$xyz.convert(x = df.0)
  # mit segments() eintragen:
  lapply( X = (1:N),
          FUN = function(i){segments(x0 = df.0.proj$x[i], x1 = df.proj$x[i], 
                                     y0 = df.0.proj$y[i], y1 = df.proj$y[i], 
                                     col = clrs[i], lwd = 1);
            return()} )
  
  # Ebene z = 0
  sc$plane3d(Intercept = 0, 
             x.coef = 0, y.coef = 0, lty = "dotted", col = "blue")
  
  return(df)
}

Einige Besonderheiten des Skriptes sollen näher erklärt werden:

Zeile 1: Der Import der Bibliothek scatterplot3d. In folgenden Skripten wird er nicht mehr gezeigt.

Zeile 4 bis 6: Die Winkel von 0° bis 359° werden ins Bogenmaß umgerechnet und die x- und y-Koordinaten der Punkte auf dem Einheitskreis berechnet.

Zeile 8: Als Farbpalette zum Färben des Histogramms wird rainbow() verwendet.

Zeile 10: Da die quadratische Form als Eingabewert einen Vektor v erwartet, müssen die x- und y-Koordinaten jeweils zu einem Vektor zusammengefasst werden. Dies geschieht mit Hilfe von mapply(), wodurch eine Liste mit 360 Komponenten entsteht. Da in mapply() angegeben wird FUN = c , werden die x- und y-Koordinaten tatsächlich zu einem Vektor zusammengefasst und ergeben je eine Komponente der Liste vs. Dadurch dass in mapply() SIMPLIFY = FALSE gesetzt wurde, entsteht eine Liste; hätte man SIMPLIFY = TRUE gesetzt, wird eine Matrix erzeugt – dann muss man aber bei der Weiterverarbeitung darauf achten, dass man nicht Zeilen und Spalten vertauscht.

Zeile 11: Auf die 360 Komponenten der Liste vs wird die quadratische Form qf angewendet. Dies geschieht mis sapply(), wodurch die berechneten Werte zu einem Vektor zusammengefasst werden.

Man könnte Zeile 10 und 11 auch in eine Anweisung zusammenfassen, die dann aber schwer verständlich ist. Das Objekt vs wird nicht benötigt und wurde hier nur berechnet, um den Quelltext klarer zu formulieren.

Zeile 13: Die x- und y-Koordinaten der Punkte auf dem Einheitskreis sowie die von der quadratischen Form berechneten z-Werte werden in einem Dataframe (mit drei Spalten) abgespeichert. Auch die Berechnung des Dataframes kann man abkürzen, indem man die Anweisungen verschachtelt. Das Dataframe wird nicht mehr verändert und am ist Ende der Rückgabewert.

Zeile 14: Zur Skalierung der z-Achse des Diagramms benötigt man den Minimal- und Maximalwert der z-Werte; sie werden in z.lim abgespeichert. Allerdings kann es sein, dass beide Werte positiv beziehungsweise negativ sind. In diesem Fall wäre die Ebene z = 0 nicht sichtbar. Damit z = 0 immer dargestellt wird, beginnt die z-Achse bei min(df$z, 0) und endet bei max(df$z, 0)

Zeile 16 bis 23: Die Funktion scatterplot3d() wird vorerst nur verwendet, um zur Orientierung den Ursprung in das Koordinatensystem einzutragen. Ansonsten wird hier lediglich der Plot konfiguriert: es wird die Skalierung der Achsen festgelegt und das Diagramm beschriftet.

Zeile 26 bis 35: Das eigentliche Histogramm wird gezeichnet. Die Aufgabe ist deswegen so aufwendig, weil ein Histogramm in scatterplot3d immer vom "Boden" des Diagramms aus gezeichnet wird – aber der muss nicht mit der z = 0 Ebene übereinstimmen. Um die Balken des Histogramms immer bei z = 0 beginnen zu lassen, sind einige Verrenkungen nötig:

Das Dataframe df enthält die dreidimensionalen Koordinaten der zu zeichnenden Funktionswerte, also die Endpunkte der Balken. Die Anfangspunkte erhält man, indem man in df jeweils die z-Koordinate gleich null setzt; dies geschieht im weiteren Dataframe df.0 (siehe Zeile 28).

Um jetzt diese Anfangs- und Endpunkte mit Hilfe von segments() verbinden zu können, muss man sie erst in das zweidimensionale Koordinatensystem des Plots projizieren; dies geschieht mit Hilfe der Funktion xyz.convert() (Zeile 26 und 29). Die Funktion xyz.convert() muss nicht die drei Koordinaten als Eingabewerte erhalten, sondern kann sofort mit dem entsprechenden Dataframe aufgerufen werden.

Mit Hilfe von lapply() iteriert man jetzt über alle 360 Winkel und verbindet die Anfangs- und Endpunkte (der Projektionen), siehe Zeile 31 bis 35. Als Farbe des jeweiligen Balkens wird dann eine Farbe aus der rainbow-Farbpalette genommen.

Zeile 38 und 39: Zur besseren Orientierung wird die xy-Ebene blau und gestrichelt eingetragen.

Zeile 41: Rückgabewert ist das Dataframe mit den Koordinaten der gezeichneten Balken.

Das folgende Skript zeigt ein Beispiel für den Aufruf der Funktion plot_QF_Hist():

df <- plot_QF_Hist(qf = function(v){A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE);
                                return(quadrForm(A = A, v = v))})

Man beachte, dass an das Argument qf von plot_QF_Hist() nicht einfach eine quadratische Form wie quadrForm(A = A, v = v) übergeben werden kann, da die beiden Argumente A und v eine unterschiedliche Rolle spielen:

• Die Matrix A konfiguriert die quadratische Form, die von quadrForm() erzeugt wird.
• Der Eingabewert v ist der eigentliche Eingabewert der quadratischen Form, der hier nicht mit einem speziellen Vektor gesetzt werden darf. An das Argument qf soll ja eine quadratische Form weitergegeben werden, die v als ihren Eingabewert besitzt. Daher wird im Argument qf eine anonyme Funktion implementiert, die v als Argument besitzt. Dagegen wird die Matrix A in quadrForm() auf den speziellen Wert der Matrix matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE) gesetzt. Damit wird an die Funktion plot_QF_Hist() eine quadratische Form übergeben, die einen Vektor als Eingabewert besitzt; die Matrix kann nicht mehr variiert werden.

Das Diagramm, das durch den Aufruf erzeugt wird stimmt nahezu mit Abbildung 3 überein: es fehlen lediglich die beiden Eigenvektoren der Matrix A in der z = 0 Ebene.

Die Implementierung einer quadratischen Form mit Hilfe einer Funktionsfabrik

Die Erklärung des Aufrufs von plot_QF_Hist() aus dem letzten Skript liefert zugleich den Zugang, wie die Implementierung einer Funktionsfabrik aussehen muss, die die Funktion quadrForm() ersetzt. Die folgende Abbildung 4 soll die Vorgehensweise zeigen.

Dazu ist links nochmals die Funktion quadrForm() dargestellt: sie erhält die Matrix A und den Vektor v als Eingabewerte und berechnet als Rückgabewert den Funktionswert der quadratischen Form. Was an dieser Funktion nicht klar getrennt wird, ist der unterschiedliche Zweck, den die beiden Eingabewerte haben:

• Die Matrix A konfiguriert die quadratische Form.
• Der Vektor v ist dann für diese spezielle quadratische Form der eigentliche Eingabewert.

Mit Hilfe einer Funktionsfabrik lassen sich diese beiden Aufgaben der Eingabewerte voneinander trennen (siehe rechter Teil der Abbildung 4):

• Die Funktionsfabrik, hier mit Namen factory_quadrForm() erhält als Eingabewert eine Matrix A.
• Der Rückgabewert der Funktionsfabrik ist die von A erzeugte quadratische Form. Und zwar als eine Funktion, die als Eingabewert einen Vektor v erhält.

Abbildung 4: Um die Arbeitsweise der Funktionsfabrik besser zu verstehen, ist links nochmals die herkömmliche Funktion dargestellt, die eine quadratische Form realisiert: Sie erhält als Eingabewert eine Matrix A und einen Vektor v. Rechts ist die entsprechende Funktionsfabrik gezeigt. Sie erhält nur die Matrix und erzeugt damit eine Funktion qf(). Diese Funktion wird dann mit einem Vektor v aufgerufen.

Das folgende Skript zeigt eine Möglichkeit, wie sich die Fabrik factory_quadrForm() implementieren lässt:

factory_quadrForm <- function(A){
  stopifnot(is.matrix(A))
  N <- nrow(A)
  stopifnot( identical(N, ncol(A)) )
  
  qf <- function(v){
    stopifnot( identical(N, length(v)) )
    return( sum( v * (A %*% v) ) )
  }
  
  return(qf)
}

Zeile 2 bis 4: Die Funktion factory_quadrForm() enthält jetzt nur die Prüfung von Eingabewerten, die sich auf die Matrix A beziehen.

Zeile 6 bis 9: Es wird eine interne Funktion qf implementiert, die

• prüft, ob der Eingabe-Vektor v die geeignete Länge zur gegebenen Matrix A besitzt (Zeile 7) und
• genau diejenigen Berechnungen durchführt, die in der Funktion quadrForm() enthalten waren (Zeile 8).

Zeile 11: Die interne Funktion qf wird zurückgegeben.

Soll jetzt eine spezielle quadratische Form ausgewertet werden – etwa indem sie an die Funktion plot_QF_Hist() übergeben wird –, so wird sie zuerst von der Fabrik erzeugt und an die Auswertung weitergereicht. Das folgende Skript zeigt die Anweisungen, die auch wieder verschachtelt als eine Anweisung geschrieben werden können.

A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE)
qf <- factory_quadrForm(A)
plot_QF_Hist(qf)

Wann soll man Funktionsfabriken einsetzen?

Die Argumente

Man kann sich natürlich auf den Standpunkt stellen, dass eine Funktionsfabrik niemals nötig ist, da man sie immer durch eine "herkömmliche" Funktion ersetzen kann. Das bisher vorgestellte Beispiel ist nicht komplex genug, dass eine Funktionsfabrik zu bevorzugen wäre.

Die drei Argumente, die dennoch für eine Funktionsfabrik sprechen, die hier diskutiert werden, kann man kurz wie folgt beschreiben:

1. Eine Funktionsfabrik erlaubt eine bessere Trennung der Aufgaben der eingesetzten Funktionen.
2. Das Beispiel der quadratischen Form zeigt eine typische Eigenschaft, wie Funktionen meist eingesetzt werden: sie werden einmal erzeugt und dann oft aufgerufen. Und kann man beim Erzeugen der Funktion schon Vorbereitungen für ihren Einsatz treffen, die nur einmalig geschehen müssen, so kann man diese Vorbereitungen in die Funktionsfabrik verlegen.
3. Müssen sehr viele, unterschiedlich konfigurierte Funktionen erzeugt werden, kann die Funktionsfabrik – zusammen mit der Möglichkeit, Funktionen in Listen abzuspeichern, – eine sehr schlanke Lösung gegenüber einer herkömmlichen Funktion liefern.

Diese drei Argumente sollen etwas ausführlicher besprochen werden.

Trennung von Aufgaben

Um für mehr Klarheit in den Quelltexten zu sorgen, sollte man aber das Prinzip anwenden, dass eine Funktion immer nur eine klar definierte Aufgabe erfüllen soll. Oder umgekehrt formuliert: Wenn man feststellt, dass eine Funktion mehrere Aufgaben erfüllt, sollte man sich überlegen, ob es nicht ratsam ist, sie in mehrere Funktionen aufzuspalten. Und genau dies geschieht hier in der Trennung der Verarbeitung der Matrix und des Vektors bei der Berechnung einer quadratischen Form.

So konnte die Prüfung der Eingabewerte, die sich in quadrForm() sowohl auf die Matrix A als auch den Vektor v bezogen, mit der Funktionsfabrik factory_quadrForm() klar in zwei Gruppen getrennt werden.

Zugegeben das Beispiel mag wenig überzeugend wirken, da es hier keinen großen Unterschied macht, ob man mit der Fabrik factory_quadrForm() oder der Funktion quadrForm() arbeitet. Sind allerdings, die Eingabewerte, die zur Konfiguration nötig sind – wie hier die Matrix A –, sehr umfangreich, so gelingt es mit der Fabrik, diese komplexe Konfiguration abzuspalten und kann dann mit der schlanken (von der Fabrik erzeugten) Funktion arbeiten.

Vorbereitungen für den Einsatz der Funktion

Im Beispiel oben mit factory_quadrForm() wird die Fabrik mit der Matrix A konfiguriert und die erzeugte Funktion besitzt einen Vektor v als Eingabewert. Oft gibt es weitere Aufgaben, die ausschließlich mit den Objekten der Konfiguration erledigt werden können. Sie sollten dann vor dem Erzeugen der Funktion schon innerhalb der Fabrik geschehen und nicht erst in der erzeugten Funktion. Denn der typische Fall ist, dass eine einmal erzeugte Funktion sehr oft eingesetzt wird (im Beispiel oben: zum Plotten der Graphen, wozu sehr viele Funktionswerte berechnet werden müssen). Die vorbereitenden Aufgaben werden dann nur einmal, nämlich in der Funktionsfabrik, erledigt.

(Nebenbei: Dies widerspricht zwar der soeben besprochenen "Trennung von Aufgaben", aber man kann die Trennung wieder herstellen, indem man die vorbereitenden Aufgaben in eine eigene Funktion auslagert, die dann von der Fabrik aufgerufen wird.)

Im Beispiel von factory_quadrForm() könnte es vielleicht wichtig sein, zusätzliche Informationen über die Matrix A bereitzustellen, etwa ihre Eigenwerte und Eigenvektoren, auf die dann beim Einsatz der erzeugten quadratischen Form zugegriffen wird.

Im nächsten Abschnitt "Die Umgebung (environment) einer Funktion" wird dazu ein Beispiel ausführlich besprochen. Dazu ist es nötig zu zeigen, wie Funktionen auf ihre Umgebung zugreifen können und was mit Umgebung einer Funktion (der sogenannten environment) in R überhaupt gemeint ist.

Erzeugen vieler Funktionen

Solange nur wenige Funktionen von der Fabrik erzeugt werden, ist es tatsächlich schwer zu argumentieren, warum sie implementiert werden soll. Dazu wird im letzten Abschnitt "Weiteres Beispiel: Erzeugen von mehreren quadratischen Formen" eine Simulation vorgestellt, in der mehrere quadratische Formen (zufällig) erzeugt werden, die dann ausgewertet werden (hier wieder mit Plots, für die zahlreiche Funktionsaufrufe stattfinden müssen).

In diesem Beispiel sollte dann besser verständlich werden, warum es so wichtig ist, die Aufgaben zwischen der Fabrik und der erzeugten Funktion zu trennen. Zusammen mit der Möglichkeit, Funktionen in Listen abzuspeichern, und zur Auswertung über diese Listen zu iterieren, lassen sich schlankere und leicht veränderbare Quelltexte erstellen.

Die Umgebung (environment) einer Funktion

Eine Erweiterung der Funktionsfabrik factory_quadrForm()

Das folgende Skript zeigt eine kleine Erweiterung der Funktionsfabrik factory_quadrForm() gegenüber der Implementierung, die oben vorgestellt wurde:

• Die Fabrik besitzt jetzt einen weiteren Eingabewert eigen (mit default-Wert FALSE).
• Falls eigen gesetzt ist, werden die Eigenwerte und Eigenvektoren der Matrix A berechnet und ausgegeben. (Die Ausgabe dient natürlich nur zum Kennenlernen des Quelltextes und sollte später abgeschaltet werden.)

factory_quadrForm <- function(A, eigen = FALSE){
  stopifnot(is.matrix(A))
  N <- nrow(A)
  stopifnot( identical(N, ncol(A)) )
  
  qf <- function(v){
    stopifnot( identical(N, length(v)) )
    return( sum( v * (A %*% v) ) )
  }
  
  if(eigen){
    e <- eigen(A)
    ew <- e$values
    ev <- e$vectors
    cat("Eigenschaften der Matrix A:\n")
    print(A)
    cat("\nEW und EV:\n")
    print(ew)
    print(ev)
  }
  
  return(qf)
}

Neu gegenüber der alten Version sind somit nur der Eingabewert eigen und die Zeilen 11 bis 20; dabei steht ew für Eigenwerte und ev für Eigenvektoren.

Einen Test zeigt das folgende Skript:

A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE)
qf <- factory_quadrForm(A, eigen = TRUE)

# Eigenschaften der Matrix A:
#   [,1] [,2]
# [1,]    1   -1
# [2,]   -1    1
# 
# EW und EV:
#   [1] 2 0
# [,1]       [,2]
# [1,] -0.7071068 -0.7071068
# [2,]  0.7071068 -0.7071068

Die Matrix A aus Zeile 1 wurde schon verwendet, um die quadratische Form für Abbildung 3 zu erzeugen. Man erkennt, dass dort die Eigenvektoren (in der x₁-x₂-Ebene eingetragen sind). Und wer mit Linearer Algebra vertraut ist, kann die Eigenwerte verwenden, um die Funktionswerte der quadratischen Form auf dem Schnittpunkten des Einheitskreises mit den Eigenvektoren zu berechnen.

Die Bedeutung der Umgebung einer Funktion

Um wie in Abbildung 3 die Eigenvektoren in das Diagramm einzutragen, muss die Funktion plot_QF_Hist() erweitert werden. Und dazu soll auf die zuletzt in factory_quadrForm() berechneten Eigenwerte und Eigenschaften zugegriffen werden.

Wer mit lokalen Variablen und allgemeiner dem Gültigkeitsbereich von Variablen vertraut ist, wird sich jetzt wundern, wie es möglich sein soll in der Funktion plot_QF_Hist() auf die Variablen in der Funktion factory_quadrForm() zuzugreifen. Eine Lösung mit globalen Variablen soll hier nicht präsentiert werden – derartige Lösungen sind meist nur in kleinen Anwendungen zu empfehlen, in großen Projekten erzeugen sie meist mehr Schaden als Nutzen.

Aber hier kommt die Umgebung einer Funktion zum Einsatz: Die quadratische Form, die von der Fabrik factory_quadrForm() erzeugt wird, hat über ihre Umgebung (environment) Zugriff auf die lokalen Variablen der Funktion factory_quadrForm(). Man muss also nur noch verstehen, wie dieser Zugriff syntaktisch richtig realisiert wird.

(In Diagnose-Funktionen für Funktionen in R wurde bereits erwähnt, dass jede Funktion eine Umgebung besitzt. Dort wurde lediglich diskutiert, dass sich damit Namenskonflikte vermeiden lassen, wenn es Funktionen mit identischen Namen in unterschiedlichen Paketen gibt. Ansonsten konnte dort noch kein relevantes Anwendungsbeispiel der Umgebung gegeben werden.)

Das folgende Skript soll die Syntax für den Zugriff auf die Variablen in der Fabrik erklären, die praktische Anwendung erfolgt dann im nächsten Unterabschnitt, der die Funktion plot_QF_Hist() erweitert.

Die oben erweiterte Fabrik factory_quadrForm() (erweitert um die Eigenwerte und Eigenvektoren der Matrix A) wird jetzt eingesetzt, um von der erzeugten quadratischen Form aus, auf die Variablen ew und ev zuzugreifen:

A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE)
qf <- factory_quadrForm(A, eigen = TRUE)

env <- environment(fun = qf)
str(env)
# <environment: 0x00000000074ec1c8> 

env.list <- as.list(env)
str(env.list)
# List of 7
# $ ev   : num [1:2, 1:2] -0.707 0.707 -0.707 -0.707
# $ ew   : num [1:2] 2 0
# $ e    :List of 2
# ..$ values : num [1:2] 2 0
# ..$ vectors: num [1:2, 1:2] -0.707 0.707 -0.707 -0.707
# ..- attr(*, "class")= chr "eigen"
# $ qf   :function (v)  
#   ..- attr(*, "srcref")= 'srcref' int [1:8] 6 9 9 3 9 3 6 9
# .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x0000000004da0c98> 
#   $ N    : int 2
# $ A    : num [1:2, 1:2] 1 -1 -1 1
# $ eigen: logi TRUE

env.list$ew
# [1] 2 0

Zeile 1 und 2: Es wird wieder die bekannte Matrix erzeugt und mit Hilfe der Fabrik eine quadratische Form qf erzeugt.

Zeile 4 bis 6: Man kann zwar die environment der Funktion qf aufrufen und untersuchen, aber damit ist nichts gewonnen.

Zeile 8 bis 22: Verwandelt man dagegen die environment in eine Liste (mit as.list() ), so erkennt man an der Struktur dieser Liste, dass alle Variablen enthalten sind, die innerhalb der Fabrik factory_quadrForm() erzeugt wurden. Das heißt die quadratische Form qf hat über ihre Umgebung Zugang zu den Variablen in der Fabrik, in der sie hergestellt wurde. Und jetzt kann man (mit dem Zugriffsoperator für Listen-Komponenten $ ) auf die einzelnen Komponenten zugreifen. Man erkennt an den letzten beiden Einträgen der Liste, dass auch die Eingabewerte der Fabrik in der Umgebung enthalten sind.

Zeile 24 und 25: Als Beispiel ist der Zugriff auf die beiden Eigenwerte der Matrix A gezeigt.

Erweiterung der Funktion plot_QF_Hist()

Das folgende Skript erweitert jetzt die Implementierung der Funktion plot_QF_Hist():

• Über die Umgebung der Funktion qf (also dem Eingabewert von plot_QF_Hist()) wird auf die Eigenvektoren der Matrix A zugegriffen.
• Die beiden Eigenvektoren werden in der z = 0 Ebene des Plots eingetragen.

plot_QF_Hist <- function(qf){
  phis = (0:359) * pi / 180
  xs <- cos(phis)
  ys = sin(phis)
  N <- length(phis)
  clrs <- rainbow(n = N)
  
  vs <- mapply(FUN = c, xs, ys, SIMPLIFY = FALSE)
  zs <- sapply(X = vs, FUN = qf)
  
  df <- data.frame(x = xs, y = ys, z = zs)
  z.lim <- c( min(df$z, 0), max(df$z, 0) )
  
  sc <- scatterplot3d(x = 0, y = 0, z = 0,
                      color = "black",
                      type = "p",
                      xlim = c(-2, 2), ylim = c(-2, 2),
                      zlim = z.lim,
                      xlab = "x1", ylab = "x2", zlab = "z",
                      main = "z = Q(x)",
                      lwd = 2, angle = 72, box = TRUE)
  # Histogramm:
  # Alle Punkte aus df projizieren:
  df.proj <- sc$xyz.convert(x = df)
  # Alle Punkte auf dem Einheitskreis projizieren:
  df.0 <- cbind(df$x, df$y, rep(x = 0, times = N))
  df.0.proj <- sc$xyz.convert(x = df.0)
  # mit segments() eintragen:
  lapply( X = (1:N),
          FUN = function(i){segments(x0 = df.0.proj$x[i], x1 = df.proj$x[i], 
                                     y0 = df.0.proj$y[i], y1 = df.proj$y[i], 
                                     col = clrs[i], lwd = 1);
            return()} )
  
  # Ebene z = 0
  sc$plane3d(Intercept = 0, 
             x.coef = 0, y.coef = 0, lty = "dotted", col = "blue")
  
  # Eintragen der EV von A in die xy-Ebene:
  ev <- 2 * as.list(environment(fun = qf))$ev
  
  ev1.proj <- sc$xyz.convert(x = ev[ 1, 1], y =  ev[ 2, 1], z = 0)
  ev1.proj.neg <- sc$xyz.convert(x = -ev[ 1, 1], y =  -ev[ 2, 1], z = 0)
  # str(ev1.proj)    # List of 2: $x, $y
  ev2.proj <- sc$xyz.convert(x = ev[ 1, 2], y =  ev[ 2, 2], z = 0)
  ev2.proj.neg <- sc$xyz.convert(x = -ev[ 1, 2], y =  -ev[ 2, 2], z = 0)
  
  # EV1:
  segments(x0 = ev1.proj$x, x1 = ev1.proj.neg$x, y0 = ev1.proj$y, y1 = ev1.proj.neg$y, lwd = 2, col = "blue")
  # EV2:
  segments(x0 = ev2.proj$x, x1 = ev2.proj.neg$x, y0 = ev2.proj$y, y1 = ev2.proj.neg$y, lwd = 2, col = "darkcyan")
  
  return(df)
}

Zeile 39 bis 51: Nur dieser Teil ist neu gegenüber der Implementierung von plot_QF_Hist(), die oben gezeigt wurde.

Zeile 40: So wie es oben gezeigt wurde, wird hier die Umgebung der quadratischen Form qf in eine Liste verwandelt und daraus die Variable ev geladen und wiederum in einer Variable ev abgespeichert. In ev (einer 2x2-Matrix) sind die Eigenvektoren der Matrix A als Spalten abgespeichert. Die Multiplikation mit dem Faktor 2 bewirkt, dass sie vergrößert dargestellt werden (mit eigen() werden normierte Eigenvektoren berechnet).

Zeile 42 bis 51: Hier werden die beiden Eigenvektoren verarbeitet, um sie in der xy-Ebene des Diagramms darzustellen. Es handelt sich hier um spezielle Befehle, die die Bibliothek scatterplot3d betreffen und die hier nicht erklärt werden. (Wie man xyz.convert() und segments() zum Zeichnen von Linien einsetzt wird in Erzeugen von dreidimensionalen Graphiken mit scatterplot3d erklärt.)

Der Aufruf von plot_QF_Hist() geschieht jetzt mit folgendem Skript, wodurch Abbildung 3 erzeugt wird:

A = matrix(data = c(1, -1, -1, 1), nrow = 2, byrow = TRUE)
qf <- factory_quadrForm(A, eigen = TRUE)
df <- plot_Q_Hist(qf)

Weiteres Beispiel: Erzeugen von mehreren quadratischen Formen

Zuletzt soll eine Anwendung gezeigt werden, die simuliert, dass mehrere, aber unterschiedliche quadratische Formen ausgewertet werden. Die Simulation enthält folgende Schritte:

1. Es werden 2x2-Matrizen mit zufälligen Eintragen erzeugt; die Einträge sind ganzzahlig. Allerdings werden ausschließlich symmetrische Matrizen erzeugt. Denn nur bei ihnen ist sichergestellt, dass die Eigenvektoren reell sind.
2. Zu jeder Matrix wird von der Funktionsfabrik die zugehörige quadratische Form erzeugt. Die quadratischen Formen werden in einer Liste abgespeichert.
3. An die Stelle der Auswertung tritt wieder der Aufruf der Funktion plot_QF_Hist().

Für den ersten genannten Punkt wird zunächst eine geeignete Funktion randomMatrix() implementiert. Dass man auch Funktionen wie "herkömmliche" Objekte behandeln kann, wurde in Spezielle selbstdefinierte Funktionen in R gezeigt. Es wird sich jetzt zeigen, dass sich Funktionsfabriken und Listen von Funktionen hervorragend ergänzen, da man die Auswertung mit einer einzigen lapply()-Anweisung erledigen kann.

Zunächst die Implementierung von randomMatrix(). Die Eingabewerte sind (siehe Zeile 1):

1. Die Anzahl nrow der Zeilen der Matrix; dies ist zugleich die Anzahl der Spalten, da eine quadratische Form eine quadratische Matrix erfordert.
2. Die ganzen Zahlen min = -1L und max = -min , für den Bereich, aus dem die Komponenten der Matrix ausgewählt werden.
3. Eine Zahl seed , mit der der Zufallsgenerator initialisiert werden kann. Gerade zum Testen ist es hilfreich, wenn die Ergebnisse reproduzierbar sind.

Der Rückgabewert ist dann eine symmetrische Matrix; die Symmetrie wird dadurch garantiert, dass zuerst eine Zufallsmatrix m gebildet wird und zuletzt die Matrix m und ihre transponierte Matrix t(m) addiert werden (siehe Zeile 9).

randomMatrix <- function(nrow = 2L, min = -1L, max = -min, seed = NULL){
  if(is.null(seed)){
    seed <- 42
  }
  set.seed(seed)
  
  m <- matrix(data = sample(x = seq(from = min, to = max, by = 1), size = nrow * nrow, replace = TRUE),
              nrow = nrow, byrow = TRUE)
  return( m + t(m) )
}

Zeile 2 bis 5: Durch seed kann der Zufallsgenerator entweder selbst initialisiert werden, oder er wird mit 42 initialisiert.

Zeile 7 und 8: Die Zufallsmatrix m wird mit Hilfe von sample() erzeugt. Dazu wird aus den Zahlen von min bis max zunächst ein Vektor der Länge nrow * nrow gebildet, der anschließend in eine Matrix verwandelt wird.

Zeile 9: Dadurch dass nicht die Matrix m sondern m + t(m) zurückgegeben wird, liegen die Einträge der Matrix jetzt zwischen 2 * min und 2 * max .

Damit lässt sich die oben beschriebene Simulation mit wenigen Anweisungen erzeugen:

QFs <- replicate(n = 4,
                expr = factory_quadrForm(randomMatrix(min = -2, max = 2, 
                seed = sample(x = (1:999), size = 1)),
                eigen = TRUE))

par(mfrow = c(2, 2))

dfs <- lapply(X = QFs, FUN = plot_QF_Hist)

Zeile 1 bis 4: Die Funktion replicate() sorgt dafür, dass die Anweisung im Argument expr viermal ausgeführt wird. In expr wird die Funktionsfabrik aufgerufen, der wiederum eine zufällig erzeugte Matrix übergeben wird. Die Funktionsfabrik wird mit dem Argument eigen = TRUE aufgerufen, da die Plots später die Eigenvektoren enthalten sollen. An der Konfiguration der Funktion randomMatrix() kann man leicht ablesen, welche Zufallsmatrizen erzeugt werden können. Die erzeugten quadratischen Formen werden in der Liste QFs abgespeichert.

Zeile 6: Die 4 Diagramme sollen in einem Quadrat angeordnet werden.

Zeile 8: Mit lapply() wird jetzt die Liste QFs abgearbeitet, wodurch 4 Plots erzeugt werden. In der Liste dfs sind jetzt die 4 Dataframes enthalten, die die Koordinaten der jeweils 360 ausgewerteten Punkte enthalten.

Abbildung 5 zeigt ein mögliches Ergebnis; die Matrizen und ihre Eigenwerte sind im Skript darunter angegeben.

Abbildung 5: Darstellung von vier zufällig erzeugten quadratischen Formen (Darstellung wie in Abbildung 3).

Die 4 Matrizen und ihre Eigenwerte lauten:

#      [,1] [,2]
# [1,]    4    1
# [2,]    1    4
# mit EW: 5 3
# 
#      [,1] [,2]
# [1,]    4    2
# [2,]    2    0
# mit EW: 4.8284271 -0.8284271
# 
#      [,1] [,2]
# [1,]   -4   -1
# [2,]   -1   -2
# mit EW: -1.585786 -4.414214
# 
#      [,1] [,2]
# [1,]   -4    2
# [2,]    2    0
# mit EW: 0.8284271 -4.8284271

Abbildung 6 zeigt die vier quadratischen Formen als "Gebirge" über der x₁-x₂-Ebene.

Abbildung 6: Die vier zufällig erzeugten quadratischen Formen aus Abbildung 5 als Gebirge über der x₁-x₂-Ebene. Die Skalierung der Achsen ist jetzt an die Flächen angepasst, so dass Abbildung 5 und 6 nur schwer zu vergleichen sind.