Taille du pas de descente de gradient adaptatif lorsque vous ne pouvez pas faire de recherche de ligne

9

J'ai une fonction objective E dépendante d'une valeur ϕ(x,t=1.0) , où ϕ(x,t) est la solution d'un PDE. J'optimise E par descente de gradient sur la condition initiale de la PDE: ϕ(x,t=0.0) . Autrement dit, je mets à jour ϕ(x,t=0.0)puis dois intégrer le PDE pour calculer mon résidu. Cela signifie que si je faisais une recherche de ligne pour la taille du pas de descente de gradient (appelez-la α ), pour chaque valeur potentielle de α je devrais réintégrer la PDE.

Dans mon cas, ce serait trop cher. Existe-t-il une autre option pour la taille de pas de descente de gradient adaptatif?

Je ne recherche pas seulement des schémas de principes mathématiques ici (bien que ce soit mieux si quelque chose existe), mais je serais satisfait de tout ce qui est généralement mieux qu'une taille de pas statique.

Merci!

NLi10Me
la source
Je ne pense pas que je veux modifier la façon dont j'intègre le PDE pour le moment, car pour moi ce serait une réécriture majeure du code. De plus, ce n'est pas tant que le PDE est délicat, car je dois le résoudre sur une grille très dense dans l'espace-temps car j'ai besoin d'une très grande précision numérique.
NLi10Me
D'un autre côté, la méthode BB (que je ne connaissais pas) semble assez bonne; tout ce que j'ai à faire, c'est de suivre l'état et le gradient de l'itération précédente et j'obtiens une approximation de second ordre ... cela semble très bien. Cependant, la dérivation suppose un quadratique convexe et mon problème ne l'est certainement pas. Cependant, je trouve aussi (et je suis satisfait) des minima locaux plutôt que globaux. Savez-vous à quel point BB a performé sur des problèmes dimensionnels très élevés?
NLi10Me
Je suppose que ce que je voulais dire sur les minima locaux est que, au voisinage d'un minimum local, aucune fonction n'est-elle approximativement quadratique? Je pense que mon état initial est suffisamment proche d'un minimum, car dans de nombreux cas, j'obtiens une convergence régulière même avec la taille de pas statique. Donc, même si sa dimension est très élevée, et en général, si vous considérez tout l'espace de recherche, le problème est non convexe / non quadratique, BB pourrait-il toujours être un bon choix sans recherche de ligne? ϕ(0)(x,t=0.0)
NLi10Me
Les autres «ingrédients» de sont des données d'images expérimentales. ϕ ( x , t = 1.0 ) tente de déformer une image pour qu'elle "corresponde" à l'autre (mesurée par des fonctions de correspondance comme la norme L2 intégrée sur les voxels). Pour certaines paires d'images, j'obtiens une convergence fluide avec (mon choix actuel) la taille de pas statique. Pour les autres paires d'images, j'obtiens beaucoup d'oscillation. Le système doit être entièrement automatisé, donc je ne peux pas revenir en arrière et modifier manuellement la taille des pas pour les paires d'images gênantes. Eϕ(x,t=1.0)
NLi10Me
À droite, je dois résoudre le système adjoint pour obtenir le gradient (qui est un système plus méchant et prend plus de temps). Ok, je pense que je vais essayer BB avec la recherche de ligne de retour en arrière. Merci très bien pour les conseils; mes conseillers sont souvent difficiles à trouver et beaucoup d'entre eux ne s'intéressent pas tant à la mise en œuvre qu'au modèle. Je trouve que les méthodes numériques sont l' élément crucial pour démontrer si un modèle est bon ou non en premier lieu, donc merci encore, je l'apprécie vraiment.
NLi10Me

Réponses:

15

