JuliaのBLAS.hemv()だけじゃない！線形代数計算のベストプラクティス

LinearAlgebra.BLAS.hemv() とは

LinearAlgebra.BLAS.hemv() は、Juliaの標準ライブラリである LinearAlgebra モジュールの一部であり、BLAS (Basic Linear Algebra Subprograms) の中の hemv (Hermitian Matrix-Vector product) ルーチンへのラッパー関数です。

BLASは、ベクトルと行列の基本的な線形代数演算（ベクトルの内積、行列とベクトルの積、行列の積など）を効率的に実行するための標準化されたルーチンの集合です。これらのルーチンは、CやFortranで最適化された実装（例えばOpenBLASやIntel MKLなど）が利用され、非常に高速な計算が可能です。

hemv は特に、エルミート行列とベクトルの積を計算するために使用されます。エルミート行列は、自身の共役転置と等しい複素行列のことです。実数の場合は対称行列に相当します。

関数の目的

LinearAlgebra.BLAS.hemv!() (感嘆符!が付いている場合) は、以下のような演算を実行します。

y=αAx+βy

ここで、

• α と β はスカラー値です。
• y は出力ベクトルであり、この関数によって更新されます。
• x は入力ベクトルです。
• A はエルミート行列です。

感嘆符 ! が付いている hemv! は、入力引数 y を直接変更（上書き）します。もし y を変更せずに新しいベクトルを返したい場合は、感嘆符なしの hemv を使用します。

使用例 (概念)

具体的な引数の詳細はJuliaのドキュメントを参照する必要がありますが、一般的には以下のような引数を取ります。

using LinearAlgebra

# エルミート行列の例 (実際には複素数の場合が多い)
A = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
x = [1.0+0.0im, 2.0+0.0im]
y = [0.0+0.0im, 0.0+0.0im] # y を初期化

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

# y = alpha * A * x + beta * y を計算し、y を更新
LinearAlgebra.BLAS.hemv!('U', A, x, y, alpha, beta)
# または
LinearAlgebra.BLAS.hemv!('L', A, x, y, alpha, beta)
# 'U' または 'L' は、行列 A の上三角または下三角部分がどのように格納されているかを示します。
# エルミート行列は対称なので、通常は片方の三角部分のみを渡します。

# 新しいベクトルを返す場合
y_new = LinearAlgebra.BLAS.hemv('U', A, x) # alpha=1, beta=0 がデフォルト

Juliaで線形代数演算を行う際、通常はより高レベルな関数（例: A * x や mul!(y, A, x)）を使用します。これらの高レベルな関数は、内部的に最適なBLASルーチンを自動的に呼び出します。

しかし、LinearAlgebra.BLAS.hemv() のような低レベルのBLASラッパーを直接使用するメリットは、特定の非常にパフォーマンスが重要なシナリオや、BLASルーチンの持つ細かな制御（例えば、行列の特定の三角部分のみを使う、特定の定数をかけるなど）が必要な場合にあります。

ほとんどのユーザーは、Juliaの標準的な線形代数演算子や関数を使うことで、十分なパフォーマンスを得られます。BLAS.hemv() は、より詳細な最適化や、BLASの動作について深く理解したい開発者向けと言えるでしょう。

• 通常は高レベルの線形代数関数を使うが、特定の最適化や制御が必要な場合に直接使用される。
• CやFortranで最適化されたBLASライブラリ（OpenBLAS, MKLなど）を利用するため、非常に高速。
• hemv! は結果を既存のベクトル y に上書きし、hemv は新しいベクトルを返す。
• エルミート行列とベクトルの積 y=αAx+βy を計算する。
• LinearAlgebra.BLAS.hemv() はJuliaのBLASラッパーの一つ。

---

型の不一致 (Type Mismatches)

BLAS.hemv() はBLASルーチンへの直接のラッパーであるため、引数の型には厳密な要件があります。

