système linéaire le plus rapide à résoudre pour les petites matrices carrées (10x10)

Je suis très intéressé par l'optimisation de la résolution de systèmes linéaires pour les petites matrices (10x10), parfois appelées minuscules matrices. Existe-t-il une solution prête à cela? La matrice peut être supposée non singulière.

Ce solveur doit être exécuté plus de 1 000 000 de fois en microsecondes sur un processeur Intel. Je parle du niveau d'optimisation utilisé dans les jeux informatiques. Peu importe si je le code dans un assemblage et une architecture spécifiques, ou si j'étudie des réductions de compromis de précision ou de fiabilité et que j'utilise des hacks à virgule flottante (j'utilise l'indicateur de compilation -ffast-math, pas de problème). La résolution peut même échouer pendant environ 20% du temps!

Le partialPivLu d'Eigen est le plus rapide de mon benchmark actuel, surpassant LAPACK lorsqu'il est optimisé avec -O3 et un bon compilateur. Mais maintenant, je suis sur le point de fabriquer un solveur linéaire personnalisé. Tout avis serait grandement apprécié. Je vais rendre ma solution open source et je vais avoir des informations clés dans les publications, etc.

Connexes: Vitesse de résolution d'un système linéaire avec une matrice diagonale de bloc Quelle est la méthode la plus rapide pour inverser des millions de matrices? https://stackoverflow.com/q/50909385/1489510

L'utilisation d'un type de matrice propre où le nombre de lignes et de colonnes est codé dans le type au moment de la compilation vous donne un avantage sur LAPACK, où la taille de la matrice n'est connue qu'au moment de l'exécution. Ces informations supplémentaires permettent au compilateur d'effectuer un déroulement complet ou partiel de la boucle, éliminant ainsi de nombreuses instructions de branchement. Si vous envisagez d'utiliser une bibliothèque existante plutôt que d'écrire vos propres noyaux, avoir un type de données où la taille de la matrice peut être incluse en tant que paramètres de modèle C ++ sera probablement essentiel. La seule autre bibliothèque que je connaisse qui le fasse est Blaze , donc cela pourrait valoir la peine d'être comparé à Eigen.

Si vous décidez de lancer votre propre implémentation, vous trouverez peut-être que ce que PETSc fait pour son format CSR de bloc est un exemple utile, bien que PETSc lui-même ne soit probablement pas le bon outil pour ce que vous avez en tête. Plutôt que d'écrire une boucle, ils écrivent explicitement chaque opération pour les petites multiplications matricielles (voir ce fichier dans leur référentiel). Cela garantit qu'il n'y a pas d'instructions de branchement comme vous pourriez obtenir avec une boucle. Les versions du code avec des instructions AVX sont un bon exemple de la façon d'utiliser réellement les extensions vectorielles. Par exemple, cette fonction utilise le__m256dtype de données pour fonctionner simultanément sur quatre doubles en même temps. Vous pouvez obtenir une amélioration sensible des performances en écrivant explicitement toutes les opérations à l'aide d'extensions vectorielles, uniquement pour la factorisation LU au lieu de la multiplication matrice-vecteur. Plutôt que d'écrire le code C à la main, vous feriez mieux d'utiliser un script pour le générer. Il peut également être amusant de voir s'il y a une différence de performance appréciable lorsque vous réorganisez certaines opérations pour mieux tirer parti du pipelining des instructions.

Vous pouvez également obtenir un certain kilométrage de l'outil STOKE , qui explorera au hasard l'espace des transformations de programme possibles pour trouver une version plus rapide.

Une autre idée pourrait être d'utiliser une approche générative (un programme écrivant un programme). Créez un (méta) programme qui crache la séquence d'instructions C / C ++ pour exécuter ** LU non pivoté sur un système 10x10 .. en prenant essentiellement le nid de boucle k / i / j et en l'aplatissant en O (1000) ou deux lignes de arithmétique scalaire. Ensuite, alimentez ce programme généré dans le compilateur d'optimisation. Ce que je pense est en quelque sorte intéressant ici, la suppression des boucles expose chaque dépendance aux données et sous-expression redondante, et donne au compilateur la possibilité maximale de réorganiser les instructions afin qu'elles correspondent bien au matériel réel (par exemple, le nombre d'unités d'exécution, les dangers / blocages, donc sur).

