ベクトルと行列(linalg)
linalg パッケージは、ベクタと行列のための numpy スタイルの API を提供します。コンストラクタ、形状操作、要素ごとの算術演算、積、集約、離散微積分(差分と数値微分)、線形代数(行列式、逆行列、連立一次方程式の求解)が含まれます。
JSON ライブラリと同様に、linalg は Lisp ソース(linalg.lisp)として一度だけ実装されています。インタプリタは linalg: 関数が最初に使われたときに定義を遅延ロードし、コンパイルパスはパッケージが参照されたときに定義をプログラムへ継ぎ足します。バックエンドごとのコードは存在しないため、すべての関数はインタプリタ、JVM コンパイラ、WASM Preview 1、WASM コンポーネントで同一に振る舞います。
データ表現
linalg のコンストラクタは packed float 配列 を作ります。これは #d(...) リテラルと同じ、アンボックスな (array double-float) です。ベクタはランク 1 の配列で #d(1.0 2.0 ...) と印字され、行列はランク 2 の配列でネストした #d((...) ...) 形式で印字されます。#d はアンボックスな packed 表現を表すため、その印字結果を読み戻すと packed 配列になります。個々の要素は aref で読み書きでき、プログラムの他の場所で構築された配列 (packed でも一般のボックス配列でも) も linalg 関数に渡せます。より高いランクの配列も扱えます。要素ごとの演算、リダクション、reshape/flatten、array-equal はフラットな行優先順で要素を走査するため任意のランクを受け付けます。一方 matmul は rank 3 以上を最後の 2 軸でスタックし (numpy の np.matmul と同じ規則)、dot/outer/det/inv/solve/trace/transpose は numpy の専用ルーチンと同様、ベクタと行列 (ランク 2 以下) に対して定義されたままです。linalg:from-list / linalg:to-list は配列とリストを相互に変換します。
linalg は速度を優先して浮動小数点で計算します。すべてのコンストラクタと配列を生成する演算はデフォルトで packed double-float 配列を返し(単精度も利用できます。単精度浮動小数点を参照)、linalg:det、linalg:inv、linalg:solve は (numpy と同様に) 浮動小数点で計算されるため、一般の逆行列には通常の丸めが生じ、ほぼ特異な行列式は厳密な 0 ではなく微小値になることがあります。リダクションは numpy と同様に要素の型に従います。packed または float 配列に対するリダクションは double を、素の整数配列 (#(1 2 3) のようなリテラル) に対するリダクションは整数または厳密な比を返します。linalg:norm は sqrt が浮動小数点数を返すため常に浮動小数点数です。クロスバックエンドの注意点が 1 つあります。WASM バックエンドは非終端の浮動小数点数をインタプリタや JVM より少ない有効桁数で印字するため、丸めのある逆行列や無理数のノルムは、内部の double が同一でもバックエンド間で見た目が異なることがあります。
実例
上記の inv と solve の行列は、float の結果が厳密になりどのバックエンドでも同一に印字されるよう選んでいます。(linalg:inv #2A((1 2) (3 4))) のような一般の逆行列は同じ値を計算しますが浮動小数点の丸めを伴います。
la は linalg の組み込みニックネームです。すべての linalg: 呼び出しは、より短い la: プレフィックスでも記述できます。
要素ごとの算術演算とブロードキャスト
linalg:add、linalg:sub、linalg:mul、linalg:div は要素ごとに演算し、numpy の規則でブロードキャストします: どちらか一方のスカラーのオペランドはもう一方のオペランドの形状にブロードキャストされ、形状の異なる 2 つの配列は末尾の軸から揃えられます -- 揃えた各軸の長さは等しいか、どちらかが 1 でなければならず(先頭側の欠けている軸は 1 と見なされます)、長さ 1 の軸がもう一方の長さに引き伸ばされます。どちらにも当てはまらないペアは shape mismatch エラーを通知します。結果は最初の配列オペランドの要素型を保持します(混合幅の規則と同じ)。mul はアダマール積(要素ごとの積)であることに注意してください。行列積は linalg:matmul(またはランクに応じてディスパッチする linalg:dot)です。要素ごとの任意の変換には linalg:emap を使います。
この 4 つには CL 演算子スペル linalg:+・linalg:-・linalg:*・linalg:/ もあります。これらは add / sub / mul / div の可変長引数の左畳み込みなので、(linalg:+ a b c) は段階的にブロードキャストし、各段階は加速されたカーネルのままです。引数が退化した場合は CL に従います: 引数なしは単位元(0 / 1)、+ と * の引数 1 つはその引数自身、- と / の引数 1 つは符号反転 / 逆数です。(vec: パッケージの演算子別名は代わりに厳密に 2 引数です。あちらのカーネルは設計上すべて固定アリティだからです。)
よく使う要素ごとの演算は、numpy の ufunc 名でも用意されています: linalg:exp・linalg:log・linalg:tanh・linalg:sin・linalg:cos・linalg:tan・linalg:asin・linalg:acos・linalg:atan・linalg:sinh・linalg:cosh・linalg:sqrt・linalg:abs・linalg:square・linalg:negative・linalg:sign・linalg:reciprocal、2 引数の linalg:power、さらに比較セレクトの linalg:maximum・linalg:minimum・linalg:clip・linalg:relu(厳密比較 (if (> x y) x y) とその鏡像で定義され、比較が偽なら第 2 被演算子または境界が選ばれます — タイや NaN を含め、すべてのバックエンドで同一の規則です)。いずれも対応する emap(または mul / div / maximum / minimum の呼び出し)と等価ですが、名前付き関数なので --simd で加速されます。任意のコールバックを取る emap は決して加速されません。linalg:softmax・linalg:log-softmax・linalg:erf が relu と同じ理由でここにあります。numpy 本体にはない (scipy / torch のものです) ものの、活性化層が必要とする配列レベルのプリミティブだからです。2 つの softmax は最大値を引いた数値的に安定な形で、erf は厳密な GELU を組み立てる素です。交代級数の Maclaurin 展開では |x| ~ 3 付近で有効数字を失うため、全項が正の級数を加算して double の下位数 ulp まで正確に保っています。
軸に沿った還元
還元関数 linalg:sum、linalg:mean、linalg:amax、linalg:amin は numpy と同じキーワード引数を取ります。整数の :axis(負は numpy 流に末尾から数える)を渡すと配列全体ではなくその軸に沿って還元し、結果からその軸は除去されますが、:keepdims が非 nil のときは長さ 1 の軸として保持されます — 入力にそのままブロードキャストで戻せる形状で、バッチ softmax が行ごとの最大値を引くのに使うのはこの形です。linalg:argmax と linalg:argmin も同じ :axis キーワードを取り、スライスごとの index を返します(行列に対しては packed double 配列 — linalg 配列に整数幅はありません)。linalg:var と linalg:std は同じ :axis / :keepdims に加えて除数の補正 :ddof を取ります (デフォルトは 0 で numpy の np.var、torch の unbiased=False に相当し、1 は標本分散です)。同じ軸の mean と組み合わせると LayerNorm の正規化になります。linalg:reshape は 1 つの -1 extent を受け付け、要素数から推論します。
rank N の形状操作: 連結・スライス・スタックされた積
ここまでの演算はすべて rank を問いませんが、その rank 自体を作り変える操作もいくつかあります。linalg:expand-dims は extent 1 の軸を挿入し、linalg:squeeze はそれを取り除きます (numpy の expand_dims / squeeze、torch の unsqueeze / squeeze)。linalg:concatenate は配列のリストを既存の軸に沿って、linalg:stack は新しい軸に沿って連結します。linalg:slice は numpy の基本スライシングで、linalg:triu / linalg:tril は行列の片側の三角だけを残します。
rontolisp には x[:, :n] という構文がないため、slice は軸ごとに 1 つの spec を並べたリストで表現します。nil はその軸をそのまま残し、(start end) や (start end step) はその軸から選択します。負のインデックスは末尾から数え、start や end の位置の nil は「先頭から」「末尾まで」を意味し、末尾の spec を省略すると残りの軸はそのままです。numpy の x[:, 0:3] が両方の軸を残すのと同様に、すべての軸が保持されます。軸を落とすのは linalg:row の役目です。
linalg:matmul も rank を問いません。いずれかの rank が 3 以上ならスタックされた積になり (torch の bmm / matmul)、最後の 2 軸が行列で先頭側のすべての軸がブロードキャストされます。これはまさにバッチ化されたアテンションスコアの形状で、(batch heads n d) と (batch heads d n) の積が 1 回の呼び出しで (batch heads n n) になります。(linalg:dot は意図的に rank 2 以下のままです。numpy の np.dot は高い rank では別の軸と縮約するため、スタックを渡すと誤った答えを返すのではなくこちらを指すエラーを通知します。)
インデックス操作・選択・マスク
linalg:take-rows は index ベクタで axis-0 スライスを選択し(numpy の x[mask]、任意 rank)、axis 0 を残します。一方 linalg:row は整数で 1 スライスだけ取り出し、axis 0 を落とします(numpy の x[i]。バッチから 1 枚の画像を取り出すと、そのままベクタとして順伝播に渡せます)。linalg:gather は行ごとに 1 要素を取り出し(y[np.arange(n), t])、linalg:one-hot はラベル行列を作ります。要素ごとの比較 linalg:equal、linalg:greater、linalg:greater-equal、linalg:less、linalg:less-equal は 0.0/1.0 のマスクを返します(スカラー被演算子とブロードキャスト対応)。numpy が boolean index する場面では、マスクを掛け算するか、より良い方法として linalg:where に渡してください。非ゼロのマスクに応じて 2 つのオペランドから選択します (numpy の np.where)。掛け算ではなく選択することで、-infinity のマスクがちょうど 0 の重みとして linalg:softmax に届きます。無限大にゼロを掛けると NaN になってしまうからです。linalg:zeros-like は同じ形状・同じ幅のゼロ配列を確保します。
乱数
np.random に相当する乱数はシード可能で、バックエンド間で決定的です。linalg:seed は Wichmann-Hill 生成器をリセットし、その draw は正確な整数演算と IEEE double 演算だけで構成されるため、シード済みの linalg:rand、linalg:randn、linalg:uniform、linalg:choice、linalg:permutation の列は interpreter・JVM・両 WASM ターゲットで bit-identical です — 重み初期化とミニバッチ抽出がどこでも正確に再現されます。randn は Box-Muller ではなく Irwin-Hall(一様乱数 12 個の和)を使います(Box-Muller の log/cos は WASM で発散するため)。そのため裾は 6σ でクリップされます。初期化には十分ですが、分布が np.random.randn と厳密に一致するわけではありません。
離散微積分
linalg:diff と linalg:gradient は numpy の離散微積分ペア(np.diff / np.gradient)です。diff は :axis(デフォルトは最後の軸)に沿った :n 階の離散差分(デフォルト 1)を取ります。1 ステップごとにその軸が 1 つ短くなり、行列はデフォルトでは各行内で、:axis 0 では各列に沿って差分されます。gradient はサンプル値のベクタの微分を、2 次精度の中心差分(両端は 1 次精度の片側差分)で推定するため、結果は入力と同じ長さになります。省略可能な第 2 引数には、一様なサンプル間隔(数値、デフォルト 1)か、非一様なサンプルのための同じ長さの座標ベクタを渡せます。どちらも他の linalg 変換と同様に入力の幅を保持します。算術は通常どおり浮動小数点ですが、厳密に微分できるサンプル値 — 以下の例のような、整数座標で読んだ多項式 — はすべてのバックエンドで同一に印字されます。
#(0 1 4 9 16) の gradient — 放物線 y = x^2 を x = 0..4 でサンプリングしたもの — は、内部の点で真の導関数 2x を厳密に復元します(中心差分は 2 次関数に対して厳密で、両端は 1 次精度の推定です)。座標ベクタ形式なら、最後の行のような不等間隔のサンプルでも厳密なままです。examples/ml/numerical-calculus.lisp は、これらの考え方を投射体の運動の例で一通り示します。
単精度浮動小数点
linalg はデフォルトで double-float で計算しますが、幅多相 (width-polymorphic) です。メモリが半分で SIMD レーン数が倍になる packed 単精度浮動小数点 (#f) 配列を受け付け、その幅を保持します。すべてのコンストラクタは :element-type キーワードを取り(デフォルトは 'double-float。#f の結果が欲しければ :element-type 'single-float を渡します)、すべての変換 -- add/sub/mul/div/emap、transpose/reshape、dot/matmul/outer、inv/solve -- は入力の幅を保持します。したがって単精度の値は最後まで単精度のまま流れます。関数的な重み更新 (linalg:sub W grad) は W の幅を保持し、暗黙に double へ戻す(JVM の --simd パスでは、続く vec:matvec で幅不一致エラーを強制する)ことはありません。f32 の速度とメモリが欲しく精度の低下を許容できるときは単精度を使い、det/inv/solve のような精度が重要な処理にはデフォルトの倍精度を使ってください。
SIMD アクセラレーション
linalg はどこでもフラグなしで正しく動きますが、--simd フラグで加速されます。32 の関数 — add・sub・mul・div・sum・norm・amax・amin・argmax・argmin・trace・transpose・reshape・dot・outer、単項 ufunc の exp・log・tanh・sin・cos・tan・asin・acos・atan・sinh・cosh・sqrt・abs・negative・sign、および比較セレクトの maximum・minimum — がネイティブなベクトルカーネル(インタプリタと JVM では jdk.incubator.vector、wasm-GC では WebAssembly の v128)にルーティングされ、それらを使って書かれている mean・matmul・flatten・solve・square・reciprocal・clip・relu も一緒に加速されます。プログラムが何を受け付け何を拒否するかは一切変わりません。カーネルが扱えない入力(一般の boxed 配列、幅の混在、素の数値)は移植可能な linalg.lisp の定義がそのまま実行され、同じ結果と同じエラーメッセージになります。観測可能な違いは単精度の精度規則だけで、これはリダクションと行列積にかかります。要素ごとの結果はビット一致のままです。最後の一文は --simd に限った話です。--gpu は要素ごとの超越関数も受け取り、デバイスは自前のライブラリでそれらを計算するので、このフラグの下では移植可能な定義と等しくはなく、近いだけになります。
速度のためにパッケージを乗り換える理由はありません。--simd の下では linalg と vec は同じカーネルに行き着きます。vec と linalg の使い分けを参照してください — 要約すると: デフォルトでは linalg に対して書き、vec に手を伸ばすのは -into の書き込み先渡しループ・--no-gc ターゲット(linalg はそこではコンパイルできません)・fail-fast な幅の厳格さが必要なときだけです。
第一級関数
linalg の関数は通常の defun なので、#'linalg:norm などは関数が期待されるあらゆる場所で第一級の値として動作します。
配列は同一性 (eq) でしか比較されないため、結果の比較には形状と数値の等価性を検査する linalg:array-equal を使います(1 と 1.0 は等しいと判定されます)。
パッケージ
linalg はそれ自体が独立したパッケージであり、cl を use していません。(in-package linalg) の中では標準関数は非修飾の名前では見えず、cl: 修飾(cl:print、cl:mapcar など)が必要になります。したがってほとんどのプログラムは、このページのすべての例がそうしているように、デフォルトの cl-user パッケージにとどまり、修飾された linalg: 名で呼び出すべきです。