• トラブルシューティング
  
  • すべての数値引数が浮動小数点型（Float64、ComplexF64 が一般的）であることを確認してください。
  • A, x, y がAbstractMatrix および AbstractVector のサブタイプであることを確認してください。
  • 例:

    using LinearAlgebra
    
    # エラーになる例 (Int型)
    # A = [1 2; 3 4]
    # x = [1, 2]
    # LinearAlgebra.BLAS.hemv!('U', A, x, similar(x), 1, 0)
    
    # 正しい例 (Float64型)
    A_ok = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
    x_ok = [1.0+0.0im, 2.0+0.0im]
    y_ok = similar(x_ok) # x_ok と同じ型とサイズで初期化
    LinearAlgebra.BLAS.hemv!('U', A_ok, x_ok, y_ok, 1.0+0.0im, 0.0+0.0im)
    

• 原因
  行列 A、ベクトル x、y、およびスカラー alpha、beta の型がBLASルーチンが期待する型（Float32, Float64, ComplexF32, ComplexF64 など）と一致しない場合。特に、整数型や他の非数値型を渡そうとすると発生します。
• エラーの例
  MethodError: no method matching hemv!(...) または ArgumentError: no BLAS method for type ...

行列の形状と次元の不一致 (Dimension Mismatches)

線形代数演算では、行列とベクトルの次元が一致する必要があります。

• トラブルシューティング
  
  • size(A)、length(x)、length(y) を確認し、適切な次元であることを確認してください。
  • A は N×N の行列である必要があります。
  • x と y は N の長さのベクトルである必要があります。
• 原因
  
  • 行列 A が正方行列でない場合（エルミート行列は定義により正方行列である必要があります）。
  • ベクトル x のサイズが行列 A の列数（または行数）と一致しない場合。
  • ベクトル y のサイズが行列 A の行数（または列数）と一致しない場合。
• エラーの例
  DimensionMismatch("...")

エルミート行列の指定 ('U' または 'L') の誤り

hemv() 関数は、エルミート行列 A の上三角部分 ('U') または下三角部分 ('L') のどちらを使用するかを指定する引数を持ちます。これは、BLASルーチンが行列全体ではなく、効率のために半分のみを読み取るためです。

• トラブルシューティング
  
  • 行列 A を定義する際に、BLASに渡す上三角または下三角部分に適切なデータがあることを確認してください。
  • Juliaの標準的な行列は列優先で格納されるため、通常は上三角部分 'U' を指定するのが直感的かもしれませんが、データがどのように用意されているかに依存します。
• 原因
  実際には上三角部分にデータが格納されているにも関わらず 'L' を指定したり、その逆を行ったりする場合。
• エラーの例
  明示的なエラーメッセージは出にくいですが、間違った結果や予期しない動作につながることがあります。

出力ベクトル y の初期化と上書き (hemv!)

hemv! (感嘆符付き) 関数は、引数として渡された y ベクトルを直接変更します。

• トラブルシューティング
  
  • y を適切なサイズと型のゼロベクトルで初期化するか、既存の値を上書きしても問題ないことを確認してください。
  • similar(x) を使用して、x と同じ型とサイズのベクトルを簡単に作成できます。
  • もし y を変更したくない場合は、感嘆符なしの LinearAlgebra.BLAS.hemv() を使用してください。これは新しいベクトルを返します。
• 原因
  y を適切に初期化せずに渡す、または y の初期値が結果に影響することを忘れる。
• エラーの例
  初期化されていない、または意図しない値が入っている y を渡すと、予期しない結果になります。

BLASライブラリに関する問題 (OpenBLAS/MKLなど)

稀に、Juliaが使用しているBLASライブラリ（デフォルトはOpenBLAS）自体に問題がある場合があります。特にWindows環境や特定のCPUアーキテクチャで報告されることがあります。

• トラブルシューティング
  
  • Juliaのアップデート
    最新のJuliaバージョンでは、BLASライブラリも更新されている可能性があります。
  • BLASバックエンドの切り替え
    Intel MKL (MKL.jl パッケージ) など、異なるBLASライブラリを試すことで問題が解決する場合があります。
  • スレッド数の設定
    BLASのスレッド数が多すぎるとパフォーマンス問題やデッドロックを引き起こすことがあります。環境変数 OPENBLAS_NUM_THREADS または JULIA_NUM_THREADS を設定して試してみてください。
  • 問題の特定
    問題が再現可能であれば、より高レベルの A * x や mul!(y, A, x) で同じ問題が起きるか確認することで、BLASレベルの問題なのかJuliaのラッパーの問題なのかを切り分けられます。
