Solve * escaso * sistema triangular superior

Question

Si quiero resolver un sistema triangular superior completo, puedo llamar al linsolve(A,b,'UT'). Sin embargo, esto actualmente no es compatible con matrices dispersas. ¿Cómo puedo superar esto?

Answer 1

Editar Dado que lo que necesita es un procedimiento resolver triangular, también llamada hacia atrás/hacia delante sustitución, puede utilizar MATLAB ordinaria barra invertida \ operador para que:

x = U\b 

Como se mencionó en la respuesta original, MATLAB reconocerá el hecho de que su matriz es triangular. Para estar seguro de eso, puede comparar el rendimiento con el procedimiento cs_usolve que se encuentra en SuiteSparse. Es una función mex implementada en C que calcula la resolución triangular dispersa para la matriz escasa triangular superior (también hay funciones similares: cs_lsolve, cs_utsolve y cs_ltsolve).

Puede echar un vistazo a performance comparison de MATLAB nativo y cs_l(t)solve en el contexto de la factorización de Cholesky dispersa. Esencialmente, el rendimiento de MATLAB es bueno. El único escollo es que si quieres resolver un sistema de transposición

x = U'\b 

MATLAB no no reconocen que y crea explícitamente una transpuesta de U. En ese caso, debe llamar al cs_utsolve explícitamente.

Respuesta original Si el sistema es simétrico y sólo se almacena la parte de la matriz triangular superior (así es como entendí completa en su pregunta), y si la descomposición de Cholesky es adecuado para usted, chol maneja matrices simétricas, si su matriz es positiva definida. Para matrices indefinidas, puede usar ldl. Ambos manejan el almacenamiento escaso y trabajan en las partes simétricas de la matriz.

Las versiones de matlab más nuevas usan cholmod and suitesparse para eso. Esa es de lejos la factorización Cholesky de mejor rendimiento que conozco. En matlab también se paraleliza usándonos BALS paralelos.

El factor que obtenga de las funciones anteriores es matriz triangular superior L tal que

A=LL' 

único que tiene que hacer ahora es realizar el avance y la sustitución hacia atrás, que es simple y barato. En Matlab esto se realiza automáticamente en el operador tha backslash

x=L'\(L\b) 

la matriz puede ser escasa, y MATLAB reconocerá que es superior/triangular inferior. También usaría esta llamada junto con la sustitución directa de los factores obtenidos mediante la factorización de cholesky.

Answer 2

puede utilizar MLDIVIDE (\) o MRDIVIDE (/) los operadores en sus matrices dispersas ...

Answer 3

Los sistemas UT y LT se encuentran entre los sistemas más fáciles de resolver. Lea on the wiki al respecto. Sabiendo esto, es fácil para escribir su propio UT o LT solucionador:

%# some example data 
A = sparse(triu(rand(100))); 
b = rand(100,1); 

%# solve UT system by back substitution  
x = zeros(size(b)); 
for n = size(A,1):-1:1  
    x(n) = (b(n) - A(n,n+1:end)*x(n+1:end))/A(n,n);  
end 

El procedimiento es bastante similar para los sistemas de LT.

Dicho esto, en general es mucho más fácil y rápido de usar el operador barra invertida de Matlab:

x = A\b 

que también trabaja para los repuestos matrices, como ya se ha indicado Nate.

Tenga en cuenta que este operador también resuelve los sistemas UT que tienen A no cuadrado o si A tiene algunos elementos iguales a cero (o < eps) en la diagonal principal. Resuelve estos casos en un sentido de mínimos cuadrados, lo que puede o no ser deseable para ti. Usted puede comprobar para estos casos antes de llevar a cabo la solución:

if size(A,1)==size(A,2) && all(abs(diag(A)) > eps) 
    x = A\b; 
else 
    %# error, warning, whatever you want 
end 

Lea más sobre la (posterior) slash operador escribiendo

>> help \ 

o

>> help slash 

en la línea de comandos de Matlab .