[PYTHON] Seien Sie vorsichtig, wenn Sie die Eigenvektoren einer Matrix unterscheiden
Zusammenfassung
Bei der Differenzierung der Eigenwerte einer Matrix divergiert sie, wenn die Eigenwerte entartet sind. Gehen wir angemessen damit um.
Auslösen
Vorheriger Artikel
https://qiita.com/sage-git/items/1afa0bb8b3a7ee36600d
Ich dachte, wenn die Eigenwerte differenziert und optimiert werden könnten, könnten die Eigenvektoren auf die gleiche Weise erzeugt werden, und eine Matrix, die die gewünschten Eigenvektoren ergeben würde, würde erhalten, also habe ich es versucht. Ich habe es tatsächlich verstanden, aber es gab ein kleines Problem, also habe ich mich damit befasst.
Mathematik
Informationen zur Ableitung finden Sie in diesem Forum. https://mathoverflow.net/questions/229425/derivative-of-eigenvectors-of-a-matrix-with-respect-to-its-components
Betrachten Sie den Eigenwert $ \ lambda_i $ für eine Matrix $ A $ und den entsprechenden Eigenvektor $ \ vec {n} _i $ und differenzieren Sie diesen.
Differenzierung von Eigenwerten
\frac{\partial \lambda_i}{\partial A_{kl}} = \left(\frac{\partial A}{\partial A_{kl}}\vec{n}_i\right)\cdot\vec{n}_i
Hier sieht $ \ frac {\ partielles A} {\ partielles A_ {kl}} \ vec {n} \ _i $ wie ein Projektionsoperator aus. Bringen Sie die $ k $ -te Komponente von $ \ vec {n} \ _i $ zur $ l $ -ten Komponente, und der Rest ist ein Vektor von 0. Dann ist diese Unterscheidung der Wert von $ \ vec {n} _i $ multipliziert mit dem $ k $ th und dem $ l $ th. Es ist einfacher als ich erwartet hatte.
Differenzierung des Eigenvektors
Für den $ i $ -ten Eigenvektor $ \ vec {n} \ _i $
\frac{\partial \vec{n}_i}{\partial A_{kl}} = \sum_{j\neq i}\left[\frac{1}{\lambda_i - \lambda_j}\left(\frac{\partial A}{\partial A_{kl}}\vec{n}_i\right)\cdot\vec{n}_j\right]\vec{n}_j
Kann geschrieben werden. Kurz gesagt, es fühlt sich an, als würde man anderen Eigenvektoren geeignete Gewichte hinzufügen. Der hier zu beachtende Punkt ist der Abschnitt $ \ frac {1} {\ lambda_i- \ lambda_j} $. Dies führt dazu, dass dieses Differential divergiert, wenn mehrere identische Eigenwerte vorhanden sind (physikalisch entsprechend der Regression).
Überprüfung der Konten
Vorbereitung
Bestimmen Sie die geeignete Matrix "X" und überprüfen Sie die Eigenwerte und Eigenvektoren.
import tensorflow as tf
X = tf.Variable(initial_value=[[1.0, 0.0, 0.12975], [0.0, 1.0, 0.0], [0.12975, 0.0, 1.1373545]])
eigval, eigvec = tf.linalg.eigh(X)
print(eigval)
print(eigvec)
eigval (eindeutiger Wert)
tf.Tensor([0.9218725 1. 1.2154819], shape=(3,), dtype=float32)
eigvec (eindeutiger Vektor)
tf.Tensor(
[[-0.8566836 0. 0.51584214]
[-0. -1. 0. ]
[ 0.51584214 0. 0.8566836 ]], shape=(3, 3), dtype=float32)
Differenzieren Sie den Eigenwert
Versuchen Sie, den minimalen Eigenwert mit "GradientTape" zu berechnen.
with tf.GradientTape() as g:
g.watch(X)
eigval, eigvec = tf.linalg.eigh(X)
Y = eigval[0]
dYdX = g.gradient(Y, X)
print(dYdX)
Differenzierung des Eigenwertes 0
tf.Tensor(
[[ 0.7339068 0. 0. ]
[ 0. 0. 0. ]
[-0.88382703 0. 0.2660931 ]], shape=(3, 3), dtype=float32)
0.7339068 \approx 0.8566836^2 0.2660931 \approx 0.51584214^2 0.88382703 \approx 2 \times 0.8566836 \times 0.51584214
Es ist also ein vernünftiges Ergebnis. Es scheint, dass $ 2 \ times $ darauf zurückzuführen ist, dass die symmetrische Matrix nur die untere Hälfte verwendet.
Differenzieren Sie den Eigenvektor
eigvec [i, j] ist die $ i $ -te Komponente des Eigenvektors für den $ j $ -ten Eigenwert.
with tf.GradientTape() as g:
g.watch(X)
eigval, eigvec = tf.linalg.eigh(X)
Y = eigvec[0, 1]
dYdX = g.gradient(Y, X)
print(dYdX)
Erste Komponente des Eigenvektors 1
tf.Tensor(
[[ 0. 0. 0. ]
[-8.158832 0. 0. ]
[ 0. 7.707127 0. ]], shape=(3, 3), dtype=float32)
Das ist ärgerlich, deshalb werde ich den Scheck überspringen.
Bis zu diesem Punkt ist normal.
Schrumpfen
Wenn Sie den Wert von "X" wie folgt ändern, sind die beiden eindeutigen Werte "1".
X = tf.Variable(initial_value=[[1.1225665, 0.0, 0.12975], [0.0, 1.0, 0.0], [0.12975, 0.0, 1.1373545]])
Ich habe den Code aus Vorheriger Artikel verwendet, um ihn zu finden.
eigval
tf.Tensor([1. 1. 1.2599211], shape=(3,), dtype=float32)
eigvec
tf.Tensor(
[[-0.7269436 0. 0.6866972]
[-0. -1. 0. ]
[ 0.6866972 0. 0.7269436]], shape=(3, 3), dtype=float32)
Beides differenzieren
dYdX
tf.Tensor(
[[nan 0. 0.]
[nan nan 0.]
[nan nan nan]], shape=(3, 3), dtype=float32)
ist geworden. Es ist unverständlich, dass der Eigenwert "nan" wird, aber die Differenzierung des Eigenvektors wird "nan" gemäß der theoretischen Formel.
Es ist zu beachten, dass der zu differenzierende Eigenvektor auch der dritte ist, dh die Differenzierung des Eigenvektors bei $ i $, wobei $ j $ nicht existiert, so dass $ \ lambda_i- \ lambda_j = 0 $ ebenfalls nan ist. Mit anderen Worten, in der Tensorflow-Implementierung scheinen alle Differentiale in einer Matrix mit Regression gnadenlos "nan" zu sein.
Gegenmaßnahmen
Es gibt mehrere mögliche Problemumgehungen.
- Transformieren Sie die Matrix im Bereich gleicher Werte im Problem
- Verschieben Sie den Wert mit der Erwartung, dass sich der Wert geringfügig ändert
- Schreiben Sie den Algorithmus in eine Richtung um, in der die Differenzierung von Eigenwerten und Eigenvektoren überhaupt nicht verwendet wird
Hier werde ich zufällig eine Störung geben. Da diese Differenzierung für die Gradientenmethode verwendet wird, kann sie wie eine Art Tempern gestört werden, ohne das Endergebnis zu beeinflussen. Natürlich muss man besonders überlegen, ob das Endergebnis ist, dass die Eigenwerte dekrementiert werden.
Entnahme finden
In jedem Fall sollten Sie herausfinden, ob es denselben Eigenwert gibt, bevor Sie differenzieren.
eigval[1:] - eigval[:-1]
Dadurch erhalten Sie ein Array, das die Unterschiede zwischen den eindeutigen Werten nebeneinander enthält. Wenn wir die Tatsache ausnutzen, dass die von tf.linalg.eigh zurückgegebenen eindeutigen Werte bereits in aufsteigender Reihenfolge sortiert sind, können wir sehen, dass sie $ 0 $ oder mehr sind, ohne absolute Werte zu verwenden. Und wenn auch nur eine von ihnen fast 0 Komponenten hat, wird angenommen, dass sie schrumpft.
tf.math.reduce_any((eigval[1:] - eigval[:-1]) < 1e-20)
Führen Sie danach unter dieser Bedingung eine Schleife durch, bis sie nicht mehr erfüllt ist. "A" ist eine Matrix aus "N" Zeilen und "N" Spalten.
eigval, eigvec = tf.linalg.eigh(A)
while tf.math.reduce_any((eigval[1:] - eigval[:-1]) < 1e-20):
Ap = A + tf.linalg.diag(tf.random.uniform(shape=[N])*0.001)
eigval, eigvec = tf.linalg.eigh(Ap)
Ich denke, dass sich das Kriterium "1e-20" und die Stärke der Störung "0,001" je nach Problem ändern werden. Dies löste vorerst das, was ich jetzt tun wollte.
Bonus
Berechnen wir eine eindimensionale Quantenmulde. Die physikalische Erklärung ist
https://qiita.com/sage-git/items/e5ced4c0f555e2b4d77b
Geben Sie es auf eine andere Seite.
Betrachtet man das Potential $ U (x) $ so, dass der Bereich von $ x \ in \ left [- \ pi, \ pi \ right] $ endlich ist, sonst ist es $ + \ infty $. Wenn Sie $ U (x) $ setzen, ist die Wellenfunktion im Basiszustand
\psi(x) = A\exp\left(-2x^2\right)
Finden Sie numerisch heraus, ob Die Methode besteht darin, Hamiltonian $ H $ in eine Matrix zu schreiben, diesen Eigenwert / Eigenvektor zu finden und ihn mit der Gradientenmethode zu erreichen, so dass der Eigenvektor für den kleinsten Eigenwert zu diesem $ \ psi (x) $ wird. Bei der numerischen Berechnung kann die Funktion jedoch nicht als Funktion behandelt werden, sodass der Bereich $ \ left [- \ pi, \ pi \ right] $ durch $ N $ Punkte geteilt wird. Wenn die Koordinaten des $ i $ -ten Punktes $ x_i $ sind, ist die Wellenfunktion der Vektor $ \ vec {\ psi} der $ i $ -ten Komponente von $ \ vec {\ psi} $ = \ psi (x_i) $. Sie können es mit $ schreiben.
Eine Sache, die zu beachten ist, ist, dass das Vorzeichen des Eigenvektors einen gewissen Freiheitsgrad aufweist, also die Verlustfunktion
L_+ = \sum_i^N(n_i - \psi(x_i))^2
L_- = \sum_i^N(n_i + \psi(x_i))^2
Beides kann berücksichtigt werden. Es ist üblich, dass die Iteration beim Drehen umgedreht wird. Diesmal die kleinste davon
L = \min(L_+, L_-)
Dann hat es funktioniert.
Basierend auf diesen habe ich den folgenden Code geschrieben.
import tensorflow as tf
import matplotlib.pyplot as plt
import numpy as np
def main():
max_x = np.pi
N = 150 + 1
dx = max_x*2/N
x = np.arange(-max_x, max_x, dx)
D = np.eye(N, k=1)
D += -1 * np.eye(N)
D += D.T
D = D/(dx**2)
m = 1.0
D_tf = tf.constant(D/(2.0*m), dtype=tf.float32)
V0_np = np.exp( - x**2/0.5)
V0_np = V0_np/np.linalg.norm(V0_np)
V0_target = tf.constant(V0_np, dtype=tf.float32)
U0 = np.zeros(shape=[N])
U = tf.Variable(initial_value=U0, dtype=tf.float32)
def calc_V(n):
H = - D_tf + tf.linalg.diag(U)
eigval, eigvec = tf.linalg.eigh(H)
while tf.math.reduce_any((eigval[1:] - eigval[:-1]) < 1e-20):
H = - D_tf + tf.linalg.diag(U) + tf.linalg.diag(tf.random.uniform(shape=[N])*0.001)
eigval, eigvec = tf.linalg.eigh(H)
print("found lambda_i+1 - lambda_i = 0")
v_raw = eigvec[:, n]
return v_raw
def calc_loss():
v0 = calc_V(0)
dplus = tf.reduce_sum((v0 - V0_target)**2)
dminus = tf.reduce_sum((v0 + V0_target)**2)
return tf.minimum(dplus, dminus)
opt = tf.keras.optimizers.Adam(learning_rate=0.001)
L = calc_loss()
v0_current = calc_V(0)
print(L)
while L > 1e-11:
opt.minimize(calc_loss, var_list=[U])
v0_current = calc_V(0)
L = (tf.abs(tf.reduce_sum(v0_current*V0_target)) - 1)**2
print(L)
plt.plot(x, U.numpy())
plt.show()
if __name__ == "__main__":
main()
Nachdem Sie dies getan und einige Minuten lang auf einer mit Ryzen 5 und GTX 1060 beladenen Maschine belassen hatten, wurde das unten gezeigte Diagramm erhalten.
Mit anderen Worten wurde gefunden, dass, wenn ein solches Potential eingestellt wird, die Wellenfunktion des Grundzustands in der Quantentopfwelle zu einem Gaußschen Wellenfluss wird.
Leider ist dies eine numerische Lösung und kann nicht analytisch bestätigt werden. Ich möchte es mit einer guten Funktion anpassen und lösen. Ich weiß nicht, ob Menschen das können.