Bioinfo-skript_2007

modlab

Grundlagen der Bioinformatik

Ergänzende Unterlagen zu ausgewählten Kapiteln

Prof. Dr. G. Schneider erstellt u.M.v. Markus Hartenfeller, Jan A. Hiß, Michael Meissner, Ewgenij Proschak, Michael Schmuker, Yusuf Tanrikulu, Martin Weisel, Matthias Wirth

© modlab® 2007

© modlab® 2007

Inhaltsverzeichnis INHALTSVERZEICHNIS...................................................................................................... 3 SEQUENZALIGNMENT........................................................................................................ 5 PAARWEISES ALIGNMENT ....................................................................................................... 7 HEURISTIKEN FÜR LOKALES ALIGNMENT .............................................................................. 10 SIGNIFIKANZ LOKALER ALIGNMENTS .................................................................................... 14 MULTIPLES ALIGNMENT........................................................................................................ 15 Heuristiken für multiple Alignments ................................................................................ 17 LITERATUR ............................................................................................................................ 20 AUTOKORRELATIONSVEKTOREN............................................................................... 21 MATHEMATISCHE DEFINITION .............................................................................................. 22 TRANSFER IN DIE CHEMIEINFORMATIK .................................................................................. 23 Beispiel ............................................................................................................................. 25 ANWENDUNG VON KORRELATIONSVEKTOREN ...................................................................... 26 LITERATUR ............................................................................................................................ 27 NIPALS-ALGORITHMUS ................................................................................................... 29 PRINZIP UND SCHRITTE DES ALGORITHMUS .......................................................................... 30 LITERATUR ............................................................................................................................ 31 GRUNDLAGEN GRAPHENTHEORIE ............................................................................. 33 Graphen............................................................................................................................ 33 Subgraphen, Pfade, Zyklen, Cliquen ................................................................................ 35 GRAPHENTHEORETISCHE PROBLEME UND IHRE LÖSUNGEN IN DER CHEMINFORMATIK ......... 36 Subgraph-Isomorphismus................................................................................................. 36 Maximum Common Substructure ..................................................................................... 40 Cliquendetektion auf dem Assoziationsgraphen .............................................................. 41 LITERATUR ............................................................................................................................ 45 KÜNSTLICHE NEURONALE NETZE .............................................................................. 47 Genereller Aufbau von Feed-Forward-Netzen ................................................................ 47 DAS PERZEPTRON .................................................................................................................. 48 Das Single Layer Perzeptron ........................................................................................... 48 Die Lösung des XOR-Problems........................................................................................ 51 MEHRLAGIGE NETZE ............................................................................................................. 52 NETZTRAINING ...................................................................................................................... 53 PROBLEME UND LÖSUNGSANSÄTZE....................................................................................... 58 LITERATUR ............................................................................................................................ 61 PARTIKELSCHWARM-OPTIMIERUNG (PSO)............................................................. 63 DER ALGORITHMUS............................................................................................................... 64 ZUSAMMENFASSUNG ............................................................................................................. 66 LITERATUR ............................................................................................................................ 67 AMEISENALGORITHMEN................................................................................................ 69 AMEISENALGORITHMEN IN DER ANWENDUNG ...................................................................... 70 LITERATUR ............................................................................................................................ 72

© modlab® 2007

© modlab® 2007

Sequenzalignment

5

Sequenzalignment Die ersten Algorithmen zur Sequenzanalyse wurden in den 50er Jahren des 20. Jahrhunderts erdacht, als die ersten Proteinsequenzen bekannt wurden. Mit der von Fred Sanger entwickelten Methode zur Sequenzierung von DNA (Kettenabbruchreaktion, 1975) wurden zudem schnell große Anzahlen an Nukleotidsequenzen verfügbar. Die hier vorgestellten Alignment-Algorithmen dienen dem Vergleich und der Analyse von Sequenzalignments und werden verwendet, um die strukturelle, funktionelle oder evolutionäre Verwandtschaft von DNA- oder Proteinsequenzen aufzuklären. Das

Ausrichten

zweier

Sequenzen

gelingt

beim

paarweisen

Alignment

durch

Gegenüberstellung von gleichartigen Positionen (Matches) oder möglichst ähnlichen Positionen (Mismatches) und dem Einfügen von Lücken (Gaps) an den Stellen, an denen sich die Sequenzen durch Punktmutationen (Substitutionen, Insertionen oder Deletionen) unterscheiden. Dabei finden Scoring-Matrizen Anwendung, die überwiegend aus multiplen Alignments bekannter homologer Sequenzen abgeleitet wurden und konservierten oder sehr ähnlichen Positionen positive Werte (Scores) zuweisen, während unähnliche Positionen meist negative Werte erhalten (es gibt auch Scoring-Matrizen, die auf physikochemischen Eigenschaften beruhen und nicht aus Alignments abgeleitet wurden). Die Bestrafung von Lücken erfolgt durch die Gap-Funktion γ (Formel 1), die negativ in den Score eines Alignments eingeht. Dabei finden meist „affine“ Gap-Penalties Verwendung, bei denen ein langes Gap weniger schlecht bewertet wird, als mehrere kurze. Dies gelingt dadurch, dass die erste Stelle eines Gaps der Länge g mit einer Gap-Open-Penalty d stärker bestraft wird als die restlichen Positionen der gleichen Lücke, die mit der Gap-Extension-Penalty (Formel 1) e bestraft werden (d > e). © modlab® 2007

Sequenzalignment

γ(g) = d + (g-1) · e.

6

(1)

(Die lineare Gap-Penalty kann als Sonderfall der affinen Gap Penalty angesehen werden , wobei d = e.) Bewertung eines Alignments am Beispiel: Für die beiden DNA-Sequenzen TAAACAC und TGCATAC schlagen wir folgendes Alignment vor: TAA__ACAC T__GCATAC Wenn wir jeden Match mit +1 belohnen, jeden Mismatch mit –m und jeden Gap mit obiger affiner Gap-Penalty γ bestrafen, erhält unser Alignment einen Score von 4 – m –2d – 2e. Wie bereits erwähnt, werden üblicherweise Scoring-Matrizen verwendet, die Matches positive Scores zuweisen, während Mismatches meist negative Werte erhalten. Dabei sorgt eine Konvention dafür, dass ein Mismatch nicht strenger bestraft wird als zwei einzelne GapPositionen, durch die man den Mismatch vermeiden könnte. Während das optimale paarweise Alignment in Laufzeit O(n²) relativ schnell (aber oft immer noch zu langsam für eine effiziente Datenbanksuche) berechnet werden kann, ist das optimale multiple Alignment mit O(nk) (k = Anzahl der Sequenzen) nicht mehr effizient zu berechnen. Hier finden Heuristiken (CLUSTAL W, Feng & Doolittle Ansatz) Verwendung, die paarweises Alignment mit phylogenetischen Methoden kombinieren. So werden von den Sequenzen ausgehend phylogenetische Bäume erstellt und die einander ähnlichsten Sequenzen in paarweisen Alignments vereinigt, die wiederum zum fertigen multiplen Alignment zusammengesetzt werden. Im nun folgenden Abschnitt sollen Algorithmen zum paarweisen Alignment betrachtet werden.

© modlab® 2007

Sequenzalignment

7

Paarweises Alignment Globales Alignment Wir haben zwei Sequenzen x,y der Länge n,m gegeben und vermuten, dass diese miteinander verwandt sind. Um nicht alle exponentiell viele Alignments berechnen zu müssen, benutzen wir einen dynamischen Programmierungsansatz: Der Needleman-Wunsch-Algorithmus (1970) liefert uns optimale globale Alignment in Laufzeit O(n · m). Dazu erstellen wir uns einen (n × m) Alignmentgraphen, in dem wir in allen Feldern je drei Kanten beschriften (Abb. 1). Die mit Di,j beschrifteten diagonalen Kanten repräsentieren einen Match bzw. einen Mismatch zwischen den jeweiligen Positionen i,j in den Sequenzen x,y und werden mit den zugehörigen Scores S beschriftet. Die vertikalen bzw. horizontalen Kanten (V, H) repräsentieren jeweils einen Gap in der ersten bzw. zweiten Sequenz. Indem wir uns bereits errechnete Scores der besten Sequenzanfänge merken und mit dynamischen Programmieren weitergeben, erhalten wir die o.g. Laufzeitbeschleunigung. So trägt jede mit Di,j beschriftete Kante den Score des besten Alignments bis zu dieser Stelle (Formel 2, hier ausgehend von Di+1,j+1). Di+1,j+1 = S(xi+1, yj+1) + max {Di,j + Vi,j + Hi,j}

x1

x2

(2)

xn

H1,0

y1 V0,1 D1,1 xi y2

yj

Di,j Hi,j

yj+1

ym

xi+1

Vi,j

Hi+1,j Vi,j+1 Di+1,j+1

Dn,m Vn,m Hn,m

Abbildung 1: Die Kanten des Alignmentgraphen tragen den Score des besten Alignments bis zu dieser Stelle.

© modlab® 2007

Sequenzalignment

8

Bei den Rekursionen für horizontale und vertikale Kanten werden keine Scores verteilt, da Gaps keine Beiträge zum Alignment liefern. An dieser Stelle fließen jedoch die Gap-Openoder Gap-Extension-Penalties (d,e) ein (Formel 3). Hi+1,j = max {Di,j - d, Vi,j - d, Hi,j - e}

(3)

Nachdem wir alle Kanten beschriftet haben, erhalten wir den Score des optimalen globalen Alignments (Formel 4). S* = max {Dn,m, Vn,m, Hn,m}

(4)

Den dazugehörigen Pfad erhalten wir per Traceback über die maximalen Kanten. Das bedeutet wir laufen die Kanten rückwärts bis zum Start und entscheiden uns jeweils für die Kante, die den maximalen Beitrag zum Score geliefert hat (Abb. 2).

S*

Abbildung 2: Der Traceback über die Kanten mit den maximalen Beschriftungen liefert den Pfad zum optimalen Alignment. Ein Wort zur Laufzeit: Wir beschriften insgesamt 3·n·m Kanten und benötigen für den Traceback maximal n+m Kanten und bleiben mit 3·n·m +n+m asymptotisch in O(n·m).

Lokales Alignment Wenn wir vermuten, dass nur ein Teil (xk,...,xi) einer Sequenz x mit einem Teil (yl,...,yj) einer Sequenz y übereinstimmt, verwenden wir den Smith-Waterman-Algorithmus (1981), der das optimale lokale Alignment für die beiden Sequenzen in O(n·m) findet. Diese Suchstrategie wird z.B. verwendet, wenn wir eine Sequenz gegen eine ganze Datenbank alignieren möchten oder überprüfen wollen, ob z.B. ein Hefeprotein irgendwo im Humangenom vorkommt. © modlab® 2007

Sequenzalignment

9

Wir benutzen wieder die affine Gap-Penalty und beschriften die Kanten Di,j, Hi,j, Vi,j wie beim Needleman-Wunsch-Algorithmus, jedoch mit folgendem Unterschied: Wir geben jeder Kante die Möglichkeit, ein neues lokales Alignment mit Score ≥ 0 zu starten. Das bedeutet, dass eine Kante nicht ein möglicherweise negatives Erbe von schlechten Sequenzanfängen weiterführen muss, sondern stets ein neues lokales Alignment starten kann. Mit anderen Worten: Die Teile die außerhalb des optimalen lokalen Alignments liegen, gehen nicht negativ in den Score ein (Abb. 3).

xk...

...xi

y.l ..

... yj

Abbildung 3: Das optimale lokale Alignment – kein anderes lokales Alignment erreicht einen größeren Score. Die Rekursionen ändern sich wie folgt (Formeln 5, 6, 7): Di+1,j+1 = max {0, Di,j + S(xi+1, yj+1), Hi,j + S(xi+1, yj+1), Vi,j + S(xi+1, yj+1)} (5) Hi+1,j = max {0, Di,j - d, Hi,j - e, Vi,j - d}

(6)

Vi,j+1 = max {0, Di,j - d, Hi,j - d, Vi,j - e}

(7)

Das optimale lokale Alignment erhalten wir per Traceback über die Kanten, die den maximalen Beitrag geliefert hatten. Dabei haben wir uns während des Beschriftens der Kanten den maximalen Score gemerkt und gehen von dieser Kante (nicht von Dn,m, Hn,m oder Vn,m wie beim Needleman-Wunsch-Algorithmus) nach links oben vor, bis wir eine Kante erreichen, die den Score 0 trägt. Der Teil zwischen dem Score 0 und dem maximalen Score ist unser lokal optimales Alignment. Laufzeit und Speicherplatz betragen wie beim Needleman-Wunsch-Algorithmus O(n·m), was deutlich besser ist, als exponentielle Laufzeit, aber für die Suche in großen Datenbanken zu © modlab® 2007

Sequenzalignment

10

langsam ist. Hier finden die Heuristiken BLAST oder FASTA Anwendung, die zwar nicht unbedingt das optimale lokale Alignment finden, dafür aber deutlich schneller laufen und meist gute Ergebnisse liefern.

Heuristiken für lokales Alignment BLAST (Basic Local Alignment Search Tool) BLAST ist eine Heuristik für lokale Alignments, die mit einer Laufzeit von O(n+m) deutlich schneller ist als der Needleman-Wunsch-Algorithmus, jedoch nicht in jedem Fall das optimale lokale Alignment findet. Gegeben seien eine kürzere Sequenz S der Länge m und eine lange Sequenz T mit Länge n, gegen die S in einem (oder mehreren) lokalen Alignment(s) aligniert werden soll. Dazu wird S in kleine Teilwörter zerteilt (Default: Wortlänge |w| = 3 bei Proteinen, |w| = 11 bei Nukleinsäuren). Hierbei fügt BLAST auch ähnliche Wörter hinzu, deren Score einen bestimmten Schwellenwert t überschreiten (diese leicht veränderten Suchwörter stehen für Mutationen, die sich während des Evolutionsprozesses eingeschlichen haben könnten). Aus diesen Teilwörtern lässt sich ein Schlüsselwortbaum so konstruieren, dass die lange Sequenz T fast in linearer Zeit durch den Suchbaum aus den Fragmenten des Suchstrings S geschoben werden kann. Dies gelingt dadurch, dass die Suchmuster so in den Baum eingefügt werden, dass sie sich erst aufspalten, wenn sie sich unterscheiden. Z.B. für Sequenz S = KLLGV (Abb. 4). L

G

V

L

L

G

K

L

Abbildung 4: Schlüsselwortbaum für Sequenz S.

© modlab® 2007

Sequenzalignment

11

So kann man im Fall eines Mismatches bei der Mustersuche direkt an die richtige Stelle eines anderen

Musters

mit

gleichem

Präfix

springen.

Dies

bringt

die

gewünschte

