Loading [MathJax]/jax/output/CommonHTML/jax.js

nAG Library for Pythonを使用した三角行列の最適化乗算

nAG Library for Python Example集

このページは、nAGライブラリのJupyterノートブックExampleの日本語翻訳版です。オリジナルのノートブックはインタラクティブに操作することができます。

nAG Library for Pythonを使用した三角行列の最適化乗算

nAGブログ記事「線形システムの解法における行列構造の活用」では、行列構造に注目することで計算時間の最適化において大きな利益が得られることを示しています。適切なタイミングで特殊なソルバーを使用することで、計算ワークフローを大幅に高速化できます。

当社のクライアントの1人が興味を持っていたケースは、両方の入力行列が三角行列である行列-行列乗算でした。標準BLASの仕様にはdtrmmと呼ばれるルーチン(nAG Library for Pythonにも実装されています)があり、入力行列の1つが三角行列である場合の最適化が含まれていますが、両方の入力行列が三角行列である状況を活用するものはありませんでした。

nAGライブラリのMark 27の2つの新しいルーチンは、両方の入力行列が上三角または下三角構造である場合の行列-行列乗算を考慮しています。2つの新しいルーチン、real_tri_matmul_inplacecomplex_tri_matmul_inplaceは、それぞれこの計算の実数版と複素数版を実装しています。

新しいルーチンをできるだけ「BLAS風」にしようと試みたので、dtrmmの使用に慣れている方なら、新しい素材の使用に問題はないでしょう。

このデモンストレーションでは、nAG Library for Pythonからこのルーチンを呼び出すことで、Intel MKLのような最適化されたBLASを使用している場合でも、デフォルトの行列-行列乗算や、より特殊なdtrmmと比較して、大規模行列の乗算がより高速になることを示します。

real_tri_matmul_inplaceの使用方法

まず、このルーチンの使用方法を簡単に示すデモンストレーションから始めましょう。これはFortranドキュメントの例に基づいており、以下の行列積を計算します。

C=αATB

ここで、AとBは両方とも上三角行列です。

from naginterfaces.library.matop import real_tri_matmul_inplace
import numpy as np

A = np.array([
    [1.5,2.3,6.7,1.9],
    [0.0,3.4,5.4,8.6],
    [0.0,0.0,8.1,2.0],
    [0.0,0.0,0.0,5.9]]
)
    
B = np.array([
    [3.5,2.1,4.0,2.1],
    [0.0,5.6,2.1,2.5],
    [0.0,0.0,1.7,1.1],
    [0.0,0.0,0.0,7.4]]
)

# B は左 'L' から操作されるのか、右 'R' から操作されるのか
side   = 'L'
# AとBは上三角行列'U'ですか、それとも下三角行列'L'ですか
uplo   = 'U' 
# 行列 A を転置しない 'N' か、転置する 'T' かを扱っているのか
transa = 'T'
# スカラー alpha の値は何ですか
alpha = 0.4
real_tri_matmul_inplace(side, uplo, transa, alpha, A, B)

この関数はインプレース関数であり、戻り値はありません。代わりに、入力引数Bがその場で変更され、結果を含むようになります。

# B には計算結果 alpha*transpose(A)*B が格納されています
B
array([[ 2.1  ,  1.26 ,  2.4  ,  1.26 ],
       [ 3.22 ,  9.548,  6.536,  5.332],
       [ 9.38 , 17.724, 20.764, 14.592],
       [ 2.66 , 20.86 , 11.624, 28.54 ]])

標準的なPythonでの同等の計算方法ですが、行列構造を全く活用しようとしていません。

A = np.array([
    [1.5,2.3,6.7,1.9],
    [0.0,3.4,5.4,8.6],
    [0.0,0.0,8.1,2.0],
    [0.0,0.0,0.0,5.9]]
)
    
# Bを再定義します。なぜならreal_tri_matmul_inplaceがBを変更したからです
B = np.array([
    [3.5,2.1,4.0,2.1],
    [0.0,5.6,2.1,2.5],
    [0.0,0.0,1.7,1.1],
    [0.0,0.0,0.0,7.4]]
)

product = alpha*(np.transpose(A) @ B)
product

array([[ 2.1 , 1.26 , 2.4 , 1.26 ], [ 3.22 , 9.548, 6.536, 5.332], [ 9.38 , 17.724, 20.764, 14.592], [ 2.66 , 20.86 , 11.624, 28.54 ]])

大規模行列に対する速度の利点

行列-行列乗算は計算科学において最も重要な演算の1つです。そのため、CPUとGPUのメーカーはこの問題の最適化のためにハードウェアの一部を専用化しています。例えば、NVIDIAはTensor Coresを開発し、比較的最近のCPU命令の進歩であるAVX (Advanced Vector Extensions)FMA (Fused Multiply-Add)により、コンピューターはこれまでにない速さで行列-行列乗算を実行できるようになりました。

その重要性から、行列-行列乗算は現代のソフトウェアライブラリで大幅に最適化されています。このデモンストレーションでは、Intel’s MKLにリンクされたPythonのNumpyのバージョンが使用されており、標準の行列乗算演算子がIntelハードウェア上で可能な限り高速に動作します。

三角行列とnAGの新しいルーチンを使用すれば、さらに良い結果を得ることができます!

