この投稿はもともと Quansight Labs ブログで公開されたもので、Quansight の許可を得て変更され、ここに再公開されています。
Web アプリケーションは、高性能科学計算と AI 対応のエンドユーザー エクスペリエンスの新たなフロンティアとして急速に台頭しています。 ML/AI 革命を支えているのは、線形方程式とそのベクトル空間および行列による表現に関する数学の一分野である線形代数です。 LAPACK ("Linear Algebra Package") は数値線形代数の基本的なソフトウェア ライブラリであり、一般的な行列演算の堅牢で実戦テストされた実装を提供します。 。 LAPACK はほとんどの数値計算プログラミング言語およびライブラリの基本コンポーネントであるにもかかわらず、Web 固有の制約に合わせた包括的で高品質な LAPACK 実装はまだ実現していません。それは…今までです。
今年の初め、私は幸運なことに、Quansight の公益部門であり、科学的な Python エコシステムのリーダーである Quansight Labs で夏季インターンとして働くことができました。インターンシップ中に、私は stdlib に初期 LAPACK サポートを追加することに取り組みました。stdlib は、C と JavaScript で記述され、Web ブラウザーや他の Web ネイティブ環境 (Node.js や Deno など) で使用するために最適化された科学計算用の基本ライブラリです。このブログ投稿では、私のこれまでの道のり、予期されたおよび予期せぬ (!) 課題、そして今後の道のりについて説明します。私の希望は、少しの幸運があれば、この成果が Web ブラウザーを数値計算と機械学習のための一流の環境にするための重要な構成要素を提供し、より強力な AI 対応 Web アプリケーションの未来の前兆となることです。
面白いと思いませんか？行きましょう!
LAPACK に詳しいこのブログの読者は、Web テクノロジーの荒々しい世界には詳しくない可能性があります。数値計算および科学計算の世界に属し、科学 Python エコシステムに精通している人にとって、stdlib を NumPy と SciPy の型をとったオープンソースの科学計算ライブラリとして考えるのが最も簡単です。多次元配列データ構造と、数学、統計、線形代数に関連するルーチンを提供しますが、主なスクリプト言語として Python ではなく JavaScript を使用します。そのため、stdlib は Web エコシステムとそのアプリケーション開発パラダイムに重点を置いています。この焦点には、いくつかの興味深い設計とプロジェクト アーキテクチャの決定が必要であり、数値計算用に設計された従来のライブラリと比較した場合、stdlib がかなりユニークになります。
NumPy を例に挙げると、NumPy は単一のモノリシック ライブラリであり、OpenBLAS などのオプションのサードパーティ依存関係を除くすべてのコンポーネントが単一の分割不可能なユニットを形成します。すべての NumPy をインストールせずに、単純に配列操作用の NumPy ルーチンをインストールすることはできません。 NumPy の ndarray オブジェクトとその操作ルーチンのいくつかだけを必要とするアプリケーションをデプロイする場合、NumPy をすべてインストールしてバンドルすることは、かなりの量の「デッド コード」を含めることを意味します。 Web 開発用語で言うと、NumPy は「ツリーシェイク可能」ではないと言えます。通常の NumPy インストールの場合、これは少なくとも 30 MB のディスク容量を意味し、すべてのデバッグ ステートメントを除外するカスタマイズされたビルドの場合は少なくとも 15 MB のディスク容量を意味します。 SciPy の場合、これらの数値はそれぞれ 130MB と 50MB にまで膨れ上がる可能性があります。言うまでもなく、Web アプリケーションにほんの少数の機能を提供するために 15MB のライブラリを出荷することは、特にネットワーク接続が不十分なデバイスやメモリの制約があるデバイスに Web アプリケーションを展開する必要がある開発者にとっては、初心者には適していません。
Web アプリケーション開発特有の制約を考慮して、stdlib は設計にボトムアップ アプローチを採用しており、コードベースの無関係な未使用部分とは独立して機能のすべてのユニットをインストールして使用できます。 stdlib は、分解可能なソフトウェア アーキテクチャと徹底的なモジュール性を採用することで、必要な API セットとその明示的な依存関係を超える余分なコードをほとんどまたはまったく使用せずに、必要なものを正確にインストールして使用できる機能をユーザーに提供します。これにより、メモリ フットプリント、バンドルの削減が保証されます。サイズを拡大し、展開を高速化します。
例として、2 つの行列スタック (つまり、3 次元立方体の 2 次元スライス) を操作していて、1 つおきのスライスを選択して、一般的な BLAS 演算 y = a * x を実行するとします。ここで、x と y は ndarray で、a はスカラー定数です。 NumPy を使用してこれを行うには、まず NumPy
をすべてインストールします。
pip install numpy
その後、さまざまな操作を実行します
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
stdlib を使用すると、プロジェクトをモノリシック ライブラリとしてインストールできることに加えて、さまざまな機能ユニットを個別のパッケージとしてインストールできます
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
その後、さまざまな操作を実行します
// 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,:,:']);
重要なのは、stdlib の 4,000 を超えるパッケージのいずれかを個別にインストールできるだけでなく、関連する GitHub リポジトリをフォークすることで、それらのパッケージのいずれかを修正、改善、リミックスすることもできることです (例: @stdlib/ndarray-fancy を参照) ）。 stdlib では、抽象化層と依存関係ツリーを明示的に定義することにより、アプリケーションに適切な抽象化層を自由に選択できます。ある意味、これは単純なアイデアであり、従来の科学ソフトウェア ライブラリの設計に慣れている人にとっては、おそらく型破りなアイデアですが、Web プラットフォームと緊密に統合されると、強力な結果をもたらし、刺激的な新しい可能性が生まれます。
なるほど、興味が湧いたのかもしれません。 stdlib は興味深いようです。しかし、これは Web ブラウザの LAPACK とどのような関係があるのでしょうか?さて、この夏の私たちの目標の 1 つは、LAPACK を Web に導入する際に、stdlib の精神 (1 つのことを実行し、1 つのことをうまく実行する、小規模で範囲が狭いパッケージ) を適用することでした。
でも、ちょっと待ってください！それは極端な取り組みです。 LAPACK には約 1,700 のルーチンがあり、その 10% を妥当な期間内に実装することは大きな課題です。 LAPACK を、C、Go、Rust などのプログラミング言語用の移植可能なコンパイル ターゲットである WebAssembly にコンパイルするだけで、Web 上でのデプロイメントが可能になり、それで終わりにする方がよいのではないでしょうか?
残念ながら、このアプローチにはいくつかの問題があります。
最後の点を説明するために、BLAS ルーチン daxpy に戻りましょう。このルーチンは、演算 y = a*x y を実行します。ここで、x と y はストライドされたベクトル、a はスカラー定数です。 C で実装された場合、基本的な実装は次のコード スニペットのようになります。
pip install numpy
WebAssembly にコンパイルし、WebAssembly バイナリを Web アプリケーションにロードした後、JavaScript から c_daxpy ルーチンを呼び出す前に、一連の手順を実行する必要があります。まず、新しい 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,:,:]
次に、モジュール メモリを定義し、新しい WebAssembly モジュール インスタンスを作成する必要があります。
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
モジュール インスタンスを作成した後、エクスポートされた BLAS ルーチンを呼び出すことができます。ただし、データがモジュール メモリの外部で定義されている場合は、まずそのデータをメモリ インスタンスにコピーする必要があり、常にリトル エンディアンのバイト オーダーでコピーする必要があります。
// 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,:,:']);
データがモジュール メモリに書き込まれたので、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;
}
そして最後に、視覚化やさらなる分析のために、WebAssembly メモリ「ポインター」(つまり、バイト オフセット) をサポートしていないダウンストリーム ライブラリ (D3 など) に結果を渡す必要がある場合は、モジュールからデータをコピーする必要があります。メモリを元の出力配列に戻します。
const binary = new UintArray([...]); const mod = new WebAssembly.Module(binary);
y = a*x y を計算するだけでも大変な作業です。対照的に、次のコード スニペットのように見える単純な JavaScript 実装と比較してください。
// 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
}
});
JavaScript 実装を使用すると、上記の WebAssembly の例で必要なデータ移動を行わずに、外部定義されたデータを使用して daxpy を直接呼び出すことができます。
// 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);
}
少なくともこの場合、WebAssembly のアプローチは人間工学的ではないだけでなく、必要なデータ移動を考えると予想されるとおり、次の図に示すようにパフォーマンスにも悪影響があります。
図 1: 配列の長さ (x 軸) を増加させた場合の、BLAS ルーチン daxpy の stdlib の C、JavaScript、および WebAssembly (Wasm) 実装のパフォーマンスの比較。 Wasm (コピー) ベンチマークでは、入出力データが Wasm メモリとの間でコピーされるため、パフォーマンスが低下します。
上の図では、配列の長さを増加させた場合の BLAS ルーチン daxpy の stdlib の C、JavaScript、および WebAssembly (Wasm) 実装のパフォーマンス比較を x 軸に沿って列挙して表示しています。 y 軸は、ベースライン C 実装に対する正規化されたレートを示します。 Wasm ベンチマークでは、入出力データは WebAssembly モジュール メモリ内で直接割り当ておよび操作されます。Wasm (コピー) ベンチマークでは、上で説明したように、入出力データは WebAssembly モジュール メモリとの間でコピーされます。チャートから次のことがわかります:
全体的に、WebAssembly はパフォーマンスを向上させることができます。ただし、このテクノロジーは特効薬ではないため、望ましい利益を実現するには慎重に使用する必要があります。また、優れたパフォーマンスを提供する場合でも、そのような利益は、複雑さの増加、潜在的にバンドル サイズが大きくなる可能性、およびより複雑なツールチェーンによるコストとバランスを取る必要があります。多くのアプリケーションでは、プレーンな JavaScript 実装で問題なく機能します。
LAPACK 全体を WebAssembly にコンパイルするだけで訴訟を起こして終わりにしたので、あとはどうなるでしょうか?そうですね、stdlib の理念を受け入れるつもりなら、抜本的なモジュール化が必要になります。
ラディカルなモジュール性を採用するということは、何が最善であるかは高度にコンテキストに依存しており、ユーザー アプリケーションのニーズと制約に応じて、開発者は適切な抽象化を選択する柔軟性が必要であることを認識することです。開発者が Node.js アプリケーションを作成している場合、優れたパフォーマンスを実現するために、OpenBLAS、Intel MKL、Apple Accelerate などのハードウェアに最適化されたライブラリにバインドする必要がある場合があります。開発者が少数の数値ルーチンを必要とする Web アプリケーションをデプロイしている場合、JavaScript がその作業に適したツールである可能性があります。また、開発者が大規模でリソースを大量に消費する WebAssembly アプリケーション (画像編集やゲーム エンジンなど) に取り組んでいる場合、より大きなアプリケーションの一部として個々のルーチンを簡単にコンパイルできることが最も重要になります。つまり、根本的にモジュール化された LAPACK が必要なのです。
私の使命は、そのような取り組みの基礎を築き、問題点を解決してギャップを見つけ、できればウェブ上の高性能線形代数に数歩近づけることでした。しかし、根本的なモジュール化とはどのようなものでしょうか?すべては、機能の基本単位である パッケージ から始まります。
stdlib 内のすべてのパッケージは、独自のスタンドアロンのものであり、同時ローカライズされたテスト、ベンチマーク、サンプル、ドキュメント、ビルド ファイル、および関連するメタデータ (依存関係の列挙を含む) を含み、外部との明確な API サーフェスを定義します。 。 LAPACK サポートを stdlib に追加するには、次の構造を持つ LAPACK ルーチンごとに個別のスタンドアロン パッケージを作成することを意味します。
pip install numpy
簡単に言うと、
stdlib の厳しいドキュメントとテスト要件を考慮すると、各ルーチンのサポートを追加するのはかなりの作業量ですが、最終結果は堅牢で高品質で、最も重要なことに、科学計算の基盤として機能するのに適したモジュール式コードになります。現代のウェブ上で。しかし、前置きはこれで十分です！本題に入りましょう!
BLAS サポートを stdlib に追加した以前の取り組みに基づいて、LAPACK サポートを追加するときにも同様の多段階アプローチに従うことにしました。このアプローチでは、最初に JavaScript 実装とそれに関連するテストとドキュメントを優先し、次にテストとドキュメントが存在するようになります。 、C および Fortran 実装、およびハードウェア最適化ライブラリへの関連するネイティブ バインディングをバックフィルします。このアプローチにより、構築ツールチェーンとパフォーマンスの雑草に飛び込む前に、いわばボード上でいくつかの初期のポイントを設定し、ユーザーの前に API を迅速に提供し、堅牢なテスト手順とベンチマークを確立し、ツールと自動化の潜在的な手段を調査することができます。最適化。しかし、どこから始めればよいのでしょうか?
どの LAPACK ルーチンを最初にターゲットにするかを決定するために、LAPACK の Fortran ソース コードを解析してコール グラフを生成しました。これにより、各 LAPACK ルーチンの依存関係ツリーを推測できるようになりました。次に、グラフを手にしてトポロジカル ソートを実行しました。これにより、依存関係のないルーチン、および必然的に他のルーチンの構成要素となるルーチンを特定するのに役立ちました。特定の高レベルのルーチンを選択して逆方向に作業する深さ優先のアプローチでは、特定の機能を実現できますが、そのようなアプローチでは、ますます複雑になるルーチンを実装しようとして行き詰まってしまう可能性があります。グラフの「葉」に焦点を当てることで、一般的に使用されるルーチン (つまり、度数 が高いルーチン) に優先順位を付けることができ、後から複数の高レベルのルーチンを提供できるようにすることで効果を最大化できます。私の取り組み、または他の貢献者による取り組み。
計画を立てたので、仕事に取り掛かることに興奮していました。最初のルーチンとして、dalaswp を選択しました。これは、提供されたピボット インデックスのリストに従って一般的な長方形行列上で一連の行交換を実行し、LAPACK の LU 分解ルーチンの重要な構成要素です。そしてそこから私の挑戦が始まりました...
Quansight Labs でのインターンシップに参加する前、私は、LLVM 上に構築された最新の対話型 Fortran コンパイラーである LFortran に定期的に貢献しており (そして今でも!)、自分の Fortran スキルにはかなり自信を持っていました。しかし、私の最初の課題の 1 つは、現在「レガシー」と考えられている Fortran コードを単純に理解することでした。以下に最初の 3 つのハードルを取り上げます。
LAPACK はもともと FORTRAN 77 (F77) で書かれました。ライブラリはバージョン 3.2 (2008) で Fortran 90 に移行されましたが、従来の規約がリファレンス実装にまだ残っています。これらの規則の中で最も顕著なものの 1 つは書式設定です。
F77 プログラムを作成する開発者は、パンチカードから継承した固定形式のレイアウトを使用してプログラムを作成しました。このレイアウトには、文字列の使用に関する厳しい要件がありました:
Fortran 90 は、列と行の長さの制限を取り除き、 ! に落ち着いたフリーフォームレイアウトを導入しました。コメント文字として。次のコード スニペットは、LAPACK ルーチン dlacpy のリファレンス実装を示しています。
pip install numpy
次のコード スニペットは同じルーチンを示していますが、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,:,:]
ご覧のとおり、列制限を削除し、指定子をすべて大文字で記述するという F77 規則から離れることにより、最新の Fortran コードはより一貫性が高く、読みやすくなりました。
LAPACK ルーチンのもう 1 つの一般的な方法は、ラベル付きの制御構造の使用です。たとえば、ラベル 10 が対応する CONTINUE.
と一致する必要がある次のコード スニペットを考えてみましょう。
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
Fortran 90 では、do ループを終了するために end do を使用できるようにすることで、この慣行の必要性がなくなり、コードの可読性が向上しました。この変更は、上記で提供されている dlacpy のフリーフォーム バージョンに示されています。
さまざまなサイズの配列を柔軟に処理できるようにするために、LAPACK ルーチンは通常、想定されたサイズの配列を操作します。上記の dlacpy ルーチンでは、入力行列 A は、式 A(LDA, *) に従って仮定されたサイズを持つ 2 次元配列として宣言されます。この式は、A に LDA 行数があることを宣言し、* をプレースホルダーとして使用して、2 番目の次元のサイズが呼び出し側プログラムによって決定されることを示します。
サイズ引継ぎ配列を使用すると、コンパイラが未指定の次元で境界チェックを実行できなくなることが起こります。したがって、現在のベスト プラクティスは、境界外のメモリ アクセスを防ぐために、明示的なインターフェイスと想定形状配列 (例: A(LDA,:)) を使用することです。これは、部分行列を他の関数に渡す必要がある場合、形状を仮定した配列の使用は問題となる可能性があることを示しています。そうするためにはスライスが必要となり、コンパイラーが配列データの内部コピーを作成することがよくあります。
言うまでもなく、LAPACK の慣例に適応し、LAPACK の考え方を取り入れるまでには時間がかかりました。しかし、純粋主義者の私としては、どうせルーチンを移植するのであれば、コードの可読性と将来のメンテナンスを改善するために、移植に成功したルーチンを少なくともより現代的な時代に導入したいと考えていました。そこで、stdlib のメンテナと話し合った後、ルーチンを Fortran 95 に移行することにしました。これは、最新かつ最高の Fortran バージョンではありませんが、元の実装のルック アンド フィールの維持と、(これで十分) 下位互換性があり、新しい構文機能を活用できます。
LAPACK サポートを追加するボトムアップ アプローチを追求する場合の問題の 1 つは、下位レベルのユーティリティ ルーチンの明示的な単体テストが LAPACK に存在しないことが多いことです。 LAPACK のテスト スイートは主に階層的なテスト哲学を採用しており、上位レベルのルーチンをテストすることで、依存する下位レベルのルーチンがワークフロー全体の一部として正しく機能していることを確認します。すべてのルーチンにテストを追加すると、保守の負担と LAPACK のテスト フレームワークの複雑さが増大する可能性があるため、下位レベルのルーチンの単体テストよりも統合テストに重点を置くのが合理的であると主張する人もいます。これは、以前のテストに簡単に依存できないことを意味します。単体テストには技術が必要であり、下位レベルのルーチンごとに包括的なスタンドアロン単体テストを自分たちで考え出す必要があります。
テスト カバレッジと同様に、LAPACK 自体の外で、下位レベルのルーチンの使用を示す実際の文書化された例を見つけることは困難でした。 LAPACK ルーチンの前には常に、入力引数と考えられる戻り値の説明を提供するドキュメント コメントが付いていますが、コード例がなければ、特に特殊な行列を扱う場合、予想される入力値と出力値を視覚化して把握するのは困難な場合があります。単体テストや文書化されたサンプルが存在しないからといって世界が終わるわけではありませんが、stdlib に LAPACK サポートを追加するのは予想以上に大変な作業になることを意味します。ベンチマーク、テスト、サンプル、ドキュメントを作成するには、より多くの時間と労力が必要となり、インターンシップ中に実装できるルーチンの数が制限される可能性がありました。
行列要素を線形メモリに格納する場合、列を連続して格納するか、行を連続して格納するかの 2 つの選択肢があります (図 2 を参照)。前者のメモリ レイアウトは 列優先 順序と呼ばれ、後者は 行優先 順序と呼ばれます。
図 2: (a) 列優先 (Fortran スタイル) または (b) 行優先 (C スタイル) のいずれかの順序で行列要素を線形メモリに保存する様子を示す回路図。どのレイアウトを使用するかの選択は、主に慣例によるものです。
どのレイアウトを使用するかの選択は、主に慣例の問題です。たとえば、Fortran は要素を列優先の順序で格納し、C は要素を行優先の順序で格納します。 NumPy や stdlib などの高レベルのライブラリは、列優先と行優先の両方をサポートしているため、配列の作成中に多次元配列のレイアウトを構成できます。
pip install numpy
どちらのメモリ レイアウトも本質的に他方より優れているわけではありませんが、最適なパフォーマンスを確保するには、基礎となるストレージ モデルの規則に従ってシーケンシャル アクセスを保証するようにデータを配置することが重要です。最新の CPU は、シーケンシャル データを非シーケンシャル データよりも効率的に処理できます。これは主に CPU キャッシュによるもので、これにより参照の空間的局所性が利用されます。
順次要素アクセスと非順次要素アクセスのパフォーマンスへの影響を示すために、MxN 行列 A から別の MxN 行列 B にすべての要素をコピーし、行列要素が列優先で格納されていると仮定してコピーする次の関数を考えてみましょう。注文します。
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
A と B を次の 3x2 行列とします:
A と B の両方が列優先の順序で格納されている場合、次のようにコピー ルーチンを呼び出すことができます。
pip install numpy
ただし、A と B が両方とも行優先の順序で格納されている場合、呼び出し署名は次のように変更されます。
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
後者のシナリオでは、da0 が 2、da1 が -5 であり、db0 と db1 についても同様であるため、最も内側のループ内の要素に順番にアクセスできないことに注意してください。代わりに、配列インデックスの「ポインター」は、ia = {0, 2, 4, 1, 3, 5} と ib が同じ場合に、線形メモリ内の以前の要素に戻る前に繰り返し先へスキップします。図 3 は、非順次アクセスによるパフォーマンスへの影響を示しています。
図 3: copy が列優先の順序に従って要素に順次アクセスすることを想定している場合、正方列優先行列と行優先行列を copy に提供した場合のパフォーマンスの比較。 X 軸は、増加する行列のサイズ (つまり、要素の数) を列挙します。すべてのレートは、対応する行列サイズの列優先の結果に対して正規化されます。
図から、要素数が 1e5 (M = N = ~316) を超える正方行列を操作するまでは、列優先のパフォーマンスと行優先のパフォーマンスがほぼ同等であることがわかります。 1e6 要素 (M = N = ~1000) の場合、コピーする行優先の行列を指定すると、パフォーマンスが 25% 以上低下します。 1e7 要素 (M = N = ~3160) では、85% を超えるパフォーマンスの低下が観察されます。パフォーマンスに大きな影響を与えるのは、行サイズが大きい行優先の行列を操作する場合の参照の局所性の低下に起因する可能性があります。
LAPACK は Fortran で記述されているため、列優先のアクセス順序を想定し、それに応じてアルゴリズムを実装します。これは、行優先の順序をサポートするだけでなく、それをデフォルトのメモリ レイアウトにする stdlib などのライブラリに問題を引き起こします。 LAPACK の Fortran 実装を JavaScript に単純に移植した場合、行優先の行列を提供するユーザーは、非順次アクセスに起因するパフォーマンスへの悪影響を経験することになります。
パフォーマンスへの悪影響を軽減するために、BLAS ルーチンで行優先メモリ レイアウトと列優先メモリ レイアウトの両方をサポートする BLAS に似たライブラリである BLIS からアイデアを借用し、ルーチンを Fortran から JavaScript に移植するときに修正された LAPACK 実装を作成することにしました。 C では、次元ごとに個別のストライド パラメーターを使用して、列優先と行優先の両方のメモリ レイアウトに明示的に対応します。上記で定義したコピー関数に似た dlacpy などの一部の実装では、別個の独立したストライドを組み込むのは簡単で、多くの場合、ストライド トリックやループ交換が必要になりますが、他の実装では、次の理由により、変更がはるかに簡単ではないことが判明しました。特殊な行列処理、さまざまなアクセス パターン、および組み合わせパラメータ化。
LAPACK ルーチンは主に線形メモリに格納された行列を操作し、その要素は指定された次元と先頭 (つまり、最初の) 次元のストライドに従ってアクセスされます。次元は、それぞれの行と列の要素の数を指定します。ストライドは、行の次の要素にアクセスするためにリニア メモリ内の要素をいくつスキップする必要があるかを指定します。 LAPACK は、同じ列に属する要素が常に連続している (つまり、線形メモリ内で隣接している) ことを前提としています。図 4 は、LAPACK 規約を視覚的に表したものです (具体的には、図 (a) および (b))。
図 4: LAPACK ストライド配列規則を非連続ストライド配列に一般化した図。 a) 列優先の順序で格納された 5 行 5 列の連続行列。 b) 列優先の順序で格納された 3 行 3 列の不連続部分行列。 LAPACK では、最初のインデックス付き要素へのポインターを提供し、先頭 (つまり、最初の) 次元のストライドを指定することで部分行列を操作できます。この場合、列ごとに要素が 3 つしかないにもかかわらず、先頭次元のストライドは 5 になります。これは、より大きな行列の一部として格納された場合、線形メモリ内のサブ行列要素が連続していないためです。 LAPACK では、後続 (つまり 2 番目) 次元のストライドは常に 1 であると想定されます。 c) 非単位ストライドを持ち、先頭次元と末尾次元の両方に LAPACK ストライド規約を一般化した、列優先の順序で格納された 3 行 3 列の不連続部分行列。この一般化は、stdlib の多次元配列 (「ndarray」とも呼ばれます) を支えています。
NumPy や stdlib などのライブラリは、LAPACK のストライド配列規則を一般化してサポートします
最後の次元での非ユニット ストライドのサポートにより、明示的なデータ移動を必要とせずに、線形メモリの不連続ビューの O(1) 作成が確実にサポートされます。これらのビューは、多くの場合「スライス」と呼ばれます。例として、stdlib.
によって提供される API を使用してそのようなビューを作成する次のコード スニペットを考えてみましょう。
pip install numpy
最後の次元での非ユニット ストライドのサポートがなければ、選択した要素を新しい線形オブジェクトにコピーする必要があるため、式 x['::2,::2'] からビューを返すことはできません。連続性を確保するためにメモリバッファを使用します。
図 5: ストライド操作を使用して線形メモリに格納された行列要素の反転および回転ビューを作成する様子を示す概略図。すべてのサブ回路図について、ストライドは [trailing_dimension, leading_dimension] としてリストされます。各回路図には暗黙的な「オフセット」があり、これは、行列 A の場合、要素 Aij が次のようになる、最初のインデックス付き要素のインデックスを示します。 i⋅strides[1] j⋅strides[0] オフセットに従って解決されます。 a) 列優先の順序で格納された 3 行 3 列の行列を指定すると、先頭次元と末尾の次元のストライドを操作して、1 つ以上の軸に沿った行列要素が逆の順序でアクセスされるビューを作成できます。 b) 同様のストライド操作を使用して、線形メモリ内の配置を基準にして行列要素の回転ビューを作成できます。
負のストライドのサポートにより、1 つ以上の次元に沿った要素の O(1) 反転および回転が可能になります (図 5 を参照)。たとえば、行列を上から下、左から右に反転するには、ストライドを無効にするだけで済みます。前のコード スニペットを基にして、次のコード スニペットは 1 つ以上の軸を中心に要素を反転する方法を示します。
pip install numpy
負のストライドの議論には、線形メモリ内の最初のインデックス付き要素のインデックスを示す「オフセット」パラメータの必要性が暗黙的に含まれています。ストライド多次元配列 A およびストライドのリスト s の場合、要素 Aij⋅⋅⋅n に対応するインデックス方程式に従って解くことができます
ここで、N は配列の次元数であり、sk は k 番目のストライドに対応します。
負のストライドをサポートする BLAS および LAPACK ルーチン (ストライド ベクトルを操作する場合にのみサポートされるもの (たとえば、上記の daxpy を参照)) では、インデックス オフセットは、次のコード スニペットのようなロジックを使用して計算されます。
pip install numpy
ここで、M はベクトル要素の数です。これは、提供されたデータ ポインターがベクトルの線形メモリの先頭を指していることを暗黙的に前提としています。 C などのポインタをサポートする言語では、線形メモリの異なる領域で操作するために、関数呼び出しの前にポインタ演算を使用してポインタを調整するのが一般的です。これは、少なくとも 1 次元の場合には比較的安価で簡単です。
たとえば、上で定義した c_daxpy に戻ると、ポインター演算を使用して、入出力配列の 11 番目と 16 番目の要素 (注: ゼロベースのインデックス付け) から始まる線形メモリ内の 5 つの要素への要素アクセスを制限できます。次のコード スニペットに示すように、それぞれ。
# Import all of NumPy: import numpy as np # Define arrays: x = np.asarray(...) y = np.asarray(...) # Perform operation: y[::2,:,:] += 5.0 * x[::2,:,:]
ただし、JavaScript ではバイナリ バッファに対する明示的なポインタ演算をサポートしていないため、必要なバイト オフセットを持つ新しい型付き配列オブジェクトを明示的にインスタンス化する必要があります。次のコード スニペットでは、上記の C の例と同じ結果を得るために、型付き配列コンストラクターを解決し、新しいバイト オフセットを計算し、新しい型付き配列の長さを計算して、新しい型付き配列インスタンスを作成する必要があります。
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
配列サイズが大きい場合、型付き配列のインスタンス化のコストは、個々の配列要素へのアクセスと操作に費やす時間と比較すると無視できます。ただし、配列サイズが小さい場合、オブジェクトのインスタンス化がパフォーマンスに大きな影響を与える可能性があります。
したがって、オブジェクトのインスタンス化のパフォーマンスへの悪影響を回避するために、stdlib は ndarray ビューの先頭に対応するバッファー要素の位置から ndarray のデータ バッファーを切り離します。これにより、次のコード スニペットに示すように、スライス式 x[2:,3:] および x[3:,1:] は、新しいバッファー インスタンスをインスタンス化する必要がなくせずに、新しい ndarray ビューを返すことができます。
// 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,:,:']);
ndarray ビューの先頭からデータ バッファーを切り離した結果、同様に、ndarray データを使用して LAPACK ルーチンを呼び出すときに、新しい型付き配列インスタンスをインスタンス化する必要を回避しようとしました。これは、すべてのストライドされたベクトルと行列の明示的なオフセット パラメーターをサポートする、修正された LAPACK API シグネチャを作成することを意味します。
わかりやすくするために、上で定義した daxpy の JavaScript 実装に戻りましょう。
pip install numpy
次のコード スニペットに示すように、最初のインデックス付き要素を解決する責任が 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,:,:]
ndarray の場合、解決は ndarray のインスタンス化中に行われるため、ndarray データを使用した daxpy_ndarray の呼び出しは、関連付けられた ndarray メタ データを直接渡すことになります。これは、次のコード スニペットで示されています。
npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy
BLIS と同様に、私たちは従来の LAPACK API シグネチャ (例: 下位互換性) と修正された API シグネチャ (例: パフォーマンスへの悪影響を最小限に抑えるため) の両方に価値があると考え、従来の LAPACK API シグネチャと修正された API シグネチャの両方を提供する計画に落ち着きました。各 LAPACK ルーチンの API を変更しました。コードの重複を最小限に抑えるために、私たちは、上位レベルの API でラップできる共通の下位レベルの「ベース」実装を実装することを目指しました。上記の BLAS ルーチン daxpy の変更は比較的簡単に見えるかもしれませんが、従来の LAPACK ルーチンとその期待される動作の一般化された実装への変換は、それほど簡単ではないことがよくあります。
課題はもうたくさんです!最終的な製品はどのようなものになるでしょうか?!
一周して、これを dlswp に戻しましょう。これは、ピボット インデックスのリストに従って入力行列に対して一連の行交換を実行する LAPACK ルーチンです。次のコード スニペットは、リファレンス 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,:,:']);
C からの Fortran 実装とのインターフェイスを容易にするために、LAPACK は LAPACKE と呼ばれる 2 レベルの C インターフェイスを提供します。これは Fortran 実装をラップし、行と列の両方の入力行列と出力行列に対応します。 dlswp の中間レベルのインターフェイスを次のコード スニペットに示します。
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;
}
列優先行列 a を指定して呼び出された場合、ラッパー LAPACKE_dlaswp_work は、指定された引数を Fortran 実装に渡すだけです。ただし、行優先行列 a を使用して呼び出された場合、ラッパーはメモリを割り当て、明示的に a を転置して一時行列 a_t にコピーし、先頭の次元のストライドを再計算し、a_t で dlaswp を呼び出し、a_t に格納された結果を転置してコピーする必要があります。に変換し、最後に割り当てられたメモリを解放します。これはかなりの作業量であり、ほとんどの LAPACK ルーチンで共通です。
次のコード スニペットは、JavaScript に移植された参照 LAPACK 実装を示しており、先頭と末尾の次元ストライド、インデックス オフセット、ピボット インデックスを含むストライド ベクトルをサポートしています。
const binary = new UintArray([...]); const mod = new WebAssembly.Module(binary);
従来の LAPACK と一貫した動作を持つ API を提供するために、次のコード スニペットに示すように、上記の実装をラップし、入力引数を「基本」実装に適合させました。
pip install numpy
その後、私は別の同様のラッパーを作成しました。これは、stdlib の多次元配列に直接マッピングする 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,:,:]
いくつかの注意点:
OpenBLAS などのハードウェアに最適化されたライブラリを利用したいサーバー側アプリケーションの場合、一般化された署名引数を最適化された同等の API に適合させる個別のラッパーを提供します。このコンテキストでは、少なくとも十分に大きな配列の場合、一時コピーを作成することはオーバーヘッドの価値がある可能性があります。
課題、予期せぬ挫折、設計の複数回の反復にもかかわらず、上記の dlswp に加えて、さまざまな LAPACK ルーチンと関連ユーティリティのサポートを追加する 35 の PR を開くことができたことをうれしく報告します。明らかに 1,700 のルーチンには達していませんが、良いスタートです。 :)
それでも、未来は明るく、私たちはこの仕事にとても興奮しています。改善や追加の研究開発の余地はまだたくさんあります。特に私たちは
に熱心に取り組んでいます。Quansight Labs でのインターンシップは終了しましたが、私の計画は引き続きパッケージを追加し、この取り組みを推進することです。 LAPACK の計り知れない可能性と基本的な重要性を考慮すると、LAPACK を Web に導入するというこの取り組みが今後も成長し続けることを願っています。そのため、この推進を支援することに興味がある場合は、遠慮なくご連絡ください。開発のスポンサーにご興味がございましたら、Quansight の担当者が喜んでお話しさせていただきます。
それでは、このインターンシップの機会を提供してくださったQuansightに感謝したいと思います。とても多くのことを学べて本当に幸運だと感じています。クアンサイトでインターンになることは私の長年の夢であり、それを実現できたことにとても感謝しています。 Athan Reines と、素晴らしい指導者であり、すべてにおいて素晴らしい人である Melissa Mendonça に特別な感謝の意を表したいと思います。そして、途中で大なり小なり私を助けてくれたすべての stdlib コア開発者と Quansight の他の全員に感謝します。
乾杯!
stdlib は、プロジェクトの開発を加速し、専門家が作成した高品質のソフトウェアに依存しているという安心感を与える、堅牢で高性能のライブラリの包括的なスイートを提供することに特化したオープン ソース ソフトウェア プロジェクトです。
この投稿が気に入ったら、スターを付けてください。 GitHub でプロジェクトを財政的に支援することを検討してください。あなたの貢献と継続的なサポートは、プロジェクトの長期的な成功を確実にするのに役立ち、非常に感謝しています!
以上がWebブラウザのLAPACKの詳細内容です。詳細については、PHP 中国語 Web サイトの他の関連記事を参照してください。
![[Web フロントエンド] Node.js クイック スタート](https://img.php.cn/upload/course/000/000/067/662b5d34ba7c0227.png)