• 原因
  BLASライブラリのバグ、特定のCPUの最適化との非互換性、BLASスレッド設定の問題など。
• エラーの例
  計算結果が誤っている、クラッシュする、パフォーマンスが著しく悪い。

複素数と実数の取り扱い

hemv はエルミート行列（複素数）を扱うために設計されていますが、実対称行列（実数）もエルミート行列の特殊なケースです。

• トラブルシューティング
  
  • 複素数を扱う場合は、ComplexF32 または ComplexF64 を使用していることを確認してください。
  • 実数を扱う場合は、Float32 または Float64 を使用していることを確認してください。BLAS.hemv は実数に対しても有効ですが、実際には内部で対称行列-ベクトル積のルーチン（symv）が使われる可能性があります。
• 原因
  引数の数値型が、想定しているBLASルーチンに合致しない場合。
• エラーの例
  実数しか扱えないBLASルーチンを複素数に適用しようとする、またはその逆。

---

感嘆符 ! が付く hemv! は、結果が出力引数 y に上書きされます。感嘆符なしの hemv は新しいベクトルを返します。

LinearAlgebra.BLAS.hemv!() の基本的な使い方

hemv! の基本的な呼び出し形式は以下のようになります。

hemv!(uplo, A, x, y, alpha, beta)

• beta：スカラー値。y に乗算される係数。
• alpha：スカラー値。A*x に乗算される係数。
• y：出力ベクトル（Vector{ComplexF64} または Vector{Float64} など）。このベクトルに結果が上書きされます。x と同じサイズである必要があります。
• x：入力ベクトル（Vector{ComplexF64} または Vector{Float64} など）。
• A：エルミート行列（Matrix{ComplexF64} または Matrix{Float64} など）。
• uplo：Char 型で、行列 A の上三角部分を使う場合は 'U'、下三角部分を使う場合は 'L' を指定します。

計算される内容は y=αAx+βy です。

例1: 複素エルミート行列の場合 (hemv!)

エルミート行列は、共役転置が自身と等しい複素行列です（A=AH）。対角成分は実数、非対角成分は共役複素数の関係にあります。

using LinearAlgebra

println("--- 例1: 複素エルミート行列とベクトル積 (hemv!) ---")

# エルミート行列 A を定義
# A = [ 2.0+0.0im  1.0-2.0im ]
#     [ 1.0+2.0im  3.0+0.0im ]
# ここで A[1,2] の共役が A[2,1] となっていることに注意
A = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
println("行列 A:\n", A)
println("行列 A の型: ", typeof(A))

# 入力ベクトル x
x = [1.0+0.0im, 2.0+0.0im]
println("\n入力ベクトル x:\n", x)
println("ベクトル x の型: ", typeof(x))

# 出力ベクトル y を初期化（x と同じ型とサイズで）
y = similar(x) # [0.0+0.0im, 0.0+0.0im] となる
println("\n初期化された出力ベクトル y:\n", y)
println("ベクトル y の型: ", typeof(y))

# スカラー alpha と beta
alpha = 1.0 + 0.0im # Ax をそのまま加算
beta = 0.0 + 0.0im  # y をゼロクリアして結果を格納

# hemv! を呼び出して y = alpha * A * x + beta * y を計算
# 'U' は A の上三角部分（A[1,1], A[1,2], A[2,2]）を使用することを意味する
LinearAlgebra.BLAS.hemv!('U', A, x, y, alpha, beta)

println("\n計算後の出力ベクトル y:\n", y)

# 結果の検証 (Juliaの標準演算子と比較)
expected_y = A * x
println("A * x (Julia標準演算子による検証):\n", expected_y)
println("結果の一致: ", isapprox(y, expected_y))

