Problème d'optimisation contraint dans l'entropie matricielle

J'ai un problème d'optimisation contraint dans l'entropie de la matrice (Shannon) . La matrice peut être écrite comme la somme des matrices de rang 1 de la forme où est un vecteur normalisé donné. Les coefficients des matrices de rang 1 sont les inconnues dans lesquelles nous optimisons et ils doivent être supérieurs à zéro et totaliser 1. A [ v $\mathtt{(sum(entr(eig(A))))}$$A$v $[v_i\,v_i^T]$$v_i$

Dans une syntaxe de type CVX, le problème se présente comme suit: variable donnée$\mathtt{c(n)}$

$$\text{minimize} \qquad \mathtt{sum(entr(eig(A)))}$$

$$\begin{align} \text{subject to} \qquad A &= \sum c_i v_i v_i^T\\ \sum c_i &= 1\\ c_i &\ge 0\end{align}$$ .

Quelqu'un at-il une idée de la façon de résoudre ce problème efficacement? Je sais déjà qu'il ne peut probablement pas être converti en problème de programmation semi-définie (SDP).

Edit: Un collègue m'a informé que ma méthode ci-dessous est une instance de la méthode générale dans l'article suivant, lorsqu'elle est spécialisée dans la fonction d'entropie,

Overton, Michael L. et Robert S. Womersley. "Dérivées secondes pour optimiser les valeurs propres des matrices symétriques." SIAM Journal on Matrix Analysis and Applications 16.3 (1995): 697-718. http://ftp.cs.nyu.edu/cs/faculty/overton/papers/pdffiles/eighess.pdf

Aperçu

Dans cet article, je montre que le problème d'optimisation est bien posé et que les contraintes d'inégalité sont inactives à la solution, puis calcule les première et deuxième dérivées de Frechet de la fonction d'entropie, puis propose la méthode de Newton sur le problème avec la contrainte d'égalité éliminée. Enfin, le code Matlab et les résultats numériques sont présentés.

Bien posé du problème d'optimisation

Premièrement, la somme des matrices définies positives est définie positive, donc pour , la somme des matrices de rang 1 est définie positive. Si l'ensemble de est de rang complet, les valeurs propres de sont positives, de sorte que les logarithmes des valeurs propres peuvent être pris. Ainsi, la fonction objectif est bien définie à l'intérieur de l'ensemble réalisable.A ( c ) : = N ∑ i = 1 c i v i v T i v i$c_i > 0$

$$A(c):=\sum_{i=1}^N c_i v_i v_i^T$$ $v_i$$A$

Deuxièmement, comme tout , perd le rang, donc la plus petite valeur propre de passe à zéro. C'est-à-dire, comme . Puisque la dérivée de explose en tant que , on ne peut pas avoir une séquence de points successivement meilleurs et meilleurs approchant la frontière de l'ensemble faisable. Ainsi le problème est bien défini et de plus les contraintes d'inégalité sont inactives.A A σ m i n ( A ( c ) ) → 0 c i → 0 - σ log ( σ ) σ → 0 c i ≥ $c_i \rightarrow 0$$A$$A$$\sigma_{min}(A(c)) \rightarrow 0$$c_i \rightarrow 0$$-\sigma \log(\sigma)$$\sigma \rightarrow 0$$c_i \ge 0$

Dérivés de Frechet de la fonction d'entropie

À l'intérieur de la région réalisable, la fonction d'entropie est Frechet différenciable partout, et deux fois Frechet différenciable partout où les valeurs propres ne sont pas répétées. Pour faire la méthode de Newton, nous devons calculer les dérivées de l'entropie matricielle, qui dépend des valeurs propres de la matrice. Cela nécessite de calculer les sensibilités de la décomposition des valeurs propres d'une matrice par rapport aux changements dans la matrice.

Rappelons que pour une matrice avec décomposition de valeurs propres , la dérivée de la matrice de valeurs propres par rapport aux changements dans la matrice d'origine est, et la dérivée de la matrice de vecteur propre est, où est le produit de Hadamard , avec le coefficient matrice A = U Λ U T d Λ = I ∘ ( U T d A U ) , d U = U C ( d A ) , ∘ C = { u T i d A u $A$$A = U \Lambda U^T$

$$d\Lambda = I \circ (U^T dA U),$$ $$dU = UC(dA),$$ $\circ$ $$C = \begin{cases} \frac{u_i^T dA u_j}{\lambda_j - \lambda_i}, & i=j \\ 0, &i=j \end{cases}$$

Ces formules sont dérivées en différenciant l'équation des valeurs propres , et les formules sont valables chaque fois que les valeurs propres sont distinctes. Lorsqu'il y a des valeurs propres répétées, la formule de a une discontinuité amovible qui peut être étendue tant que les vecteurs propres non uniques sont choisis avec soin. Pour plus de détails à ce sujet, consultez la présentation et le document suivants .d $AU=\Lambda U$$d\Lambda$

