o sistema linear mais rápido resolve pequenas matrizes quadradas (10x10)

9

Estou muito interessado em otimizar a solução do sistema linear para matrizes pequenas (10x10), às vezes chamadas de matrizes minúsculas . Existe uma solução pronta para isso? A matriz pode ser considerada não singular.

Esse solucionador deve ser executado em excesso de 1 000 000 vezes em microssegundos em uma CPU Intel. Estou falando do nível de otimização usado em jogos de computador. Independentemente de codificá-lo em montagem e arquitetura específica, ou estudar reduções de trocas de precisão ou confiabilidade e usar hacks de ponto flutuante (eu uso o sinalizador de compilação -ffast-math, não há problema). A solução pode até falhar por cerca de 20% do tempo!

O parcialPivLu da Eigen é o mais rápido no meu benchmark atual, superando o LAPACK quando otimizado com -O3 e um bom compilador. Mas agora estou no ponto de criar manualmente um solucionador linear personalizado. Qualquer conselho seria muito apreciado. Farei minha solução de código aberto e conheço informações importantes em publicações, etc.

Relacionado: Velocidade de resolução de sistemas lineares com matriz diagonal de blocos Qual é o método mais rápido para inverter milhões de matrizes? https://stackoverflow.com/q/50909385/1489510

rfabbri
fonte
7
Parece uma meta esticada. Vamos supor que usamos o Skylake-X Xeon Platinum 8180 mais rápido, com um pico de rendimento teórico de 4 TFLOPs de precisão única, e que um sistema 10x10 requer cerca de 700 operações de ponto flutuante (aproximadamente 2n ** 3/3) para serem resolvidas. Em seguida, um lote de 1 milhão desses sistemas poderia ser teoricamente resolvido em 175 microssegundos. Esse é um número de velocidade da luz que não pode exceder. Você pode compartilhar o desempenho que está alcançando atualmente com o código existente mais rápido? BTW, os dados são de precisão única ou dupla precisão?
Njuffa 02/03/19
@njuffa sim, eu pretendia atingir perto de 1ms, mas micro é outra história. Para micro, considerei explorar a estrutura inversa incremental no lote, detectando matrizes semelhantes, que ocorrem com frequência. O Perf é atualmente na faixa de 10 a 500ms, dependendo do processador. A precisão é dupla ou mesmo complexa. A precisão única é mais lenta.
precisa saber é o seguinte
@njuffa posso reduzir ou até de precisão para a velocidade
rfabbri
2
Parece que a precisão não é sua prioridade. Para seu objetivo, talvez seja útil um método iterativo truncado em um número relativamente pequeno de avaliações? Especialmente se você tiver um palpite inicial razoável.
Spencer Bryngelson 03/03/19
11
Você gira? Você poderia fazer uma fatoração QR em vez da eliminação gaussiana. Você intercala seus sistemas para poder usar as instruções do SIMD e executar vários sistemas ao mesmo tempo? Você escreve programas lineares sem loops e sem endereçamento indireto? Qual precisão você deseja e como condicionarei seu sistema? Eles têm alguma estrutura que possa ser explorada.
314 Carl Carl Christian

Respostas:

7

O uso de um tipo de matriz Eigen em que o número de linhas e colunas é codificado no tipo em tempo de compilação fornece uma vantagem sobre o LAPACK, onde o tamanho da matriz é conhecido apenas em tempo de execução. Essas informações extras permitem ao compilador desenrolar total ou parcialmente o loop, eliminando muitas instruções de ramificação. Se você estiver olhando para usar uma biblioteca existente em vez de escrever seus próprios kernels, provavelmente será essencial ter um tipo de dados em que o tamanho da matriz possa ser incluído como parâmetros do modelo C ++. A única outra biblioteca que conheço que faz isso é a chama , de modo que talvez valha a pena comparar o Eigen.

Se você decidir implementar sua própria implementação, poderá descobrir o que o PETSc faz para o formato CSR de bloco ser um exemplo útil, embora o próprio PETSc provavelmente não seja a ferramenta certa para o que você tem em mente. Em vez de escrever um loop, eles escrevem todas as operações para pequenos vetores matriciais explicitamente (veja este arquivo em seu repositório). Isso garante que não há instruções de ramificação, como você pode obter com um loop. As versões do código com instruções AVX são um bom exemplo de como realmente usar extensões de vetor. Por exemplo, esta função usa o__m256dtipo de dados para operar simultaneamente em quatro duplos ao mesmo tempo. Você pode obter um aumento considerável no desempenho escrevendo explicitamente todas as operações usando extensões de vetor, apenas para fatoração LU em vez de multiplicar vetor matriz. Em vez de realmente escrever o código C manualmente, seria melhor usar um script para gerá-lo. Também pode ser divertido verificar se há uma diferença significativa de desempenho quando você reordena algumas das operações para aproveitar melhor o pipelining de instruções.

Você também pode obter alguma milhagem da ferramenta STOKE , que explorará aleatoriamente o espaço de possíveis transformações de programas para encontrar uma versão mais rápida.

Daniel Shapero
fonte
tx. Eu já uso Eigen como Map <const Matrix <complex, 10, 10>> AA (A) com sucesso. irá verificar as outras coisas.
rfabbri
O Eigen também possui o AVX e até um cabeçalho complex.h. Por que PETSc para isso? É difícil competir com Eigen neste caso. Eu especializei o Eigen ainda mais para o meu problema e com uma estratégia de pivô aproximada que, em vez de assumir o máximo de uma coluna, troca imediatamente o pivô quando encontra outro que é três ordens de magnitude maior.
rfabbri
11
@rfabbri Eu não estava sugerindo que você usasse o PETSc para isso, apenas que o que eles fazem nessa instância específica pode ser instrutivo. Eu editei a resposta para deixar isso mais claro.
Daniel Shapero
4

