Cálculo do determinante ao resolver

11

Estou resolvendo para uma enorme matriz definida positiva esparsa A usando o método do gradiente conjugado (CG). É possível calcular o determinante de A usando as informações produzidas durante a resolução?Ax=bAA

Manuel Schmidt
fonte
Por que você deseja calcular o determinante? Esse resultado certamente será um underflow ou um estouro para uma matriz enorme de qualquer maneira. Seria mais caridoso se você pedisse para calcular o número da condição, mas não perca seu tempo com o determinante!
Você provavelmente já sabe disso, mas os valores de Ritz durante o processo de gradiente conjugado convergem para os valores próprios da matriz e é possível derivar estimativas simples para o determinante.
shuhalo

Respostas:

10

Computar o determinante de uma matriz esparsa é tipicamente tão caro quanto uma solução direta, e sou cético quanto ao fato de que o CG seria de grande ajuda para computá-lo. Seria possível executar o CG para iterações (onde A é n × n ) para gerar informações para todo o espectro de A e depois calcular o determinante como o produto dos autovalores, mas isso seria lento e numericamente instável.nAn×nA

Seria uma idéia melhor calcular a fatoração de Cholesky esparsa direta de sua matriz, digamos , onde L é triangular inferior. Então det ( A ) = det ( L ) det ( L H ) = | det ( L ) | 2 , em que det ( L ) é simplesmente o produto das entradas diagonais da matriz triangular inferior L, uma vez que os valores próprios de uma matriz triangular estão ao longo de sua diagonal.A=LLHL

det(A)=det(L)det(LH)=|det(L)|2,
det(L)L

No caso de uma matriz geral não singular, uma decomposição LU pivotada deve ser usada, digamos, , onde P é uma matriz de permutação, de modo que det ( A ) = det ( P - 1 ) det ( L ) det ( U ) . Como P é uma matriz de permutação, det ( P ) = ± 1 e, por construção, LPA=LUP

det(A)=det(P1)det(L)det(U).
Pdet(P)=±1Lnormalmente terá uma diagonal de todos, o que implica que . Assim, você pode calcular det ( A ) como ± det ( U ) e novamente reconhecer que o determinante de uma matriz triangular é simplesmente o produto de suas entradas diagonais. Assim, o custo da computação do determinante é essencialmente apenas o de uma fatoração.det(L)=1det(A)±det(U)
Jack Poulson
fonte
A106x106
@ManuelSchmidt Matrizes esparsas desse tamanho, resultantes de discretizações do tipo de elementos finitos, geralmente podem ser facilmente fatoradas com (por exemplo) métodos multifrontais. Concordo que a fatoração de Cholesky deve ser usada se sua matriz for HPD (e a generalização do meu argumento acima é óbvia).
Jack Poulson
Obrigado pela sua resposta e resposta rápidas. Infelizmente, a matriz não possui estrutura espezial (o que permitiria uma fatoração fácil).
Manuel Schmidt
2
Estou curioso para saber por que você precisa calcular o determinante da matriz. Os valores próprios mais altos e mais baixos não são suficientes?
precisa
Faz parte de uma função complexa de distribuição de probabilidade e não apenas uma constante de normalização. Eu sei que as distribuições podem ser fatoradas (e é isso que estamos fazendo no momento), no entanto, temos toneladas de dados para modelar e cada um dos fatores se torna enorme.
Manuel Schmidt
6

ABdimAdimBdimB=

BABABdetAdetBAB

detA=j=1dimAλi(A)j=1dimAλi(B)j=dimA+1dimBλi(B)
BAdimB=detAdetB
Wolfgang Bangerth
fonte
Acontece que existem alguns algoritmos verdadeiramente bonitos e práticos que envolvem o cálculo de determinantes consideráveis. Verifique www-m3.ma.tum.de/foswiki/pub/M3/Allgemeines/…
Matt Knepley 15/14
2

Sem entrar (novamente) no porquê e como os determinantes são maus, vamos supor que seu operador não seja facilmente fatorável ou simplesmente não esteja disponível como uma matriz e que você realmente precise estimar seu determinante.

AA

Provavelmente, você pode fazer engenharia reversa de como essa estimativa do determinante ocorre na implementação padrão do CG, seguindo atentamente a Seção 6.7.3 do livro.

Dominique
fonte
2

det(A)=i=1nαk1,
αk=rkTrkpkTApkrk0k=1,,nRrkPpk
pk=rk+i=1k1γiri.
det(P)=(1)ndet(R)rkpkA
k=1nαk=k=1nrkTrkpkTApk=det(RTR)det(PTAP)=det(RTR)det(A)det(PTP)=(det(A))1.
Yermek
fonte