# N x N サイズのランダムな上三角行列を生成する
N = 10000
# 再現性のため
np.random.seed(1)
# ランダムな上三角行列を構築する
A = np.triu(np.random.random((N, N)))
B = np.triu(np.random.random((N, N)))

# B は左 'L' から操作されるのか、右 'R' から操作されるのか
side   = 'L'
# AとBは上三角行列'U'または下三角行列'L'ですか
uplo   = 'U' 
# 行列 A を転置しない 'N' か、転置する 'T' かを扱っているのか
transa = 'N'
# スカラー alpha の値は何ですか
alpha = 1.0

まず、標準的なPython演算子を使用してこれにどのくらい時間がかかるか見てみましょう。今回は、nAGルーチンのtransaalpha引数を変更して、単純に標準的な行列-行列乗算 C=AB を実行するようにしました。

%%timeit N = 10000;np.random.seed(1);A = np.triu(np.random.random((N,N)));B = np.triu(np.random.random((N,N)))
# 上記のコード行は、timeitの各実行のためのセットアップです。これは計測されません。以下の本体のみが計測されます
C = A @ B
18.5 s ± 1.3 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

BLAS ルーチン dtrmm には、入力行列の1つが三角行列である場合の最適化が含まれています。2つ目の行列も三角行列の場合、この最適化は全く活用されません。しかし、ここではそれでも有用であり、標準的な行列-行列乗算と比較してかなりの速度向上を示しています。

# 以下のdtrmmのインポートはtimeitブロック内で使用されています:
from naginterfaces.library.blas import dtrmm # pylint: disable=unused-import 
%%timeit N = 10000;np.random.seed(1);A = np.triu(np.random.random((N,N)));B = np.triu(np.random.random((N,N)));
C = dtrmm(side='L', uplo='U', transa='N', diag='N', m=N, alpha=alpha, a=A, b=B)
12.9 s ± 1.05 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

最後に、BLASに相当するものがない新しいnAGルーチンを試してみましょう

%%timeit N = 10000;np.random.seed(1);A = np.triu(np.random.random((N,N)));B = np.triu(np.random.random((N,N)))
# 上記のコード行は、timeitの各実行のためのセットアップです。これは計測されません。以下の本体のみが計測されます
# ここでセットアップが必要なのは、real_tri_matmul_inplaceの各呼び出しが入力行列Bを変更するためです
real_tri_matmul_inplace(side, uplo, transa, alpha, A, B)
7.88 s ± 617 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

BLASとLAPACKの詳細

線形代数ルーチンを比較する際、使用されたBLASとLAPACKの実装の詳細と、いくつかのシステム詳細を提供することが非常に重要です。

# Numpyの設定詳細。IntelのPythonディストリビューションによって提供され、Intel MKLにリンクされたNumpyのバージョンを使用しました。
np.show_config()
mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\include']
blas_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\include']
blas_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\include']
lapack_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\include']
lapack_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/michael.croucher/AppData/Local/Continuum/anaconda3/envs/intel\\Library\\include']
# nAGライブラリの詳細
from naginterfaces.library.info import impl_details
impl_details()
{'title': 'Microsoft Windows x64, 64-bit, Intel C/C++ or Microsoft C/C++ or Intel Fortran',
 'integer_bits': 32,
 'precision': 'double',
 'product_code': 'NLW6I27DE',
 'mark': (27, 0, 0),
 'vendor_lib': 'MKL 2019.0.2',
 'build_id': 188890,
 'hardware': 'x86_64',
 'os': 'Windows 64-bit',
 'compilers': {'C': {'title': 'Intel(R) C++ Intel(R) 64 Compiler for Intel(R) 64, Version 19.0.2.190 Build 20190117',
   'flags': '-O3 -QaxCORE-AVX2,AVX -Qfma- -fp:precise -Qfp-speculation:safe -MD -nologo -w /EHsc -Z7'},
  'C++': {'title': 'Intel(R) C++ Intel(R) 64 Compiler for Intel(R) 64, Version 19.0.2.190 Build 20190117',
   'flags': '-O3 -QaxCORE-AVX2,AVX -Qfma- -fp:precise -Qfp-speculation:safe -MD -nologo -w /EHsc -Z7'},
  'Fortran': {'title': 'Intel(R) Visual Fortran Intel(R) 64 Compiler for Intel(R) 64, Version 19.0.2.190 Build 20190117',
   'flags': '-O3 -QaxCORE-AVX2,AVX -Qfma- -fp:precise -Qfp-speculation:safe -auto -MD -nologo -w -Z7'}}}
# オペレーティングシステムとCPUの詳細
import platform
import psutil
print(platform.architecture())
print(platform.processor())
print("Number of Physical CPUs = {0}".format(psutil.cpu_count(logical=False)))
print("Number of Logical CPUs = {0}".format(psutil.cpu_count(logical=True)))
      
('64bit', 'WindowsPE')
Intel64 Family 6 Model 142 Stepping 10, GenuineIntel
Number of Physical CPUs = 4
Number of Logical CPUs = 8
関連情報
MENU
© 譌・譛ャ繝九Η繝シ繝。繝ェ繧ォ繝ォ繧「繝ォ繧エ繝ェ繧コ繝�繧コ繧ー繝ォ繝シ繝玲�ェ蠑丈シ夂、セ 2024
Privacy Policy  /  Trademarks