Laufzeiteinsparung (Vergleiche hierzu auch weiterführende Literatur. Stichwort: „AhoCorasick-Algorithmus“). Wir schieben nun den langen Text T durch den Suchbaum und geben alle Matches aus (alle Fragmente mit Score ≥ t). Diese lassen sich in einem Alignmentgraphen (Dotplot) darstellen und repräsentieren Stücke, die in beiden Sequenzen identisch (bzw. ähnlich im Falle der zusätzlich eingefügten Strings mit Score ≥ t) sind. In einem weiteren Schritt verlängern wir diese gefundenen Alignments entlang der Diagonalen zu sogenannten HSPs (High Scoring Segment Pairs). Diese sind so definiert, dass sie sich durch Verlängern oder Verkürzen in ihrem Score nicht mehr verbessern lassen. Wir können uns nun eine Reihe guter lokaler Alignments ausgeben lassen, die einen bestimmten Mindestscore überschreiten oder einfach das MSP (Maximum Scoring Segment Pair) betrachten. Wie schon angedeutet findet BLAST nicht unbedingt das optimale lokale Alignment, liefert aber meist sehr gute Ergebnisse. Die Laufzeit ist mit O(n+m) deutlich schneller als beim Smith-Waterman-Algorithmus. Anmerkung: In der klassischen BLAST Version werden HSPs nur entlang der Diagonalen verlängert. Beim sog. GAPPED BLAST werden nur solche HSPs betrachtet, die (innerhalb eines begrenzten Abstands zueinander) auf einer gemeinsamen Diagonalen liegen, dass heißt einem gemeinsamen gaplosen Alignment angehören. Zwischen diesen HSPs suchen wir nun nach örtlichen Exkursionen, die unsere beiden HSPs mit einem weiteren HSP über Gaps verbinden kann (Abb. 5). Sequenz 1

Abbildung 5: Gapped BLAST: einzelne HSPs werden über Gaps verbunden, wenn sich dadurch der Score verbessern lässt. © modlab® 2007

Sequenzalignment

12

Um weiter entfernte Sequenzen zu finden empfiehlt sich PSI-BLAST (Position Sensitive Iterated-BLAST). Hierbei wird aus signifikanten Treffern einer BLAST-Suche ein multiples Alignment erstellt, aus dem ein Profil abgeleitet werden kann. Im Profil sind die jeweiligen Wahrscheinlichkeiten für jedes Nukleotid bzw. jede Aminosäure für die jeweilige Position im Alignment vermerkt. Dieses Profil wird nun zur Datenbanksuche eingesetzt, die resultierenden signifikanten Treffer wieder aligniert und ein neues (noch variableres) Profil erstellt. Nach mehreren Iterationen endet das Verfahren. PSI-BLAST findet auch weiter entfernte Sequenzen (auch < 30% Identität). Die wichtigsten Datenbanken (Genbank, EMBL und DDBJ) lassen sich allesamt mit BLAST durchsuchen. Die zahlreichen BLAST-Varianten sind in der BLAST-Suite zusammengefasst:

http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#pstab Die wichtigsten Varianten: blastn: zur Suche in Nukleotidsequenzdatenbanken. blastp: zur Suche in Proteinsequenzdatenbanken. blastx: Vergleicht eine Nukleotidsequenz (in allen Leserastern translatiert) gegen eine Proteindatenbank. tblastn: Vergleicht eine Proteinsequenz gegen eine Nukleotiddatenbank (dynamisch in allen Leserastern translatiert) tblastx: Vergleicht die six-frame-Translation einer Nukleotidsequenz gegen die six-frame-Translationen einer Nukleotidsequenzdatenbank.

FASTA Auch hier soll ein kurzer Text S in einem langen Text T gefunden werden. Wir gehen ähnlich wie bei BLAST vor und zerteilen Text S in kleinere Bruchstücke. Der Parameter „ktup“ bestimmt hier die Länge der Bruchstücke, die gematcht werden sollen (Default: ktup = 6 für DNA, 2 für Proteine). Da die Bruchstücke hier kleiner sind als bei BLAST, suchen wir nicht mit Pattern-Matching-Algorithmen nach Übereinstimmungen, sondern mit dem sog. 'lookup © modlab® 2007

Sequenzalignment

13

table', einer Hashing-Tabelle. Im Lookup-Table wird das Auftreten aller Aminosäuren bzw. Nukleotide in der Sequenz T aufgeführt, z.B. für die Sequenz T= A1G2T3A4C5A6 (Tab. 1). Tabelle 1: Lookup-Table für Sequenz T Symbol A G T C

Position 1,4,6 2 3 5

Beginnt ein Suchstring der Länge ktup aus S mit einem 'A', erkennen wir im lookup-table, dass im Text T ein 'A' an den Positionen 1,4 und 6 auftritt. Steht in S an zweiter Postion ein 'G', prüfen wir, ob es ein 'G' in T an den Positionen 2,5 oder 7 gibt (also ein Zeichen nach den schon gefundenen 'A's.), usw. Man kann so sehr schnell alle "hotspots" aufspüren. Das sind genau die Stellen, an denen T und S an mindestens ktup Positionen übereinstimmen. Die Übereinstimmungen tragen wir in einen Dotplot ein. Wir können nun die Init1-Scores berechnen, d.h. wir verlängern alle hotspots zu HSPs (s.o.) entlang der Diagonalen. Danach berechnen wir die Initn-Scores, dass heißt, wir versuchen lokal die besten n HSPs durch Gaps zu verbinden, s.d. ein noch höherer Alignment-Score resultiert. Die Breite der Nachbarschaftssuche wird dabei auf einen bestimmten Rahmen beschränkt: Ausgehend von den Init1's betrachten wir die umliegenden 32 Diagonalen (16 höher, 16 tiefer). Ein besonderes Merkmal von FASTA ist die anschließende "Exakte Scorebestimmung": Für den Bereich in dem das Alignment mit dem höchsten Initn-Score wird nochmal ein exaktes lokales Alignment mit dem Smith-Waterman-Algorithmus berechnet. Dabei ist der betrachtete Rahmen noch etwas größer als der Bereich in dem das beste Initn-Alignment liegt, um eventuelle Fehler von FASTA auszugleichen (Abb. 6). Fazit: BLAST vs. FASTA BLAST läuft schneller, weil die Suchkriterien hier strenger sind als bei FASTA (FASTA findet mehr Treffer, wg. Default ktup = 2, während |w| = 3 bei BLAST (für Proteine)). Dafür findet FASTA aber mit höherer Wahrscheinlichkeit das optimale Alignment, durch Anwendung des Smith-Waterman-Algorithmus. © modlab® 2007

Sequenzalignment

14

Sequenz 1 Suchraum um Init1's Smith-Waterman

Abbildung 6: Exakte Scorebestimmung von FASTA.

Signifikanz lokaler Alignments Die Signifikanz eines Alignments drückt aus, wie aussagekräftig dessen Score tatsächlich ist: Ist der Score das Ergebnis eines zufälligen Prozesses oder ist der Score so hoch, dass wir ausschließen können, dass die Ähnlichkeit der Sequenzen zufällig ist? In diesem Fall gehen wir davon aus, dass die Sequenzen homolog sind. Hier sollen mit dem Z-Wert, P-Wert und dem E-Wert drei Maße für Signifikanz vorgestellt werden.

Z-Wert Ist der Score S unseres Alignments besser als der einer zufällig (aus den gleichen Bausteinen) zusammengewürfelten Sequenz? Dazu berechnen wir üblicherweise 100 neue Alignments, dazu den Mittelwert ihrer Scores und die Standardabweichung (Formel 8). Z-Score von S = (S – Mittelwert) / Standardabweichung

(8)

Ein Z-Wert von ≥ 5 gilt als signifikant. Dies tritt ein, wenn der beobachtete Score deutlich höher ist, als der Mittelwert der zufällig erstellten Sequenzen. Je höher ein Z-Score, desto geringer ist dabei die Wahrscheinlichkeit, dass der beobachtete Alignmentscore nur Zufall ist. Ein Z-Score von 0 bedeutet, dass der beobachtete Score nicht größer ist als der einer „zufälligen Sequenz“. © modlab® 2007

Sequenzalignment

15

P-Wert Hier berechnen wir die Wahrscheinlichkeit für das Ereignis, dass ein so hoher Score wie der beobachtete Score S zufällig zustande kam (Nullhypothese). Fällt die Wahrscheinlichkeit für alle Alignments der gleichen Länge, die einen Score S erreichen oder übertreffen unter das Signifikanzniveau von 5%, verwerfen wir die Nullhypothese und gehen von Homologie der Sequenzen aus. Das heißt, die Ähnlichkeit der betrachteten Sequenzen ergibt sich durch gemeinsame Abstammung und nicht durch Zufall.

E-Wert Der E-Wert (Formel 9) gibt den Erwartungswert für ein unabhängiges Alignment an, dass dieses einen Score ≥ S erreicht. E-Wert für S = n · m · k · e-λS

.

(9)

λ und k sind hier Skalierungsfaktoren für die Größe des Suchraums und die gewählte Scoringmatrix. Der E-Wert wird besonders häufig zur Bestimmung der Signifikanz von homologen Alignments verwendet, da er (anders als Z- und P-Wert) Sequenzlängen berücksichtigt. Das bedeutet, liefern zwei Sequenzen A,B den gleichen Score S und A ist kürzer als B, so hat B einen deutlich höheren E-Wert (ist also weniger signifikant), da es ja mehr Positionen benötigt hat, um den Score S zu erreichen. Die unterschiedlichen Sequenzlängen gehen mit (n·m) in die Formel des E-Wertes ein.

Multiples Alignment Das multiple Alignment liefert im Gegensatz zum paarweisen Alignment genauere Informationen über Aminosäureverteilungen an bestimmten Positionen, da sich hieraus ableiten lässt, wie stark konserviert eine Position tatsächlich ist. Darüberhinaus geben sie Aufschluss darüber, wie wahrscheinlich der Austausch einer bestimmten Aminosäure durch © modlab® 2007

Sequenzalignment

16

eine andere ist. Hierauf beruhen z.B. die PAM-Scoring-Matrizen, die sich aus den Austauschwahrscheinlichkeiten von Aminosäuren aus Alignmentblöcken errechnen lassen. Zudem lassen sich aus multiplen Alignments von Proteinen einer Familie Profile ableiten, mit denen man nach weiteren Vertretern dieser Familie suchen kann. Während für paarweises Alignment effiziente Algorithmen vorhanden sind, ist die Berechnung des optimalen Alignments für mehr als zwei Sequenzen nicht mehr effizient möglich. Hier müssten wir den Weg mit optimalem Score durch eine n-dimensionale Matrix finden, was exponentielle Laufzeit benötigt. Abhilfe schaffen hier Heuristiken für multiple Alignments.

Scores für multiples Alignment Der Gesamtscore eines multiplen Alignments m errechnet sich als die Summe der Scores seiner Spalten i. Zusätzlich gehen hier die Gap-Penalties G negativ ein (Formel 10).

S(m) = G + ∑ S(mi )

(10)

i

Die Scores der einzelnen Spalten lassen sich auf verschiedene Weise berechnen:

Sum of Pair Scores Hierzu berechnen wir die Scores aller n! möglichen Paare einer Alignmentposition (z.B. mit PAM- oder BLOSUM-Matrizen). Ein Problem beim Sum of Pair Score ist, dass kleine evolutionäre Ereignisse überbewertet werden. Ist nämlich nur eine Position einer sonst konservierten Spalte mutiert, so geht der (meist negative) Score für den Mismatch mit allen anderen (n-1) Positionen n-1 mal in den Spaltenscore ein.

Entropie-Scores •

Shannon-Entropie (Formel 11) | A|

H = −∑ p k log 2 p k , wobei k=1

© modlab® 2007

(11)

Sequenzalignment

17

|A| = Die Länge des Alphabets (z.B. 20 Aminosäuren) pk = Die relative Häufigkeit des Elementes k in der betrachteten Spalte Die Shannon-Entropie ist ein Minimum-Entropie-Score, d. h. eine durchweg konservierte Spalte (Optimum) erhält den Score 0, während der Score bei immer größer werdender Unordnung (Entropie) immer größere (positive!) Werte annimmt (das Maximum ist abhängig von der Größe des Alphabets). Der sog. Information-Content I (Formel 12) kann als die „reverse“ Shannon-Entropie aufgefasst werden. Er drückt aus, wie aussagekräftig eine Spalte im multiplen Alignment ist (je höher der Wert, umso besser). | A|

I = H max − H act = log 2 | A | +∑ p k log 2 p k .

(12)

k =1

•

Kullbach-Leibler-Distanz (Formel 13)

H(P || Q) = ∑ p k log k

pk qk

(13)

Die Kullbach-Leibler-Distanz beschreibt die „Distanz“ zwischen der beobachteten Häufigkeit pk (also der relativen Häufigkeit einer Aminosäure oder eines Nukleotids k in der Spalte) und der Hintergrundverteilung qk (der erwarteten Wahrscheinlichkeit für das Auftreten einer Aminosäure oder eines Nukleotids). Mit anderen Worten: Die durch einen Mismatch in der Spalte hervorgerufene Entropie wird noch erhöht, wenn die betreffende Position (AS oder Nukleotid) besonders selten auftritt.

Heuristiken für multiple Alignments Die Berechnung des optimalen multiplen Alignments von n Sequenzen ist zeitaufwendig. So ließe sich zwar in einem n-dimensionalen Ansatz ähnlich zum Needleman-WunschAlgorithmus durch die Berechnung der jeweils besten Scores in einer n-dimensionalen Matrix durch dynamisches Programmieren das optimale multiple Alignment berechnen, jedoch wäre © modlab® 2007

Sequenzalignment

18

die Laufzeit immer noch exponentiell. Das Problem ist, das jede der n Positionen in einer Spalte entweder ein Gap oder eine tatsächliche Position (AS oder Nukleotid) sein könnte. Wenn wir verbieten, dass alle Positionen einer Spalte Gaps sein können, beschert uns dies immer noch 2n-1 mögliche Kombinationen, deren Scores es über m Spalten hinweg zu maximieren gilt. Es gibt Branch & Bound-Ansätze, die nicht die Sequenzanfänge aller Alignments berechnen (wie im Brute-Force-Ansatz), sondern unwahrscheinliche Alignments von der Berechnung ausschließen. Diese Methoden benötigen jedoch immer noch exponentielle Laufzeit, auch wenn sie um eine Konstante besser sind. Nachfolgend sollen zwei Heuristiken vorgestellt werden, die zwar einen Verlust der Genauigkeit in Kauf nehmen, dabei aber deutlich schneller sind als Brute-Force-Ansätze und dennoch gute Ergebnisse bringen.

Feng-Doolittle-Algorithmus Vorgehen: 1.

Aligniere alle Sequenzen paarweise mit dem Needleman-Wunsch-Algorithmus.

2.

Weise allen n-1 Sequenzpaaren Distanzen zu. Der Feng-Doolittle-Algorithmus berechnet hierzu den sog. „difference score“. Dieser entspricht in etwa dem Reversal des erhaltenen Scores eines paarweisen Alignments, d.h. je höher der Score des Alignments, desto niedriger der difference score. So wird eine hohe Ähnlichkeit zwischen zwei Sequenzen in eine niedrige Distanz übersetzt.

3.

Berechne aus den Distanzen einen phylogenetischen Baum mit UPGMA.

4.

Aligniere längs der Kanten (ausgehend von den Blättern) Sequenz-Paare, SequenzBlock-Paare und Block-Paare (Abb. 7).

Seq 1

Seq 2

Seq 3

Seq 4

Seq 5

Block (1,2) Block (1,2,3) Block (4,5)

Block (1,2,3,4,5)

