Mejor manera de resolver exactamente un sistema lineal (300 x 300) con coeficientes enteros

Question

Quiero resolver el sistema de ecuaciones lineales de tamaño 300 x 300 (300 variables y 300 ecuaciones) exactamente (no con punto flotante, aka dgesl, pero con fracciones en la respuesta). Todos los coeficientes de este sistema son entero (por ejemplo, 32 bits), sólo hay 1 solución a la misma. Hay no-cero términos constantes en la columna de la derecha (b).

A*x = b

• donde A es la matriz de coeficientes, b es el vector de términos constantes.
• la respuesta es x vector, dado que en los números racionales (fracciones de pares de muy larga de números enteros).

La matriz a es denso (caso general), pero puede tener hasta un 60 % de cero de los coeficientes.

¿Cuál es la mejor manera (algoritmo más rápido para x86/x86_64) para resolver este tipo de sistema?

Gracias!

PS típica respuesta de dichos sistemas enteros en la fracción hasta 50-80 dígitos, así que, por favor, no sugieren nada basa sólo en float/double. Ellos no tienen la necesaria precisión.

Answer 1

No sé si es el mejor, pero funciona linalg.solve de Scipy . Usted podría mirar el código fuente...

Editar para los casos que traté de scipy.linalg.solve(A,b) * scipy.linalg.det(A) da un resultado exacto, por ejemplo, da una matriz correspondiente a (1/det(A)) * x , por lo que el resultado fraccionario real es cada elemento dividida por ese denominador común.

Answer 2

Esto estaba destinado a ser el comentario -- en Mathematica me tarda 0,5 segundos para 300x300 matriz con alrededor del 60% de los ceros y el 40% de los enteros positivos con muestreo uniforme en una escala de 1-10.

Para integrar Mathematica biblioteca con C, es necesario utilizar MathLink de la interfaz. No estoy seguro de si se puede enlazar estáticamente, se puede requerir el lanzamiento de Mathematica kernel junto con su programa, con el fin de utilizar el interfaz.

La función relevante es LinearSolve. Aquí está el código completo que he usado para ejecutar sincronización de comparación

mat1 = RandomVariate[BernoulliDistribution[.4], {300, 300}];
mat2 = RandomInteger[{1, 10}, {300, 300}];
mat3 = mat1*mat2;
answer = RandomInteger[{1, 10}, 300];
result = LinearSolve[mat3, answer]

No sé el algoritmo que se utiliza. Dicen que la versión 8 es aproximadamente 1000 veces más rápido que la versión anterior, así que es probable que haya algo especializado y de propiedad

Answer 3

Este tipo de cosas se puede hacer por paquetes de cálculo simbólico Maple, Mathematica, Maxima, etcetera.

Si usted necesita una subrutina para hacer esto que te llaman de un programa más grande, entonces la respuesta va a depender mucho el lenguaje de programación que está usando y qué tipo de licencia te está dispuesto a trabajar con.

Answer 4

El Dixon del algoritmo. Se hace una inversión de la matriz en la aritmética modular, basada en el primer módulo P (este paso es muy sencillo y exacto). A continuación, se resuelve el trivial del sistema para cada x para obtener un p-ádico número de la representación de la solución. Puede ser convertido a un número racional. Método se describe en "Dixon, J. D. solución Exacta de las ecuaciones lineales usando P-ádico expansiones // Numer. Math., 40, 1982, P. 137-141.", pero este artículo está fuera de línea.

Hay un artículo acerca de Dixon del algoritmo y sus modificaciones en línea: http://www2.isye.gatech.edu/~dsteffy/papers/OSLifting.pdf

El autor también tiene todos describen los algoritmos implementados y publicado en http://www2.isye.gatech.edu/~dsteffy/racional/

Otro enfoque interesante es http://en.wikipedia.org/wiki/Montante's_method>Montante del método, que es una variante de gauss para un entero coeficiente de matrices. Esta variante no crear números racionales en la matriz, mientras que la conversión a triangular. Echa un vistazo a la última foto en el link http://en.wikipedia.org/wiki/Montante's_method

Answer 5

Para el algoritmo, sencillo eliminación Gaussiana debe estar bien. Esta parte puede ser en realidad más simple que cuando se trata con números de punto flotante, ya que usted no tiene que preocuparse acerca de la estabilidad numérica.

Dependiendo de su matriz, usted podría ser capaz de hacer un poco mejor que la eliminación Gaussiana, por ejemplo, si la matriz es simétrica, se puede utilizar una factorización de Cholesky. Pero el 60% de escasez no es mucho en el esquema de las cosas. Si se tratara de tridiagonal, bandas, etc. entonces usted puede intentar algunos métodos especializados.

Usted necesitará un buen número racional de la biblioteca para manejar la real de las operaciones aritméticas. GNU MP debe ajustarse a la ley, y las reclamaciones se pueden hacer números racionales, pero nunca la he utilizado.