Resolviendo un sistema escaso y altamente mal acondicionado

Tengo la intención de resolver Ax = b donde A es matriz cuadrada o rectangular compleja, escasa, asimétrica y altamente mal acondicionada (número de condición ~ 1E + 20). He podido resolver el sistema con ZGELSS en LAPACK con precisión. Pero a medida que crecen los grados de libertad en mi sistema, lleva mucho tiempo resolver el sistema en una PC con ZGELSS ya que no se explota la escasez. Recientemente probé SuperLU (usando el almacenamiento Harwell-Boeing) para el mismo sistema, pero los resultados fueron inexactos para el número de condición> 1E + 12 (no estoy seguro de si este es un problema numérico con el pivote).

Estoy más inclinado a usar solucionadores ya desarrollados. ¿Existe un solucionador robusto que pueda resolver el sistema que mencioné rápidamente (es decir, explotar la escasez) y de manera confiable (en vista de los números de condición)?

Cuando usa ZGELSS para resolver este problema, está usando la descomposición de valores singulares truncados para regularizar este problema extremadamente mal condicionado. Es importante comprender que esta rutina de biblioteca no intenta encontrar una solución de mínimos cuadrados para , sino que intenta equilibrar la búsqueda de una solución que minimicecontra minimizar. $Ax=b$$\| x \|$$\| Ax-b \|$

Tenga en cuenta que el parámetro RCOND pasado a ZGELSS puede usarse para especificar qué valores singulares deben incluirse y excluirse del cálculo de la solución. Cualquier valor singular menor que RCOND * S (1) (S (1) es el valor singular más grande) será ignorado. No nos ha dicho cómo ha configurado el parámetro RCOND en ZGELSS, y no sabemos nada sobre el nivel de ruido de los coeficientes en su matriz o en el lado derecho , por lo que es difícil decir si ha utilizado Una cantidad adecuada de regularización. $A$$b$

Parece estar contento con las soluciones regularizadas que está obteniendo con ZGELSS, por lo que parece que la regularización efectuada por el método SVD truncado (que encuentra una solución mínima entre las soluciones de mínimos cuadrados que minimizan sobre el espacio de soluciones abarcado por los vectores singulares asociados con los valores singulares mayores que RCOND * S (1)) es satisfactorio para usted. $\| x \|$$\| Ax-b \|$

Su pregunta podría reformularse como "¿Cómo puedo obtener eficientemente soluciones de mínimos cuadrados regularizadas para este gran problema de mínimos cuadrados lineales muy dispersos y muy mal condicionados?"

Mi recomendación sería utilizar un método iterativo (como CGLS o LSQR) para minimizar el problema de mínimos cuadrados explícitamente regularizado

$\min \| Ax-b \|^{2} + \alpha^{2} \| x \|^{2}$

donde el parámetro de regularización se ajusta para que el problema de mínimos cuadrados amortiguados esté bien condicionado y para que esté satisfecho con las soluciones regularizadas resultantes. $\alpha$

Jed Brown ya ha señalado esto en los comentarios a la pregunta, pero en realidad no hay mucho que pueda hacer con la doble precisión habitual si su número de condición es grande: en la mayoría de los casos, es probable que no obtenga un solo dígito de precisión en su solución y, lo que es peor, ni siquiera puede saberlo porque no puede evaluar con precisión el residuo correspondiente a su vector de solución. En otras palabras: no se trata de qué solucionador lineal debe elegir, ningún solucionador lineal puede hacer algo útil para tales matrices.

$1,x,x^2,x^3,...$

La forma más simple / rápida de resolver problemas mal condicionados es aumentar la precisión de los cálculos (por fuerza bruta). Otra forma (aunque no siempre posible) es reformular su problema.

Es posible que deba usar precisión cuádruple (34 dígitos decimales). Aunque se perderán 20 dígitos en un curso (debido al número de condición), obtendrá 14 dígitos correctos.

Si le interesa, ahora los solucionadores dispersos de precisión cuádruple también están disponibles en MATLAB.

(Soy el autor de la caja de herramientas mencionada).