Abbildung 7: Multiples Alignment mit Feng-Doolittle: Aligniere zunächst Blätter untereinander, dann Blätter gegen Blöcke und schließlich Blöcke untereinander. © modlab® 2007

Sequenzalignment

19

Zunächst werden je zwei Sequenzen zu einem Block zusammengefügt, und gegen diesen Block können später weiter Sequenzen oder andere Blöcke aligniert werden. Ein einmal in einem Block fixiertes Alignment bleibt für das multiple Alignment mit all seinen Gaps nach der Regel „once gap, always gap“ erhalten. Mit einer Ausnahme: In den fixierten Block kann ein globales Gap eingefügt werden, um eine weitere Sequenz gegen den Block zu alignieren. Eine Sequenz, die neu zu einem Block hinzugefügt werden soll, wird zunächst mit allen Sequenzen des Blocks paarweise aligniert (n-1 viele paarweise Alignments). Die Sequenz wird dann gemäß des besten paarweisen Alignments in den Block eingefügt. Als Ausgabe erhalten wir einen Block, der das Alignment aller unserer Sequenzen beinhaltet, also das multiple Alignment.

CLUSTAL W Clustal W verfolgt eine ähnliche Strategie wie Feng-Doolittle. Deshalb sollen hier nur die wesentlichen Unterschiede vorgestellt werden. Zunächst werden auch hier alle paarweisen Alignments berechnet und Distanzen zugewiesen. Dabei kennt CLUSTAL W neben dem difference score aus dem Alignment Score wahlweise noch einen heuristischen Ansatz zur Berechung der Distanzen. So kann die Distanz in einem beschleunigten Verfahren über die Anzahl der identischen Treffer bei einer ktup-Suche (siehe FASTA) berechnet werden. Aus den Distanzen wird wieder ein phylogenetischer Baum aufgebaut, hier aber nach der Neighbor-Joining-Methode. Ausgehend von den Blättern werden nun wieder zunächst Sequenzen miteinander aligniert. Diese werden aber nicht zu fixierten Blöcken zusammengefasst, sondern an den Knoten zu Profilen geclustert. Ein Profil ist eine Matrix, die für jede der 20 Aminosäuren (oder 4 Nukleotide) die relative Häufigkeit ihres Auftretens an dieser Stelle im multiplen Alignment beinhaltet. Um zu verhindern, dass eng verwandte Sequenzen das Profil in ihre Richtung verzerren, wird ihr Einfluss auf das Profil rechnerisch herabgesetzt. Als Ausgabe erhalten wir ein Profil, aus dem sich ablesen lässt, welche Stellen im Alignment hoch konserviert sind (evtl. funktionstragend) und welche eher nicht.

© modlab® 2007

Sequenzalignment

20

Literatur Needleman SB, Wunsch CD (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. J Mol Biol. 48(3): 443-53. Smith TF, Waterman MS (1981). Identification of common molecular subsequences. J Mol Biol. 147(1): 195-7. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990). Basic local alignment search tool. J Mol Biol. 215(3): 403-10. Pearson WR (1994). Using the FASTA program to search protein and DNA sequence databases. Methods Mol Biol. 24: 307-31. Lesk A. (2002). Introduction to Bioinformatics. Oxford University Press, Oxford. Feng DF, Doolittle RF (1990). Progressive alignment and phylogenetic tree construction of protein sequences. Methods Enzymol. 183: 375-87. Thompson JD, Higgins DG, Gibson TJ (1994). CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22(22): 4673-80.

© modlab® 2007

Autokorrelationsvektoren

21

Autokorrelationsvektoren Molekulare Deskriptoren sind im Wirkstoffdesign und für QSAR Modelle unentbehrlich, um physikochemische Eigenschaften von Molekülen mit mathematischen Methoden zu beschreiben. Sie ermöglichen es, umfangreiche Datenbanken mit Millionen von Verbindungen digital abzuspeichern und chemoinformatische Methoden anzuwenden. Autokorrelationsvektoren stellen nur eine Klasse der molekularen Deskriptoren dar. Durch die Zusammenfassung von Korrelationswerten zu Vektoren lässt sich ein molekularer Deskriptor erstellen, der Ähnlichkeitsberechnungen von Substanzen auf Basis der Verteilung von Atomeigenschaften und deren Abständen zueinander zulässt. Hierbei erfolgt eine Abstraktion der chemischen Bindungen zwischen den Atomen eines Moleküls und der Molekülstruktur im kartesischen Raum. Grundlegend gibt es zwei Anforderungen an einen Deskriptor. Zum einen soll er den Vergleich von Molekülen zulassen, und zum anderen soll er praktisch anwendbar sein und eindeutig aus den bekannten Daten berechnet werden können. Beides trifft auf Korrelationsvektor-basierte Deskriptoren zu: a) die Darstellung als Vektor ermöglicht eine einfache mathematische Quantifizierung der Ähnlichkeit von Molekülen (durch Vektoralgebra), und b) die minimale Anforderung zur Erstellung eines Korrelationsvektors ist die zweidimensionale Molekülstruktur, die meist vorhanden ist oder leicht bestimmt werden kann. Bei der Verwendung von Autokorrelationsvektoren als molekularen Deskriptor sollten auch die Nachteile stets bedacht werden. Erstens ist hier die Vernachlässigung der vollständigen Molekülstruktur inklusive wichtiger dreidimensionaler Informationen zu nennen. Außerdem © modlab® 2007

Autokorrelationsvektoren

22

ist keine Rekonstruktion der Molekülstruktur anhand des zugehörigen Korrelationsvektors möglich. Der Abstraktionsgrad des Deskriptors ist zu groß, als dass eine bijektive Darstellung möglich wäre.

Mathematische Definition Die Wurzeln der Autokorrelationsfunktionen liegen in der Funktionsanalyse von physikalischen Methoden. Im Allgemeinen hat sie die Form:

b

AC t = ∫ f ( x) ⋅ f ( x + t )dx .

(1)

a

Dabei stellt die Funktion f(x) eine beliebige Funktion dar, die auf dem Intervall [a, b+t] definiert ist. Die Berechnung des Korrelationswertes an der Stelle t lässt sich folgendermaßen beschreiben (Abbildung 1): Gegeben seien zwei Punkte in dem Intervall [a, b+t], deren Abstand zueinander t beträgt. Beide Punkte laufen nun unter Beibehaltung des Abstands über das Intervall. Für jeden Zeitpunkt berechnet die Autokorrelationsfunktion das Produkt der Funktionswerte an beiden Punkten und bildet das Integral über das gesamte Intervall. Als Ergebnis geht aus der Rechnung ein Wert hervor, der nicht von einem Funktionswert an einer bestimmten Stelle x abhängt, sondern vielmehr von der Verteilung der Funktionswerte im Intervall zu einem gegebenen Zeitpunkt.

Abbildung 1: Darstellung von zwei möglichen Punkten im Abstand t einer gegebenen Funktion f(x), die autokorreliert werden können. © modlab® 2007

Autokorrelationsvektoren

23

Diese Tatsache ermöglicht vor allem die Anwendung von Autokorrelationsfunktionen in der Analyse von physikalischen Spektren. Beispiele hierfür sind zum Beispiel NMRSpektrenauswertungen und Spracherkennungsprozesse. In den meisten Fällen ist zu einem gegebenen Spektrum keine Funktion vorhanden, sondern es liegt nur der Graph vor. Die Stärke der Autokorrelationsfunktion ist die Erkennung von sich periodisch wiederholenden

Signalen. Da die Produkte aus den Spitzen mit einem Abstand t in das Integral mit einfließen, machen sich wiederkehrende Elemente mit einer Periodenlänge t durch Maxima in der Autokorrelationsfunktion bemerkbar. Bei einer unterschiedlichen Periode wird der Funktionswert an der Spitze der Funktion mit einem möglicherweise kleineren Wert multipliziert, woraufhin die Wertigkeit des gesamten Integrals kleiner wird. Diese Eigenschaft ist es, die die Verwendung von Autokorrelationsfunktionen insbesondere für das Filtern von stark verrauschten Signalen auszeichnet.

Transfer in die Chemieinformatik Moreau und Broto waren 1980 die ersten, die die Autokorrelationsfunktionen in die Welt der Chemieinformatik einführten. Die Hauptintention war es, Aussagen über die Verteilung von Atomeigenschaften, die als numerische Werte vorlagen, im Molekülgraphen machen zu können. Der Transfer der klassischen Autokorrelations-funktionen auf ein als Strukturformel vorliegendes Molekül wurde in folgender Weise vollzogen:

•

Anstatt der Stelle x auf einem Funktionsgraphen wurden die Atome im Molekül als zu betrachtende Punkte verwendet. Die Wasserstoffatome wurden hierbei nicht in Betracht gezogen. Alle anderen Atome wurden beliebig durchnummeriert, um eine eindeutige Identifizierung der Punkte zu ermöglichen. Das Intervall [a, b+t] wurde also durch das gesamte Molekül dargestellt. Aufgrund der Tatsache, dass es nur eine endliche Anzahl von Punkten (Atome im Molekül) gibt, und keine kontinuierliche Funktion vorliegt, konnte das Integral der Funktion durch die Aufsummierung der Korrelationswerte ersetzt werden.

•

Als Funktionswerte f(x) wurden die Atomeigenschaften des Atoms x betrachtet. Meist wurde

die

Funktion

als

p(x),

wie

property,

bezeichnet.

In

Autokorrelationsfunktion wirdnur eine Atomeigenschaft mit sich selbst korreliert.

© modlab® 2007

jeder

Autokorrelationsvektoren •

24

Schließlich wurde der Abstand t als die minimale Anzahl der Bindungen zwischen zwei Atomen definiert (topologische Distanz). Hierbei wird nicht zwischen Doppelund Einfachbindungen unterschieden.

• Die auf diese Weise modifizierte Autokorrelationsfunktion entspricht der Formel: AC d =

∑ p(i) ⋅ p( j )

(2)

i , j∈M ( d )

wobei M(d) als die Menge der Atomtupel (i,j) mit einer topologischen Distanz d zueinander definiert ist. Gleichung (2) lässt nun die Berechnung von Korrelationswerten für jede gegebene Distanz und jede bekannte Atomeigenschaft zu. Es wird lediglich die Strukturformel des zu kodierenden Moleküls benötigt. Betrachtet man zum Beispiel Distanzen von 0 bis 10 Bindungen und die Autokorrelation einer Atomeigenschaft, so hat der als Deskriptor verwendete Korrelationsvektor folgende Datenstruktur: AV = (AC0, AC1, AC2, AC3, AC4, AC5, AC6, AC7, AC8, AC9, AC10) Der Deskriptor entspricht der Aneinanderreihung der Korrelationswerte zu den auftretenden Atomeigenschaft bei einer Distanz d. Erfahrungsgemäß kann (2) aufgrund der Aufsummierung zu besonders großen Korrelationswerten führen, wenn eine große Anzahl von Atompaaren mit einer Distanz d vorliegt. Dies kann sich gerade bei dem Vergleich von unterschiedlich großen Molekülen negativ auswirken, weshalb dies meist durch eine Mittelung über die Anzahl der Atompaare eliminiert wird. Die erweiterte Formel lautet nun:

AC d =

1 ∑ p(i) ⋅ p( j ) , n i , j∈M ( d )

wobei n der Anzahl der Atompaare mit einer Distanz d entspricht.

© modlab® 2007

(3)

Autokorrelationsvektoren

25

Beispiel In der Tabelle 1 soll die Berechnung des Autokorrelationsvektors für zwei Moleküle A und B beispielhaft gezeigt werden. Berechnet werden die Autokorrelationswerte für die Elektronegativität der vorliegenden Atome (ENC=2,5; ENO = 3,5) mit einer Distanz von 1 bis 4 Bindungen unter Benutzung der Formel (3).

Tabelle 1: Berechnung der Korrelationswerte der Elektronegativität für zwei Moleküle A und B. Verbindung A B O

OH

Molekülstruktur

OH

OH

Butan-1,3-diol

Name

Vereinfachte

Struktur

6 O

mit

willkürlichem Atomindex

Buttersäure

2

1

6 O 3

5 4 O

1

2

AC1

7,25

7,25

AC2

7,75

8,45

AC3

4,75

4,75

AC4

3,7

2,5

3

5 4 O

Als Ergebnis erhält man folgende Vektoren: VerbindungA

7,25

7,75

4,75

3,7

VerbindungB

7,25

8,45

4,75

2,5

Die Deskriptoren unterscheiden sich nur in den Funktionswerten für die Distanzen 2 und 4, weil immer bei derjenigen Verbindung ein größerer Wert auftritt, bei dem die Sauerstoffe in dem entsprechenden Abstand zueinander vorlagen. Eine Interpretation der resultierenden Autokorrelationsvektoren ist in folgender Weise möglich: In der Verbindung A liegen mehr

© modlab® 2007

Autokorrelationsvektoren

26

Atome mit der betrachteten Eigenschaft in einem Abstand von 4 Bindungen vor, als in der Verbindung B. Hier tritt dies bei einer Distanz von 2 Bindungen auf. Die betrachteten Eigenschaften müssen nicht direkt miteinander korreliert werden. In bestimmten Fällen ist es ratsam modifizierte Werte der numerischen Eigenschaften zu multiplizieren. Hier können unter anderem der Betrag, der Logarithmus u.ä. verwendet werden. Oft ist zudem die Verteilung von mehreren Atomeigenschaften auf den Molekülen zur Beschreibung

der

biologischen

Aktivität

wichtig.

Hier

bedient

man

sich

der

Kreuzkorrelation, bei der verschiedene Eigenschaften miteinander korreliert werden. Die verwendete Formel lässt sich in folgender Weise zusammenfassen:

AC d =

1 n

∑ ∑ p(i) ⋅ q( j ) ,

(4)

p , q∈P i , j∈M ( d )

wobei die Funktionen p und q den Elementen aus einer Funktionsmenge P mit allen betrachteten Eigenschaften entspricht. Autokorrelationsvektoren werden nicht nur zur Kodierung von zweidimensionalen Molekülstrukturen verwendet. Eine Erweiterung der Methode zur Benutzung von 3-D Molekülinformationen ist durch ein simples Austauschen der topologischen Distanz mit der euklidischen Distanz möglich. Die Korrelationsvektoren kodieren dann für eine Beschreibung der Eigenschaftsverteilungen im dreidimensionalen Raum, die durch das Molekül gegeben ist, unter Berücksichtigung der Atomdistanzen in Ångstrom.

Anwendung von Korrelationsvektoren Die Kodierung von Korrelationsvektor-basierten molekularen Deskriptoren liefert gleichgroße Vektoren zur Beschreibung der vorliegenden Substanzen. Bei der Betrachtung von k Atomeigenschaften und l Distanzen erhält man einen (k ⋅ l ) -dimensionalen Deskriptor für jedes Molekül. Zur Berechnung der Ähnlichkeit von zwei Molekülen bedient man sich einfachen mathematischen Vektorfunktionen, wie der euklidischen Norm zweier Vektoren:

v v a −b =

n

∑ (a − b ) i

i

2

,

i =1

wobei n der Dimensionalität des Deskriptors entspricht.

© modlab® 2007

(5)

Autokorrelationsvektoren

27

