[PYTHON] Soyez prudent lors de la différenciation des vecteurs propres d'une matrice
Résumé
Lors de la différenciation des valeurs propres d'une matrice, il diverge si les valeurs propres sont dégénérées. Traitons-le de manière appropriée.
Déclencheur
Article précédent
https://qiita.com/sage-git/items/1afa0bb8b3a7ee36600d
Je pensais que si les valeurs propres pouvaient être différenciées et optimisées, alors les vecteurs propres pourraient être créés de la même manière, et une matrice qui donnerait les vecteurs propres souhaités serait obtenue, alors je l'ai essayé. Je l'ai effectivement compris, mais il y avait un petit problème, alors je l'ai réglé.
Math
Référez-vous à ce forum pour la dérivation. https://mathoverflow.net/questions/229425/derivative-of-eigenvectors-of-a-matrix-with-respect-to-its-components
Pour une matrice $ A $, considérez la valeur propre $ \ lambda_i $ et le vecteur propre correspondant $ \ vec {n} _i $, et envisagez de la différencier.
Différenciation des valeurs propres
\frac{\partial \lambda_i}{\partial A_{kl}} = \left(\frac{\partial A}{\partial A_{kl}}\vec{n}_i\right)\cdot\vec{n}_i
Ici, $ \ frac {\ partial A} {\ partial A_ {kl}} \ vec {n} \ _i $ ressemble à un opérateur de projection. Apportez le composant $ k $ th de $ \ vec {n} \ _i $ au composant $ l $ th, et le reste sera un vecteur de 0. Alors cette différenciation sera la valeur de $ \ vec {n} _i $ multipliée par les $ k $ th et les $ l $ th. C'est plus simple que ce à quoi je m'attendais.
Différenciation du vecteur propre
Pour le $ i $ ème vecteur propre $ \ 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
Peut être écrit. En bref, c'est comme ajouter des poids appropriés à d'autres vecteurs propres. Le point à noter ici est la section $ \ frac {1} {\ lambda_i- \ lambda_j} $. Cela fait diverger cette différenciation lorsqu'il existe plusieurs valeurs propres identiques (correspondant physiquement à la régression).
vérification des comptes
Préparation
Déterminez la matrice appropriée «X» et vérifiez les valeurs propres et les vecteurs propres.
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 (valeur unique)
tf.Tensor([0.9218725 1. 1.2154819], shape=(3,), dtype=float32)
eigvec (vecteur unique)
tf.Tensor(
[[-0.8566836 0. 0.51584214]
[-0. -1. 0. ]
[ 0.51584214 0. 0.8566836 ]], shape=(3, 3), dtype=float32)
Différencier la valeur propre
Essayez de calculer la valeur propre minimale en utilisant GradientTape.
with tf.GradientTape() as g:
g.watch(X)
eigval, eigvec = tf.linalg.eigh(X)
Y = eigval[0]
dYdX = g.gradient(Y, X)
print(dYdX)
Différenciation de la valeur propre 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
C'est donc un résultat raisonnable. Il semble que $ 2 \ times $ est dû au fait que la matrice symétrique n'utilise que la moitié inférieure.
Différencier le vecteur propre
ʻEigvec [i, j] `est le $ i $ ème composant du vecteur propre pour la $ j $ ème valeur propre.
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)
Premier composant du vecteur propre 1
tf.Tensor(
[[ 0. 0. 0. ]
[-8.158832 0. 0. ]
[ 0. 7.707127 0. ]], shape=(3, 3), dtype=float32)
C'est ennuyeux, alors je vais sauter le chèque.
Jusqu'à ce point, c'est normal.
Rétrécir
Si vous modifiez la valeur de «X» comme suit, les deux valeurs uniques seront «1».
X = tf.Variable(initial_value=[[1.1225665, 0.0, 0.12975], [0.0, 1.0, 0.0], [0.12975, 0.0, 1.1373545]])
J'ai utilisé le code de Article précédent pour le trouver.
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)
Différencier cela, à la fois
dYdX
tf.Tensor(
[[nan 0. 0.]
[nan nan 0.]
[nan nan nan]], shape=(3, 3), dtype=float32)
est devenu. Il est incompréhensible que la valeur propre devienne "nan", mais la différenciation du vecteur propre devient "nan" selon la formule théorique.
Notez que le vecteur propre à différencier est aussi le troisième, c'est-à-dire que la différenciation du vecteur propre en $ i $ où $ j $ n'existe pas de telle sorte que $ \ lambda_i- \ lambda_j = 0 $ est aussi nan. En d'autres termes, dans l'implémentation de Tensorflow, il semble que tous les différentiels soient impitoyablement «nan» dans une matrice avec régression.
Contre-mesures
Il existe plusieurs solutions de contournement possibles.
- Transformez la matrice dans la plage de valeurs égales dans le problème
- Déplacez la valeur avec l'espoir que la valeur changera légèrement
- Réécrire l'algorithme dans une direction qui n'utilise pas la différenciation des valeurs propres et des vecteurs propres en premier lieu
Ici, je vais donner une perturbation au hasard. Comme cette différenciation est utilisée pour la méthode du gradient, elle peut être perturbée comme une sorte de recuit sans affecter le résultat final. Bien sûr, vous devez penser spécialement si le résultat final est que les valeurs propres sont décrémentées.
Trouver un retrait
Dans tous les cas, envisagez de savoir s'il existe la même valeur propre avant de différencier.
eigval[1:] - eigval[:-1]
Cela vous donnera un tableau qui contient les différences entre les valeurs uniques les unes à côté des autres. Profitant du fait que les valeurs uniques renvoyées par tf.linalg.eigh sont déjà triées par ordre croissant, nous pouvons voir qu'elles sont $ 0 $ ou plus sans utiliser de valeurs absolues. Et si même l'un d'entre eux a presque 0 composant, on suppose qu'il rétrécit.
tf.math.reduce_any((eigval[1:] - eigval[:-1]) < 1e-20)
Après cela, avec ceci comme condition, bouclez jusqu'à ce qu'il ne soit pas satisfait. «A» est une matrice de «N» lignes et «N» colonnes.
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)
Je pense que le critère «1e-20» et l'ampleur de la perturbation «0,001» changeront en fonction du problème. Pour le moment, cela a résolu ce que je voulais faire maintenant.
prime
Calculons un puits quantique unidimensionnel. L'explication physique est
https://qiita.com/sage-git/items/e5ced4c0f555e2b4d77b
Etc., donnez-le à une autre page.
Considérant le potentiel $ U (x) $ tel que la plage de $ x \ in \ left [- \ pi, \ pi \ right] $ soit finie, sinon c'est $ + \ infty $. Si vous définissez $ U (x) $, la fonction d'onde dans l'état de base sera
\psi(x) = A\exp\left(-2x^2\right)
Découvrez numériquement si La méthode consiste à écrire hamiltonien $ H $ dans une matrice, trouver cette valeur propre / vecteur propre et l'approcher par la méthode du gradient de sorte que le vecteur propre pour la plus petite valeur propre devienne ce $ \ psi (x) $. Cependant, dans le calcul numérique, la fonction ne peut pas être traitée comme une fonction, donc la plage $ \ left [- \ pi, \ pi \ right] $ est divisée par $ N $ points. Si les coordonnées du $ i $ ème point sont $ x_i $, la fonction d'onde est le vecteur $ \ vec {\ psi} du $ i $ ème composant de $ \ vec {\ psi} $ = \ psi (x_i) $ Vous pouvez l'écrire avec $.
Une chose à garder à l'esprit est qu'il y a un degré de liberté dans le signe du vecteur propre, donc la fonction de perte est
L_+ = \sum_i^N(n_i - \psi(x_i))^2
L_- = \sum_i^N(n_i + \psi(x_i))^2
Les deux peuvent être envisagés. Il est courant que l'itération se retourne lorsqu'elle est tournée. Cette fois, le plus petit d'entre eux
L = \min(L_+, L_-)
Ensuite, cela a fonctionné.
Sur cette base, j'ai écrit le code suivant.
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()
Après avoir fait cela et l'avoir laissé sur une machine avec Ryzen 5 et GTX 1060 pendant des dizaines de minutes, le graphique ci-dessous a été obtenu.
En d'autres termes, il a été constaté que si un tel potentiel est défini, la fonction d'onde de l'état fondamental dans le puits quantique deviendra un flux d'onde gaussien.
Malheureusement, il s'agit d'une solution numérique et ne peut être confirmée analytiquement. Je veux l'adapter et le résoudre avec une bonne fonction. Je ne sais pas si les humains peuvent le faire.