# beta を非ゼロにする例
y_initial = [10.0+1.0im, 20.0+2.0im] # y の初期値を設定
y_modified = copy(y_initial) # 上書きされるのでコピーを作成
alpha_val = 0.5 + 0.0im
beta_val = 2.0 + 0.0im
println("\n--- beta を非ゼロにした例 ---")
println("初期 y_modified:\n", y_modified)
LinearAlgebra.BLAS.hemv!('U', A, x, y_modified, alpha_val, beta_val)
println("計算後の y_modified:\n", y_modified)

# 検証: y_modified = 0.5 * A * x + 2.0 * y_initial
expected_y_modified = alpha_val * (A * x) + beta_val * y_initial
println("検証 (0.5 * A * x + 2.0 * y_initial):\n", expected_y_modified)
println("結果の一致: ", isapprox(y_modified, expected_y_modified))

解説

• hemv!('U', A, x, y, alpha, beta) は、行列 A の上三角部分のみを見て計算を行います。エルミート行列の定義により、下三角部分は上三角部分から自動的に導出されます。
• similar(x) は x と同じ要素型と次元を持つ新しい配列を割り当て、初期値は設定されません（数値型の場合は通常ゼロ）。
• x と y も同じく ComplexF64 型のベクトルです。
• A は ComplexF64 型のエルミート行列として定義されています。A[2,1] は A[1,2] の共役になっています。

例2: 実対称行列の場合 (hemv! / symv!)

実対称行列は、転置が自身と等しい実数行列です（A=AT）。これはエルミート行列の特殊なケースであり、BLASの hemv ルーチンは実数に対しても機能します。内部的には、より特化した symv ルーチンが呼び出されることもあります。

using LinearAlgebra

println("\n--- 例2: 実対称行列とベクトル積 (hemv! / symv!) ---")

# 対称行列 A を定義
# A = [ 2.0  1.0 ]
#     [ 1.0  3.0 ]
A_real = [2.0 1.0; 1.0 3.0]
println("行列 A_real:\n", A_real)
println("行列 A_real の型: ", typeof(A_real))

# 入力ベクトル x
x_real = [1.0, 2.0]
println("\n入力ベクトル x_real:\n", x_real)
println("ベクトル x_real の型: ", typeof(x_real))

# 出力ベクトル y を初期化
y_real = similar(x_real)
println("\n初期化された出力ベクトル y_real:\n", y_real)
println("ベクトル y_real の型: ", typeof(y_real))

# スカラー alpha と beta
alpha_real = 1.0
beta_real = 0.0

# hemv! を呼び出して y = alpha * A * x + beta * y を計算
# 実対称行列に対しても hemv! を使用できる
LinearAlgebra.BLAS.hemv!('L', A_real, x_real, y_real, alpha_real, beta_real)
# または BLAS.symv! を明示的に使うことも可能 (同じ引数)
# LinearAlgebra.BLAS.symv!('L', A_real, x_real, y_real, alpha_real, beta_real)

println("\n計算後の出力ベクトル y_real:\n", y_real)

# 結果の検証
expected_y_real = A_real * x_real
println("A_real * x_real (Julia標準演算子による検証):\n", expected_y_real)
println("結果の一致: ", isapprox(y_real, expected_y_real))

解説

• この例では 'L' を指定しています。これは、行列 A_real の下三角部分（A_real[1,1], A_real[2,1], A_real[2,2]）が使用されることを意味します。
• 実数の場合でも hemv! を使用できます。これは、実対称行列がエルミート行列の特殊なケースだからです。
• A_real, x_real, y_real は Float64 型です。

hemv() (感嘆符なし) は、新しいベクトルを返します。この場合、y を出力引数として渡す必要はなく、alpha と beta のデフォルト値はそれぞれ 1.0 と 0.0 です。

using LinearAlgebra

println("\n--- 例3: 感嘆符なしの hemv() ---")

A_complex = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
x_complex = [1.0+0.0im, 2.0+0.0im]

# hemv() は新しいベクトルを返す
# デフォルトで alpha=1.0, beta=0.0
result_vector = LinearAlgebra.BLAS.hemv('U', A_complex, x_complex)