Falls man mehrere Eigenschaften bei der Korrelation in Betracht zieht, können Fälle auftreten, die sich aufgrund der Werte je Eigenschaft in (5) stärker auf die Ähnlichkeitsanalysen auswirken. Ein einfaches Beispiel bildet zum Beispiel die Korrelation der Elektronegativität von Atomen im Vergleich zu Korrelationen der Atommasse, welche wesentlich höher ausfallen kann. Hier bietet es sich an, nach der Berechnung des Korrelationsvektors auf eine Normierung des Vektors zurückzugreifen (z.B. eine Skalierung der Werte von 0 bis 1). Neben Ähnlichkeitsanalysen werden Korrelationsvektor-basierte Deskriptoren für viele weitere chemieinformatische Methoden verwendet.

Literatur Moreau, G., Broto, P., 1980, Autocorrelation of Molecular Structures, Application to SAR Studies. Nouv. J. Chim., 4, 757. Wagener, M., Sadowski, J., Gasteiger, J., 1995, Autocorrelation of Molecular Surface Properties for Modelung Corticosteroid Binding Globulin and Cytosolic Ah Receptor Activity by Neural Networks. J. Am. Chem. Soc., 117, 7769. Schneider, G., Clément-Chomienne, O., Hilfiger, L., Schneider, P., Kirsch, S., Böhm, H.-J., Neidhart, W., 2000, Evolutionäres De-novo-Design bioaktiver Moleküle: ein Ansatz zum virtuellen Screening. Angew. Chem., 112, 4305. Schneider, G., Wrede, P., 1998, Artificial neural networks for computer-based molecular design. Prog. Biophys. Mol. Biol., 70, 175.

© modlab® 2007

Autokorrelationsvektoren

© modlab® 2007

28

NIPALS-Algorithmus

29

NIPALS-Algorithmus Die Hauptkomponentenanalyse (Principal Component Analysis, im Folgenden mit „PCA“ abgekürzt) ist eine statistische Methode zur Analyse von multidimensionalen Daten. Mit Hilfe dieser Technik ist es unter anderem möglich, die Dimensionalität eines Datensatzes mit einem definiertem Informationsverlust zu reduzieren. Die Basis der PCA ist die Berechung von

Eigenvektoren. Diese haben als Eigenschaft, dass sie alle orthogonal zu einander sind. Sie treten immer paarweise mit dem dazugehörigen Eigenwert auf, anhand welcher man die Eigenvektoren sortieren kann. Der Eigenvektor mit dem größten Eigenwert ist als die erste

Hauptkomponente definiert. Mit ihr ist es möglich, die stärkste Beziehung in der Datenmatrix aufzuzeigen. Aufgrund dieser Eigenschaften ist es möglich, die Daten einer Matrix durch Ihre Projektion auf die Hauptkomponenten (Eigenvektoren mit den größten Eigenwerten) darzustellen. Bei einer geschickten Auswahl von Eigenvektoren können so dieselben Daten mit einer geringeren Dimensionalität wiedergegeben werden. In der Regel wählt man hierzu alle Eigenvektoren mit einem Eigenwert > 1 aus. Zur Berechnung der Eigenvektoren greifen viele Anwender auf die integrierten Lösungen der verwendeten Programmpakete, wie Java oder Matlab, zu. Hier werden meist numerisch exakte, mathematische Verfahren wie der QZ-Algorithmus, die Potenzmethode, die Inverse Iteration, das Lanczos-, Arnoldi- oder Jacobi-Davidson-Verfahren verwendet. Alle numerischen Verfahren haben gemeinsam, dass die Komplexität bei der Anwendung auf großen Matrizen enorm steigt. Dadurch wird die Berechnung der Eigenvektoren sehr aufwändig und die Methoden nicht praktikabel in der Anwendung. Eine

Alternative

bilden

approximative

Lösungen.

Ein

populäres

Beispiel

zur

Eigenvektorabschätzung ist der 1966 von Herman Wold entwickelte NIPALS-Algorithmus (nonlinear iterative partial least squares), welcher im Folgenden erklärt wird. Der NIPALS© modlab® 2007

NIPALS-Algorithmus Algorithmus

wird

30 nicht

nur

in

der

Hauptkomponentenanalyse

angewendet.

Weiterentwicklungen und Derivate des Algorithmus finden sich auch in der multiplen Regression, der Faktorzerlegung und der kanonischen Korrelationsanalyse wieder.

Prinzip und Schritte des Algorithmus Die Hauptkomponenten werden nacheinander durch die iterative Anwendung der Schritte 1-9 (Tabelle 1) gewonnen.

Tabelle 1: Die Arbeitsschritte des NIPALS-Algorithmus. Schritt

Funktion

Erklärung

1

s := xi

Wähle eine beliebige Spalte der Matrix X als den initialen Scoringvektor s.

2 3 4 5 6

l=

XTs sT s

Projiziere X auf s, um den entsprechenden

l=

l l

Normalisiere l auf die Länge 1.

s old := s

Loadingvektor zu erhalten.

Zwischenspeichern des aktuellen Scores.

Xl lT l

Projiziere X auf l, um den entsprechenden

d := s old − s

Bestimme die Differenz d des neuen und

s=

neuen Scoringvektor zu erhalten. alten Scores.

7

d < t?

Falls d größer als ein vorher definierter Schwellwert t (z.B. 10-6) ist, gehe zurück zu Schritt 2. Ansonsten fahre fort bei Schritt 8.

8

E := X − sl T

Entferne

den

Anteil

der

gefundenen

Hauptkomponente aus der Datenmatrix X. 9

X := E

Um weitere Hauptkomponenten zu finden, wiederhole den Algorithmus ab Schritt 1 mit der neuen Datenmatrix.

© modlab® 2007

NIPALS-Algorithmus

31

Der NIPALS-Algorithmus approximiert eine Datenmatrix X durch die kleinere Scoringmatrix T S und die Loadingmatrix LT mit Hilfe der Matrixmultiplikation X = SL (Abb. 1). Das

hochgestellte „T“ an einer Matrix verdeutlicht, dass die jeweilige Matrix transponiert wird.

Abbildung 2: Prinzip des NIPALS-Algorithmus. Ziel ist die Reproduktion der Datenmenge X durch eine Matrixmultiplikation von S und LT, die der Scoringmatrix und der transponierten

Loadingmatrix entsprechen. Die Rechtecke skizzieren die jeweiligen Matrizen mit den an den Ecken angegebenen Dimensionen.

Literatur Wold, H., Estimation of principal components and related models by iterative least squares, Multivariate Analysis (Ed., Krishnaiah, P. R.), Academic Press, NY, pp. 391-420 (1966). Wold S., 1974, A theoretical foundation of extrathermodynamic relationships (linear free energy rela-tionships), Chemica Scripta, 5, 97-106. Trygg, J., 2003, Chemometrics, Umeå University, Sweden, URL: http://www.chemometrics.se/editorial/sep2003.html#4

© modlab® 2007

NIPALS-Algorithmus

© modlab® 2007

32

Graphentheorie

33

Grundlagen Graphentheorie Graphen Der Graph als mathematisches Modell und Datenstruktur hat in der Cheminformatik eine zentrale Bedeutung zur Repräsentation von chemischen. Seien FE die Menge der Kanteneigenschaften und FV die Menge der Knoteneigenschaften. Ein Graph ist ein 4-Tupel G=(V, E, µV, µE), wobei •

V eine endliche Menge der Knoten (engl. vertices),

•

E ⊆ V×V die Menge der Kanten (engl. edges),

•

µV: V→ FV eine Funktion ist, die Knoten Eigenschaften zuweist,

•

µE: E→ FE eine Funktion ist, die Kanten Eigenschaften zuweist.

Man unterscheidet zwischen gerichteten und ungerichteten Graphen. Für ungerichtete Graphen gilt die Bedingung: •

∀ v,w ∈V: (v,w) ∈ E ⇔ (w,v) ∈ E

Die Kanten eines gerichteten Graphen haben also einen Anfangs- und einen Endpunkt, daher werden sie häufig als Pfeile dargestellt (Abbildung ). Im ungerichteten Graphen ist es gleich, welcher Knoten der Anfangs- und welcher der Endpunkt einer Kante ist.

© modlab® 2007

Graphentheorie

34

Abbildung 1: Beispiel a) für einen ungerichteten und b) für einen gerichteten Graphen.

Stellt ein Graph das Modell eines Moleküls dar, spricht man vom molekularen Graphen. Dabei repräsentieren die Knoten die Atome und die Kanten die Bindungen zwischen den Atomen. Die Knoten- und Kanteneigenschaften sind in diesem Fall Atom- und Bindungseigenschaften, z.B. das Elementsymbol oder die Bindungsordnung. Die graphische Repräsentation eines Moleküls ist in der organischen Chemie ein Standard (Abbildung ).

O

OH

Abbildung 2: Molekularer Graph von Tetrahydrocannabinol.

Da es meistens nicht von Bedeutung ist, welches Atom der Ausgangspunkt und welches der Endpunkt einer Bindung ist, handelt es sich beim molekularen Graphen um einen ungerichteten Graphen. Es gibt mehrere Möglichkeiten, einen Graphen im Rechner zu repräsentieren:

•

die Adjazenzliste

•

die Adjazenzmatrix,

•

und die Inzidenzmatrix.

Die Datenstrukturen mit den darauf operierenden effizienten Algorithmen werden in Java von der Graphenbibliothek jgrapht (http://jgrapht.sourceforge.net/, Version 0.8.3) bereitgestellt.

© modlab® 2007

Graphentheorie

35

Subgraphen, Pfade, Zyklen, Cliquen Sei G=(V, E, µv, µE) ein Graph. Dann heißt G’ = (V’, E’, µv’, µE’) ein Subgraph von G, wenn •

V’⊆V

•

E’⊆E

•

µV’(v) = µV(v) wenn v ∈ V’

•

µE’(e) = µE(E) wenn e ∈ E’.

Ein Pfad in einem Graphen ist definiert als eine Reihe von Kanten e1, e2,..., en für die gilt: •

ei und ei+1 haben einen gemeinsamen Knoten.

•

wenn für alle i ∉ {1, n} gilt: ei hat einen gemeinsamen Knoten mit ei-1 und den anderen mit ei+1.

Ein Zykel (Kreis) ist ein Pfad, für den gilt •

e1 und en haben einen gemeinsamen Knoten.

ein ungerichteter Graph

ein Pfad (fett gekennzeichnet)

ein Kreis (fett gekennzeichnet)

eine Clique (fett gekennzeichnet)

Abbildung 3: Übersicht über die Definitionen zur Graphentheorie. Ein Graph heißt vollständig, wenn jeder Knoten mit jedem anderen Knoten verbunden ist. Eine Clique ist definiert als ein vollständiger Subgraph. Die Größe der Clique ist definiert als die Anzahl der enthaltenen Knoten. Eine maximale Clique ist echt nicht in einer größeren © modlab® 2007

Graphentheorie

36

Clique enthalten. Eine größtmögliche Clique ist eine Clique der Größe s, wobei gilt: der Graph enthält keine weiteren Cliquen größer s. Alle Begriffe sind mit Beispielen in der Abbildung zusammengefasst.

Graphentheoretische Probleme und ihre Lösungen in der Cheminformatik Subgraph-Isomorphismus Die Suche von Substrukturen in Molekülen kann auf graphentheoretische Methoden zurückgeführt werden. In diesem Fall ist sie komplett innerhalb eines anderen gut untersuchten Problems enthalten: dem Problem des Subgraph-Isomorphismus. Eine bijektive Funktion f: V→ V’ ist ein Graph-Isomorphismus vom Graphen G=(V, E, µV,

µE ) zum Graphen G’ = (V’, E’, µv’, µE’), wenn gilt: •

µV(v) = µV’(f(v)) für alle v∈V.

•

für jede Kante e = (v,w) ∈ E existiert eine Kante e’ = (f(v),f(w)) ∈ E’, für die gilt µE(e) =µE(e’), und für jede Kante e’ = (v’,w’) ∈ E’ existiert eine Kante e = (f-1(v’),f-1(w’)) ∈ E für die gilt µE(e’) =µE(e).

Eine injektive Funktion f: V→ V’ ist ein Subgraph-Isomorphismus von G zu G’ wenn gilt: •

es existiert ein Subgraph S ⊆ G’ so dass f ein Graph-Isomorphismus von G zu S ist.

Informell kann man das Problem des Subgraph-Isomorphismus so formulieren: •

„Gegeben seien die Graphen G1 und G2. Enthält G1 eine Kopie von G2 als Subgraph? Können wir also eine Teilmenge der Knoten von G1 finden, die zusammen mit ihren Kanten in G1 eine exakte Kopie von G2 bilden, wenn wir die Zuordnung zwischen den Knoten und Kanten von G2 und den Knoten und Kanten des Teilgraphen von G1 entsprechend wählen?“

© modlab® 2007

Graphentheorie

37

Das Subgraph-Isomorphismus-Problem ist NP-vollständig. Folglich konnte kein Algorithmus geschrieben werden, der den Subgraph-Isomorphismus in polynomialer Zeit findet. Dennoch wurden Algorithmen gefunden, die eine akzeptable mittlere Laufzeit besitzen. In der Chemieinformatik hat der Ullman-Algorithmus breiten Zuspruch gefunden. Der Ullman-Algorithmus, der in 1976 von J.D. Ullman publiziert wurde, ist einer der schnellsten Algorithmen für das Problem des Subgraph-Isomorphismus. Insbesondere eignet er sich zum Auffinden von Subgraph-Isomorphismen auf Graphen, die in ihren topologischen Eigenschaften den Molekulargraphen ähneln. Da der Original-Algorithmus für nicht markierte Graphen geschrieben wurde, wird im Folgenden die Version für markierte Graphen aus vorgestellt. Der Ullman-Algorithmus vereinigt die Backtracking-Suche mit der „forward-checking“Technik. Der Pseudocode gibt Aufschluss über die genaue Funktionsweise des UllmanAlgorithmus.

function ullman(G=(V,E,µV,µE), G’=(V’,E’,µ’V,µ’E)) //finds all subgraph isomorphisms of an arbitrary graph G in a graph G’// P=(pij): a |V|×|V’| matrix F: set of vertex pairs, partial matching 01:for i = 1 to |V| 02:

for j = 1 to |V’|

03:

if µV(vi)= µ’V(v’j)

04:

then pij = 1

05:

else

06: 07:

pij = 0 end;

08:end; 09:backtrack(P,1,∅) end function

© modlab® 2007

Graphentheorie

38

function backtrack(P,i,F) 10:

if i > |V|

11: 12:

then F is an isomorphism else

13:

for all pij = 1 and j = 1 to |V’|

14:

F ← F ∪ {vi,v’j};

15:

P’ ← P;

16:

for all k > i

17:

p’ij = 0;

18:

end;

19:

if forwardchecking(P’,i,F) = true

20:

then backtrack(P’,i+1,F); F ← F\{vi,v’j};

21: 22:

end;

end function

function forwardchecking(P,i,F) 23:

for k = i+1 to |V|

24:

for l = 1 to |V’|

25:

if pkl = 1

26:

if consistent({vk,v’l},F)= false

27:

pkl = 0

28:

end;

29:

end;

30:

if there exists a row k in P such that pkl = 0 for l = 1.. |V’|

31: 32: 33:

then return false else return true

end function

© modlab® 2007

Graphentheorie

39

function consistent({vk,v’l},F) 34:

if for each {v,v’} ∈ F

35:

if there is an edge (vk,v) ∈ E then there must be an edge (v’,v’l) ∈ E’ with µE(vk,v)= µ’E(v’,v’l)

36:

if there is an edge (v’,v’l) ∈ E’ then there must be an edge (vk,v) ∈ E with µE(vk,v)= µ’E(v’,v’l)

37:

then return true

38: 39:

else return false

end function

Seien G=(V, E, µV, µE) und G’ = (V’, E’, µv’, µE’) zwei Graphen. Zunächst wird eine |V|×|V’| Matrix P = (pij) erschaffen (Zeilen 1-8), für die gilt: •

wenn µV(vi)= µV’(v’i), dann ist pij= 1,

•

sonst ist pij= 0.

In P steht also immer dann eine 1 in der i-ten Spalte und j-ten Zeile, wenn der Knoten vi den knoten v’i matcht (Abbildung ).

Abbildung 4: a) Querygraph, b) Molekülgraph, c) daraus resultierende Matrix P.

