herk!()だけじゃない!Juliaで行列積を最適化する代替手法

2025-05-26

JuliaにおけるLinearAlgebra.BLAS.herk!()は、線形代数演算の高速化のために使用される関数で、特にエルミート行列のランクk更新を行います。

以下に詳しく説明します。

BLASとは?

まず、BLASは「Basic Linear Algebra Subprograms」の略で、基本的な線形代数演算(ベクトルとベクトルの積、行列とベクトルの積、行列と行列の積など)を高速に実行するための標準的なライブラリです。Juliaは内部的にBLASライブラリ(通常はOpenBLASなど)を利用しており、これらの関数は非常に最適化されています。

!の記号の意味

Juliaの関数名に付いている!は、その関数が引数として渡されたオブジェクトを「インプレースで(その場で)変更する」ことを示します。つまり、新しいメモリを割り当てるのではなく、既存のメモリの内容を直接書き換えます。herk!()も例外ではなく、計算結果を引数として渡された行列Cに直接書き込みます。

herkとは?

herkは「Hermitian Rank-k update」の略です。エルミート行列のランクk更新とは、以下のような計算を指します。

C:=αAAH+βC

または

C:=αAHA+βC

ここで、

  • \alpha, \beta: スカラー値(複素数)
  • A^H: 行列Aの共役転置(エルミート転置)
  • A: 一般の行列
  • C: エルミート行列(更新される行列)

この計算は、行列Aとその共役転置A^Hの積にスカラー\alphaを掛けたものと、既存の行列Cにスカラー\betaを掛けたものを足し合わせ、その結果をCに格納するというものです。

LinearAlgebra.BLAS.herk!()の引数

herk!()関数は通常、以下のような引数を取ります。

herk!(uplo, trans, alpha, A, beta, C)

  • C: エルミート行列。計算結果が格納される行列です。この引数はインプレースで変更されます。
  • beta: スカラー値。既存の行列Cに掛けられる係数です。
  • A: 行列。ランクk更新に使用される行列です。
  • alpha: スカラー値(通常はFloat64, Float32, Complex128, Complex64)。A A^HまたはA^H Aに掛けられる係数です。
  • trans: Char型。'N'または'C'を指定します。
    • 'N' (No transpose): C:=αAAH+βC の形式で計算します。
    • 'C' (Conjugate transpose): C:=αAHA+βC の形式で計算します。
  • uplo: Char型。'U'または'L'を指定します。Cがエルミート行列であるため、その上半部分('U')または下半部分('L')のみを計算・更新することを示します。これにより、計算効率が向上します。
  • 標準化: BLASは線形代数計算の業界標準であり、様々な高性能ライブラリ(OpenBLAS, Intel MKLなど)がそのAPIを実装しています。Juliaはこれらのライブラリをバックエンドとして利用することで、高いパフォーマンスを保証します。
  • 専門性: エルミート行列の特性(H=HH)を活かした効率的な計算が行われます。全ての要素を計算するのではなく、上半部分または下半部分のみを計算することで、計算量とメモリ使用量を削減します。
  • パフォーマンス: BLAS関数は高度に最適化されており、手動でループを記述するよりもはるかに高速です。特に大規模な行列演算においてその効果は顕著です。


引数の型(データ型)の不一致

herk!()はBLASルーチンを呼び出すため、引数の型が厳密に期待される型と一致している必要があります。特に、スカラー値(alpha, beta)と行列(A, C)の要素型に注意が必要です。

よくあるエラーの例

  • herkは「Hermitian」とあるように、複素数行列を対象としています。実数行列に対してはLinearAlgebra.BLAS.syrk!()(symmetric rank-k update)を使用します。間違って実数行列にherk!()を使おうとすると、型エラーになるか、意図しない結果を招く可能性があります。
  • Cが行列ではない、あるいはAの要素型とCの要素型が異なる場合。
  • alphabetaに整数型(Int)を渡してしまう。BLASは浮動小数点数(Float32, Float64)または複素数(ComplexF32, ComplexF64)を期待します。