La dérivée seconde est alors trouvée en différenciant à nouveau,

$$\begin{align} d^2 \Lambda &= d(I \circ (U^T dA_1U)) \\ &= I \circ (dU_2^T dA_1 U + U^T dA_1 dU_2) \\ &= 2 I \circ (dU_2^T dA_1 U). \end{align}$$

Alors que la première dérivée de la matrice des valeurs propres peut être rendue continue à des valeurs propres répétées, la deuxième dérivée ne peut pas puisque dépend de , qui dépend de , qui explose lorsque les valeurs propres dégénèrent l'une vers l'autre. Cependant, tant que la vraie solution n'a pas de valeurs propres répétées, alors c'est OK. Des expériences numériques suggèrent que c'est le cas pour le générique , bien que je n'ai pas de preuve à ce stade. Ceci est vraiment important à comprendre, car la maximisation de l'entropie tenterait généralement de rapprocher les valeurs propres si possible.d U 2 C v $d^2 \Lambda$$dU_2$$C$$v_i$

Éliminer la contrainte d'égalité

Nous pouvons éliminer la contrainte en travaillant uniquement sur les premiers coefficients et en réglant le dernier sur $\sum_{i=1}^N c_i = 1$$N-1$

$$c_N = 1-\sum_{i=1}^{N-1} c_i.$$

Globalement, après environ 4 pages de calculs matriciels, les dérivées premières et secondes réduites de la fonction objectif par rapport aux changements dans les premiers coefficients sont données par, où $N-1$

$$df = dC_1^T M^T [I \circ (V^T U B U^T V)]$$ $$ddf = dC_1^T M^T [I \circ (V^T[2dU_2 B_a U^T + U B_b U^T]V)],$$ $$M = \begin{bmatrix} 1 & \\ & 1 & \\ &&\ddots& \\ &&&1\\ -1 & -1 & \dots & -1 \end{bmatrix},$$

$$B_a = \mathrm{diag}(1+\log \lambda_1, 1 + \log \lambda_2, \ldots, 1 + \log \lambda_N),$$

$$B_b = \mathrm{diag}(\frac{d_2\lambda_1}{\lambda_1},\ldots,\frac{d_2\lambda_N}{\lambda_N}).$$

La méthode de Newton après élimination de la contrainte

Étant donné que les contraintes d'inégalité sont inactives, nous commençons simplement dans l'ensemble faisable et exécutons la région de confiance ou la recherche en ligne newton-CG inexacte pour la convergence quadratique vers les maxima intérieurs.

La méthode est la suivante (sans inclure les détails de la recherche par région de confiance / ligne)

1. Commencez par .$\tilde{c} = [1/N,1/N,\ldots,1/N]$
2. Construisez le dernier coefficient, .$c = [\tilde{c},1 - \sum_{i=1}^{N-1} c_i]$
3. Construct .$A = \sum_i c_i v_i v_i^T$
4. Trouvez les vecteurs propres et les valeurs propres de .$U$$\Lambda$$A$
5. Construire le gradient .$G = M^T [I \circ (V^T U B U^T V)]$
6. Résoudre pour via le gradient conjugué (seule la capacité d'appliquer est nécessaire, pas les entrées réelles). est appliqué au vecteur en trouvant , et puis en se à la formule, p H H δ ˜ c d U 2 B a B b M T [ I ∘ ( V T [ 2 d U 2 B a U T + U B b U T ] V ) $H G = p$$p$$H$$H$$\delta \tilde{c}$$dU_2$$B_a$$B_b$ $$M^T [I \circ (V^T[2dU_2 B_a U^T + U B_b U^T]V)]$$
7. Définissez .$\tilde{c} \leftarrow \tilde{c} - p$
8. Aller à 2.

Résultats

Pour aléatoire , avec la recherche de ligne pour la longueur de pas, la méthode converge très rapidement. Par exemple, les résultats suivants avec (100 ) sont typiques - la méthode converge quadratique. N = 100 v $v_i$$N=100$$v_i$

>> N = 100;
>> V = randn (N, N);
>> pour k = 1: NV (:, k) = V (:, k) / norme (V (:, k)); fin
>> maxEntropyMatrix (V);
Itération de Newton = 1, norme (grad f) = 0,67748
Itération de Newton = 2, norme (grad f) = 0,03644
Itération de Newton = 3, norme (grad f) = 0,0012167
Itération de Newton = 4, norme (grad f) = 1,3239e-06
Itération de Newton = 5, norme (grad f) = 7.7114e-13

Pour voir que le point optimal calculé est en fait le maximum, voici un graphique de la façon dont l'entropie change lorsque le point optimal est perturbé de manière aléatoire. Toutes les perturbations font diminuer l'entropie.

Code Matlab

Fonction tout en 1 pour minimiser l'entropie (nouvellement ajouté à ce message): https://github.com/NickAlger/various_scripts/blob/master/maxEntropyMatrix.m