© modlab® 2007

Graphentheorie

40

Danach wird die rekursive Prozedur backtrack aufgerufen (Zeile 09). Als Initialparameter werden die Matrix P, Spaltenzähler i=1 und die Menge der erfolgreichen Matchings F=∅ übergeben. Die Prozedur backtrack funktioniert wie folgt: •

ist der Spaltenzähler i > |V|, so ist F ein Isomorphismus (Zeilen 10 und 11). Dies ist dann der Fall, wenn durch die rekursiven Aufrufe jedem Knoten aus V ein Knoten aus V’ zugewiesen wurde.

•

ist i ≤ |V|, dann ist der komplette Subgraph-Isomorphismus noch nicht gefunden. In diesem Fall wird für jede 1 in der i-ten Spalte rekursiv backtrack aufgerufen. Als Parameter werden der inkrementierte Spaltenzähler (i+1; Zeile 20), der bereits gefundene Teil des Subgraphisomorphismus in Form der Menge F und die durch die Prozedur forwardchecking bearbeitete Matrix P’.

Die Prozedur forwardchecking minimiert die Anzahl der Rekursionsschritte. Dazu wird jede 1 in der Matrix P, also jedes „matchende“ Knotenpaar (v,v’) auf Konformität mit jedem Element aus der Menge der bereits gefundenen „matchenden“ Knotenpaare F überprüft. Zwei Knotenpaare (v,v’) und (w,w’) sind konform, wenn gilt: •

wenn es eine Kante (v,w) gibt, muss es auch eine Kante (v’,w’) geben für die gilt: µ E(v,w) = µ’E(v’,w’), und

•

wenn es eine Kante (v’,w’) gibt, muss es auch eine Kante (v,w) geben für die gilt: µ E(v,w) = µ’E(v’,w’).

Somit erfolgt in der Prozedur forwardchecking nicht nur die Minimierung der Backtrackschritte, sondern auch die Überprüfung, ob F auch tatsächlich ein SubgraphIsomorphismus von G zu G’ ist. Außerdem wird die Rekursion unterbrochen, wenn in P eine Zeile vorhanden ist, die nur 0er enthält, also kein Subgraph-Isomorphismus in diesem Ast des Rekursionsbaums gefunden werden kann.

Maximum Common Substructure Eine weitere graphentheoretische Anwendung in der Cheminformatik ist die Berechnung der Isomorphie größtmöglicher gemeinsamer Teilgraphen (Maximum Common Substructure, MCS). Die MCS zweier Molekülgraphen ist in Abbildung 5 dargestellt. © modlab® 2007

Graphentheorie

41

N HO

OH OH

N

O N

Abbildung 5: MCS zweier Molekülgraphen

Die MCS kann zur Berechnung der Ähnlichkeit zweier Moleküle benutzt werden und ist wie folgt definiert: Seien G1 und G2 zwei Graphen, dann ist MCS eine Struktur (S1, S2), so dass gilt

•

S1 ist Subgraph von G1

•

S2 ist Subgraph von G2

•

S1 und S2 sind isomorph

•

es gibt kein S1’ für das gilt: S1 ist Subgraph von S1’

•

es gibt kein S2’ für das gilt: S2 ist Subgraph von S2’

Um alle MCS zweier Graphen aufzuzählen, bedient man sich u.a. der Cliquendetektion auf dem Assoziationsgraphen.

Cliquendetektion auf dem Assoziationsgraphen Die Lösung des MCS-Problems durch Cliquendetektion auf dem Assoziationsgraphen besteht im wesentlichen aus zwei Schritten: •

Erstellung des Assoziationsgraphen,

•

Cliquendetektion auf dem Assoziationsgraphen.

Für zwei Graphen G=(V, E, µV, µE) und G’ = (V’, E’, µv’, µE’) ist der Assoziationsgraph (Abbildung 6) definiert als ein unmarkierter Graph GA = (VA, EA) mit VA ⊆ V × V’ und EA ⊆ VA × VA. vA ∈ VA heißt a-Knoten (Assoziationsknoten), eA ∈ EA heißt a-Kante (Assoziationskante) von GA. Für GA muß gelten: © modlab® 2007

Graphentheorie •

VA = {(v,v’)| v ∈ V, v’ ∈ V’, µV(v) = µV’(v’) }

•

EA besteht aus a-Kanten eA= (vA,wA) mit vA=(v,v’), wA=(w,w’), so dass gilt

42

o v ≠ w und v’ ≠ w’ o wenn eine Kante e = (v,w) ∈ E existiert, muss eine Kante e’ = (v’,w’) ∈ E’ existieren und µE (e) = µE’(e’) muss gelten.

o wenn keine Kante e = (v,w) ∈ E existiert, darf auch keine Kante e’ = (v’,w’) ∈ E’ existieren.

Abbildung 6: a) Querygraph. b) Molekülgraph. c) Assoziationsgraph aus a) und b), grau eingezeichnet: Clique der Größe 4. Wenn also die a-Knoten (v1,v’1), (v2,v’2)… (vn,v’n) für 1≤ n ≤ |V|⋅|V’| im Assoziationsgraphen paarweise miteinander verbunden sind, ist der Subgraph in V, der durch die Knoten v1, v2,…,vn induziert ist, isomorph zum Subgraphen in V’, der durch die Knoten v’1, v’2,…,v’n induziert ist. Der Isomorphismus ist durch (v1,v’1), (v2,v’2)… (vn,v’n) explizit gegeben. Abbildung 6 zeigt die Erstellung eines Assoziationsgraphen aus dem Molekülgraphen g1 (Alanin) und einem Querygraphen g2 (Carbonylfunktion verbunden mit einem Stickstoff über eine Kohlenstoffbrücke). Die Atomtypen entsprechen dabei den Knoteneigenschaften, die Bindungsordnungen den Kanteneigenschaften. Nun wird z.B. ein a-Knoten (3,1’) erstellt, weil der Knoten 3 dieselbe Eigenschaft trägt wie der Knoten 1’, nämlich O (Sauerstoff). Einen aKnoten (3,4’) gibt es nicht, da der Knoten 3 die Eigenschaft O (Sauerstoff) und der Knoten 4’ © modlab® 2007

Graphentheorie

43

die Eigenschaft N (Stickstoff) trägt. Ähnlich verfährt man bei der Erstellung der a-Kanten. Die a-Knoten (3,1’) und (2,2’) werden durch eine a-Kante verbunden, da sowohl in g1 die Knoten 3 und 2 als auch in g2 die Knoten 1’ und 2’ durch eine Kante mit der Eigenschaft „Doppelbindung“ verbunden sind. Die a-Knoten (3,1’) und (6,4’) sind verbunden, weil es in g1 keine Kante (3,6) und in g2 keine Kante (1’,4’) existiert. Der komplette a-Graph enthält eine Clique der Größe 4, nämlich {(3,1’), (4,3’),(2,2’),(6,4’)} (in Abbildung , rot hervorgehoben). Offensichtlich entspricht sie dem MCS von g2 zu g1. Der zweite Schritt besteht in der Suche der Cliquen auf dem a-Graphen. In der Literatur gibt es vielfältige Ansätze zur Lösung dieses Problems. Der schnellste vollkombinatorische Ansatz ist der Bron-Kerbosch-Algorithmus. Es handelt sich dabei um einen rekursiven BacktrackingAlgorithmus, der alle maximalen Cliquen in einem Graphen aufzählt. Im Wesentlichen operiert der BK-Algorithmus auf drei Mengen C, P und S. Die Menge C (Clique) ist die Menge der Knoten, die bereits als komplett verbunden gefunden wurden, also eine Clique bilden. Die Menge P (Potential) enthält diejenigen Knoten, die zu dem zuletzt zu C hinzugekommenen Knoten adjazent sind, folglich also für die Erweiterung der Clique benutzt werden können. Die Menge S enthält all die Knoten, die nicht mehr für die Erweiterung von C benutzt werden können, weil alle Cliquen, die diese Knoten enthalten, bereits gefunden wurden. Der Pseudocode gibt die Funktionsweise des BK-Algorithmus wieder. Der Algorithmus wird mit C,S = ∅ und P=V initialisiert. Wenn die Mengen S und P leer sind, ist C eine maximale Clique und kann ausgegeben oder gespeichert werden (Zeile 03). Ist dies nicht der Fall, so wird jeder Knoten aus P in der Schleife (Zeilen 04-11) bearbeitet. Ein ausgesuchter Knoten ui wird aus P entfernt (Zeile 05), P wird nach P’ (Zeile 06) und S nach S’ (Zeile 07) kopiert, um für den rekursiven Aufruf verwendet zu werden. Knoten, die mit ui verbundenen sind, werden in N gespeichert (Zeile 08). Der rekursive Aufruf erfolgt mit der Menge C∪{ui} (die bereits gefundene Clique wird um ui erweitert), P∩N (nur die Knoten in P, die mit ui verbunden sind, dürfen für die potentielle Erweiterung der Clique verwendet werden) und S∩N (alle Knoten, die bereits betrachtet wurden und mit ui verbunden sind) (Zeile 09). Anschließend wird ui der Menge S zugefügt (Zeile 10).

© modlab® 2007

Graphentheorie

44

function find_all_cliques(C,P,S) //enumerates all cliques in an arbitrary graph G // C: set of vertices belonging to the current clique P: set of vertices which can be added to C S: set of vertices which are not allowed to be added to C N[u]: set of vertices which are adjacent to vertex u in G 01:

Let P be the set {u1,...,uk}⊆ V;

02:

if P = ∅ and S = ∅

03:

then C is a maximal clique

04:

else for i ← 1 to k

05:

do P ← P\{ui};

06:

P’ ← P;

07:

S’ ← S;

08:

N ← {v∈V | {ui,v}∈E};

09:

find_all_cliques(C∪{ui},P’∩N,S’∩N);

10:

S ← S ∪ {ui};

11:

end;

12:end; end function Der BK-Algorithmus wurde bereits in vielen bioinformatischen Zusammenhängen verwendet und modifiziert. Die besondere Geschwindigkeit der sog. „Version2“ des BK-Algorithmus liegt in einer Heuristik bei der Wahl des Knotens ui aus P. Diese motiviert sich aus der Verwendung der Menge S, da man die Abbruchbedingung für die Ausgabe der Clique auch so formulieren kann: Es gibt einen Knoten s in S, der mit allen Knoten in P verbunden ist. Ein solcher Knoten s kann niemals aus S entfernt werden. Somit haben wir eine Möglichkeit, vorzeitig vorherzusagen, ob ein Ast im Backtrackingbaum zum Erfolg führt oder nicht. Allgemein nennt man diese Backtracking-Technik die „branch-and-bound“ Methode, das Abbruchkriterium ist das „branch-and-bound“ Kriterium. Die „Version2“ des BKAlgorithmus sucht ui aus P so aus, dass das „branch-and-bound“ Kriterium möglichst schnell eintritt. Dies geschieht in folgender Weise: alle Knoten in S sind mit einem Zähler assoziiert, welcher die Anzahl der Kanten zu den Knoten in P wiedergibt. Verschiebt man nun einen Knoten p aus P nach S, so vergrößern sich die Zähler aller Knoten in S höchstens um 1 und p bekommt einen neuen Zähler. Das „branch-and-bound“ Kriterium kann nun folgendermaßen © modlab® 2007

Graphentheorie

45

formuliert werden: Es gibt einen Knoten s in S mit Zähler z =|P|. Bron und Kerbosch zeigen, dass man dieses Kriterium am schnellsten erreicht, wenn man ui aus P so auswählt, dass ui mit demjenigen Knoten s in S durch eine Kante verbunden ist, welcher den höchsten Zähler z hat. Dadurch erhöht man in jedem Schritt den Zähler von s um 1, das „branch-and-bound“ Kriterium wird auf dem schnellst möglichen Wege erreicht.

Literatur

Graphentheorie, Komplexizität, grundlegende Algorithmen G. Valiente (2002), Algorithms on Trees and Graphs, Springer Berlin Heidelberg. Ottmann T, Widmayer P (1996). Algorithmen und Datenstrukturen, Spektrum Akademischer Verlag GmbH, Heidelberg, 537ff. Garey R, Johnson DS (1979). Comuters and Intractability: a Guide to the Theory of NPCompleteness, W. H. Freeman and Company, San Francisco. Schöning U (2001), Algorithmik, Spektrum Akademischer Verlag GmbH, Heidelberg.

Algorithmen zur Substruktursuche Barnard J (1993). "Substructure Searching Methods - Old and New." J. Chem. Inf. Comp. Sci. 33, 532-538. Messmer B (1995). Efficient Graph Matching Algorithms, Universität Bern, Bern, 11ff. Ullmann J (1976). "Algorithm for Subgraph Isomorphism." Journal of the ACM 23, 31-42.

Cliquendetectionsalgorithmen und MCS Bron C, Kerbosch J (1973). "Agorithm 457: Finding all Cliques of an Undirected Graph." Communications of the ACM 16, 575-577. Koch I (2001). "Enumerating all connected maximal common subgraphs in two graphs." Theoretical Computer Science 250, 1-30. Raymond J, Gardiner E, Willett P (2002). "RASCAL: Calculation of graph similarity using maximum common edge subgraphs." Computer Journal 45, 631-644.

© modlab® 2007

Graphentheorie

© modlab® 2007

46

Künstliche neuronale Netze

47