トラブルシューティング

  • 実数行列の場合
    syrk!()を使用します。
    A_real = rand(Float64, 5, 3)
    C_real = zeros(Float64, 5, 5)
    LinearAlgebra.BLAS.syrk!('U', 'N', 1.0, A_real, 0.0, C_real)
    
  • 型を変換する
    必要に応じて、Float64(value)ComplexF64(value)のように型変換を行います。
    using LinearAlgebra
    A = rand(ComplexF64, 5, 3)
    C = zeros(ComplexF64, 5, 5)
    alpha = 1.0 # Float64をComplexF64に変換
    beta = 0.0 # Float64をComplexF64に変換
    
    # 正しい例
    LinearAlgebra.BLAS.herk!('U', 'N', ComplexF64(alpha), A, ComplexF64(beta), C)
    
    # 誤った例 (alpha, betaがInt64の場合)
    # alpha_int = 1
    # beta_int = 0
    # LinearAlgebra.BLAS.herk!('U', 'N', alpha_int, A, beta_int, C) # MethodErrorまたは予期せぬ動作
    
  • 型を確認する
    typeof(variable) を使用して、各引数の型を確認します。

行列の次元の不一致

herk!()は特定の次元要件を満たす必要があります。これらの要件が満たされない場合、BLASライブラリが内部でエラーを報告し、Julia側でDimensionMismatchなどのエラーが発生します。

よくあるエラーの例

  • Cが正方行列ではない場合。
  • C:=αAHA+βC の場合、Aの列数とCの行数/列数が一致しない。
  • C:=αAAH+βC の場合、Aの行数とCの行数/列数が一致しない。

トラブルシューティング

  • 行列のサイズを確認する
    size(matrix)で各行列の次元を確認します。
    • Cはエルミート行列なので、正方行列である必要があります(size(C, 1) == size(C, 2))。
    • trans = 'N' (C:=αAAH+βC): size(C, 1)size(A, 1)が同じである必要があります。
    • trans = 'C' (C:=αAHA+βC): size(C, 1)size(A, 2)が同じである必要があります。
    • また、ランクkA A^HまたはA^H Aの計算で消える次元)も考慮する必要がありますが、BLASルーチンは通常、この部分を自動的に処理します。
  • ドキュメントを確認する
    ?LinearAlgebra.BLAS.herk!でヘルプを確認し、引数の次元要件を理解します。

uploとtransの指定ミス

uplotransは文字型であり、特定の値を指定する必要があります。

よくあるエラーの例

  • 文字型ではなく文字列型("U"など)で渡してしまう。
  • 'U''L''N''C'以外の文字を指定してしまう。

トラブルシューティング

  • 文字型を使用する
    uplo = 'U'のようにシングルクォーテーションで囲んで文字型として指定します。

パフォーマンスの問題(エラーではないが重要)

herk!()は高性能な線形代数演算を目的としていますが、使い方によっては期待通りのパフォーマンスが得られないことがあります。

よくある問題の例

  • スレッド数
    BLASは通常マルチスレッドで動作しますが、スレッド数の設定が不適切だと、パフォーマンスが低下することがあります。
  • BLAS実装の選択
    Juliaは複数のBLASライブラリを切り替えることができます(OpenBLAS, MKLなど)。環境によっては、特定のBLAS実装が最適なパフォーマンスを発揮しない場合があります。
  • キャッシュミス
    行列のメモリレイアウトやアクセスパターンが非効率な場合、キャッシュミスが発生し、パフォーマンスが低下します。Juliaの配列は通常、列優先順序(Column-Major)で格納されます。
  • 型の不安定性 (Type Instability)
    Juliaコード全体で型の不安定な部分があると、BLAS呼び出しのオーバーヘッドが大きくなり、最適化の恩恵を十分に受けられないことがあります。