Je vais commencer par une remarque générale: les informations de premier ordre (c'est-à-dire, en utilisant uniquement des gradients, qui codent la pente) ne peuvent que vous donner des informations directionnelles: elles peuvent vous dire que la valeur de la fonction diminue dans le sens de la recherche, mais pas pendant combien de temps . Pour décider jusqu'où aller dans la direction de la recherche, vous avez besoin d'informations supplémentaires (la descente de gradient avec des longueurs de pas constantes peut échouer même pour des problèmes quadratiques convexes). Pour cela, vous avez essentiellement deux choix:

  1. 1
  2. Essais et erreurs (j'entends par là bien sûr utiliser une recherche de ligne appropriée comme Armijo).

O(1)

α0>0g0:=f(x0)k=0,...

  1. sk=αk1gkxk+1=xk+sk
  2. gk+1=f(xk+1)yk=gk+1gk
  3. αk+1=(yk)Tyk(yk)Tsk

f(xk+1)f(xk)σk(0,αk1)

f(xkσkgk)maxmax(kM,1)jkf(xj)γσk(gk)Tgk,
γ(0,1)MM=10

Une approche alternative (et, à mon avis, bien meilleure) consisterait à utiliser déjà cette approximation aux différences finies dans le calcul de la direction de recherche; c'est ce qu'on appelle une méthode quasi-Newton . L'idée est de construire progressivement une approximation de la Hesse en utilisant des différences de gradients. Par exemple, vous pouvez prendre (la matrice d'identité) et pour résoudre et définissez avec comme ci-dessus et . (Cela s'appelle la mise à jour Broyden2f(xk)H0=Idk=0,

(1)Hksk=gk,
Hk+1=Hk+(ykHksk)T(sk)T(sk)Tsk
ykxk+1=xk+sket est rarement utilisé dans la pratique; une mise à jour meilleure mais légèrement plus compliquée est la mise à jour BFGS , pour laquelle - et plus d'informations - je me réfère au livre de Nocedal et Wright Numerical Optimization .) L'inconvénient est que a) cela nécessiterait la résolution d'un système linéaire à chaque étape (mais seulement de la taille de l'inconnu qui dans votre cas est une condition initiale, donc l'effort doit être dominé par la résolution des PDE pour obtenir le gradient; il existe également des règles de mise à jour pour les approximations de la Hesse inverse , qui ne nécessitent que le calcul d'une seule matrice -produit vectoriel) et b) vous avez encore besoin d'une recherche de ligne pour garantir la convergence ...

Heureusement, dans ce contexte, il existe une approche alternative qui utilise chaque évaluation de fonction. L'idée est que pour symétrique et positif défini (qui est garanti pour la mise à jour BFGS), la résolution de équivaut à minimiser le modèle quadratique Dans une méthode de région de confiance , vous le feriez avec la contrainte supplémentaire que , où est un rayon de région de confiance choisi de manière appropriée (qui joue le rôle de la longueur de pas ). L'idée clé est maintenant de choisir ce rayon de manière adaptative, en fonction de l'étape calculée. Plus précisément, vous regardez le rapport Hk(1)

qk(s)=12sTHks+sTgk.
sΔkΔkσk
ρk:=f(xk)f(xk+sk)f(xk)qk(sk)
de la réduction réelle et prévue de la valeur de la fonction. Si est très petit, votre modèle était mauvais et vous jetez et réessayez avec . Si est proche de , votre modèle est bon et vous définissez et augmentez . Sinon, vous définissez simplement et laissez seul. Pour calculer le minimiseur réel deρkskΔk+1<Δkρk1xk+1=xk+skΔk+1>Δkxk+1=xk+skΔkskminsΔkqk(s), il existe plusieurs stratégies pour éviter d'avoir à résoudre le problème d'optimisation entièrement contraint; mon préféré est la méthode CG tronquée de Steihaug . Pour plus de détails, je me réfère à nouveau à Nocedal et Wright.
Christian Clason
la source
Je suis en train de revoir cela et je me rends compte que j'ai une question. À l'étape trois pour la méthode BB, vous avez ; où et . Le numérateur et le dénominateur dans l'expression de ressemblent à des produits intérieurs. Dans mon cas, , où est un espace vectoriel avec une métrique riemannienne non triviale: K. Autrement dit, . Est-ce que cela affecte la définition de ? αk+1=(yk)Tyk(yk)Tskyk=gk+1gksk=αk1gkαk+1gkVVgk,gkV=gk,KgkL2αk+1
NLi10Me
Oui, si vous avez une structure d'espace vectoriel non triviale, vous devez respecter cela dans les algorithmes. En particulier, vous devez distinguer les produits internes de deux fonctions dans le même espace (par exemple, et ) et les produits de dualité entre une fonction dans l'espace et une dans le double espace (par exemple, et ) - pour ce dernier, vous devez d'abord inclure la cartographie Riesz pour la transformer en produit intérieur. (Cela peut être interprété comme un préconditionnement.)ykykskyk
Christian Clason
Dr.Clason, je soumettrai un article à l'ISBI 2017 détaillant certaines expériences que j'ai faites en utilisant la méthode de recherche de ligne BB + pour une tâche d'enregistrement d'images difféomorphes. Souhaitez-vous être inclus comme auteur sur le manuscrit? Je ne l'ai pas encore écrit, mais j'ai la plupart des expériences terminées ou en cours. S'il vous plaît, faites-moi savoir.
NLi10Me
@ NLi10Me Merci pour cette aimable offre, mais je n'ai rien fait qui mérite d'être coauteur - tout ce que j'ai écrit est du matériel standard. Si vous y tenez beaucoup, vous pouvez me remercier pour "des remarques utiles sur (tout ce qui a aidé)", mais même cela ne serait pas nécessaire. Savoir que ce que j'ai écrit était utile suffit!
Christian Clason
1
Désolé, vous avez raison, c'est une faute de frappe - corrigé! (La condition Armijo est souvent écrite comme , où est le sens de la recherche - pas nécessairement le négatif gradient - et la taille du pas, ce qui devrait clarifier ce qui se passe.)f(x+σs)f(x)γf(x)T(σs)sσ
Christian Clason