Künstliche neuronale Netze Künstliche neuronale Netze (engl.: artificial neural net, ANN) bilden eine wichtige Kategorie im Bereich des maschinellen Lernens und der künstlichen Intelligenz. Ihre Aufgabe ist es, anhand von Trainingsbeispielen eine Funktion (allgemein: ein Muster) zu erlernen und anschließend das Erlernte auf unbekannte Beispiele anzuwenden. Ihre generelle Fähigkeit zur Mustererkennung und Klassifizierung macht ANNs in großer Breite in der Praxis einsetzbar: Von der Spracherkennung über Komprimierung und Rauschunterdrückung in der Bildverarbeitung bis hin zur Vorhersage von Aktienkursverläufen wird auf künstliche neuronale Netze zurückgegriffen. ANNs sind in ihrem Aufbau und in ihrer Funktion biologischen neuronalen Netzen nachempfunden. Das Zusammenspiel relativ simpler Untereinheiten, die in Anlehnung an die biologische Terminologie als Neuronen bezeichnet werden, ergibt durch deren Vernetzung ein komplexes System zur Informationsverarbeitung. Die Art der Vernetzung der Neuronen ist ein Merkmal zur Unterteilung in verschiedene Klassen künstlicher neuronaler Netze. So unterscheidet man im Wesentlichen zwischen vollständig und teilweise verbundenen Netzen. Als Vertreter der vollständig verbundenen Netze sei hier das Hopfield-Neural-Network genannt. In ihm ist jedes Neuron mit allen anderen Neuronen verbunden. Dieser Teil des Skriptes konzentriert sich hauptsächlich auf eine bestimmte Klasse der teilweise verbundenen Netzwerke – den Feed-Forward-Netzen.

Genereller Aufbau von Feed-Forward-Netzen Feed-Forward-Netze folgen einem schichtartigen Aufbau. Die oberste Schicht wird als InputLayer bezeichnet. Darauf folgt eine (theoretisch) beliebige Anzahl von sog. Hidden-Layer. Die

abschließende Schicht wird als Output-Layer

© modlab® 2007

bezeichnet. Abbildung 1 zeigt den

Künstliche neuronale Netze

48

allgemeinen Aufbau eines vierlagigen Feed-Forward-Netzes. Jede Schicht wird also bei der Nomenklatur gezählt.

Abbildung 1: Aufbau eines vierlagigen Feed-Forward-Netzes.

Die

namensgebende

Eigenschaft

der

Feed-Forward-Netze

ist,

dass

die

Informationsverarbeitung nur in eine Richtung (vom Input-Layer in Richtung des OutputLayer) und zwischen direkt benachbarten Schichten erfolgt. Daraus resultiert die in Abbildung 1 gezeigte Vernetzung zwischen den benachbarten Schichten: Jedes Neuron eines Hidden-Layers ist mit allen Neuronen der darüber liegenden Schicht und mit allen Neuronen der darunter liegenden Schicht verbunden. Im Folgenden soll zunächst auf die einfachste Feed-Forward-Architektur - das Perzeptron eingegangen werden, um anschließend komplexere Beispiele zu betrachten.

Das Perzeptron Das Single Layer Perzeptron Das Single Layer Perzeptron, kurz SLP, stellt die einfachste Art eines Feed-Forward-Netzes dar. Entwickelt wurde es 1958 von Frank Rosenblatt. Es besteht aus einem Input-Layer, dessen Neuronen direkt mit einem Output-Neuron verbunden sind (Abbildung 2). Zur Vermeidung von Missverständnissen: Nach obiger Definition handelt es sich beim SLP um ein zweilagiges Netz, welches jedoch nur eine aktive Schicht besitzt und deshalb als Single Layer bezeichnet wird. © modlab® 2007

Künstliche neuronale Netze

49

Die Verbindungen zwischen den Neuronen sind mit Gewichten versehen. Diese beschreiben die Intensität, mit der der jeweilige Input auf das Ergebnis des Outputs einwirkt. In dem Ouput-Neuron werden die Eingaben verarbeitet. Es wird eine gewichtete Summe, nämlich die Summe der Produkte der Gewichte und der Inputwerte, berechnet. Weiterhin wird ein zusätzlicher Wert, das sog. Bias, zur Summe hinzuaddiert. Das Bias-Neuron ist ein zusätzliches Neuron, welches stets den Wert 1 ausgibt. Wie alle Neuronen besitzt es ein Gewicht. In Abbildung 2 ist dieses mit Ө bezeichnet.

 n  =  ∑ xi wi  + θ  i =1 

Abbildung 2: Schematische Darstellung eines Perzeptrons

Das

Ergebnis

dieser

Berechnung

wird

in

der

sog.

Aktivierungs-

oder

Schwellenwertfunktion verarbeitet. Der Wert, den diese Funktion annimmt, bezeichnet den Zustand des Neurons. Für eine binäre Schwellenwertfunktion gibt er an, ob ein Neuron aktiviert oder deaktiviert ist (Abbildung 3), bzw. bei anderen Aktivierungsfunktionen, wie stark die Aktivierung ist (Abbildung 6). Im Falle des SLPs, das mit einer binären Schwellenwertfunktion arbeitet, gibt das Neuron als Output den Wert 1 aus, wenn das Neuron aktiviert ist. Sollte das Ergebnis der gewichteten Summe unter dem gegebenen Schwellenwert liegen, so geht das Neuron in einen deaktivierten Zustand über. Die Ausgabe beträgt in diesem Fall -1. Der Schwellenwert liegt bei der in Abbildung 3 gezeigten binären Schwellenwertfunktion bei 0.

© modlab® 2007

Künstliche neuronale Netze

50

Abbildung 3: Binäre Schwellenwertfunktion des SLP. Der Schwellenwert der Funktion ist hier 0. Mit Hilfe des Bias ist eine implizite Verschiebung des Schwellenwertes möglich. Wie wir später sehen werden, ist das Training eines neuronalen Netzes durch die Veränderung der Gewichte charakterisiert. Die Veränderung des Bias zieht eine Veränderung des Schwellenwertes des Neurons nach sich. Sinkt das Bias-Gewicht, wirkt sich dies durch ein Sinken der gewichteten Summe aus. Folglich erhöht sich dadurch implizit der Schwellenwert für das Neuron, d.h. Eingaben und Gewichte müssen größer sein, um mit gleicher Stärke auf die Aktivität des Neurons einzuwirken. Der Vorteil, die Veränderung des Schwellenwertes über das Bias darzustellen, liegt darin, dass in Netzen mit mehr als einem Neuron für alle Neuronen die gleiche Aktivierungsfunktion benutzt werden kann. Explizit bleibt der Schwellenwert der Aktivierungsfunktion unverändert für alle Neuronen gleich. Durch das Verändern des Bias-Gewichts während des Trainings kann er jedoch implizit für jedes Neuron individuell verschoben werden. Ohne das zusätzliche Bias-Neuron benötigte man dafür für jedes Neuron eine individuelle Aktivierungsfunktion, die zusätzlich während des Trainings verändert werden müsste. SLPs sind in der Lage, Klassifizierungsprobleme zu lösen, die linear separierbar sind. Die Klassifizierung erfolgt hierbei implizit durch das Legen einer trennenden Hyperebene durch die Punktmenge. Diese Hyperebene ist in ihrer Dimension stets um eins geringer als die Dimension der zu erlernenden Funktion. Explizit erfolgt die Trennung durch die vom Perzeptron erlernte Funktion, die orthogonal zu dieser trennenden Hyperebene steht.

© modlab® 2007

Künstliche neuronale Netze

51

Abbildung 4 illustriert die lineare Trennbarkeit. Hierbei ist die eingezeichnete Gerade die trennende Hyperebene. Wie man sieht ist die zu erlernende Funktion zweidimensional, während die trennende Hyperebene eine Gerade ist, die sich durch die allgemeine eindimensionale Geradengleichung darstellen lässt. Sie ist in (1) gezeigt. f ( x) = ax + b

(1)

Ist das Klassifizierungsproblem nicht linear separierbar, kann es durch das SLP nicht gelöst werden. Als populäres Beispiel kann hier die XOR-Funktion angeführt werden. Es kann durch ein einzelnes Perzeptron nicht gelöst werden.

Abbildung 4: Ein linear separierbares Klassifizierungsproblem im zweidimensionalen Deskriptorraum.

Die Lösung des XOR-Problems Für logische Operatoren wie AND oder OR ist das SLP in der Lage, durch das Legen einer Geraden, die Punkte so zu separieren, dass eine Lösung möglich ist. Betrachtet man aber Abbildung 5, so sieht man, dass dies für den XOR-Fall nicht möglich ist. Egal, wie eine Gerade durch den Raum gelegt wird, ist es nicht möglich, dass diese die Punkte {(0,0),(1,1)} von den Punkten {(0,1),(1,0)} trennt. Gezeigt wurde dies 1969 von Minsky und Papert. Die Separierung des Raumes durch zwei Geraden würde zu der Lösung des Problems führen. Dies ist allerdings nur mit einem mehrschichtigen Perzeptron lösbar.

© modlab® 2007

Künstliche neuronale Netze

52

Abbildung 5: Das XOR-Problem. Die umrahmten Eingaben sollen in eine Klasse zusammengefasst werden, die anderen beiden in eine zweite Klasse.

