Testando se duas matrizes 12x12 têm o mesmo determinante

$12 \times 12$ $Q$

det (Q) = det (12 I - Q - J) (1)

$\det(Q) = \det(12I-Q-J) \; \; (1)$

J

$J$

Atualmente, estou fazendo isso com a biblioteca do tatu , mas ela acaba sendo muito lenta. O problema é que eu preciso fazer isso por um trilhão de matrizes e acontece que calcular os dois determinantes é o gargalo do meu programa. Portanto, eu tenho duas perguntas

Existe algum truque que eu possa usar para calcular o determinante mais rapidamente, desde que eu saiba o tamanho deles? Talvez uma expansão confusa para $12 \times12$ matrizes que poderiam funcionar nesse caso?
Existe alguma outra maneira eficiente de testar a igualdade $(1)$

Editar. Para responder aos comentários. Eu preciso calcular todos os gráficos não auto-complementares conectados $G$ da ordem $13$ , modo que $G$ e $\overline{G}$ tenham o mesmo número de árvores de abrangência. A motivação para isso pode ser encontrada nesta postagem do mathoverflow . Quanto à máquina, eu estou executando isso em uma máquina de 3,4 GHh de 8 núcleos em paralelo.

Editar. Consegui reduzir o tempo de execução esperado em 50% criando um programa C para calcular especificamente o determinante de uma matriz $12 \times 12$ . Sugestões ainda são bem-vindas.

linear-algebra matrix c++ graph-theory c Jernej
fonte

Quão lento é muito lento? Quanto tempo leva em que hardware? Os trilhões desses são independentes para que você possa calcular muitos desses determinantes em paralelo? Se sim, qual o tamanho de uma máquina que você pode executar? O que levou a esse problema? Tem certeza de que precisa calcular os determinantes?

Q

$Q$

Bill Barth

Com que frequência (para qual fração dos casos) os determinantes são iguais / diferentes? Se eles são diferentes na maioria das vezes, pode haver um teste mais barato para determinar se eles podem ser diferentes e você verificaria se eles são os mesmos somente se o primeiro teste falhar. O contrário, se eles são os mesmos na maioria das vezes.

Wolfgang Bangerth

Como já perguntado: você poderia fornecer alguns detalhes sobre de onde veio? Talvez exista uma abordagem melhor do que os determinantes da computação às cegas.

Q

$\mathbf Q$

A noção de que essa condição deve ser testada "para um trilhão de matrizes" sugere 1) que é conhecido a priori por ter alguma estrutura especial (caso contrário, a expectativa de que a condição seja aleatória é pequena) e 2) que uma abordagem melhor pode ser para caracterizar todas as matrizes com essa propriedade (com uma formulação eficientemente verificável).

Q

$Q$

Q

$Q$

hardmath

@hardmath Sim, é uma matriz inteira com entradas diagonais que variam de a e

Q

$Q$

1

$1$

12

$12$

- 1

$-1$ como fora elementos diagonais

Jernej

Respostas:

Como você já está usando C ++ e suas matrizes são definidas positivamente simétricas, eu realizaria um processo não dinâmico $LDL^T$ fatoração dinâmica de e também de . Aqui, estou assumindo que o também é definido positivamente, caso contrário, o exigirá um pivô para estabilidade numérica (também é possível que, embora não seja definido positivamente, o pivotamento não seja necessário, mas você deve experimentá-lo). $Q$ $12I-Q-J$ $12I-Q-J$ $LDL^T$

Isso é mais rápido que uma fatoração de LU e também mais rápido que Cholesky, porque raízes quadradas são evitadas. O determinante é simplesmente o produto dos elementos da matriz diagonal . O código para executar uma fatoração LDL é tão simples que você pode escrevê-lo em menos de 50 linhas de C. A página da Wikipedia descreve o algoritmo, e eu tenho um código de modelo simples para executar Cholesky aqui . Você pode simplificá-lo bastante e modificá-lo para evitar a raiz quadrada para implementar a fatoração $D$ $LDL^T$

Como você também pode controlar o formato de armazenamento, é possível otimizar ainda mais a rotina para armazenar apenas metade da matriz e agrupá-la em uma matriz linear para maximizar a localidade da memória. Eu também escreveria rotinas simples de produtos personalizados e atualizações de classificação 1, já que os tamanhos dos problemas são tão pequenos que você deve deixar o compilador alinhar as rotinas para reduzir a sobrecarga de chamadas. Como é um loop de tamanho fixo, o compilador deve ser capaz de alinhar e desenrolar automaticamente as coisas quando apropriado.

