Algoritmo de equilibrio de matriz

He estado escribiendo una caja de herramientas del sistema de control desde cero y puramente en Python3 (plug descarado:) harold. De mi investigación anterior, siempre tengo quejas sobre el solucionador de Riccati care.mpor razones técnicas / irrelevantes.

Por lo tanto, he estado escribiendo mi propio conjunto de rutinas. Una cosa que no puedo encontrar es evitar un algoritmo de equilibrio de alto rendimiento, al menos tan bueno como balance.m. Antes de mencionarlo, la xGEBALfamilia está expuesta en Scipy y básicamente puede llamar desde Scipy de la siguiente manera, supongamos que tiene una matriz 2D de tipo flotante 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 )

Ahora si uso la siguiente matriz de prueba

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.      ]])

yo obtengo

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.      ]])

Sin embargo, si le paso esto balance.m, me sale

>> 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 marca los patrones de permutación, son los mismos, sin embargo, la escala está desactivada. gebalda escaladas unidad mientras que MATLAB proporciona las siguientes potencias de 2: [-5,0,8,0,2].

Aparentemente, estos no están usando la misma maquinaria. He probado varias opciones, como Lemonnier, Van Dooren escala a doble cara, Parlett-Reinsch original y también algunos otros métodos menos conocidos en la literatura, como la versión densa SPBALANCE.

Un punto que quizás podría enfatizar es que conozco el trabajo de Benner; en particular, el equilibrio simpléctico de matrices hamiltonianas específicamente para este propósito. Sin embargo, tenga en cuenta que este tipo de tratamiento se realiza dentro de gcare.m(solucionador generalizado de Riccati) y el equilibrio se realiza directamente a través de balance.m. Por lo tanto, agradecería si alguien me puede señalar la implementación real.

Divulgación: Realmente no estoy tratando de revertir el código de ingeniería matemática: realmente quiero alejarme de él debido a varias razones, incluida la motivación de esta pregunta, es decir, no sé qué está haciendo lo que me costó mucho de tiempo atrás en el día. Mi intención es obtener un algoritmo de equilibrio satisfactorio que me permita pasar ejemplos de CAREX de modo que pueda implementar métodos de iteración de Newton sobre el solucionador habitual.

Me tomó bastante tiempo resolver esto y, como de costumbre, se vuelve obvio después de encontrar al culpable.

Después de verificar los casos problemáticos reportados en David S. Watkins. Un caso donde el equilibrio es dañino. Electrón. Trans. Numer. Anal, 23: 1–4, 2006 y también la discusión aquí (ambas citadas en arXiv: 1401.5766v1 ), resulta que matlab usa el equilibrio separando primero los elementos diagonales.

Mi pensamiento inicial fue, según la documentación limitada clásica sobre las funciones de LAPACK, GEBAL realizó esto automáticamente. Sin embargo, supongo que lo que los autores quieren decir al ignorar los elementos diagonales no es eliminarlos de las sumas de fila / columna.

De hecho, si elimino manualmente la diagonal de la matriz, ambos resultados coinciden, es decir

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 )  

da el mismo resultado que balance.m(sin las entradas diagonales, por supuesto).

Si algún usuario de Fortran-savy puede confirmar esto marcando dgebal.f , estaría agradecido.

EDITAR: El resultado anterior no implica que esta sea la única diferencia. También construí diferentes matrices donde GEBAL y balance.m producen resultados diferentes incluso después de separar las diagonales.

Tengo curiosidad por saber cuál puede ser la diferencia, pero parece que no hay forma de saberlo, ya que es un código incorporado de Matlab y, por lo tanto, cerrado.

EDIT2 : Resulta que matlab estaba usando una versión anterior de LAPACK (probablemente anterior a 3.5.0) y para 2016b parecen haberse actualizado a la versión más nueva. Ahora los resultados están de acuerdo en lo que puedo probar. Así que creo que eso resuelve el problema. Debería haberlo probado con versiones anteriores de LAPACK.