Mehrlagige Netze Ein einzelnes Neuron erstellt aus seinen Inputwerten und den zugehörigen Gewichten eine gewichtete Summe, also eine lineare Funktion. Das SLP, welches im Wesentlichen ein Neuron darstellt, kann also nur linear trennbare Klassifizierungsprobleme lösen. Wie am Beispiel der XOR-Funktion gesehen, sind für das Perzeptron damit viele (sogar der überwiegende Teil) der Klassifizierungsprobleme nicht lösbar. Schon das XOR-Problem zeigt, dass der Übergang vom SLP zum Multi Layer Perzeptron (abgekürzt MLP) die Menge der lösbaren Probleme vergrößert: MLPs sind generelle Funktionsapproximatoren. Der Übergang zu mehrlagigen Netzen mit mehreren Neuronen pro Hidden-Schicht ergibt eine Kombination linearer Gleichungen. Damit können beliebige Funktionen approximiert werden. MLPs enthalten nichtlineare Aktivierungsfunktionen, z.B. die Sigmoidalfunktion und der hyperbolische Tangens. Neben ihrem nichtlinearen Verlauf zeichnet diese Funktionen aus, dass sie stetig und differenzierbar sind. Dies ist, wie wir noch sehen werden, für das NetzTraining wichtig. Weiterhin bilden sie alle möglichen Ergebnisse der gewichteten Summe auf das Intervall ]0,1[ bzw. ]-1,1[ ab. Abbildung 6 zeigt die entsprechenden Funktionsverläufe. Achtung: Die Fähigkeit eines MLPs beliebige, auch nichtlineare Funktionen zu approximieren, ist nicht durch die Wahl der Aktivierungsfunktion begründet, sondern durch die Architektur des Netzes.

© modlab® 2007

Künstliche neuronale Netze

53

Abbildung 6: Links die Sigmoidalfunktion, rechts der hyperbolische Tangens

Die Neuronen der Hidden-Layer geben also nicht direkt das Ergebnis der gewichteten Summe an die nächste Schicht weiter, sondern ein zuvor durch die Aktivierungsfunktion transformierten Wert. Eine zentrale Frage beim Einsatz von ANNs zu Klassifizierung ist die der zu verwendenden Netztopologie. Wie viele verdeckte Schichten werden benötigt? Wie viele Neuronen sind in der jeweiligen Schicht optimal? Es gibt kein deterministisches mathematisches Verfahren, das die optimale Netztopologie für ein gegebenes Problem bestimmt. Bei der Wahl der Netztopologie stehen sich wie so oft zwei gegensätzliche Anforderungen gegenüber: Einerseits soll das Netz die Möglichkeit haben, sich gut an die Trainingsbeispiele anzupassen. Hierfür ist ein großes Netz von Vorteil. Andererseits soll das trainierte Netz möglichst gut generalisieren, also ihm unbekannte Eingaben korrekt klassifizieren, und nicht einfach nur den Trainingsdatensatz auswendig lernen. Diesem als Overfitting bezeichneten Problem entgeht ein kompaktes Netz eher. Deshalb bleibt hier oft nur der trial-and-error-Ansatz. Dennoch lässt sich generell sagen, dass in aller Regel mehr als zwei Hidden-Layer keinen Sinn machen. Bereits Netze mit dieser Anzahl an verdeckten Schichten können potenziell alle in der Praxis auftretenden Trennungsprobleme lösen. Eine weitere Vergrößerung der Netze würde also keinen Vorteil für deren Mächtigkeit bringen, gleichzeitig aber die Gefahr des Overfittings erhöhen und das Training verlangsamen.

Netztraining Das Training eines ANN ist nichts anderes als die Anpassung der Gewichte zwischen den Neuronen. Ziel ist es, die Gewichte so anzupassen, dass der Fehler des Netzes minimal wird. Mit Fehler ist hier die Summe der Quadrate der Differenzen zwischen den vom Netz berechneten Werten für die Trainingsmuster und deren gegebenen Zielwerten gemeint. Die © modlab® 2007

Künstliche neuronale Netze

54

Existenz der Zielwerte ist es, die das Training von ANNs zu einem Vertreter des überwachten

Lernens macht. Zusammengefasst führt dies zu Gleichung (2), wobei E den Fehlerwert des Netzes summiert über alle k Trainingsmuster, ti den Zielwert des Trainingsmusters i und oi dessen vom Netz berechneten Output-Wert darstellt.

k v E ( w) = ∑ (oi − t i ) 2

min.

(2)

i =1

Dies ist die während des Trainings zu minimierende Funktion. Macht man sich klar, dass in die Berechnung der Output-Werte oi die Netzgewichte maßgeblich einfließen, ist leicht v ersichtlich, dass E eine Funktion der Gewichtswerte des Netzes ist (in (2) ist w der Vektor aller Gewichtswerte des Netzes). Die Anzahl der Gewichte gibt also die Dimension des Raumes an, in dem die Fehlerfunktion definiert ist. Abbildung 7 zeigt eine mögliche Fehlerfunktion eines Netzes mit zwei Gewichten. Auf die Gefahren in sehr hochdimensionalen Räumen und die unter anderem daraus folgende Bevorzugung möglichst kompakter Netze soll später noch eingegangen werden. Der populärste Trainingsalgorithmus für ANNs - der Backpropagtion of Errors (kurz

Backprop) - basiert auf der Idee des Gradientenabstieges. Dabei wird stets in Richtung des steilsten Gefälles des Funktionsgraphen gelaufen. Bildlich kann man sich dies als einen zufällig im Gebirge abgesetzten Skifahrer (das entspricht der zufällige Gewichtsinitialisierung beim Erzeugen des Netzes) vorstellen, der anschließend eine Abfahrt immer entlang des stärksten Gefälles startet, bis er im nächstgelegenen Tal landet.

Abbildung 7: Fehlerlandschaft, Quelle: [1]

© modlab® 2007

Künstliche neuronale Netze

55

Der Bewegung des Skifahrers entspricht in unserem Problem eine Veränderung der Gewichte, sodass sich der Fehler des Netzes verringert. Das Gradientenabstiegverfahren ermöglicht es nun, für jedes Gewicht (jede Dimension) die Richtung (Vorzeichen) und die Intensität (numerischer Wert) der Veränderung zu finden, um den Fehler zu verkleinern. Mathematisch wird dies durch eine partielle Ableitung der Fehlerfunktion nach jedem Gewicht erreicht. Der dabei entstehende Vektor von partiellen Ableitungen wird als Gradient oder Nabla-Operator bezeichnet und ist in (3) definiert.

 ∂ ∂ ∂   g =  , , ... , ∂x n   ∂x1 ∂x 2

(3)

Backpropagation nimmt nun die Modifizierung der Gewichte proportional zum Gradienten vor, wobei man den Proportionalitätfaktor als Lernrate η bezeichnet. Sie dient als Faktor der linearen Modulierung des Einflusses des Gradienten. In unserem Bild des Skifahrers wäre dies so etwas wie die Qualität seiner Skier: Je besser das Material (je höher η), desto schneller fährt er bei gleichem Gefälle. Zusammengefasst ist dies in (4).

v ∆w = −η ⋅ g

(4)

v ∆w stellt die Veränderung der einzelnen Gewichte zusammengefasst als Vektor dar, während

g den Gradienten der Fehlerfunktion E bezeichnet. Am hypothetischen Beispiel eines Netzes mit nur einem Gewichtswert soll das Vorgehen des Gradientenabstieges verdeutlicht werden (Abbildung 8). Der Gradient ist die erste Ableitung der gezeigten Funktion nach dem Gewicht w. Hier wird auch deutlich, warum in (4) der negative Gradient zur Gewichtsänderung herangezogen wird: Der Funktionswert der ersten Ableitung ist bei w1 negativ, die Änderung des Gewichtes muss jedoch positiv sein, um sich dem Minimum der Funktion zu nähern. Genau umgekehrt liegt der Fall bei w2. Außerdem ist in Abbildung 8 ersichtlich, dass sich die Neigung der Geraden, die ja die partielle Ableitung der Funktion an einer Position darstellt, auf die Stärke der Gewichtsänderung auswirkt. Je stärker die Steigung (je höher der Wert der Ableitung), desto größer ist auch die Veränderung des Gewichtes hin zum Optimum.

© modlab® 2007

Künstliche neuronale Netze

56

Abbildung 8: Gradientenabstieg. w1 und w2 stellen zwei unterschiedliche Werte des Gewichtes w dar. Die Steigungen der Geraden entsprechen den Ableitungen der Funktion E(w) an diesen Positionen dar. Das zu findende Minimum der Funktion liegt bei w = 0. Backpropagation stellt ein iteratives Verfahren zur Verfügung, um einen Gradientenabstieg auf der Fehlerfunktion zu ermöglichen. Die Veränderung der Gewichtwerte wird dabei wie in (5) bestimmt.

wijnew = wijold + ∆wij .

(5)

wij stellt die Veränderung des Gewichtwertes wij in einem Trainingsschritt dar. Die Berechnung von wij ist in (6) gezeigt und geht auf (4) und damit auf den Gradientenabstieg zurück, wobei hier die Fehlerfunktion E partiell nach dem zu aktualisierenden Gewicht wij abgeleitet wird. (Zur Wiederholung: Dies entspricht einer Komponente des Gradienten g aus (4), der ja den Vektor aller partiellen Ableitungen darstellt).

∆wij = −η

∂E = ηδ j xi . ∂wij

(6)

wij stellt hier das Gewicht zwischen Neuron i und Neuron j dar; xi ist der Ausgabewert des Neurons i. Entscheidend für den Backpropagation-Algorithmus ist nun, dass sich diese partiellen Ableitungen der Fehlerfunktion nach den Gewichten und damit δj effizient berechnen lassen. Dabei muss allerdings unterschieden werden, ob es sich bei dem zu © modlab® 2007

Künstliche neuronale Netze

57

verändernden Gewicht um ein Gewicht zwischen dem untersten Hidden-Layer und der Output-Layer handelt, oder um ein Gewicht zwischen höher liegenden Schichten. Die Berechnung von δj unterscheidet sich also je nach dem, ob Neuron j ein Output-Neuron oder ein Hidden-Neuron ist. Zusammengefasst ist dies in (7).

act ' ( net j )(t j − o j )  δ j = act ' ( net ) δ w j ∑ k jk  1≤ k ≤ n

falls j Outputneur on ist falls j Hiddenneur on ist

(7)

Hierbei entspricht tj im oberen Teil von (7) dem Zielwert und n im unteren Teil der Anzahl der Neuronen in der Schicht unterhalb des Neurons j. Abbildung 9 soll der besseren Übersicht und der Klärung der Variablenbezeichnungen dienen.

Abbildung 9: Kantengewichte sind rechts von ihrer Kante geschrieben, Ein- bzw. Ausgabewerte der Neuronen links davon. Die Bias-Neuronen sind durch „B“ gekennzeichnet, die Kästchen stellen die Aktivierungsfunktion dar. Exemplarisch sollen die Berechnungen der w-Werte für die Gewichte wij und vjk aus Abbildung 9 demonstriert werden. Nehmen wir an, das angegebene Netz benutzt die in (8) gezeigte Sigmoidalfunktion zur Aktivierung. Da wir nach (7) deren erste Ableitung benötigen, © modlab® 2007

Künstliche neuronale Netze wird

auch

ersichtlich,

58

warum

eine

Forderung

an

Aktivierungsfunktionen

ihre

Differenzierbarkeit ist. Die erste Ableitung der Sigmoidalfunktion sig(x) ist, wie in (9) ersichtlich, leicht zu errechnen.

sig ( x) =

1 . 1 + e−x

(8)

sig ' ( x) = sig ( x) ⋅ (1 − sig ( x)) .

(9)

Damit ergibt sich nach (6) für ∆vjk (k ist Output-Neuron): ∆v jk = η ⋅ (t − o ) ⋅ sig (net ) ⋅ (1 − sig (net )) ⋅ o j , k 444444 k3 1k44k44442

δ

k

und für ∆ wij (j ist Hidden-Neuron): ∆wij = η ⋅ (δ v ) ⋅ sig (net ) ⋅ (1 − sig ( net )) ⋅ inpi . k jk j j 144444424444443

δ

j

Aus (7) ist auch ersichtlich, dass zur Berechnung der δ-Werte der Neuronen einer Schicht die

δ-Werte der Neuronen der darunter liegenden Schicht benötigt werden. Einzig die δ-Werte des Output-Layers sind direkt berechenbar. Der Netzfehler wird also durch den Algorithmus vom Output-Layer über die Hidden-Layer

von unten nach

oben durchgereicht

(„backpropagation“). Dies gab dem Trainingsalgorithmus seinen Namen. Die Gewichte der Bias-Neuronen werden wie alle anderen Gewichte nach dem beschriebenen Schema während des Trainings angepasst. Für Details der mathematischen Herleitung sei auf [1] verwiesen. Dort wird genau beschrieben, wie die partielle Ableitung von E nach einem Gewicht durch Anwendung der Kettenregel das in (7) angegebene δ ergibt.

Probleme und Lösungsansätze Die Verwendung des Gradientenabstiegs bringt nicht nur Vorteile mit sich. Man muss sich auch Gedanken um etwaige Probleme machen, die bei der Wahl dieser Methode auftreten © modlab® 2007

Künstliche neuronale Netze

59

können. Ohne Zweifel ist es mit Hilfe des Gradientenabstiegs möglich, das Netz immer besser werden zu lassen. Es stellt sich nur die Frage, ob man nach Beendigung des Trainings tatsächlich die optimale Gewichtskonstellation und dadurch auch die Funktion gefunden hat, die die Datenklassen am besten voneinander trennen kann. Hierbei trifft man oft auf ein allgemeines Problem von Optimierungsverfahren: das Auftreten lokaler Minima, die nicht wieder verlassen werden. Die Anzahl potentieller lokaler Minima erhöht sich mit der Dimension des Raumes. Es macht also Sinn, das Netz nicht unbedingt mit der größtmöglichen Anzahl Parameter, die einem zur Verfügung stehen, zu trainieren, sondern im Vorfeld wichtige Parameter zu selektieren. Im Gegensatz dazu steht das Problem der flachen Plateaus. Diese bewirken, dass der Gradient sehr niedrige Werte nahe 0 oder unter Umständen sogar 0 annimmt. Dies führt zu einer Stagnation des Verfahrens. Das Risiko des Auftretens flacher Plateaus steigt mit abnehmender Anzahl Trainingsparameter. Im Wesentlichen sind zwei Methoden zu benennen, die sich mit den genannten Problemen befassen, und gute Resultate erzielen: 1. Das Verändern der Lernrate zur Laufzeit, und 2. das Einführen eines Momentum-Terms. Wie wir bereits gesehen haben, besitzt die Lernrate Einfluss auf die Veränderung der Gewichte nach jedem Trainingszyklus. Am Anfang macht es durchaus Sinn, mit einer hohen Lernrate zu starten, da der Algorithmus dadurch in der Lage ist, sich recht schnell auf dem Suchraum zu bewegen. Im späteren Verlauf des Trainings ist eine große Lernrate allerdings nicht von Vorteil. Diese stellt nämlich auch die Gefahr dar, gute Optima wieder zu verlassen oder zwischen gleich guten Werten zu oszillieren. Eine Verringerung der Lernrate zur Laufzeit des Algorithmus bietet den Vorteil, dass man zunächst recht schnell gute Positionen des Raumes erkennt und diese durch die Verringerung der Lernrate nicht sofort wieder verlässt. Weiterhin ist dadurch eine feinere Optimierung möglich, indem das Netz mit immer kleiner werdenden Schritten weiter verbessert wird. Ein einfaches Beispiel für dieses Vorgehen ist die Verwendung einer adaptiven Lernrate nach Darken und Moody. Diese ist in (10) definiert.

ηn =

η0 n 1+ r

.

© modlab® 2007

(10)

Künstliche neuronale Netze

60

Man startet zunächst mit einer hohen Lernrate η0. Der Wert der Lernrate verringert sich im Laufe des Trainings auf ungefähr η0 / (1 + n), wobei n die Anzahl bisher präsentierter Trainingsbeispiele ist. Der Parameter r gibt hierbei die Stärke der Verringerung der Lernrate in jedem Schritt an. Nach r Schritten hat sich die Lernrate halbiert. In der Praxis liefert bereits die Verwendung dieses einfachen Verfahrens gute Ergebnisse. Genauere Informationen zu dem Vorgehen nach Darken und Moody sowie eine Auswahl weiterer möglicher Verfahren zur Veränderung der Lernrate während des Trainings sind in [4] zu finden. Die zweite erwähnte Möglichkeit, das Risiko lokaler Minima und flacher Plateaus zu minimieren, ist die Einführung des sog. Momentum-Terms, der sich am ehesten durch Impuls oder Trägheitsmoment übersetzen lässt. Dieses Konzept stammt aus der Physik. Im Prinzip kann man sich den Suchraum (Fehlerfunktion) als zerklüftetes Gebirge vorstellen, mit vielen Gipfeln und Tälern und möglicherweise auch Plateaus. Durch die Veränderung der Gewichte bewegt man sich in diesem Suchraum von Gipfel zu Tal; der Fehler wird minimiert. Durch die Stärke der Gewichtsänderung ist es möglich zu erkennen, wie stark das Gefälle ist. Man stelle sich eine Kugel vor, die in dem vorhandenen Gebirge einen Hang hinunterrollt. Je stärker das Gefälle ist, desto höher ist die Beschleunigung, die die Kugel erreicht. Erreicht die Kugel das Tal, stoppt sie nicht sofort. Der Schwung bewirkt, dass die Bewegung erst allmählich abklingt. Aufgrund dessen ist es mit Hilfe des Momentum-Terms möglich, flache Plateaus zu überqueren oder zerklüftete Flächen zu passieren. Hierbei werden bei der Berechnung der neuen Gewichtsveränderungen, die Änderungen des letzten Schrittes berücksichtigt. Diese Änderungen wirken sich, mit einer gegebenen Stärke des Momentums, auf die neuen Gewichtswerte aus. Mathematisch ist dies durch (11) definiert.

∆wi , j (t ) = η δ jο i + α ⋅ ∆wi , j (t − 1) , wobei

(11)

∆wi , j (t ) die Gewichtsveränderung im aktuellen Schritt, ∆wi , j (t − 1) die des vorherigen Schrittes, δ j den Wert des unterhalb des Gewichtes liegenden Neurons j, ο i die Ausgabe des darüber liegenden Neurons i und α den Momentum-Parameter darstellen. Wie man anhand von (11) erkennen kann, ist die Stärke des Momentums durch den Parameter

α modifizierbar. (11) stellt also eine um das Momentum erweiterte Form von (6) dar. Trotz der offensichtlichen Vorteile, die der Momentum-Term mit sich bringt, sollte dennoch darauf hingewiesen werden, dass dieser auch die Gefahr birgt, durch einen hohen Schwung © modlab® 2007

Künstliche neuronale Netze

61

aus einem besseren lokalen Optimum in ein schlechteres zu gelangen. Leider gibt es keinen allgemeinen, problemunspezifischen Wert für Parameter α. Dessen Wert muss deswegen durch trial-and-error Verfahren bestimmt werden. Generell ist zu sagen, dass man beim Training neuronaler Netze eine Multistartstrategie wählen und mehrere Netze mit verschiedenen Parametereinstellungen parallel trainieren sollte. Die Selektion des besten Netzes ist dabei durch den RMSE (Abkürzung für root mean square error) möglich. Dessen Berechnung ist in (12) gezeigt.

RMSE =

1 n (t i − oi ) 2 , wobei ∑ n i =1

(12)

n die Anzahl der Trainingsbeispiele, ti der Zielwert und oi der berechnete Output-Wert des Netzes für Trainingsbeispiel i ist.

Literatur [1] W. Oberhofer, T. Zimerer 1996: „Wie Künstliche Neuronale Netze Lernen: Ein Blick in die Black Box der Backpropagation Netzwerke“ http://www.fh-ansbach.de/~thomas.zimmerer/pdfs/Literatur/DP_287.pdf [2] P. Baldi, S. Brunak 2001: „Bioinformatics - The Machine Learning Approach“, The MIT Press, Cambridge. [3] D. Kriesel 2005: „Ein kleiner Überblick über Neuronale Netze” http://www.dkriesel.com/content/studium/neuronetzefiles/neuronalenetze.pdf

[4] W. Schiffmann, M.Joost, R.Werner: „Optimization of the Backpropagation Algorithm for Training Multilayer Perceptrons” http://www.cs.bham.ac.uk/~pxt/NC/schiffmann.bp.pdf

© modlab® 2007

Künstliche neuronale Netze

© modlab® 2007

62

Partikelschwarm-Optimierung

63

Partikelschwarm-Optimierung (PSO) Die Idee für die PSO hat ihren Ursprung in der Simulation der Dynamik von biologischen Systemen, wie Vogel- oder Fischschwärmen. Relativ früh schon wurde jedoch die Eignung des

Algorithmus

als

heuristisches

Optimierungsverfahren

erkannt

und

daraufhin

weiterentwickelt. Zum ersten Mal wurden Partikelschwärme von Eberhart und Kennedy im Jahr 1995 als Optimierungsparadigma vorgeschlagen [1, 2]. Künstliche Schwärme bestehen aus Individuen, den so genannten Agenten oder Partikeln. Indem den Partikeln einfache Regeln vorgegeben werden, nach denen sie sich zu verhalten haben, können verschiedene Verhaltensmuster des Schwarmes entstehen. Trotz der Einfachheit des Regelwerkes und somit der Einfachheit des Verhaltens der einzelnen Individuen kann das Verhalten eines Schwarmes sehr komplex und vielschichtig sein. Solche Komplexität entwickelt sich aus der Fähigkeit der Partikel, miteinander zu interagieren. Indem sie untereinander Informationen über ihre Umwelt und über sich selbst austauschen, beeinflussen sie sich auch gegenseitig in ihrem Verhalten. Durch diesen Informationsaustausch können die Partikel gegenseitig ihre „Flugrichtung“ und „Fluggeschwindigkeit“

im

durch

das

Problem

aufgespannten

Suchraum

(oder

Fitnesslandschaft) beeinflussen. Jeder Partikel hat das Bestreben, sich auf optimale, im Raum gefundene Punkte hinzuzubewegen, so dass sich letztlich der Schwarm gerichtet durch den Raum bewegt und nach einiger Zeit in einer Lösung konvergiert. Eine künstliche Optimierungslandschaft ist in Abb. 1 gezeigt.

© modlab® 2007

64

nswe

n2 s io en

e

Funkt io

m Di

Dim

io ns

rt

Partikelschwarm-Optimierung

n1

Abbildung 1: Beispiel einer Fitnesslandschaft mit vielen lokalen Maxima und Minima (multimodale Fitnesslandschaft)

Der Algorithmus Zu Anfang der PSO wird in dem von einem Problem aufgespannten Lösungsraum eine bestimmte Anzahl von Partikeln generiert. Die Positionen dieser Partikel sind zunächst zufällig gewählt, ebenso ihre Bewegungsrichtung und Geschwindigkeit. Jeder Partikel des Schwarmes kann sich frei im Suchraum bewegen (Abb. 2).

Abbildung 2: Momentaufnahme: Partikel in einer Fitnesslandschaft, das globale Optimum ist durch den Schnittpunkt der schwarzen Geraden mit der Lösungsoberfläche visualisiert.

Bei einem D-dimensionalen Problem wird die Position jedes Partikels – hier des i-ten Partikels – durch einen Ortsvektor Xi = (xi1, xi2, xi3,…xiD) beschrieben. Die Geschwindigkeit des i-ten Partikels wird analog mit einem Geschwindigkeitsvektor Vi = (vi1, vi2, vi3,…viD) beschrieben. Da jedem Partikel durch den Ortsvektor eine Position im Lösungsraum © modlab® 2007

Partikelschwarm-Optimierung

65

zugewiesen ist, ist es auch möglich, die jeweilige Fitness der Partikel in Bezug auf das Problem zu bestimmen. Des Weiteren existieren zwei Arten von „Gedächtnis“ in der PSO:

•

Zum einen haben Partikel ein so genanntes „individuelles“ Gedächtnis, in dem sie die beste Position pi abspeichern können, an der sie sich jeweils selbst im Optimierungsprozess bis zum aktuellen Zeitpunkt befunden hat.

•

Zum anderen besitzen die Partikel ein „soziales“ Gedächtnis, das man sich auch als „Schwarmgedächtnis“ vorstellen kann. Dieses „soziale“ Gedächtnis speichert die Position pb ab, die bis zum aktuellen Zeitpunkt die insgesamt beste, von allen Partikeln zusammen gefundene Position ist.

Ausgehend von diesen beiden Informationen – also den Positionen aus den individuellen Gedächtnissen und dem sozialen Gedächtnis – sowie dem aktuellen Geschwindigkeitsvektor werden nun für jeden Partikel neue Geschwindigkeitsvektoren berechnet. In dem 1998 von Shi und Eberhart vorgeschlagenen PSO-Algorithmus [3] erfolgt die Berechnung nach Gleichung 1:

v i (t + 1) = w ⋅ v i (t ) + n1 ⋅ r1 ⋅ ( p i − x i (t )) + n 2 ⋅ r2 ⋅ ( p b − x i (t )) ,

(1)

wobei t für die Epochenzahl steht, n1 und n2 Konstanten zur Gewichtung des Einflusses der beiden Gedächtnisse auf das Schwarmverhalten sind, und r1 und r2 Zufallszahlen zwischen 0 und 1 sind. Das so genannte „Massenträgheitsgewicht“ w gewichtet den Einfluss des alten Geschwindigkeitsvektors vi eines Partikels und dient maßgeblich der Kontrolle der Gesamtgeschwindigkeit, mit der sich Partikel durch den Suchraum bewegen. „Vernünftige“ Werte für die Konstanten, die jedoch von Anwendung zu Anwendung unterschiedlich passend sein können, liegen für n1 und n2 zwischen 1 und 2, für w zwischen 0,5 und 1. Diese Angaben sind jedoch nur als Anhaltspunkte zu verstehen. Um einer unkontrollierten, explosionsartigen Bewegung des Schwarmes vorzubeugen, wird oftmals der maximale Wert, den der Geschwindigkeitsvektor annehmen kann, eingeschränkt. Wird dieser maximale Wert überschritten, so wird vi auf den Grenzwert zurückgesetzt. Eine weitere Möglichkeit, die Geschwindigkeit des Schwarmes zu kontrollieren, ist das Massenträgheitsgewicht mit zunehmender Epochenzahl linear kleiner werden zu lassen, was in der Regel zu einer Konvergenz des Schwarmes über die Zeit führt [4].

© modlab® 2007

Partikelschwarm-Optimierung

66

Die Positionen der einzelnen Partikel werden durch Addition des in jeder Epoche neu errechneten Geschwindigkeitsvektors zu dem Ortsvektor neu berechnet (Gleichung 2):

xi (t + 1) = xi (t ) + vi (t + 1) .

(2)

Nach dem Update der Positionen kann ein neuer Evaluationszyklus beginnen, bis ein Stoppkriterium (wie z.B. eine maximale Anzahl an Epochen) erreicht wird. In Pseudocode lässt sich die PSO folgendermaßen beschreiben:

begin Positionen und Geschwindigkeiten der Partikel initialisieren; Fitness evaluieren; Gedächtnisse initialisieren; while (Abbruchbedingung ≠ true) neue Geschwindigkeiten berechnen; Positionen neu berechnen; Fitness evaluieren; Gedächtnisse auffrischen; end end

Zusammenfassung Die Partikelschwarm-Optimierung ist ein biologisch motiviertes, leicht zu implementierendes und effizientes heuristisches Optimierungsverfahren. Im Vergleich zu anderen biologisch motivierten Optimierungsalgorithmen, wie beispielsweise genetischen Algorithmen [5], hat es nur wenige freie Parameter, die jedoch entscheidenden Einfluss auf das Verhalten der Schwärme haben. PSO wurde in der Vergangenheit bereits erfolgreich für unterschiedlichste Anwendungen eingesetzt, vom Training künstlicher neuronaler Netze bis hin zum Lösen technischer Probleme im Ingenieurswesen.

© modlab® 2007

Partikelschwarm-Optimierung

67

Literatur 1.

Kennedy J, Eberhart RC: Particle swarm optimization. In Proceedings of IEEE International Conference on Neural Networks; Piscataway, NJ. 1995: 1942-1948.

2.

Eberhart RC, Kennedy J: A new optimizer using particle swarm theory. In Proceedings of the Sixth International Symposium on Micromachine and Human Science; Nagoya, Japan. 1995: 39-43.

3.

Shi Y, Eberhart RC: Parameter selection in particle swarm optimization. In Evolutionary Programming VII: Proceedings of the Seventh Annual Conference on Evolutionary Programming; New York, USA. 1998: 591-600.

4.

Kennedy J, Eberhart RC: Swarm Intelligence. San Diego: Academic Press; 2001.

5.

Holland J.H.: Adaptation in natural and artificial system, Ann Arbor, The University of Michigan Press; 1975.

© modlab® 2007

Partikelschwarm-Optimierung

© modlab® 2007

68

Ameisenalgorithmen

69

Ameisenalgorithmen Ameisenalgorithmen stammen aus dem Gebiet der naturanalogen bzw. biologisch motivierten Algorithmen. Das Ant Colony Optimization (ACO) Konzept wurde von Marco Dorigo (Dorgio et al., 1996) als Metaheuristik eingeführt (http://www.aco-metaheuristic.org/). Die Idee

des

Ameisenalgorithmus

wurde

dem

Foraging

(Futtersuche)

bestimmter

Ameisenunterfamilien nachempfunden. Diese sind in der Lage, bestehende Hindernisse zwischen Nest und Futterquelle bei verschiedenen möglichen Wegen kurzmöglichst zu verbinden (Deneubourg et al., 1990; Abb.1).

Abbildung 1. Zeitlicher Ablauf der Wegfindung bei Ameisen. Die Graustufe der Pfeile spiegelt die Pheromonkonzentration wieder. Die Ameisen verwenden zuerst alle zur Verfügung stehenden Verbindungen zwischen Futterquelle und Nest (Abb. 1a). Über die Zeit (Abb. 1b) ist ein autokatalytischer Prozess zu beobachten, der immer mehr Arbeiterinnen auf den kürzeren Weg zieht, bis schlussendlich (Abb. 1c) alle Wege außer dem kürzesten ignoriert werden. Diese Lösung eines natürlichen Optimierungsproblems ist unter den folgenden Gesichtspunkten für die Informatik interessant: © modlab® 2007

Ameisenalgorithmen

70

•

Keine übergeordnete Intelligenz (Königin) steuert die Ameisen.

•

Einzelne Ameisen haben kein Problembewusstsein.

•

Die optimale Lösung wird gefunden.

•

Die Ameisen kooperieren, um eine Lösung zu finden.

Der Aspekt der Kooperation ist von entscheidender Bedeutung. Den Ameisen ist es möglich, über abgegebene Duftstoffe (Pheromone) über ihre Umwelt als Medium zu kommunizieren. Diese Kommunikation über eine Veränderung der Umwelt, um einem gemeinsamen Ziel zu folgen, bezeichnet man als Stigmergie, eingeführt von Grassé 1959 für Termiten (Grassé, P. P. & Noirot, C., 1959; Theraulaz, G. & Bonabeau, E., 1999). Die Ameisen lösen das in Abb.1 gezeigte Weg-Problem mit Hilfe von Pheromonen (Deneubourg et al., 1990). Auf dem kürzeren Weg wandern pro Zeit mehr Ameisen als auf dem längeren Weg, wodurch die Pheromonkonzentration auf diesem schneller anwächst (bzw. weniger abnimmt). Die aus dem Nest nachfolgenden Ameisen werden durch die Pheromonkonzentration hinsichtlich ihrer Entscheidung für einen Weg zur Futterquelle beeinflusst. Die kürzere Route wird durch die höhere Pheromonkonzentration mit einer gewissen Wahrscheinlichkeit bevorzugt. Einzelne gleichartige Individuen werden also über eine Kommunikationsschnittstelle verbunden und dadurch in die Lage versetzt, Probleme zu lösen. Es entsteht eine Lösung aus dem System, die über die Fähigkeiten der einzelnen Ameisen hinausgeht (Emergenz). Dies wird

durch

die

Ameisenalgorithmen

für

die

Informatik

im

Bereich

von

Optimierungsproblemen adaptiert. Ameisen werden formal zu „Agenten“. Der Aspekt der Pheromone wird 1:1 übernommen.

Ameisenalgorithmen in der Anwendung ACS (Ant Colony System) für das Travelling Salesman Problem nach Dorigo & Gambardella (1997). Im Problem des Handlungsreisenden (TSP) geht es darum, eine möglichst kurze Rundreise zwischen verschiedenen Städte zu organisieren, wobei jede Stadt nur einmal besucht werden darf. Ein Beispiel ist in Abb.2 zu sehen.

© modlab® 2007

Ameisenalgorithmen

71

Abbildung 2. Vollständiger verbundener Graph für das TSP Problem. Alle Kanten haben unterschiedliche Gewichte (zufällig) alle Knoten eindeutige Bezeichnungen. In Abb.2 stellen Knoten, Städte und Kanten, Straßen mit unterschiedlichem Gewicht (Länge) entsprechen. Ein Agent (k) startet nun in einer zufälligen Stadt und wählt von dieser Stadt (r) die nächste noch nicht besuchte Stadt (s) entsprechend Gleichung (1).

{

 arg max [τ (r , u )] ⋅ [η (r , u )]β  u∉M k s=  S

}

wenn q ≤ q0

, wobei

(1)

sonst

τ (r, u ) der Pheromonkonzentration auf der Kante zwischen Stadt r und u entspricht. η (r, u ) ist eine heuristische Funktion, invers zu der Distanz zwischen den Städten r und u. β gewichtet die relative Bedeutung der Pheromonspur sowie die Entfernung (Kantengewicht) zwischen den Städten. q wird mit gleichverteilter Wahrscheinlichkeit aus [0,1] gezogen, q 0 ist ein Parameter, S ist eine zufällig nach Gleichung 2 belegte Variable.

 [τ (r , s )] ⋅ [η (r , s )]β  β pk (r , s ) =  ∑ [τ (r , u )] ⋅ [η (r , u )] u ∉M k   0

wenn s ∉ M k .

(2)

sonst

Die Formel bevorzugt Kanten mit einer stärkeren Pheromonspur und geringerem Gewicht (Länge). Die Idee ist, dass der Agent mit einer Wahrscheinlichkeit von q 0 die Möglichkeit hat, auf das kollektive Gedächtnis zuzugreifen, oder mit einer Wahrscheinlichkeit von 1 − q 0 eine „eigene“ Stadt nach den Vorgaben von Gleichung 2 zu wählen. Für eine weiterführende Beschreibung des Pheromonupdates (globales, lokales) sei auf Dorigo & Gambardella (1997) verwiesen. © modlab® 2007

Ameisenalgorithmen

72

Ein Applet für TSP nach dem ACO Konzept ist hier zu finden: http://www.geocities.com/andres_perezuribe/ Es ist zu beachten, dass das hier verwendete Pheromonupdate explizit auf die Problemstellung zugeschnitten ist. Das ACO Konzept lässt sich jedoch auf eine weites Spektrum von Optimierungsproblemen anwenden. Für eine Auflistung möglicher Anwendungen sei auf Dorigo & Stützle, 2004 verwiesen. Die Herausforderung bei der Anwendung von ACO liegt darin, a) eine geeignete Methode des Pheromonupdates zu finden, und b)

das „Pfad“-Konzept an die Problemstellung anzupassen.

Naheliegende Anwendungen sind Verkehrssimulationen oder Routing Probleme. Auch die Möglichkeit der Modellierung biologischer Systeme (echter Ameisen) stellt eine mögliche Anwendung da.

Literatur Deneubourg, J-L., Aron, S., Goss, S. & Pasteels, J.M. (1990). The self-organizing exploratory pattern of the Argentine ant. Journal of Insect Behavior 3: 159-168. Dorigo, M., Maniezzo, V. & Colorni, A. (1996). Ant system: optimization by a colony of cooperating agents. IEEE Transactions on Systems, Man, und Cybernetics-Part B. 26: 29-41. Dorigo M. & Gambardella, L.M. (1997). Ant Colonies for the Traveling Salesman Problem. BioSystems 43:73-81. Also Tecnical Report TR/IRIDIA/1996-3, IRIDIA, Université Libre de Bruxelles. Dorigo, M. & Stützle, T. (2004). Ant Colony Optimization. The MIT Press, Cambridge, Massachusetts. Grassé, P. P. & Noirot, C. (1959). The development of symbiosis in Isoptera. Experientia 15: 365-72. Theraulaz, G. & Bonabeau, E. (1999). A brief history of stigmergy. Artificial Life 5: 97-116.

© modlab® 2007