Outra idéia poderia ser usar uma abordagem generativa (um programa escrevendo um programa). Crie um (meta) programa que cite a sequência de instruções C / C ++ para executar LU ** não dinâmica em um sistema 10x10. Basicamente, pegue o ninho de loop k / i / j e aplaná-lo em O (1000) ou mais linhas de aritmética escalar. Em seguida, alimente esse programa gerado para qualquer compilador otimizador. O que eu acho interessante aqui, é remover os loops que expõem todas as dependências de dados e subexpressão redundante e oferece ao compilador a oportunidade máxima de reordenar as instruções para que elas mapeiem bem o hardware real (por exemplo, número de unidades de execução, perigos / paradas, portanto em).

Se você conhecer todas as matrizes (ou apenas algumas delas), poderá melhorar a taxa de transferência chamando as funções / intrínsecas do SIMD (SSE / AVX) em vez do código escalar. Aqui você exploraria o paralelismo embaraçoso entre as instâncias, em vez de perseguir qualquer paralelismo dentro de uma única instância. Por exemplo, você pode executar quatro LUs de precisão dupla simultaneamente usando intrínsecas AVX256, agrupando 4 matrizes "em todo o registro" e executando as mesmas operações ** em todas elas.

** Daí o foco na LU não dinâmica. A articulação estraga essa abordagem de duas maneiras. Primeiro, ele introduz ramificações devido à seleção dinâmica, o que significa que suas dependências de dados não são tão conhecidas. Segundo, significa que diferentes "slots" do SIMD teriam que fazer coisas diferentes, porque a instância A pode girar de maneira diferente da instância B. Portanto, se você prosseguir com isso, sugiro rodar estaticamente suas matrizes antes do cálculo (permita a maior entrada de cada coluna para diagonal).

rchilton1980
fonte
como as matrizes são muito pequenas, talvez seja possível eliminar a rotação se elas forem pré-dimensionadas. Nem mesmo pré-girando as matrizes. Tudo o que precisamos é que as entradas estejam entre 2 e 3 ordens de magnitude uma da outra.
rfabbri
2

Sua pergunta leva a duas considerações diferentes.

Primeiro, você precisa escolher o algoritmo certo. Portanto, a questão de saber se as matrizes têm alguma estrutura deve ser considerada. Por exemplo, quando as matrizes são simétricas, uma decomposição de Cholesky é mais eficiente que a LU. Quando você precisa apenas de uma quantidade limitada de precisão, um método iterativo pode ser mais rápido.

Segundo, você precisa implementar o algoritmo com eficiência. Para fazer isso, você precisa conhecer o gargalo do seu algoritmo. A sua implementação está vinculada à velocidade da transferência de memória ou à velocidade da computação. Como você considera apenas10×10matrizes, sua matriz deve caber completamente no cache da CPU. Portanto, você deve usar as unidades SIMD (SSE, AVX etc.) e os núcleos do seu processador para fazer o maior número possível de cálculos por ciclo.

Ao todo, a resposta à sua pergunta depende muito do hardware e das matrizes que você considera. Provavelmente não existe uma resposta definitiva e você deve tentar algumas coisas para encontrar um método ideal.

H. Rittich
fonte
Até agora, o Eigen já otimiza bastante, usa SEE, AVX, etc, e tentei métodos iterativos em um teste preliminar e eles não ajudaram. Eu tentei o Intel MKL, mas não melhor que o Eigen com sinalizadores GCC otimizados. Atualmente, estou tentando criar algo melhor e mais simples que o Eigen e fazer testes mais detalhados com métodos iterativos.
precisa saber é o seguinte
1

Eu tentaria inversão em blocos.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen usa uma rotina otimizada para calcular o inverso de uma matriz 4x4, que é provavelmente o melhor que você obterá. Tente usar isso o máximo possível.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Superior esquerdo: 8x8. Superior direito: 8x2. Parte inferior esquerda: 2x8. Em baixo à direita: 2x2. Inverta o 8x8 usando o código de inversão 4x4 otimizado. O resto são produtos matriciais.

EDIT: O uso de blocos 6x6, 6x4, 4x6 e 4x4 mostrou ser um pouco mais rápido do que o que descrevi acima.

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Aqui estão os resultados de uma execução de benchmark usando um milhão de Eigen::Matrix<double,10,10>::Random()matrizes e Eigen::Matrix<double,10,1>::Random()vetores. Em todos os meus testes, meu inverso é sempre mais rápido. Minha rotina de resolução envolve calcular o inverso e depois multiplicá-lo por um vetor. Às vezes é mais rápido que Eigen, às vezes não. Meu método de marcação de bancada pode ter falhas (não desativou o turbo boost, etc.). Além disso, as funções aleatórias de Eigen podem não representar dados reais.

  • Pivô parcial do próprio inverso: 3036 milissegundos
  • Meu inverso com bloco superior de 8x8: 1638 milissegundos
  • Meu inverso com bloco superior 6x6: 1234 milissegundos
  • Solução de pivô parcial de Eigen: 1791 milissegundos
  • Minha resolução com bloco superior 8x8: 1739 milissegundos
  • Minha resolução com bloco superior 6x6: 1286 milissegundos

Estou muito interessado em ver se alguém pode otimizar isso ainda mais, pois tenho um aplicativo de elementos finitos que inverte um zilhão de matrizes 10x10 (e sim, preciso de coeficientes individuais do inverso, portanto, resolver diretamente um sistema linear nem sempre é uma opção) .

Charlie S
fonte