Cet article a été initialement publié sur le blog Quansight Labs et a été modifié et republié ici avec la permission de Quansight.
Les applications Web apparaissent rapidement comme une nouvelle frontière pour le calcul scientifique haute performance et les expériences utilisateur basées sur l'IA. À la base de la révolution ML/AI se trouve l’algèbre linéaire, une branche des mathématiques concernant les équations linéaires et leurs représentations dans les espaces vectoriels et via des matrices. LAPACK ("Linear Algebra Package") est une bibliothèque logicielle fondamentale pour l'algèbre linéaire numérique, fournissant des implémentations robustes et éprouvées d'opérations matricielles courantes. . Bien que LAPACK soit un composant fondamental de la plupart des langages et bibliothèques de programmation informatique numérique, une implémentation LAPACK complète et de haute qualité, adaptée aux contraintes uniques du Web, ne s'est pas encore matérialisée. C'est... jusqu'à maintenant.
Plus tôt cette année, j'ai eu la grande chance d'être stagiaire d'été chez Quansight Labs, la division d'intérêt public de Quansight et un leader de l'écosystème scientifique Python. Au cours de mon stage, j'ai travaillé pour ajouter le support initial de LAPACK à stdlib, une bibliothèque fondamentale pour le calcul scientifique écrite en C et JavaScript et optimisée pour une utilisation dans les navigateurs Web et d'autres environnements Web natifs, tels que Node.js et Deno. Dans cet article de blog, je discuterai de mon parcours, de certains défis attendus et inattendus (!) et du chemin à parcourir. J'espère que ce travail, avec un peu de chance, fournira un élément essentiel pour faire des navigateurs Web un environnement de premier ordre pour le calcul numérique et l'apprentissage automatique et présage un avenir d'applications Web plus puissantes basées sur l'IA.
Cela vous semble intéressant ? C'est parti !
Les lecteurs de ce blog qui connaissent LAPACK ne connaissent probablement pas intimement le monde sauvage des technologies Web. Pour ceux qui viennent du monde du calcul numérique et scientifique et qui sont familiers avec l’écosystème scientifique Python, la façon la plus simple de considérer stdlib est comme une bibliothèque de calcul scientifique open source dans le moule de NumPy et SciPy. Il fournit des structures de données de tableaux multidimensionnels et des routines associées pour les mathématiques, les statistiques et l'algèbre linéaire, mais utilise JavaScript, plutôt que Python, comme langage de script principal. En tant que tel, stdlib se concentre sur l’écosystème Web et ses paradigmes de développement d’applications. Cette orientation nécessite des décisions intéressantes en matière de conception et d'architecture de projet, qui rendent stdlib plutôt unique par rapport aux bibliothèques plus traditionnelles conçues pour le calcul numérique.
Pour prendre NumPy comme exemple, NumPy est une bibliothèque monolithique unique, où tous ses composants, en dehors des dépendances tierces facultatives telles que OpenBLAS, forment une unité unique et indivisible. On ne peut pas simplement installer des routines NumPy pour la manipulation de tableaux sans installer tout NumPy. Si vous déployez une application qui n'a besoin que de l'objet ndarray de NumPy et de quelques-unes de ses routines de manipulation, installer et regrouper l'ensemble de NumPy signifie inclure une quantité considérable de "code mort". Dans le langage du développement Web, nous dirions que NumPy n'est pas « arborescent ». Pour une installation NumPy normale, cela implique au moins 30 Mo d'espace disque et au moins 15 Mo d'espace disque pour une version personnalisée qui exclut toutes les instructions de débogage. Pour SciPy, ces chiffres peuvent atteindre respectivement 130 Mo et 50 Mo. Inutile de dire que l'envoi d'une bibliothèque de 15 Mo dans une application Web pour quelques fonctions seulement n'est pas une solution, en particulier pour les développeurs ayant besoin de déployer des applications Web sur des appareils ayant une mauvaise connectivité réseau ou des contraintes de mémoire.
Compte tenu des contraintes uniques du développement d'applications Web, stdlib adopte une approche ascendante pour sa conception, dans laquelle chaque unité de fonctionnalité peut être installée et consommée indépendamment des parties non liées et inutilisées de la base de code. En adoptant une architecture logicielle décomposable et une modularité radicale, stdlib offre aux utilisateurs la possibilité d'installer et d'utiliser exactement ce dont ils ont besoin, avec peu ou pas d'excès de code au-delà d'un ensemble souhaité d'API et de leurs dépendances explicites, garantissant ainsi des empreintes mémoire plus petites, tailles et un déploiement plus rapide.
À titre d'exemple, supposons que vous travaillez avec deux piles de matrices (c'est-à-dire des tranches bidimensionnelles de cubes tridimensionnels) et que vous souhaitiez sélectionner une tranche sur deux et effectuer l'opération BLAS commune y = a * x, où x et y sont des ndarrays et a est une constante scalaire. Pour faire cela avec NumPy, vous devez d'abord installer tout NumPy
pip install numpy
puis effectuez les différentes opérations
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Avec stdlib, en plus d'avoir la possibilité d'installer le projet comme une bibliothèque monolithique, vous pouvez installer les différentes unités de fonctionnalités sous forme de packages séparés
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
puis effectuez les différentes opérations
// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';
// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;
// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);
Il est important de noter que non seulement vous pouvez installer indépendamment l'un des plus de 4 000 packages de stdlib, mais vous pouvez également corriger, améliorer et remixer n'importe lequel de ces packages en créant un référentiel GitHub associé (par exemple, voir @stdlib/ndarray-fancy ). En définissant des couches explicites d'abstraction et des arbres de dépendances, stdlib vous offre la liberté de choisir la bonne couche d'abstraction pour votre application. À certains égards, c'est une idée simple (et peut-être peu orthodoxe si vous êtes habitué à la conception de bibliothèques de logiciels scientifiques classiques), mais, lorsqu'elle est étroitement intégrée à la plate-forme Web, elle a des conséquences puissantes et crée de nouvelles possibilités passionnantes !
D'accord, alors peut-être que votre intérêt a été éveillé ; stdlib semble intrigant. Mais qu’est-ce que cela a à voir avec LAPACK dans les navigateurs Web ? Eh bien, l'un de nos objectifs l'été dernier était d'appliquer la philosophie stdlib (des petits packages à portée étroite qui font une chose et qui le font bien) en amenant LAPACK sur le Web.
Mais attendez, dites-vous ! C'est une entreprise extrême. LAPACK est vaste, avec environ 1 700 routines, et la mise en œuvre ne serait-ce que de 10 % d’entre elles dans un délai raisonnable constitue un défi de taille. Ne serait-il pas préférable de simplement compiler LAPACK en WebAssembly, une cible de compilation portable pour les langages de programmation tels que C, Go et Rust, qui permet le déploiement sur le Web, et de mettre fin à cela ?
Malheureusement, cette approche pose plusieurs problèmes.
Pour illustrer le dernier point, revenons à la routine BLAS daxpy, qui effectue l'opération y = a*x y et où x et y sont des vecteurs striés et a une constante scalaire. Si elle est implémentée en C, une implémentation de base pourrait ressembler à l'extrait de code suivant.
pip install numpy
Après la compilation sur WebAssembly et le chargement du binaire WebAssembly dans notre application Web, nous devons effectuer une série d'étapes avant de pouvoir appeler la routine c_daxpy depuis JavaScript. Tout d’abord, nous devons instancier un nouveau module WebAssembly.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Ensuite, nous devons définir la mémoire du module et créer une nouvelle instance du module WebAssembly.
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
Après avoir créé une instance de module, nous pouvons maintenant invoquer la routine BLAS exportée. Cependant, si les données sont définies en dehors de la mémoire du module, nous devons d'abord copier ces données dans l'instance de mémoire et toujours le faire dans l'ordre des octets petit-endien.
// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';
// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;
// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);
Maintenant que les données sont écrites dans la mémoire du module, nous pouvons appeler la routine c_daxpy.
void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
int ix;
int iy;
int i;
if (N <= 0) {
return;
}
if (alpha == 0.0) {
return;
}
if (strideX < 0) {
ix = (1-N) * strideX;
} else {
ix = 0;
}
if (strideY < 0) {
iy = (1-N) * strideY;
} else {
iy = 0;
}
for (i = 0; i < N; i++) {
Y[iy] += alpha * X[ix];
ix += strideX;
iy += strideY;
}
return;
}
Et, enfin, si nous devons transmettre les résultats à une bibliothèque en aval sans prise en charge des "pointeurs" de mémoire WebAssembly (c'est-à-dire les décalages d'octets), telle que D3, pour une visualisation ou une analyse plus approfondie, nous devons copier les données du module mémoire au tableau de sortie d'origine.
const binary = new UintArray([...]); const mod = new WebAssembly.Module(binary);
Cela fait beaucoup de travail rien que de calculer y = a*x y. En revanche, comparez-le à une implémentation JavaScript simple, qui pourrait ressembler à l'extrait de code suivant.
// Initialize 10 pages of memory and allow growth to 100 pages:
const mem = new WebAssembly.Memory({
'initial': 10, // 640KiB, where each page is 64KiB
'maximum': 100 // 6.4MiB
});
// Create a new module instance:
const instance = new WebAssembly.Instance(mod, {
'env': {
'memory': mem
}
});
Avec l'implémentation JavaScript, nous pouvons alors appeler directement daxpy avec nos données définies en externe sans le mouvement de données requis dans l'exemple WebAssembly ci-dessus.
// External data:
const xdata = new Float64Array([...]);
const ydata = new Float64Array([...]);
// Specify a vector length:
const N = 5;
// Specify vector strides (in units of elements):
const strideX = 2;
const strideY = 4;
// Define pointers (i.e., byte offsets) for storing two vectors:
const xptr = 0;
const yptr = N * 8; // 8 bytes per double
// Create a DataView over module memory:
const view = new DataView(mem.buffer);
// Resolve the first indexed elements in both `xdata` and `ydata`:
let offsetX = 0;
if (strideX < 0) {
offsetX = (1-N) * strideX;
}
let offsetY = 0;
if (strideY < 0) {
offsetY = (1-N) * strideY;
}
// Write data to the memory instance:
for (let i = 0; i < N; i++) {
view.setFloat64(xptr+(i*8), xdata[offsetX+(i*strideX)], true);
view.setFloat64(yptr+(i*8), ydata[offsetY+(i*strideY)], true);
}
Au moins dans ce cas, non seulement l'approche WebAssembly est moins ergonomique, mais, comme on pouvait s'y attendre étant donné le mouvement des données requis, il y a également un impact négatif sur les performances, comme le démontre la figure suivante.
Figure 1 : Comparaison des performances des implémentations C, JavaScript et WebAssembly (Wasm) de stdlib pour la routine BLAS daxpy pour augmenter la longueur des tableaux (axe des x). Dans le benchmark Wasm (copie), les données d'entrée et de sortie sont copiées vers et depuis la mémoire Wasm, ce qui entraîne de moins bonnes performances.
Dans la figure ci-dessus, j'affiche une comparaison des performances des implémentations C, JavaScript et WebAssembly (Wasm) de stdlib pour la routine BLAS daxpy pour augmenter la longueur des tableaux, comme indiqué le long de l'axe des x. L’axe des y montre un taux normalisé par rapport à une implémentation C de base. Dans le benchmark Wasm, les données d'entrée et de sortie sont allouées et manipulées directement dans la mémoire du module WebAssembly et, dans le benchmark Wasm (copie), les données d'entrée et de sortie sont copiées vers et depuis la mémoire du module WebAssembly, comme indiqué ci-dessus. À partir du graphique, nous pouvons observer ce qui suit :
Dans l’ensemble, WebAssembly peut offrir des améliorations de performances ; cependant, la technologie n’est pas une solution miracle et doit être utilisée avec précaution afin de réaliser les gains souhaités. Et même lorsqu’ils offrent des performances supérieures, ces gains doivent être mis en balance avec les coûts liés à une complexité accrue, à des tailles de bundles potentiellement plus grandes et à des chaînes d’outils plus complexes. Pour de nombreuses applications, une simple implémentation JavaScript fera très bien l'affaire.
Maintenant que j'ai engagé des poursuites contre la simple compilation de l'intégralité de LAPACK dans WebAssembly et que je l'ai arrêté, où cela nous mène-t-il ? Eh bien, si nous voulons adopter la philosophie stdlib, cela nous laisse dans le besoin d'une modularité radicale.
Adopter une modularité radicale, c'est reconnaître que ce qui est le mieux est hautement contextuel et, en fonction des besoins et des contraintes des applications utilisateur, les développeurs ont besoin de flexibilité pour choisir la bonne abstraction. Si un développeur écrit une application Node.js, cela peut impliquer une liaison à des bibliothèques optimisées pour le matériel, telles que OpenBLAS, Intel MKL ou Apple Accelerate, afin d'obtenir des performances supérieures. Si un développeur déploie une application Web nécessitant un petit ensemble de routines numériques, JavaScript est probablement l'outil idéal pour cette tâche. Et si un développeur travaille sur une application WebAssembly volumineuse et gourmande en ressources (par exemple, pour l'édition d'images ou un moteur de jeu), il sera alors primordial de pouvoir compiler facilement des routines individuelles dans le cadre d'une application plus vaste. Bref, nous voulons un LAPACK radicalement modulaire.
Ma mission était de jeter les bases d'un tel projet, de résoudre les problèmes et de trouver les lacunes, et, espérons-le, de nous rapprocher de quelques pas de l'algèbre linéaire haute performance sur le Web. Mais à quoi ressemble une modularité radicale ? Tout commence par l'unité fondamentale de fonctionnalité, le package.
Chaque package dans stdlib est son propre élément autonome, contenant des tests co-localisés, des benchmarks, des exemples, de la documentation, des fichiers de build et des métadonnées associées (y compris l'énumération de toutes les dépendances) et définissant une surface API claire avec le monde extérieur. . Afin d'ajouter le support LAPACK à stdlib, cela signifie créer un package autonome distinct pour chaque routine LAPACK avec la structure suivante :
pip install numpy
En bref,
Compte tenu des exigences exigeantes de stdlib en matière de documentation et de tests, l'ajout de la prise en charge de chaque routine représente une quantité de travail décente, mais le résultat final est un code robuste, de haute qualité et, surtout, modulaire, adapté pour servir de base au calcul scientifique. sur le Web moderne. Mais assez de préliminaires ! Passons aux choses sérieuses !
En nous appuyant sur les efforts précédents qui ont ajouté le support BLAS à stdlib, nous avons décidé de suivre une approche similaire en plusieurs phases lors de l'ajout du support LAPACK dans laquelle nous donnons d'abord la priorité aux implémentations JavaScript et à leurs tests et documentation associés, puis, une fois les tests et la documentation présents. , remplissez les implémentations C et Fortran et toutes les liaisons natives associées aux bibliothèques optimisées pour le matériel. Cette approche nous permet de mettre quelques premiers points sur le tableau, pour ainsi dire, en présentant rapidement les API aux utilisateurs, en établissant des procédures de test et des références robustes et en étudiant les voies potentielles en matière d'outillage et d'automatisation avant de plonger dans les mauvaises herbes des chaînes d'outils de construction et des performances. optimisations. Mais par où commencer ?
Pour déterminer quelles routines LAPACK cibler en premier, j'ai analysé le code source Fortran de LAPACK pour générer un graphique d'appel. Cela m'a permis de déduire l'arbre de dépendances pour chaque routine LAPACK. Le graphe en main, j'ai ensuite effectué un tri topologique, m'aidant ainsi à identifier des routines sans dépendances et qui seront inévitablement des briques de base pour d'autres routines. Même si une approche approfondie dans laquelle je choisissais une routine particulière de haut niveau et travaillerais à rebours me permettrait d'obtenir une fonctionnalité spécifique, une telle approche pourrait m'enliser en essayant de mettre en œuvre des routines de complexité croissante. En me concentrant sur les « feuilles » du graphique, je pourrais donner la priorité aux routines couramment utilisées (c'est-à-dire les routines avec des degrés élevés) et ainsi maximiser mon impact en débloquant la possibilité de proposer plusieurs routines de niveau supérieur soit plus tard dans mes efforts ou par d'autres contributeurs.
Avec mon plan en main, j'avais hâte de me mettre au travail. Pour ma première routine, j'ai choisi dlaswp, qui effectue une série d'échanges de lignes sur une matrice rectangulaire générale selon une liste d'indices pivots fournie et qui est un élément clé des routines de décomposition LU de LAPACK. Et c'est à ce moment-là que mes défis ont commencé...
Avant mon stage Quansight Labs, j'étais (et je suis toujours !) un contributeur régulier à LFortran, un compilateur Fortran interactif moderne construit sur LLVM, et je me sentais assez confiant dans mes compétences Fortran. Cependant, l'un de mes premiers défis consistait simplement à comprendre ce qui est désormais considéré comme du code Fortran « hérité ». Je souligne ci-dessous trois obstacles initiaux.
LAPACK a été initialement écrit en FORTRAN 77 (F77). Bien que la bibliothèque ait été déplacée vers Fortran 90 dans la version 3.2 (2008), les conventions héritées persistent dans l'implémentation de référence. L'une des conventions les plus visibles est le formatage.
Les développeurs qui écrivent des programmes F77 l'ont fait en utilisant une disposition de formulaire fixe héritée des cartes perforées. Cette mise en page avait des exigences strictes concernant l'utilisation de colonnes de caractères :
Fortran 90 a introduit la mise en page de forme libre qui a supprimé les restrictions de longueur de colonne et de ligne et s'est installée ! comme caractère de commentaire. L'extrait de code suivant montre l'implémentation de référence pour la routine LAPACK dlacpy :
pip install numpy
L'extrait de code suivant montre la même routine, mais implémentée en utilisant la mise en page de forme libre introduite dans Fortran 90.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Comme on peut l'observer, en supprimant les restrictions de colonnes et en s'éloignant de la convention F77 consistant à écrire les spécificateurs en TOUTES MAJUSCULES, le code Fortran moderne est plus visiblement cohérent et donc plus lisible.
Une autre pratique courante dans les routines LAPACK est l'utilisation de structures de contrôle étiquetées. Par exemple, considérons l'extrait de code suivant dans lequel l'étiquette 10 doit correspondre à un CONTINUE correspondant.
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
Fortran 90 a évité le besoin de cette pratique et amélioré la lisibilité du code en permettant d'utiliser end do pour terminer une boucle do. Ce changement est affiché dans la version gratuite de dlacpy fournie ci-dessus.
Pour permettre une flexibilité dans la gestion de tableaux de différentes tailles, les routines LAPACK fonctionnent généralement sur des tableaux ayant une taille supposée. Dans la routine dlacpy ci-dessus, la matrice d'entrée A est déclarée comme étant un tableau bidimensionnel ayant une taille supposée selon l'expression A(LDA, *). Cette expression déclare que A a un nombre de lignes LDA et utilise * comme espace réservé pour indiquer que la taille de la deuxième dimension est déterminée par le programme appelant.
L'une des conséquences de l'utilisation de tableaux de taille supposée est que les compilateurs sont incapables d'effectuer une vérification des limites sur la dimension non spécifiée. Ainsi, la meilleure pratique actuelle consiste à utiliser des interfaces explicites et des tableaux de forme supposée (par exemple, A(LDA, :)) afin d'empêcher l'accès à la mémoire hors limites. Cela dit, l'utilisation de tableaux de forme supposée peut être problématique lorsqu'il faut transmettre des sous-matrices à d'autres fonctions, car cela nécessite un découpage qui entraîne souvent par les compilateurs la création de copies internes des données du tableau.
Inutile de dire qu'il m'a fallu un certain temps pour m'adapter aux conventions LAPACK et adopter un état d'esprit LAPACK. Cependant, étant un puriste, si je devais quand même porter des routines, je voulais au moins amener les routines que j'ai réussi à porter dans une ère plus moderne dans l'espoir d'améliorer la lisibilité du code et la maintenance future. Ainsi, après avoir discuté des choses avec les responsables de stdlib, j'ai décidé de migrer les routines vers Fortran 95, qui, bien que n'étant pas la dernière et la meilleure version de Fortran, semblait trouver le bon équilibre entre le maintien de l'apparence des implémentations d'origine, garantissant ( assez bon) compatibilité ascendante et utilisation de fonctionnalités syntaxiques plus récentes.
L'un des problèmes liés à la poursuite d'une approche ascendante pour ajouter le support de LAPACK est que les tests unitaires explicites pour les routines utilitaires de niveau inférieur sont souvent inexistants dans LAPACK. La suite de tests de LAPACK utilise en grande partie une philosophie de test hiérarchique dans laquelle le test des routines de niveau supérieur est supposé garantir que les routines de niveau inférieur qui en dépendent fonctionnent correctement dans le cadre d'un flux de travail global. Bien que l'on puisse affirmer qu'il est raisonnable de se concentrer sur les tests d'intégration plutôt que sur les tests unitaires pour les routines de niveau inférieur, dans la mesure où l'ajout de tests pour chaque routine pourrait potentiellement augmenter la charge de maintenance et la complexité du cadre de test de LAPACK, cela signifie que nous ne pouvions pas nous fier facilement aux versions antérieures. art pour les tests unitaires et nous devrions proposer nous-mêmes des tests unitaires autonomes complets pour chaque routine de niveau inférieur.
Dans la même veine que la couverture des tests, en dehors de LAPACK lui-même, trouver des exemples documentés du monde réel illustrant l'utilisation de routines de niveau inférieur était un défi. Alors que les routines LAPACK sont systématiquement précédées d'un commentaire de documentation fournissant des descriptions des arguments d'entrée et des valeurs de retour possibles, sans exemples de code, visualiser et analyser les valeurs d'entrée et de sortie attendues peut être difficile, en particulier lorsqu'il s'agit de matrices spécialisées. Et même si ni l'absence de tests unitaires ni d'exemples documentés ne sont la fin du monde, cela signifiait que l'ajout du support de LAPACK à stdlib serait plus compliqué que prévu. La rédaction de benchmarks, de tests, d'exemples et de documentation allait simplement nécessiter plus de temps et d'efforts, limitant potentiellement le nombre de routines que je pouvais mettre en œuvre pendant le stage.
Lors du stockage d'éléments matriciels en mémoire linéaire, on a deux choix : soit stocker les colonnes de manière contiguë, soit les lignes de manière contiguë (voir Figure 2). La première disposition de la mémoire est appelée ordre colonne-major et la seconde ordre ligne-major.
Figure 2 : Schéma illustrant le stockage des éléments matriciels dans la mémoire linéaire dans l'ordre (a) de colonne majeure (style Fortran) ou (b) de ligne majeure (style C). Le choix de la mise en page à utiliser est en grande partie une question de convention.
Le choix de la mise en page à utiliser est en grande partie une question de convention. Par exemple, Fortran stocke les éléments dans l’ordre principal des colonnes et C stocke les éléments dans l’ordre principal des lignes. Les bibliothèques de niveau supérieur, telles que NumPy et stdlib, prennent en charge les ordres majeurs des colonnes et des lignes, vous permettant de configurer la disposition d'un tableau multidimensionnel lors de la création du tableau.
pip install numpy
Bien qu'aucune disposition de la mémoire ne soit intrinsèquement meilleure que l'autre, l'organisation des données pour garantir un accès séquentiel conformément aux conventions du modèle de stockage sous-jacent est essentielle pour garantir des performances optimales. Les processeurs modernes sont capables de traiter les données séquentielles plus efficacement que les données non séquentielles, principalement grâce à la mise en cache du processeur qui, à son tour, exploite la localité spatiale de référence.
Pour démontrer l'impact sur les performances de l'accès aux éléments séquentiels et non séquentiels, considérons la fonction suivante qui copie tous les éléments d'une matrice MxN A vers une autre matrice MxN B et qui le fait en supposant que les éléments de la matrice sont stockés dans la colonne majeure. commande.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Soit A et B les matrices 3x2 suivantes :
Lorsque A et B sont stockés dans l'ordre des colonnes principales, nous pouvons appeler la routine de copie comme suit :
pip install numpy
Si, toutefois, A et B sont tous deux stockés dans l'ordre des lignes principales, la signature d'appel devient
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Remarquez que, dans ce dernier scénario, nous ne parvenons pas à accéder aux éléments dans un ordre séquentiel dans la boucle la plus interne, car da0 vaut 2 et da1 vaut -5 et de même pour db0 et db1. Au lieu de cela, les "pointeurs" d'index de tableau avancent à plusieurs reprises avant de revenir aux éléments précédents dans la mémoire linéaire, avec ia = {0, 2, 4, 1, 3, 5} et ib identiques. Dans la figure 3, nous montrons l'impact sur les performances de l'accès non séquentiel.
Figure 3 : Comparaison des performances lors de la fourniture de matrices carrées de colonne principale et de ligne principale à copie lorsque copie suppose un accès séquentiel aux éléments selon l'ordre des colonnes majeures. L'axe des X énumère les tailles de matrice croissantes (c'est-à-dire le nombre d'éléments). Tous les taux sont normalisés par rapport aux résultats des colonnes principales pour une taille de matrice correspondante.
À partir de la figure, nous pouvons observer que les performances des colonnes et des lignes principales sont à peu près équivalentes jusqu'à ce que nous opérions sur des matrices carrées ayant plus de 1e5 éléments (M = N = ~ 316). Pour les éléments 1e6 (M = N = ~ 1 000), la fourniture d'une matrice de lignes principales à copier entraîne une diminution des performances de plus de 25 %. Pour les éléments 1e7 (M = N = ~3160), nous observons une diminution des performances supérieure à 85 %. L'impact significatif sur les performances peut être attribué à une diminution de la localité de référence lors de l'utilisation de matrices à lignes principales ayant de grandes tailles de lignes.
Étant donné qu'il est écrit en Fortran, LAPACK assume l'ordre d'accès aux colonnes principales et implémente ses algorithmes en conséquence. Cela présente des problèmes pour les bibliothèques, telles que stdlib, qui non seulement prennent en charge l'ordre des lignes principales, mais en font également leur disposition de mémoire par défaut. Si nous devions simplement porter les implémentations Fortran de LAPACK vers JavaScript, les utilisateurs fournissant des matrices de lignes principales subiraient des impacts négatifs sur les performances dus à un accès non séquentiel.
Pour atténuer les impacts négatifs sur les performances, nous avons emprunté une idée à BLIS, une bibliothèque de type BLAS prenant en charge les dispositions de mémoire principales en lignes et en colonnes dans les routines BLAS, et avons décidé de créer des implémentations LAPACK modifiées lors du portage des routines de Fortran vers JavaScript et C qui prend explicitement en charge les dispositions de mémoire en colonnes et en lignes principales via des paramètres de foulée distincts pour chaque dimension. Pour certaines implémentations, comme dlacpy, qui est similaire à la fonction de copie définie ci-dessus, l'incorporation de foulées séparées et indépendantes est simple, impliquant souvent des astuces de foulée et un échange de boucles, mais, pour d'autres, les modifications se sont avérées beaucoup moins simples en raison de gestion matricielle spécialisée, modèles d'accès variables et paramétrage combinatoire.
Les routines LAPACK fonctionnent principalement sur des matrices stockées dans la mémoire linéaire et dont les éléments sont accessibles en fonction de dimensions spécifiées et de la foulée de la dimension principale (c'est-à-dire la première). Les dimensions spécifient respectivement le nombre d'éléments dans chaque ligne et colonne. La foulée spécifie combien d'éléments de la mémoire linéaire doivent être ignorés afin d'accéder à l'élément suivant d'une ligne. LAPACK suppose que les éléments appartenant à la même colonne sont toujours contigus (c'est-à-dire adjacents en mémoire linéaire). La figure 4 fournit une représentation visuelle des conventions LAPACK (en particulier, les schémas (a) et (b)).
Figure 4 : Schémas illustrant la généralisation des conventions des tableaux striés LAPACK aux tableaux striés non contigus. a) Une matrice contiguë 5 x 5 stockée dans l'ordre des colonnes principales. b) Une sous-matrice 3 x 3 non contiguë stockée dans l'ordre des colonnes principales. Les sous-matrices peuvent être exploitées dans LAPACK en fournissant un pointeur vers le premier élément indexé et en spécifiant la foulée de la dimension principale (c'est-à-dire la première). Dans ce cas, le pas de la dimension principale est de cinq, même s'il n'y a que trois éléments par colonne, en raison de la non-contiguïté des éléments de la sous-matrice dans la mémoire linéaire lorsqu'ils sont stockés dans le cadre d'une matrice plus grande. Dans LAPACK, la foulée de la dimension finale (c'est-à-dire la deuxième) est toujours supposée être l'unité. c) Une sous-matrice 3 x 3 non contiguë stockée dans l'ordre des colonnes principales ayant des foulées non unitaires et généralisant les conventions de foulée LAPACK aux dimensions de début et de fin. Cette généralisation sous-tend les tableaux multidimensionnels de stdlib (également appelés « ndarrays »).
Les bibliothèques, telles que NumPy et stdlib, généralisent les conventions de tableaux stridés de LAPACK pour prendre en charge
La prise en charge des foulées non unitaires dans la dernière dimension garantit la prise en charge de la création O(1) de vues non contiguës de la mémoire linéaire sans nécessiter de mouvement explicite des données. Ces vues sont souvent appelées « tranches ». À titre d'exemple, considérons l'extrait de code suivant qui crée de telles vues à l'aide des API fournies par stdlib.
pip install numpy
Sans prise en charge des foulées non unitaires dans la dernière dimension, renvoyer une vue à partir de l'expression x['::2,::2'] ne serait pas possible, car il faudrait copier les éléments sélectionnés dans un nouveau linéaire. mémoire tampon afin d'assurer la contiguïté.
Figure 5 : Schémas illustrant l'utilisation de la manipulation de la foulée pour créer des vues inversées et pivotées d'éléments matriciels stockés dans la mémoire linéaire. Pour tous les sous-schémas, les foulées sont répertoriées sous la forme [trailing_dimension, leads_dimension]. Implicite pour chaque schéma est un "offset", qui indique l'indice du premier élément indexé tel que, pour une matrice A, l'élément Aij est résolu en fonction du décalage i⋅strides[1] j⋅strides[0]. a) Étant donné une matrice 3 par 3 stockée dans l'ordre des colonnes principales, on peut manipuler les foulées des dimensions de début et de fin pour créer des vues dans lesquelles les éléments de la matrice le long d'un ou plusieurs axes sont accessibles dans l'ordre inverse. b) En utilisant une manipulation de foulée similaire, on peut créer des vues pivotées des éléments de la matrice par rapport à leur disposition dans la mémoire linéaire.
La prise en charge des foulées négatives permet l'inversion O(1) et la rotation des éléments le long d'une ou plusieurs dimensions (voir Figure 5). Par exemple, pour retourner une matrice de haut en bas et de gauche à droite, il suffit d’annuler les foulées. S'appuyant sur l'extrait de code précédent, l'extrait de code suivant montre l'inversion d'éléments autour d'un ou plusieurs axes.
pip install numpy
La discussion sur les foulées négatives implique implicitement la nécessité d'un paramètre "offset" qui indique l'index du premier élément indexé dans la mémoire linéaire. Pour un tableau multidimensionnel foulé A et une liste de foulées s, l'index correspondant à l'élément Aij⋅⋅⋅n peut être résolu selon l'équation
où N est le nombre de dimensions du tableau et sk correspond à la kème foulée.
Dans les routines BLAS et LAPACK prenant en charge les foulées négatives (ce qui n'est pris en charge que lors d'opérations sur des vecteurs foulés (par exemple, voir daxpy ci-dessus)), le décalage d'index est calculé en utilisant une logique similaire à l'extrait de code suivant :
pip install numpy
où M est le nombre d'éléments vectoriels. Cela suppose implicitement qu'un pointeur de données fourni pointe vers le début de la mémoire linéaire pour un vecteur. Dans les langages prenant en charge les pointeurs, tels que C, afin d'opérer sur une région différente de la mémoire linéaire, on ajuste généralement un pointeur en utilisant l'arithmétique du pointeur avant l'invocation de la fonction, ce qui est relativement peu coûteux et simple, du moins pour le cas unidimensionnel.
Par exemple, en revenant à c_daxpy tel que défini ci-dessus, nous pouvons utiliser l'arithmétique du pointeur pour limiter l'accès aux éléments à cinq éléments dans la mémoire linéaire en commençant aux onzième et seizième éléments (remarque : indexation de base zéro) d'un tableau d'entrée et de sortie, respectivement, comme indiqué dans l'extrait de code suivant.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Cependant, en JavaScript, qui ne prend pas en charge l'arithmétique explicite des pointeurs pour les tampons binaires, il faut explicitement instancier de nouveaux objets tableau typés ayant un décalage d'octets souhaité. Dans l'extrait de code suivant, afin d'obtenir les mêmes résultats que l'exemple C ci-dessus, nous devons résoudre un constructeur de tableau typé, calculer un nouveau décalage d'octets, calculer une nouvelle longueur de tableau typé et créer une nouvelle instance de tableau typé.
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
Pour les tableaux de grande taille, le coût de l'instanciation d'un tableau typé est négligeable par rapport au temps passé à accéder et à fonctionner sur des éléments individuels du tableau ; cependant, pour les tableaux de plus petite taille, l'instanciation d'objets peut avoir un impact significatif sur les performances.
En conséquence, afin d'éviter des impacts négatifs sur les performances d'instanciation d'objet, stdlib dissocie le tampon de données d'un ndarray de l'emplacement de l'élément tampon correspondant au début d'une vue ndarray. Cela permet aux expressions slice x[2:,3:] et x[3:,1:] de renvoyer de nouvelles vues ndarray sans avoir besoin d'instancier de nouvelles instances de tampon, comme démontré dans l'extrait de code suivant.
// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';
// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;
// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);
En conséquence du découplage d'un tampon de données du début d'une vue ndarray, nous avons également cherché à éviter d'avoir à instancier de nouvelles instances de tableau typées lors de l'appel à des routines LAPACK avec des données ndarray. Cela signifiait créer des signatures API LAPACK modifiées prenant en charge des paramètres de décalage explicites pour tous les vecteurs et matrices striés.
Pour simplifier, revenons à l'implémentation JavaScript de daxpy, qui a été précédemment définie ci-dessus.
pip install numpy
Comme démontré dans l'extrait de code suivant, nous pouvons modifier la signature et l'implémentation ci-dessus de telle sorte que la responsabilité de la résolution du premier élément indexé soit transférée au consommateur de l'API.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Pour les ndarrays, la résolution se produit lors de l'instanciation de ndarray, ce qui fait de l'invocation de daxpy_ndarray avec les données ndarray un transfert simple des métadonnées ndarray associées. Ceci est démontré dans l’extrait de code suivant.
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
Semblable à BLIS, nous avons vu la valeur des signatures API LAPACK conventionnelles (par exemple, pour la compatibilité ascendante) et des signatures API modifiées (par exemple, pour minimiser les impacts négatifs sur les performances), et nous avons donc décidé d'un plan pour fournir à la fois des signatures conventionnelles et API modifiées pour chaque routine LAPACK. Pour minimiser la duplication de code, nous avons cherché à implémenter une implémentation de « base » commune de niveau inférieur qui pourrait ensuite être encapsulée par des API de niveau supérieur. Bien que les changements apportés à la routine BLAS daxpy présentés ci-dessus puissent sembler relativement simples, la transformation d'une routine LAPACK conventionnelle et de son comportement attendu en une implémentation généralisée l'était souvent beaucoup moins.
Assez de défis ! À quoi ressemble un produit final ?!
Bouclons la boucle et revenons à dlaswp, une routine LAPACK permettant d'effectuer une série d'échanges de lignes sur une matrice d'entrée selon une liste d'indices pivots. L'extrait de code suivant montre l'implémentation de référence de LAPACK Fortran.
// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';
// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;
// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);
Pour faciliter l'interface avec l'implémentation Fortran à partir de C, LAPACK fournit une interface C à deux niveaux appelée LAPACKE, qui enveloppe les implémentations Fortran et s'adapte aux matrices d'entrée et de sortie principales en ligne et en colonne. L'interface de niveau intermédiaire pour dlaswp est présentée dans l'extrait de code suivant.
void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
int ix;
int iy;
int i;
if (N <= 0) {
return;
}
if (alpha == 0.0) {
return;
}
if (strideX < 0) {
ix = (1-N) * strideX;
} else {
ix = 0;
}
if (strideY < 0) {
iy = (1-N) * strideY;
} else {
iy = 0;
}
for (i = 0; i < N; i++) {
Y[iy] += alpha * X[ix];
ix += strideX;
iy += strideY;
}
return;
}
Lorsqu'il est appelé avec une matrice de colonne majeure a, le wrapper LAPACKE_dlaswp_work transmet simplement les arguments fournis à l'implémentation Fortran. Cependant, lorsqu'il est appelé avec une matrice de ligne majeure a, le wrapper doit allouer de la mémoire, transposer et copier explicitement a dans une matrice temporaire a_t, recalculer la foulée de la dimension principale, invoquer dlaswp avec a_t, transposer et copier les résultats stockés dans a_t. à a, et enfin libérer la mémoire allouée. Cela représente une quantité de travail considérable et est courant dans la plupart des routines LAPACK.
L'extrait de code suivant montre l'implémentation de référence LAPACK portée en JavaScript, avec prise en charge des foulées de dimension de début et de fin, des décalages d'index et d'un vecteur foulé contenant des indices de pivotement.
const binary = new UintArray([...]); const mod = new WebAssembly.Module(binary);
Pour fournir une API ayant un comportement cohérent avec le LAPACK conventionnel, j'ai ensuite enveloppé l'implémentation ci-dessus et adapté les arguments d'entrée à l'implémentation "de base", comme indiqué dans l'extrait de code suivant.
pip install numpy
J'ai ensuite écrit un wrapper séparé mais similaire qui fournit un mappage API plus directement aux tableaux multidimensionnels de stdlib et qui effectue une gestion spéciale lorsque la direction dans laquelle appliquer les pivots est négative, comme indiqué dans l'extrait de code suivant.
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
Quelques points à noter :
Pour les applications côté serveur souhaitant exploiter des bibliothèques optimisées pour le matériel, telles qu'OpenBLAS, nous fournissons des wrappers distincts qui adaptent les arguments de signature généralisés à leurs équivalents API optimisés. Dans ce contexte, au moins pour les tableaux suffisamment grands, la création de copies temporaires peut valoir la peine.
Malgré les défis, les revers imprévus et les multiples itérations de conception, je suis heureux d'annoncer qu'en plus du dlaswp ci-dessus, j'ai pu ouvrir 35 PR ajoutant la prise en charge de diverses routines LAPACK et utilitaires associés. Evidemment pas tout à fait 1 700 routines, mais un bon début ! :)
Néanmoins, l’avenir est prometteur et nous sommes très enthousiasmés par ce travail. Il y a encore beaucoup de place à l’amélioration et à des recherches et développements supplémentaires. Nous tenons particulièrement à
Bien que mon stage chez Quansight Labs soit terminé, mon plan est de continuer à ajouter des packages et à poursuivre cet effort. Compte tenu de l'immense potentiel et de l'importance fondamentale de LAPACK, nous serions ravis de voir cette initiative visant à amener LAPACK sur le Web continuer à se développer, donc si vous souhaitez contribuer à faire avancer ce projet, n'hésitez pas à nous contacter ! Et si vous souhaitez parrainer le développement, les gens de Quansight se feront un plaisir de discuter.
Et sur ce, je voudrais remercier Quansight pour m'avoir offert cette opportunité de stage. Je me sens incroyablement chanceuse d’avoir appris autant de choses. Être stagiaire chez Quansight était depuis longtemps un de mes rêves, et je suis très reconnaissant de l'avoir réalisé. Je tiens à remercier tout particulièrement Athan Reines et Melissa Mendonça, qui est un mentor extraordinaire et une personne formidable à tous points de vue ! Et merci à tous les développeurs principaux de stdlib et à tous les autres membres de Quansight de m'avoir aidé, petits et grands, tout au long du processus.
Bravo !
stdlib est un projet logiciel open source dédié à fournir une suite complète de bibliothèques robustes et hautes performances pour accélérer le développement de votre projet et vous offrir une tranquillité d'esprit en sachant que vous dépendez d'un logiciel de haute qualité conçu par des experts.
Si vous avez apprécié ce post, donnez-nous une étoile ? sur GitHub et envisagez de soutenir financièrement le projet. Vos contributions et votre soutien continu contribuent à assurer le succès à long terme du projet et sont grandement appréciés !
Ce qui précède est le contenu détaillé de. pour plus d'informations, suivez d'autres articles connexes sur le site Web de PHP en chinois!
Méthodes de cryptage du stockage des données
Comment lire des fichiers texte en HTML
Quelles sont les méthodes de transfert de fichiers Java ?
La différence entre le langage C et Python
Utilisation de base de l'instruction insert
Les systèmes logiciels informatiques comprennent
Que diriez-vous du processeur n5095
Quels sont les logiciels de dessin ?
![[Web front-end] Démarrage rapide de Node.js](https://img.php.cn/upload/course/000/000/067/662b5d34ba7c0227.png)
