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)

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.

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.

image.png

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.

Recommended Posts

Soyez prudent lors de la différenciation des vecteurs propres d'une matrice
[Mémo Python] Soyez prudent lors de la création d'un tableau à deux dimensions (liste de listes)
Faites attention au type lorsque vous créez un masque d'image avec Numpy
[Attention] Lors de la création d'une image binaire (1 bit / pixel), faites attention au format de fichier!
Soyez prudent lorsque vous attribuez une série en tant que colonne aux pandas.
Trouver les valeurs propres d'une vraie matrice symétrique en Python
Lors de l'incrémentation de la valeur d'une clé qui n'existe pas
Soyez prudent lorsque vous spécifiez la valeur d'argument par défaut dans la série Python 3
Faites attention à LANG pour UnicodeEncodeError lors de l'impression du japonais avec Python 3
[Python] Soyez prudent lorsque vous utilisez print
Trouvez le rang de la matrice dans le monde XOR (rang de la matrice sur F2)
Note addictive: max (max (list)) ne doit pas être utilisé lors de la maximisation de la valeur d'un tableau à deux dimensions
Lors de la création d'une matrice dans une liste
L'histoire de l'exportation d'un programme
Trouver l'intersection d'un cercle et d'une droite (matrice sympy)
Je voulais faire attention au comportement des arguments par défaut de Python
Choses à prendre en compte lors de la création d'un système de recommandation avec Item2Vec
Générer semi-automatiquement une description du package à enregistrer dans PyPI
Soyez prudent lorsque vous récupérez des tweets à intervalles réguliers avec l'API Twitter
Mesurer la force de l'association dans un tableau croisé
Prédire quand l'ISS sera visible
[python] [meta] Le type de python est-il un type?
Soyez prudent lors de l'ajout d'un tableau à un tableau
Un mémo expliquant la spécification de l'axe de l'axe
Obtenez le nom de fichier du répertoire (glob)
L'histoire du traitement A du blackjack (python)
Notez l'achèvement d'une commande chronophage
Soyez prudent lorsque vous exécutez CakePHP3 avec PHP7.2
Un mémorandum de problème lors du formatage des données
Tensorflow, il semble que même la valeur propre de la matrice puisse être automatiquement différenciée
[Langage C] Faites attention à la combinaison de l'API de mise en mémoire tampon et de l'appel système sans mise en mémoire tampon
L'histoire de Django créant une bibliothèque qui pourrait être un peu plus utile
[Numpy, scipy] Comment calculer la racine carrée d'une matrice Elmeet à valeur semi-régulière
Lors de l'augmentation de la mémoire disponible du nœud, elle peut être limitée par vm.max_map_count
Six mois, lorsque le système littéraire de Gorigori a appris l'IA presque par lui-même
Comment calculer la volatilité d'une marque
Récupérer l'appelant d'une fonction en Python
Visualisez la couche interne du réseau neuronal
Pourquoi la fonction d'activation doit être une fonction non linéaire
Copiez la liste en Python
Trouvez le nombre de jours dans un mois
Essayez de décomposer la matrice daimyo par valeur singulière
Écrire une note sur la version python de python virtualenv
Calculer la probabilité de valeurs aberrantes sur les moustaches de la boîte
[Python] Une compréhension approximative du module de journalisation
Sortie sous la forme d'un tableau python
L'histoire de la création d'un générateur d'icônes mel
[Python] Trouvez la matrice de translocation en notation d'inclusion
Prise en compte des forces et faiblesses de Python
Soyez prudent lorsque vous travaillez avec des fichiers texte compressés au format gzip
Lorsqu'un fichier est placé dans le dossier partagé de Raspberry Pi, le processus est exécuté.
Lors de la lecture d'un fichier csv avec read_csv de pandas, la première colonne devient index
Mesures à prendre en cas de caractères déformés lors de la tentative de redirection / canalisation du résultat de aws-cli
Une note de malentendu lors de la tentative de chargement de l'intégralité du module self-made avec Python3
Je souhaite être informé de l'environnement de connexion lorsque RaspberryPi se connecte au réseau