¿Se puede aproximar la solución de un sistema lineal de ecuaciones solo para las primeras variables?

Tengo un sistema lineal de ecuaciones de tamaño mxm, donde m es grande. Sin embargo, las variables que me interesan son solo las primeras n variables (n es pequeña en comparación con m). ¿Hay alguna manera de aproximar la solución para los primeros valores de m sin tener que resolver todo el sistema? Si es así, ¿sería esta aproximación más rápida que resolver el sistema lineal completo?

Como otros han señalado, esto es difícil de hacer con un solucionador directo. Dicho esto, no es tan difícil de hacer con solucionadores iterativos. Con este fin, tenga en cuenta que la mayoría de los solucionadores iterativos de una forma u otra minimizan el error con respecto a alguna norma. A menudo, esta norma es inducida por la matriz misma, pero a veces también es solo la norma del vector l2. Pero ese no tiene por qué ser el caso: puede elegir en qué norma quiere minimizar el error (o residual) y, por ejemplo, puede elegir una norma en la que pese los componentes que le interesan con 1 y todos los demás con 1e-12, es decir, por ejemplo, algo como (1e-24) ∑ N i = 6 x 2 i y el producto escalar correspondiente. Luego, escriba todos los pasos del solucionador iterativo con respecto a esta norma y producto escalar, y obtendrá un solucionador iterativo que presta mucha más atención a los elementos vectoriales que le interesan que a los demás.$|| x ||^2 = \sum_{i=1}^5 x_i^2 +$$\sum_{i=6}^N x_i^2$

La pregunta, por supuesto, es si necesita menos iteraciones que con el producto estándar / escalar que pesa todos los componentes por igual. Pero ese debería ser el caso: digamos que solo te importan los cinco primeros elementos vectoriales. Entonces debería necesitar como máximo cinco iteraciones para reducir el error en un factor de 1e12 ya que cinco iteraciones es lo que se necesita para el sistema 5x5 que las describe. Eso no es una prueba, pero estoy bastante seguro de que debería salirse con una cantidad mucho menor de iteraciones si el peso en la norma (1e-12 arriba) es menor que la tolerancia con la que desea resolver el sistema lineal de forma iterativa .

Formando el complemento Schur

Suponga que ha permutado y dividido su matriz en el formulario

tal que contenga sus grados de libertad de interés y sea mucho más pequeño que A 11 , entonces uno puede formar el complemento Schur$A_{22}$$A_{11}$

ya sea a través de una factorización LU parcial de aspecto correcto o la fórmula explícita, y luego puede entenderse en el siguiente sentido:$S_{22}$

donde representa la parte 'poco interesante' de la solución. Por lo tanto, dado un lado derecho que no es cero en los grados de libertad del complemento Schur S 22 , solo necesitamos resolver contra S 22 para obtener la porción de la solución correspondiente a esos grados de libertad.$\star$$S_{22}$$S_{22}$

Complejidad computacional en caso denso no estructurado

Configuración de a la altura de A y n a la altura de A 22 , a continuación, el método estándar para el cálculo de S 22 es a primera factor de L 11 U 11 : = A 11 (Ignoremos pivotante por ahora) en aproximadamente 2 / 3 ( N - n ) 3 trabajos, luego formar$N$$A$$n$$A_{22}$$S_{22}$$L_{11} U_{11} := A_{11}$$2/3 (N-n)^3$

usando dos soluciones triangulares que requieren trabajo cada una, y luego realizando la actualización a A 22 en 2 n 2 ( N - n ) trabajo.$n(N-n)^2$$A_{22}$$2n^2 (N-n)$

Por lo tanto, el trabajo total es aproximadamente . Cuando n es muy pequeño, N - n ≈ N , por lo que el coste puede ser visto a ser aproximadamente 2 / 3 N 3 , que es el costo de una factorización completa.$2/3 (N-n)^3 + 2n(N-n)^2 + 2n^2 (N-n)$$n$$N-n \approx N$$2/3 N^3$

El beneficio es que, si hay una gran cantidad de lados derechos para resolver con el mismo sistema de ecuaciones, entonces podría potencialmente reutilizarse una gran cantidad de veces, donde cada resolución solo requeriría 2 n 2 de trabajo (en lugar de 2 N 2 funciona) si S 22 está factorizado.$S_{22}$$2n^2$$2N^2$$S_{22}$

Complejidad computacional en el caso disperso (típico)

Si su sistema disperso surgió de algún tipo de aproximación por diferencias finitas o elementos finitos, entonces los solucionadores directos dispersos seguramente podrán explotar parte de la estructura; Sistemas 2D se pueden resolver con de trabajo y O ( N log N ) de almacenamiento, mientras que los sistemas 3D se pueden resolver con O ( N 2 ) de trabajo y O ( N 4 / 3 ) de almacenamiento. Los sistemas factorizados se pueden resolver con la misma cantidad de trabajo que los requisitos de almacenamiento.$O(N^{3/2})$$O(N \log N)$$O(N^2)$$O(N^{4/3})$

