スパース行列へのアクセス
非ゼロ要素
スパース行列の非ゼロ要素の上位レベル情報を取得する関数が、いくつか用意されています。
これらのいくつかを試すため、Harwell-Boeing コレクションの 1 つであるスパース行列 west0479
を読み込みます。
load west0479
whos
Name Size Bytes Class Attributes west0479 479x479 34032 double sparse
この行列は、8 段階の化学蒸留過程をモデル化したものです。
以下のコマンドを入力します。
nnz(west0479)
ans = 1887
format short e west0479
west0479 = 479×479 sparse double matrix (1887 nonzeros) (25,1) 1.0000e+00 (31,1) -3.7648e-02 (87,1) -3.4424e-01 (26,2) 1.0000e+00 (31,2) -2.4523e-02 (88,2) -3.7371e-01 (27,3) 1.0000e+00 (31,3) -3.6613e-02 . . .
nonzeros(west0479)
ans = 1.0000e+00 -3.7648e-02 -3.4424e-01 1.0000e+00 -2.4523e-02 -3.7371e-01 1.0000e+00 -3.6613e-02 -8.3694e-01 1.3000e+02 . . .
メモ
Ctrl+C を使用すると、nonzeros
の一覧表示を常時停止できます。
既定で、nnz
の出力結果は nzmax
と同じ値である点に注意してください。すなわち、非ゼロ要素数が非ゼロに割り当てられているストレージ数と等しくなっています。ただし、追加の配列要素をゼロにしても MATLAB® は動的にメモリを解放しません。いくつかの行列要素の値をゼロにすると、nnz
の値は変わりますが nzmax
の値は変わりません。
ただし、必要に応じて行列にいくつでも非ゼロ要素を追加することができます。nzmax
の元の値が制約を与えることはありません。
インデックスと値
任意の行列 (非スパースまたはスパース行列) に対して、関数 find
は非ゼロ要素のインデックスと値を返します。構文は、次のようになります。
[i,j,s] = find(S);
find
は、ベクトル i
に非ゼロ値の行インデックス、ベクトル j
に列インデックス、ベクトル s
に非ゼロ値を返します。次の例は、find
を使用してスパース行列の非ゼロ要素のインデックスと値を求めるものです。関数 sparse
は、find
出力と行列サイズを使って行列を再作成します。
S1 = west0479; [i,j,s] = find(S1); [m,n] = size(S1); S2 = sparse(i,j,s,m,n);
スパース行列の演算におけるインデックスの使用
スパース行列は圧縮されたスパース列の形式で保存されるので、スパース行列のインデックス処理にかかる時間と非スパース行列のインデックス処理にかかる時間は異なります。この処理時間は、スパース行列の要素を数個変更するだけであれば無視できるため、そのような場合は通常の配列インデックス付けを使用して値を再度割り当てることが可能です。
B = speye(4); [i,j,s] = find(B); [i,j,s]
ans = 1 1 1 2 2 1 3 3 1 4 4 1
B(3,1) = 42; [i,j,s] = find(B); [i,j,s]
ans = 1 1 1 3 1 42 2 2 1 3 3 1 4 4 1
(3,1)
に 42
を格納するため、MATLAB は非ゼロ値のベクトルと添字ベクトルに追加行を挿入して、(3,1)
の後の行列の値をすべてシフトします。線形インデックスが 2^48-1
(行列内で許可される要素数の現在の上限) を超える場合には、大規模なスパース行列要素へのアクセスや要素の代入に、線形インデックス付けを使用できません。
S = spalloc(2^30,2^30,2); S(end) = 1
Maximum variable size allowed by the program is exceeded.
関数 intmax
よりも大きな線形インデックスをもつ要素にアクセスするには、次のような配列インデックス付けを使用します。
S(2^30,2^30) = 1
S = (1073741824,1073741824) 1
スパース行列の 1 つの要素の変更だけであればそのインデックス処理にかかる時間は無視できますが、ループの状況では所要時間が長くなり、サイズが大きい行列の場合はかなり時間がかかることもあります。そのため、スパース行列の多数の要素を変更する必要がある場合、ループを使用する代わりに演算をベクトル化することを推奨します。たとえば、次のようなスパース単位行列を考えます。
n = 10000; A = 4*speye(n);
A
の要素をループ内で変更すると、同じことをベクトル化演算で行うよりも時間がかかります。tic A(1:n-1,n) = -1; A(n,1:n-1) = -1; toc
Elapsed time is 0.003344 seconds.
tic for k = 1:n-1 C(k,n) = -1; C(n,k) = -1; end toc
Elapsed time is 0.448069 seconds.
A
内の複数の要素をシフトする必要があります。スパース行列用に事前にメモリを割り当て、要素ごとにそれを満たしていく方法でも、同様にスパース配列へのインデックス付けに際してかなりのオーバーヘッドが発生します。
S1 = spalloc(1000,1000,100000); tic; for n = 1:100000 i = ceil(1000*rand(1,1)); j = ceil(1000*rand(1,1)); S1(i,j) = rand(1,1); end toc
Elapsed time is 2.577527 seconds.
インデックスと値のベクトルを構築すると、スパース配列のインデックスが不要となり、かなりの高速化が達成できます。
i = ceil(1000*rand(100000,1)); j = ceil(1000*rand(100000,1)); v = zeros(size(i)); for n = 1:100000 v(n) = rand(1,1); end tic; S2 = sparse(i,j,v,1000,1000); toc
Elapsed time is 0.017676 seconds.
そのため、関数 sparse
や spdiags
などの構築関数を使用して、一度にすべてのスパース行列を構築するのが最善の方法です。
たとえば、スパース形式の座標行列 C
が必要であると仮定します。
行の添字、列の添字、値の 3 要素の組を使用し、関数 sparse
により直接この 5 列の行列を構築します。
i = [1 5 2 5 3 5 4 5 1 2 3 4 5]'; j = [1 1 2 2 3 3 4 4 5 5 5 5 5]'; s = [4 1 4 1 4 1 4 1 -1 -1 -1 -1 4]'; C = sparse(i,j,s)
C = 5×5 sparse double matrix (13 nonzeros) (1,1) 4 (5,1) 1 (2,2) 4 (5,2) 1 (3,3) 4 (5,3) 1 (4,4) 4 (5,4) 1 (1,5) -1 (2,5) -1 (3,5) -1 (4,5) -1 (5,5) 4
スパース行列の可視化
スパース行列で非ゼロ要素分布を表示するには、グラフィカルな形式を使用すると便利な場合があります。MATLAB の関数 spy
はスパース構造のテンプレート表示を作成し、グラフ上の各点は非ゼロの配列要素の位置を表します。
以下に例を示します。
Harwell-Boeing コレクションの 1 つであるスパース行列 west0479
を指定して読み込みます。
load west0479
スパース構造を表示します。
spy(west0479)