Algorithme d'équilibrage de matrice

J'ai écrit une boîte à outils de système de contrôle à partir de zéro et uniquement en Python3 (plug sans vergogne:) harold. De mes recherches passées, je me suis toujours plaint du solveur Riccati care.mpour des raisons techniques / non pertinentes.

Par conséquent, j'ai écrit mon propre ensemble de routines. Une chose que je ne peux pas trouver est d'obtenir un algorithme d'équilibrage hautes performances, au moins aussi bon que balance.m. Avant de le mentionner, la xGEBALfamille est exposée dans Scipy et vous pouvez essentiellement appeler à partir de Scipy comme suit, supposons que vous ayez un tableau 2D de type flottant A:

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A, scale=1 , permute=1 , overwrite_a=0 )

Maintenant, si j'utilise la matrice de test suivante

array([[ 6.      ,  0.      ,  0.      ,  0.      ,  0.000002],
       [ 0.      ,  8.      ,  0.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  6.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  0.      ,  8.      ,  0.      ],
       [ 0.      ,  0.      ,  0.000002,  0.      ,  2.      ]])

Je reçois

array([[ 8.      ,  0.      ,  0.      ,  2.      ,  2.      ],
       [ 0.      ,  2.      ,  0.000002,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  6.      ,  2.      ,  2.      ],
       [ 0.      ,  0.000002,  0.      ,  6.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  8.      ]])

Cependant, si je passe cela à balance.m, je reçois

>> balance(A)

ans =

    8.0000         0         0    0.0625    2.0000
         0    2.0000    0.0001         0         0
         0         0    6.0000    0.0002    0.0078
         0    0.0003         0    6.0000         0
         0         0         0         0    8.0000

Si vous vérifiez les modèles de permutation, ils sont les mêmes, mais l'échelle est désactivée. gebaldonne l' unité battitures alors que Matlab donne les pouvoirs suivants 2: [-5,0,8,0,2].

Donc apparemment, ceux-ci n'utilisent pas les mêmes machines. J'ai essayé diverses options telles que Lemonnier, la mise à l'échelle bilatérale Van Dooren, la Parlett-Reinsch originale et également d'autres méthodes moins connues dans la littérature comme la version dense de SPBALANCE.

Je voudrais peut-être souligner un point: je connais le travail de Benner; en particulier l' équilibrage symplectique des matrices hamiltoniennes spécifiquement à cet effet. Cependant, notez que ce type de traitement se fait dans gcare.m(solveur Riccati généralisé) et l'équilibrage se fait directement via balance.m. Par conséquent, j'apprécierais que quelqu'un me montre la mise en œuvre réelle.

Divulgation: Je n'essaie vraiment pas de rétroconcevoir le code mathworks: je veux en fait m'en éloigner pour diverses raisons, notamment la motivation de cette question, c'est-à-dire que je ne sais pas ce qu'il fait qui m'a coûté beaucoup de temps dans la journée. Mon intention est d'obtenir un algorithme d'équilibrage satisfaisant qui me permette de passer des exemples CAREX de telle sorte que je puisse implémenter des méthodes d'itération Newton au-dessus du solveur standard.

Cela m'a pris un certain temps pour comprendre cela et comme d'habitude, cela devient évident après avoir trouvé le coupable.

Après avoir vérifié les cas problématiques signalés dans David S. Watkins. Un cas où l'équilibre est nuisible. Électron. Trans. Numer. Anal, 23: 1–4, 2006 et aussi la discussion ici (tous deux cités dans arXiv: 1401.5766v1 ), il s'avère que matlab utilise l'équilibrage en séparant d'abord les éléments diagonaux.

Ma pensée initiale était, selon la documentation limitée classique sur les fonctions LAPACK, GEBAL a effectué cela automatiquement. Cependant, je suppose que ce que les auteurs veulent dire en ignorant les éléments diagonaux ne les supprime pas des sommes de ligne / colonne.

En fait, si je supprime manuellement la diagonale du tableau, les deux résultats coïncident, c'est-à-dire

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A - np.diag(np.diag(A)), scale=1 , permute=1 , overwrite_a=0 )  

donne le même résultat que balance.m(sans les entrées diagonales bien sûr).

Si un utilisateur Fortran-savy peut le confirmer en vérifiant dgebal.f , je vous en serais reconnaissant.

EDIT: Le résultat ci-dessus ne signifie pas que c'est la seule différence. J'ai également construit différentes matrices où GEBAL et balance.m produisent des résultats différents même après la séparation des diagonales.

Je suis assez curieux de savoir quelle pourrait être la différence, mais il semble qu'il n'y ait aucun moyen de le savoir car il s'agit d'un code matlab intégré et donc fermé.

EDIT2 : Il s'avère que matlab utilisait une ancienne version de LAPACK (probablement avant 3.5.0) et en 2016b, ils semblent avoir été mis à niveau vers la nouvelle version. Maintenant, les résultats sont d'accord pour autant que je puisse tester. Je pense donc que cela règle le problème. J'aurais dû le tester avec d'anciennes versions de LAPACK.