Eu evitaria tentar fazer truques para aproveitar o fato de que $12I-Q-J$ conter dentro da expressão. É provável que, para tamanhos de problemas tão pequenos, esses truques acabem sendo mais lentos do que simplesmente realizar dois cálculos determinantes separados. Obviamente, a única maneira de verificar essas alegações é tentar. $Q$

Victor Liu
fonte

Eu apóio a recomendação de implementar o , também conhecido como Cholesky sem raiz, já que o Armadillo não parece ter uma maneira de tirar vantagem da definição positiva / dominância diagonal.

L D L^{T}

$LDL^T$

hardmath

Sem algumas informações sobre a construção dessas matrizes simétricas reais definidas positivas, as sugestões a serem feitas são necessariamente limitadas. $12\times 12$

Fiz o download do pacote Armadillo do Sourceforge e dei uma olhada na documentação. Tente melhorar o desempenho da computação separada e , onde é a matriz de classificação um de todos, definindo por exemplo . A documentação observa que esse é o padrão para matrizes de tamanho $\det(Q)$ $\det(12I - Q - J)$ $J$ det(Q,slow=false) $4\times 4$ , portanto, por omissão, presumo que a slow=trueopção seja um padrão para o caso . $12\times 12$

O que slow=true presumivelmente faz é uma rotação parcial ou total na obtenção de um formulário de escalão de linha, a partir do qual o determinante é facilmente encontrado. No entanto, você sabe com antecedência que a matriz é definitiva positiva, portanto, o giro é desnecessário para a estabilidade (pelo menos presumivelmente para a maior parte de seus cálculos. Não está claro se o pacote Armadillo lança uma exceção se os pivôs se tornarem indevidamente pequenos, mas isso deve ser um recurso de um pacote de álgebra linear numérica razoável.Edit : Encontrei o código do Tatu que implementa no arquivo de cabeçalho , usando modelos C ++ para funcionalidade substancial. A configuração não parece afetar como um $Q$ detinclude\armadillo_bits\auxlib_meat.hppslow=false $12\times 12$ o determinante será feito porque o cálculo é "jogado por cima de um muro" para o LAPACK (ou ATLAS) nesse ponto, sem indicação de que não é necessário girar; veja det_lapacke suas invocações nesse arquivo.

O outro ponto seria seguir a recomendação de criar o pacote Armadillo vinculado a substituições de alta velocidade para BLAS e LAPACK, se você realmente estiver usando; veja a seção 5 do arquivo README.TXT do Armadillo para obter detalhes. [O uso de uma versão dedicada de 64 bits do BLAS ou LAPACK também é recomendado para obter velocidade nas máquinas atuais de 64 bits.]

A redução de linha na forma de escalão é essencialmente uma eliminação gaussiana e tem complexidade aritmética . Para ambas as matrizes, isso equivale ao dobro do trabalho, ou . Essas operações podem muito bem ser o "gargalo" em seu processamento, mas há pouca esperança de que, sem uma estrutura especial em $\frac{2}{3} n^3 + O(n^2)$ $\frac{4}{3} n^3 + O(n^2)$ $Q$ (ou algumas relações conhecidas entre os trilhões de casos de teste que permitem a amortização), o trabalho possa ser reduzido para . $O(n^2)$

Para comparação, a expansão por co-fatores de uma matriz envolveoperações de multiplicação (e aproximadamente tantas adições / subtrações); portanto, para a comparação ( vs. ) favorece claramente a eliminação dos cofatores. $n\times n$ $n!$ $n=12$ $12! = 479001600$ $\frac{2}{3} n^3 = 1152$

Outra abordagem que requer trabalho seria reduzir para a forma tridiagonal com transformações de Householder, que também coloca na forma tridiagonal. A computação e pode ser realizada posteriormente em operações. [O efeito da atualização do ranking um $\frac{4}{3} n^3 + O(n^2)$ $Q$ $12I - Q$ $\det(Q)$ $\det(12I - Q - J)$ $O(n)$ $-J$ no segundo determinante pode ser expresso como um fator escalar dado pela resolução de um sistema tridiagonal.]

A implementação de uma computação independente pode valer a pena como uma verificação dos resultados de chamadas bem-sucedidas (ou com falha) à detfunção do Armadillo .

Caso especial: Como sugerido por um Comentário de Jernej, suponha que onde como antes é a matriz (rank 1) de todos e seja um matriz diagonal não singular (positiva). De fato, para a aplicação proposta na teoria dos grafos, essas seriam matrizes inteiras. Em seguida, uma fórmula explícita para $Q = D - J$ $J$ $D=\text{diag}(d_1,\ldots,d_n)$ $\det(Q)$ é:

det (Q) = (\prod_{i = 1}^{n} d_{i}) (1 - \sum_{i = 1}^{n} d_{i}^{- 1})

$\det(Q) = \left(\prod_{i=1}^n d_i \right)\left(1 - \sum_{i=1}^n d_i^{-1} \right)$

Um esboço de sua prova oferece uma oportunidade para ilustrar uma aplicabilidade mais ampla, ou seja, sempre que tem um determinante conhecido e o sistema $D$ $Dv = (1\ldots 1)^T$ é rapidamente resolvido. Comece fatorando:

det (D - J) = det (D) \cdot det (I - D^{- 1} J)

$\det(D - J) = \det(D) \cdot \det(I - D^{-1}J)$

Agora está novamente no ranking 1, a saber $D^{-1}J$ $(d_1^{-1}\ldots d_n^{-1})^T(1\ldots 1)$ . Observe que o segundo determinante é simplesmente:

f (1) = det (I - D^{- 1} J)

$f(1) = \det(I- D^{-1}J)$

onde é o polinómio característico de . Como uma matriz de classificação 1, deve ter (pelo menos) fatores de para explicar seu espaço nulo. O valor próprio "ausente" é $f(x)$ $D^{-1}J$ $f(x)$ $n-1$ $x$ $\sum d_i^{-1}$ , como pode ser visto no cálculo:

D^{- 1} J (d_{1}^{- 1} \dots d_{n}^{- 1})^{T} = (\sum d_{i}^{- 1}) (d_{1}^{- 1} \dots d_{n}^{- 1})^{T}

$D^{-1}J\; (d_1^{-1}\ldots d_n^{-1})^T = \left(\sum d_i^{-1} \right) (d_1^{-1}\ldots d_n^{-1})^T$

Segue-se que a característica polinomial , e é tal como mostrado acima para , . $f(x) = x^{n-1} (x - \sum d_i^{-1})$ $f(1)$ $\det(I - D^{-1}J)$ $1 - \sum d_i^{-1}$

Observe também que se , então , uma matriz diagonal cujo determinante é simplesmente o produto de suas entradas diagonais. $Q = D-J$ $12I - Q - J = 12I - D + J - J = 12I - D$

hardmath
fonte

Hm .. é de fato onde é a matriz de adjacência de então acho que esse resultado pode não estar correto. Em particular, isso implicaria que o número de árvores abrangidas de um gráfico é determinado por sua sequência de graus que não é válida.

Q

$Q$

D - A

$D-A$

A

$A$

G

$G$

G

$G$

Jernej

As entradas fora da diagonal de geralmente conterão 0 e -1. A decomposição do sugerida por Victor tira proveito da simetria e reduz o termo principal na contagem de operações de para . Existe uma abordagem inteira exata, mas isso provavelmente não é necessário para sua matriz e entradas de tamanho modesto. Se eu entendo a construção, é definido positivamente pela mesma razão que é.

Q

$Q$

L D L^{T}

$LDL^T$

\frac{2}{3} n^{3}

$\frac{2}{3}n^3$

\frac{1}{3} n^{3}

$\frac{1}{3}n^3$

12 I - Q - J

$12I-Q-J$

Q

$Q$

hardmath

@Jernej: Se você acredita que algo que afirmei está incorreto, criei uma sala de bate-papo com base nessa pergunta, na qual a discussão pode ser encadeada sem comentários desnecessários aqui.

hardmath

Se você possui uma maneira estruturada de enumerar os gráficos dos quais deseja calcular os determinantes, talvez seja possível encontrar atualizações de baixa classificação que o transferem de um gráfico para outro.

Nesse caso, você poderia usar o lema determinante da matriz para calcular de forma barata o determinante do gráfico subsequente a ser enumerado usando seu conhecimento do determinante do gráfico atual.

Ou seja, para uma matriz e vetores : Isso pode ser generalizado se U e V forem matrizes e é : $A$ $u,v$

det (A + u v^{T}) = (1 + v^{T} A^{- 1} u) det (A)

$\det(A+uv^T)=(1+v^T A^{-1} u)\det(A)$

n \times m

$n \times m$

A

$A$

n \times n

$n \times n$

det (A + U V^{T}) = det (I_{m} + V^{T} A^{- 1} U) det (A)

$\det(A+UV^T)=\det(I_m + V^TA^{-1}U)\det(A)$

Para calcular com eficiência o inverso, você pode usar a fórmula de Sherman-Morrison para obter o inverso da matriz subsequente da matriz atual:

(A + u v^{T})^{- 1} = A^{- 1} - \frac{A^{- 1} u v^{T} A^{- 1}}{1 + v^{T} A^{- 1} u}

$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1} uv^T A^{-1}}{1+v^T A^{-1} u}$

Richard
fonte