println("A_complex:\n", A_complex)
println("x_complex:\n", x_complex)
println("\nhemv('U', A_complex, x_complex) の結果:\n", result_vector)

# 結果の検証
expected_result = A_complex * x_complex
println("A_complex * x_complex (Julia標準演算子による検証):\n", expected_result)
println("結果の一致: ", isapprox(result_vector, expected_result))

解説

• これは、既存のベクトルを上書きしたくない場合に便利です。
• hemv() を呼び出すだけで、計算結果を含む新しいベクトルが返されます。

LinearAlgebra.BLAS.hemv() 関数は、エルミート行列（または実対称行列）とベクトルの積を効率的に計算するための低レベルなインターフェースを提供します。通常、Juliaの線形代数演算は内部でこれらのBLASルーチンを自動的に呼び出すため、ほとんどのユーザーは A * x のような高レベルな構文を使用するだけで十分です。しかし、特定のパフォーマンスチューニングやBLASの詳細な動作を理解する必要がある場合に、これらの低レベルなラッパー関数を直接使用することができます。

---

以下に主な代替方法を説明します。

標準の行列とベクトルの乗算演算子 (*)

最も一般的で推奨される方法です。Juliaは線形代数演算子を非常に効率的に実装しており、内部的に最適なBLASルーチンを自動的に呼び出します。

特徴

• メモリ効率
  新しいベクトルを割り当てて結果を返します。既存のベクトルを上書きしたい場合は、mul! を使用します。
• 効率的
  Juliaの線形代数実装が内部で自動的に最適なBLAS/LAPACKルーチン（この場合は hemv または symv）を選択し、呼び出します。
• 最もシンプルで直感的
  コードが読みやすく、書きやすい。

コード例

using LinearAlgebra

println("--- 代替方法1: 標準の行列とベクトルの乗算演算子 (*) ---")

# 複素エルミート行列
A_complex = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
x_complex = [1.0+0.0im, 2.0+0.0im]

# A * x を計算 (内部で hemv が呼ばれる)
result_complex = A_complex * x_complex
println("複素エルミート行列の場合 (A * x):\n", result_complex)

# 実対称行列
A_real = [2.0 1.0; 1.0 3.0]
x_real = [1.0, 2.0]

# A * x を計算 (内部で symv が呼ばれる)
result_real = A_real * x_real
println("\n実対称行列の場合 (A * x):\n", result_real)

インプレース乗算関数 (mul!)

mul! 関数は、行列とベクトルの積の結果を事前に割り当てられた出力ベクトルに「インプレース（in-place）」で書き込みます。これにより、新しいメモリ割り当てを避けることができるため、ループ内などでの頻繁な計算でパフォーマンスが向上することがあります。

特徴

• 直感的
  * 演算子と同様に、内部で最適なBLASルーチン（この場合は hemv! または symv!）を呼び出します。
• パフォーマンス
  メモリ割り当てのオーバーヘッドがないため、大規模な計算や反復計算で有利です。
• メモリ効率
  結果を既存のベクトルに直接書き込むため、新しいメモリ割り当てを最小限に抑えます。

コード例

using LinearAlgebra

println("\n--- 代替方法2: インプレース乗算関数 (mul!) ---")

# 複素エルミート行列
A_complex = [2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im]
x_complex = [1.0+0.0im, 2.0+0.0im]

# 出力ベクトルを事前に割り当てる
y_complex = similar(x_complex) # x と同じ型とサイズで初期化される

# mul!(y, A, x) で y = A * x を計算
mul!(y_complex, A_complex, x_complex)
println("複素エルミート行列の場合 (mul!(y, A, x)):\n", y_complex)


# 実対称行列
A_real = [2.0 1.0; 1.0 3.0]
x_real = [1.0, 2.0]

y_real = similar(x_real)

# mul!(y, A, x) で y = A * x を計算
mul!(y_real, A_real, x_real)
println("\n実対称行列の場合 (mul!(y, A, x)):\n", y_real)