トラブルシューティング

  • ベンチマーク
    @btime (BenchmarkTools.jl) を使用して、様々な設定や実装でのパフォーマンスを比較します。
  • システムイメージの再構築
    稀に、BLASライブラリが正しくリンクされていないなどの問題が発生することがあります。JuliaのPkg.build("LinearAlgebra")を実行したり、Julia自体を再インストールしたりすることで解決する場合があります。
  • BLAS設定の確認
    • LinearAlgebra.BLAS.vendor() で現在使用されているBLASベンダーを確認できます。
    • BLAS.set_num_threads(n) でBLASが使用するスレッド数を設定できます。最適なスレッド数はCPUコア数やワークロードによって異なります。
  • 列優先順序を意識する
    行列操作において、列方向のアクセスがメモリ上で連続しているため、効率的です。可能であれば、列方向の操作を優先するようにコードを記述します。
  • @code_warntypeの使用
    関数呼び出しの前に@code_warntypeマクロを使用して、型の不安定性がないか確認します。
    function my_herk_operation(A, C, alpha, beta)
        LinearAlgebra.BLAS.herk!('U', 'N', alpha, A, beta, C)
    end
    # @code_warntype my_herk_operation(A, C, ComplexF64(1.0), ComplexF64(0.0))
    

これはherk!()固有のエラーではありませんが、BLAS関数を直接使用する際に発生する一般的な問題です。LinearAlgebraモジュールもdotなどの一部の関数をエクスポートしており、LinearAlgebra.BLASも同じ名前の関数をエクスポートしている場合があります。

よくあるエラーの例

using LinearAlgebra
using LinearAlgebra.BLAS # これによりdotなどの名前が衝突する可能性がある

# ERROR: UndefVarError: dot not defined
# または WARNING: both BLAS and LinearAlgebra export "dot"; uses of it in module Main must be qualified
  • または、using LinearAlgebra; const LBLAS = LinearAlgebra.BLAS; のようにエイリアスを付けて、LBLAS.herk!()と呼び出すこともできます。
  • BLAS関数を頻繁に呼び出す場合は、必要な関数だけをimportまたはusing LinearAlgebra.BLAS: herk!, gemm!のように指定することもできます。
  • LinearAlgebra.BLAS.herk!()のように、常に完全修飾名でBLAS関数を呼び出すのが最も安全な方法です。


事前準備

herk!() を使用するには、LinearAlgebra モジュールをインポートする必要があります。

using LinearAlgebra
# BLAS関数を直接呼び出すので、通常は LinearAlgebra.BLAS.herk!() のように完全修飾名で呼び出します。
# もし BLAS を頻繁に使うなら、以下のようにエイリアスを付けると便利です。
# const LBLAS = LinearAlgebra.BLAS

例1: 基本的な使い方 (C:=AAH)

最も基本的なケースとして、C = A * A' (ここで A' は共役転置) を計算する例です。これは alpha = 1.0, beta = 0.0 の場合に対応します。

using LinearAlgebra

# 行列 A を定義 (複素数行列である必要がある)
A = rand(ComplexF64, 5, 3) # 5x3 の複素数行列

# 結果を格納するエルミート行列 C を定義
# C は A の行数 (5) と同じ次元の正方行列になる
C = zeros(ComplexF64, 5, 5) # 初期値はゼロ行列

# herk!() を呼び出す
# 'U': C の上三角部分を更新 (Lower 'L' も可能)
# 'N': A を転置しない (C := alpha * A * A^H + beta * C)
# 1.0 + 0.0im: alpha (複素数スカラー)
# A: 行列 A
# 0.0 + 0.0im: beta (複素数スカラー)
# C: 更新される行列 C
LinearAlgebra.BLAS.herk!('U', 'N', 1.0 + 0.0im, A, 0.0 + 0.0im, C)

println("行列 A:")
display(A)
println("\n計算結果 C (herk!):")
display(C)

# 検証: Julia の通常の行列積と比較
C_ref = A * A'
println("\n検証結果 C (A * A'):")
display(C_ref)

# C がエルミート行列であることを確認
# C の上三角部分のみを計算しているので、下三角部分はゼロのままです。
# 通常、エルミート行列の処理では、片側の三角部分を計算し、
# もう片側は共役転置で導出するか、対称性を利用してアクセスします。
# 完全に充填されたエルミート行列が必要な場合は、例えば C = C + C' - Diagonal(C) のような操作が必要です。

# この場合、herk!() は上三角部分のみを埋めるため、C == C_ref は false になることがあります。
# C の上三角部分と C_ref の上三角部分を比較します。
println("\nC の上三角部分と C_ref の上三角部分が一致するか:")
println(C[triu(trues(size(C)))] ≈ C_ref[triu(trues(size(C)))])