El punto de plantear las complejidades computacionales es que, si y usted tiene un sistema 2d, entonces dado que el complemento de Schur probablemente será denso, la complejidad de resolución dado el complemento de Schur factorizado seráO(n2)=O(N), que solo le falta un factor logarítmico en lugar de resolver el completo ¡sistema! En 3D, que requiereO(N)trabajo en lugar deO(N 4 / 3 ).$n \approx \sqrt{N}$$O(n^2) = O(N)$$O(N)$$O(N^{4/3})$

Por lo tanto, es importante tener en cuenta que, en su caso, donde , solo habrá ahorros significativos si está trabajando en varias dimensiones y tiene muchos lados derechos para resolver.$n=\sqrt{N}$

El enfoque de reducción modelo

Como Paul preguntó, hablaré sobre lo que sucede si usa métodos de reducción de modelos basados ​​en proyección para este problema. Suponga que podría llegar a un proyector tal que el rango de P , denotado R ( P ) , contenga la solución a su sistema lineal A x = b , y tenga una dimensión k , donde k es el número de incógnitas para las cuales usted desea resolver en un sistema lineal.$\mathbf{P}$$\mathbf{P}$$\mathcal{R}(\mathbf{P})$$\mathbf{Ax} = \mathbf{b}$$k$$k$

Una descomposición de valor singular de producirá la siguiente matriz particionada:$\mathbf{P}$

Las matrices oscurecidas por las estrellas son importantes para otras cosas (como error de estimación, etc.), pero por ahora, evitaremos tratar con detalles extraños. Resulta que

es una descomposición de rango completo de .$\mathbf{P}$

Esencialmente, resolverás el sistema

de una manera inteligente, ya que y W también tienen la propiedad de que W T V = I . Multiplicando ambos lados de P A x = P b por W T y dejando que y = V x ser una aproximación para x rendimientos$\mathbf{V}$$\mathbf{W}$$\mathbf{W}^{T}\mathbf{V} = \mathbf{I}$$\mathbf{PAx} = \mathbf{Pb}$$\mathbf{W}^{T}$$\mathbf{y} = \mathbf{V}\widehat{\mathbf{x}}$$\mathbf{x}$

Resuelve para x , premultiplicar por V , y usted tiene y , a su aproximación para x .$\widehat{\mathbf{x}}$$\mathbf{V}$$\mathbf{y}$$\mathbf{x}$

¿Por qué el enfoque de complemento de Schur es probablemente mejor?

Para empezar, debes elegir alguna manera. Si la solución a A x = b está en R ( P ) , entonces y = x , e y no es una aproximación. De lo contrario, y ≠ x , e introduce algún error de aproximación. Este enfoque realmente no aprovecha toda la estructura que mencionó que desea explotar. Si elegimos P de modo que su rango sea la base de la unidad estándar en las coordenadas de x que desea calcular, las coordenadas correspondientes de y tendrán errores. No está claro cómo te gustaría elegir$\mathbf{P}$$\mathbf{Ax} = \mathbf{b}$$\mathcal{R}(\mathbf{P})$$\mathbf{y} = \mathbf{x}$$\mathbf{y}$$\mathbf{y} \neq \mathbf{x}$$\mathbf{P}$$\mathbf{x}$$\mathbf{y}$ . Podría utilizar una SVD de A , por ejemplo, y seleccionar P para que sea el producto de los primeros k vectores singulares izquierdos de A y el adjunto de los primeros k vectores singulares derechos de A , suponiendo que los vectores singulares estén dispuestos en orden decreciente de valor singular Esta elección de proyector sería equivalente a realizar una descomposición ortogonal adecuada en A , y minimizaría elerror deL 2 en la solución aproximada.$\mathbf{P}$$\mathbf{A}$$\mathbf{P}$$k$$\mathbf{A}$$k$$\mathbf{A}$$\mathbf{A}$$_{2}$

$\mathbf{V}$$\mathbf{W}$$\mathbf{P}$

Los inconvenientes son muy parecidos al enfoque de JackPoulson, excepto que no estás aprovechando la estructura que mencionaste.

La respuesta larga es ... más o menos.

$k$

$k$

$n-k$$n$

Además, tenga en cuenta que la restricción de la orden en el que se va a realizar la copia de substituion puede restringir la forma de la matriz (que le quita la capacidad de las columnas de intercambio) que podría posiblemente conducir a un sistema mal condicionado, pero no estoy seguro de eso, solo algo a tener en cuenta.