# さらに、alpha と beta を指定することも可能 (y = alpha*A*x + beta*y)
y_initial = [10.0+1.0im, 20.0+2.0im]
y_modified = copy(y_initial) # y_initial を変更したくないのでコピー
alpha_val = 0.5 + 0.0im
beta_val = 2.0 + 0.0im

println("\n--- mul! で alpha と beta を指定する例 ---")
println("初期 y_modified:\n", y_modified)
mul!(y_modified, A_complex, x_complex, alpha_val, beta_val)
println("計算後の y_modified:\n", y_modified)

# 検証: y_modified = 0.5 * A_complex * x_complex + 2.0 * y_initial
expected_y_modified = alpha_val * (A_complex * x_complex) + beta_val * y_initial
println("検証 (0.5 * A * x + 2.0 * y_initial):\n", expected_y_modified)
println("結果の一致: ", isapprox(y_modified, expected_y_modified))

特殊な行列型 (Hermitian, Symmetric) を使用する

Juliaの LinearAlgebra モジュールは、行列がエルミート（または対称）であることを示すための特別な型を提供します。これらの型を使用すると、Juliaは行列の対称性（またはエルミート性）を利用して、さらに効率的なBLASルーチンを呼び出すことができます。

特徴

• 引数 uplo の自動処理
  BLAS.hemv() で手動で 'U' や 'L' を指定する必要がなくなります。
• 潜在的なパフォーマンス向上
  Juliaコンパイラがこれらの型を認識し、エルミート性（または対称性）を最大限に活用できるBLASルーチン（例えば hemv や symv）をより確実に選択できます。これにより、行列の片側しか読み込まないなどの最適化が自動的に行われます。
• 型安全性と明確性
  行列の特性を明示的に示すことで、コードの意図が明確になります。

コード例

using LinearAlgebra

println("\n--- 代替方法3: 特殊な行列型 (Hermitian, Symmetric) を使用する ---")

# 複素エルミート行列 (Hermitian 型でラップ)
# Hermitian(A, :U) は A の上三角部分を使ってエルミート行列を構築
# Hermitian(A, :L) は A の下三角部分を使ってエルミート行列を構築
A_herm = Hermitian([2.0+0.0im 1.0-2.0im; 1.0+2.0im 3.0+0.0im], :U)
x_complex = [1.0+0.0im, 2.0+0.0im]

println("Hermitian 型の行列 A_herm:\n", A_herm)

# Hermitian 型の行列とベクトルの積 (内部で hemv が最適化される)
result_herm = A_herm * x_complex
println("Hermitian(A) * x の結果:\n", result_herm)


# 実対称行列 (Symmetric 型でラップ)
A_symm = Symmetric([2.0 1.0; 1.0 3.0], :L)
x_real = [1.0, 2.0]

println("\nSymmetric 型の行列 A_symm:\n", A_symm)

# Symmetric 型の行列とベクトルの積 (内部で symv が最適化される)
result_symm = A_symm * x_real
println("Symmetric(A) * x の結果:\n", result_symm)

• これらの特殊な型を使うと、Juliaは内部的に LinearAlgebra.BLAS.hemv() や LinearAlgebra.BLAS.symv() を呼び出す際に、適切な引数（'U' や 'L'）を自動的に渡します。
• 同様に、Symmetric(A, :U) または Symmetric(A, :L) は実対称行列用です。
• Hermitian(A, :U) または Hermitian(A, :L) は、A がエルミート行列であることをJuliaに伝えます。:U または :L は、与えられた行列 A のどちらの三角部分に有効なデータが格納されているかを示します。

• 行列の特性を明示したい、またはさらに微細な最適化を期待したい場合
  Hermitian や Symmetric 型で行列をラップすることを検討してください。これはコードの可読性も高めます。
• メモリ割り当てを避けたい、または大規模な反復計算
  mul! 関数を使用します。これにより、GC (ガベージコレクション) のオーバーヘッドを減らし、パフォーマンスを向上させることができます。
• ほとんどの場合
  標準の行列とベクトルの乗算演算子 (*) を使用するのが最もシンプルで読みやすく、かつ効率的です。Juliaが最適なBLASルーチンを自動的に選択します。