解説

  • スカラー alphabetaComplexF64 型で渡す必要があります。1.00.0 を直接渡すと、Julia は Float64 と推論し、メソッドエラーになる可能性があります。1.0 + 0.0im のように明示的に複素数として指定するのが安全です。
  • 'N'A が転置されないことを意味し、計算式は C:=αAAH+βC になります。
  • 'U'C の上三角部分(Upper part)のみを計算・更新することを意味します。エルミート行列は対称性を持つため、片側の三角部分を計算すれば十分です。

例2: C:=αAHA+βC の計算

trans 引数を 'C' に設定することで、C:=αAHA+βC の形式で計算できます。

using LinearAlgebra

A = rand(ComplexF64, 3, 5) # 3x5 の複素数行列
# この場合、A^H A は 5x5 行列になる
C = zeros(ComplexF64, 5, 5)

alpha = 2.0 + 1.0im
beta = 0.5 + 0.0im

# C の初期値を設定(beta * C の影響を見るため)
C_initial = rand(ComplexF64, 5, 5)
C_initial = C_initial + C_initial' # C をエルミートにする
C .= C_initial # C に初期値をコピー

println("初期の C:")
display(C)

# herk!() を呼び出す
# 'L': C の下三角部分を更新
# 'C': A を共役転置して使用 (C := alpha * A^H * A + beta * C)
LinearAlgebra.BLAS.herk!('L', 'C', alpha, A, beta, C)

println("\n計算結果 C (herk!):")
display(C)

