Main Content

スパース行列へのアクセス

非ゼロ要素

スパース行列の非ゼロ要素の上位レベル情報を取得する関数が、いくつか用意されています。

  • nnz は、スパース行列の非ゼロ要素数を返します。

  • nonzeros は、スパース行列の非ゼロ要素すべてを含む列ベクトルを返します。

  • nzmax は、スパース行列の非ゼロ要素に割り当てられたストレージの総容量を返します。

これらのいくつかを試すため、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.
MATLAB はスパース行列を圧縮されたスパース列の形式で保存するので、1 回ループするたびに 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.

そのため、関数 sparsespdiags などの構築関数を使用して、一度にすべてのスパース行列を構築するのが最善の方法です。

たとえば、スパース形式の座標行列 C が必要であると仮定します。

C=(4000104001004010101014114)

行の添字、列の添字、値の 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 でスパース行列を保存する方法の詳細については、John R. Gilbert、Cleve Moler、および Robert Schreiber の Sparse Matrices In MATLAB: Design and Implementation, (SIAM Journal on Matrix Analysis and Applications , 13:1, 333–356 (1992)) を参照してください。

スパース行列の可視化

スパース行列で非ゼロ要素分布を表示するには、グラフィカルな形式を使用すると便利な場合があります。MATLAB の関数 spy はスパース構造のテンプレート表示を作成し、グラフ上の各点は非ゼロの配列要素の位置を表します。

以下に例を示します。

Harwell-Boeing コレクションの 1 つであるスパース行列 west0479 を指定して読み込みます。

load west0479

スパース構造を表示します。

spy(west0479)

Figure contains an axes object. The axes object with xlabel nz = 1887 contains a line object which displays its values using only markers.

参考

関連するトピック