Si vous connaissez toutes les matrices (ou même seulement quelques-unes d'entre elles), vous pouvez améliorer le débit en appelant SIMD intrinsèques / fonctions (SSE / AVX) au lieu du code scalaire. Ici, vous exploiteriez le parallélisme embarrassant entre les instances, au lieu de chasser tout parallélisme au sein d'une seule instance. Par exemple, vous pouvez effectuer 4 LU double précision simultanément en utilisant les intrinsèques AVX256, en plaçant 4 matrices "à travers" le registre et en effectuant les mêmes opérations ** sur chacune d'elles.

** D'où l'accent mis sur les LU non pivotés. Le pivot gâche cette approche de deux manières. Tout d'abord, il introduit des branches en raison de la sélection du pivot, ce qui signifie que vos dépendances de données ne sont pas aussi parfaitement connues. Deuxièmement, cela signifie que différents "emplacements" SIMD devraient faire des choses différentes, car l'instance A peut pivoter différemment de l'instance B. Donc, si vous poursuivez tout cela, je suggère de faire pivoter statiquement vos matrices avant le calcul (permuter la plus grande entrée de chaque colonne à la diagonale).

Votre question mène à deux considérations différentes.

Tout d'abord, vous devez choisir le bon algorithme. Par conséquent, la question de savoir si les matrices ont une structure doit être examinée. Par exemple, lorsque les matrices sont symétriques, une décomposition de Cholesky est plus efficace que LU. Lorsque vous n'avez besoin que d'une précision limitée, une méthode itérative peut être plus rapide.

Deuxièmement, vous devez implémenter l'algorithme efficacement. Pour ce faire, vous devez connaître le goulot d'étranglement de votre algorithme. Votre implémentation est-elle liée par la vitesse du transfert de mémoire ou par la vitesse du calcul. Puisque vous ne considérez que$10 \times 10$matrices, votre matrice doit s'intégrer complètement dans le cache du processeur. Ainsi, vous devez utiliser les unités SIMD (SSE, AVX, etc.) et les cœurs de votre processeur, pour faire autant de calculs par cycle que possible.

En tout, la réponse à votre question dépend fortement du matériel et des matrices que vous considérez. Il n'y a probablement pas de réponse définitive et vous devez essayer quelques choses pour trouver une méthode optimale.

J'essaierais l'inversion par blocs.

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

Eigen utilise une routine optimisée pour calculer l'inverse d'une matrice 4x4, qui est probablement la meilleure que vous obtiendrez. Essayez d'utiliser autant que possible.

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

En haut à gauche: 8x8. En haut à droite: 8x2. En bas à gauche: 2x8. En bas à droite: 2x2. Inversez le 8x8 en utilisant le code d'inversion 4x4 optimisé. Le reste est constitué de produits matriciels.

EDIT: L'utilisation de blocs 6x6, 6x4, 4x6 et 4x4 s'est avérée être un peu plus rapide que ce que j'ai décrit ci-dessus.

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;
}

Voici les résultats d'une analyse comparative utilisant un million de Eigen::Matrix<double,10,10>::Random()matrices et de Eigen::Matrix<double,10,1>::Random()vecteurs. Dans tous mes tests, mon inverse est toujours plus rapide. Ma routine de résolution consiste à calculer l'inverse puis à le multiplier par un vecteur. Parfois, c'est plus rapide que Eigen, parfois ce n'est pas le cas. Ma méthode de banc d'essai peut être défectueuse (n'a pas désactivé le turbo boost, etc.). De plus, les fonctions aléatoires d'Eigen peuvent ne pas être représentatives de données réelles.

• Inverse du pivot partiel propre: 3036 millisecondes
• Mon inverse avec bloc supérieur 8x8: 1638 millisecondes
• Mon inverse avec bloc supérieur 6x6: 1234 millisecondes
• Résolution de pivot partiel propre: 1791 millisecondes
• Ma résolution avec le bloc supérieur 8x8: 1739 millisecondes
• Ma résolution avec le bloc supérieur 6x6: 1286 millisecondes

Je suis très intéressé de voir si quelqu'un peut optimiser cela davantage, car j'ai une application par éléments finis qui inverse un million de matrices gazeuses 10x10 (et oui, j'ai besoin de coefficients individuels de l'inverse, donc résoudre directement un système linéaire n'est pas toujours une option) .