# 検証: Julia の通常の行列積と比較
C_ref = alpha * (A' * A) + beta * C_initial
println("\n検証結果 C (alpha * A' * A + beta * C_initial):")
display(C_ref)

# C の下三角部分と C_ref の下三角部分を比較
println("\nC の下三角部分と C_ref の下三角部分が一致するか:")
println(C[tril(trues(size(C)))] ≈ C_ref[tril(trues(size(C)))])

解説

  • beta に非ゼロの値を指定することで、既存の C の値が計算にどのように影響するかを確認できます。
  • 'C'A が共役転置されて使用されることを意味し、計算式は C:=αAHA+βC になります。
  • 'L'C の下三角部分(Lower part)のみを計算・更新することを意味します。

herk!() は「Hermitian」(エルミート)と名前にある通り、複素数行列向けです。実数行列の場合には、LinearAlgebra.BLAS.syrk!() (symmetric rank-k update) を使用します。使い方は herk!() と非常に似ていますが、スカラー引数は実数型になります。

using LinearAlgebra

A_real = rand(Float64, 5, 3) # 実数行列
C_real = zeros(Float64, 5, 5) # 結果を格納する実数対称行列

alpha_real = 2.0
beta_real = 0.0

# syrk!() を呼び出す
# 'U': C_real の上三角部分を更新
# 'N': A_real を転置しない (C := alpha * A * A' + beta * C)
# alpha_real: スカラー (Float64)
# A_real: 行列 A_real
# beta_real: スカラー (Float64)
# C_real: 更新される行列 C_real
LinearAlgebra.BLAS.syrk!('U', 'N', alpha_real, A_real, beta_real, C_real)

println("\n実数行列 A_real:")
display(A_real)
println("\n計算結果 C_real (syrk!):")
display(C_real)

# 検証
C_real_ref = alpha_real * (A_real * A_real') + beta_real * zeros(Float64, 5, 5)
println("\n検証結果 C_real (alpha_real * A_real * A_real'):")
display(C_real_ref)

println("\nC_real の上三角部分と C_real_ref の上三角部分が一致するか:")
println(C_real[triu(trues(size(C_real)))] ≈ C_real_ref[triu(trues(size(C_real)))])

解説

  • 引数のスカラーは Float64 型で問題ありません。
  • syrk!() は対称行列 (C = C') の更新に使用されます。

herk!()syrk!() のような BLAS 関数は、特に大規模な行列演算において絶大なパフォーマンスを発揮します。手動でループを記述するよりもはるかに高速です。これは、BLAS の実装が高度に最適化されており、CPU のキャッシュ効率、SIMD 命令、マルチスレッド処理などを最大限に活用するように設計されているためです。



標準的な行列乗算 (* 演算子)

Julia の最も一般的な代替手段は、行列乗算演算子 * を使用することです。これは非常に直感的で、Julia が内部的に最適化された BLAS ルーチンを呼び出すため、多くの場合に十分なパフォーマンスを提供します。

コード例

using LinearAlgebra

A = rand(ComplexF64, 500, 300) # 500x300 の複素数行列
C = zeros(ComplexF64, 500, 500) # 結果を格納する行列

alpha = 1.0 + 0.0im
beta = 0.0 + 0.0im

# herk! を使用した計算
C_herk = zeros(ComplexF64, 500, 500)
LinearAlgebra.BLAS.herk!('U', 'N', alpha, A, beta, C_herk)
# C_herk の下三角部分も埋める必要がある場合は、
# C_herk = C_herk + C_herk' - Diagonal(C_herk)
# または C_herk = Hermitian(C_herk, :U) のように Hermitian 型に変換する

# 標準的な行列乗算と加算を使用した計算
C_std = alpha * A * A' + beta * C

println("herk! と標準演算子の結果の差:")
# herk! は上三角部分のみを計算するので、比較のために C_herk の下三角部分をコピーする必要がある
C_herk_full = C_herk + C_herk' - Diagonal(C_herk)
println(norm(C_herk_full - (alpha * A * A' + beta * C)))

利点

  • インプレースでない操作
    herk!() と異なり、C = alpha * A * A' + beta * C は新しい行列 C を割り当てます。これは、元の C を変更したくない場合に便利です。
  • 自動最適化
    Julia の行列乗算は、通常、バックエンドの BLAS ライブラリ(OpenBLAS, MKL など)によって最適化されているため、多くの場合は十分なパフォーマンスが得られます。
  • 安全性
    Julia が次元チェックなどを自動的に行い、エラーを早期に検出します。
  • 直感的で読みやすい
    数学的な式に非常に近い形で記述できます。

欠点

  • エルミート性の利用不足
    herk!() はエルミート行列の対称性を利用して計算量を半減させますが、標準的な A * A' は全ての要素を計算します。これは、エルミート行列 C が最終的に必要とされる場合に効率が悪い可能性があります。
  • メモリ割り当て
    alpha * A * A' の結果として一時的な行列が作成され、その後に beta * C と加算されるため、追加のメモリ割り当てが発生する可能性があります。大規模な行列や繰り返し計算では、これがパフォーマンスのボトルネックになることがあります。

LinearAlgebra.mul! 関数

mul! (multiply-add!) 関数は、行列乗算と加算を組み合わせた操作をインプレースで行うための関数です。これにより、一時的なメモリ割り当てを削減し、パフォーマンスを向上させることができます。これは BLAS.gemm! (General Matrix-Matrix Multiplication) の高レベルラッパーと考えることができます。

コード例

using LinearAlgebra

A = rand(ComplexF64, 500, 300)
C = zeros(ComplexF64, 500, 500)

alpha = 1.0 + 0.0im
beta = 0.0 + 0.0im

# mul! を使用した計算: C = alpha * A * A' + beta * C
# C はエルミート行列として扱われないため、全てを計算する
LinearAlgebra.mul!(C, A, A', alpha, beta)

println("mul! と herk! の結果の差:")
C_herk = zeros(ComplexF64, 500, 500)
LinearAlgebra.BLAS.herk!('U', 'N', alpha, A, beta, C_herk)
C_herk_full = C_herk + C_herk' - Diagonal(C_herk) # herk! の結果を完全なエルミート行列にする

println(norm(C_herk_full - C))

利点

  • 柔軟性
    mul!(C, A, B) のように、一般的な行列 B との積も計算できます。
  • パフォーマンス
    * 演算子よりも効率的な場合があり、特に大きな行列で顕著です。
  • インプレース操作
    結果を既存の行列 C に直接書き込むため、メモリ割り当てが削減されます。

欠点

  • mul! の引数
    mul!(C, A, A') のように A' を明示的に指定する必要があります。
  • エルミート性の利用不足
    mul! はエルミート行列の特別な構造を認識しません。herk!() と異なり、全ての要素を計算し、メモリ上でのエルミート性の利用は行いません。そのため、エルミート行列 C の計算としては、herk!() の方が最終的に効率的である可能性があります。

Symmetric または Hermitian 型の利用

Julia の LinearAlgebra モジュールには、Symmetric および Hermitian という特殊な行列型が用意されています。これらの型は、行列が対称性またはエルミート性を持つことをコンパイラに伝え、特定の操作で最適化を可能にします。

コード例

using LinearAlgebra

A = rand(ComplexF64, 500, 300)
C_init = zeros(ComplexF64, 500, 500) # 初期値が必要な場合

alpha = 1.0 + 0.0im
beta = 0.0 + 0.0im

# Hermes を使用しない場合(通常の計算)
C_standard = alpha * A * A' + beta * C_init

# Hermitian 型を使用した計算
# A * A' は Hermitian 型を返す
C_herm = Hermitian(alpha * (A * A')) + beta * C_init

# または、インプレースで Hermitian(C, :U) のように直接操作
# (ただし、これは C が既にエルミート構造を持っている場合に最も効果的)
# herk!() の直接の代替というよりは、結果をエルミート行列として扱う方法

println("標準と Hermitian 型の差:")
println(norm(C_standard - C_herm)) # 数値的には同じ結果が得られるはず

利点

  • 読みやすさ
    コードの意図がより明確になります。
  • 一部の操作の最適化
    Hermitian 型に対して特定の操作(例: 固有値計算)を行う場合、内部的に最適化されたルーチンが選択されることがあります。
  • 意味論的な正確さ
    行列がエルミートであることを明示的に表現できます。

欠点

  • メモリ割り当て
    alpha * (A * A') の部分は、herk!() のようにインプレースでなく、新しい行列を割り当てます。
  • 直接的なランクk更新ではない
    これは herk!() のような直接的なランクk更新の代替というよりも、行列がエルミートであることを宣言し、その後の操作で恩恵を受けるためのものです。

低レベルな BLAS/LAPACK 関数への直接アクセス(上級者向け)

Julia は、LinearAlgebra.BLAS を通じて様々な BLAS ルーチンへのアクセスを提供します。herk!() はその一つですが、もし特殊なケースで herk!() が利用できない、またはより詳細な制御が必要な場合は、他の BLAS ルーチンを直接呼び出すことも可能です。しかし、これは非常に稀なケースであり、通常は推奨されません。


gemm! (General Matrix-Matrix Multiplication) は C:=αAB+βC を計算しますが、herk! が利用する最適化(B=AH のような関係を利用した計算量削減)は行いません。

利点

  • ニッチな最適化
    特定の非常に特殊な状況で、より高いパフォーマンスを引き出せる可能性があります。
  • 究極の制御
    BLAS ライブラリの各引数に直接アクセスできます。
  • 通常は不要
    herk!()mul! が既に高度に最適化されているため、ほとんどのユースケースでこれらの直接呼び出しは不要です。
  • 安全性の欠如
    型や次元の不一致が、分かりにくいエラーやクラッシュにつながる可能性があります。
  • 複雑性
    引数の意味を完全に理解し、正しく指定する必要があります。
  1. 最も一般的なケース
    ほとんどの場合、標準的な行列乗算 (* 演算子) を使用するのが最も簡単で、パフォーマンスも十分です。
  2. メモリ割り当てを最小限に抑えたい場合
    大規模な行列で繰り返し計算を行う場合や、メモリ使用量が懸念される場合は、LinearAlgebra.mul! が良い選択肢です。ただし、herk!() のようにエルミート性を利用した計算量削減は行いません。
  3. エルミート行列のランクk更新を高速かつ効率的に行いたい場合
    LinearAlgebra.BLAS.herk!() が最も推奨されます。これはエルミート行列の特性を最大限に利用し、計算量とメモリ使用量の両方を最適化します。
  4. コードの意図を明確にしたい場合
    エルミート行列を扱う際には、Hermitian 型を使用することで、コードの可読性を高めることができます。