ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

svds

特異値とベクトルのサブセット

構文

  • s = svds(A)
  • s = svds(A,k)
  • s = svds(A,k,sigma)
  • s = svds(A,k,sigma,opts)
  • s = svds(Afun,n,___)
  • [U,S,V] = svds(___)
  • [U,S,V,flag] = svds(___)

説明

s = svds(A) は、行列 A の最も大きい 6 個の特異値からなるベクトルを返します。

s = svds(A,k) は最も大きい k 個の特異値を返します。

s = svds(A,k,sigma) は、sigma の値に基づいて k 個の特異値を返します。たとえば、svds(A,k,'smallest') は、最も小さい k 個の特異値を返します。

s = svds(A,k,sigma,opts) はさらに、構造体を使用してオプションを指定します。

s = svds(Afun,n,___) は、行列 A ではなく関数ハンドル Afun を指定します。オプションで、ksigma または opts を追加の入力引数として指定できます。

[U,S,V] = svds(___) は、左特異ベクトル U、特異値の対角行列 S および右特異ベクトル V を返します。前述の構文にある任意の入力引数の組み合わせが使用できます。

[U,S,V,flag] = svds(___) は、収束フラグも返します。flag0 の場合、すべての特異値が収束しています。

すべて折りたたむ

スパース行列の最も大きい 6 個の特異値を計算します。

A = delsq(numgrid('C',15));
s = svds(A)
s =

    7.8666
    7.7324
    7.6531
    7.5213
    7.4480
    7.3517

2 番目の入力に、計算する最も大きい特異値の数を指定します。

s = svds(A,3)
s =

    7.8666
    7.7324
    7.6531

スパース行列の最も小さい 5 個の特異値を計算します。

A = delsq(numgrid('C',15));
s = svds(A,5,'smallest')
s =

    0.5520
    0.4787
    0.3469
    0.2676
    0.1334

スパースの 100 行 100 列のノイマン行列を作成します。

C = gallery('neumann',100);

最も小さい 10 個の特異値を計算します。

ss = svds(C,10,'smallest')
ss =

    0.9828
    0.9049
    0.5625
    0.5625
    0.4541
    0.4506
    0.2256
    0.1139
    0.1139
         0

最も小さい 10 個の非ゼロの特異値を計算します。行列にはゼロに等しい特異値があるため、'smallestnz' オプションによりその値が省略されます。

snz = svds(C,10,'smallestnz')
snz =

    0.9828
    0.9828
    0.9049
    0.5625
    0.5625
    0.4541
    0.4506
    0.2256
    0.1139
    0.1139

スパース行列の右上および左下の非ゼロのブロックを表す 2 つの行列を作成します。

n = 500;
B = rand(500);
C = rand(500);

svds で使用できるように、Afun を現在のディレクトリに保存します。

function y = Afun(x,tflag,B,C,n)
if strcmp(tflag,'notransp')
    y = [B*x(n+1:end); C*x(1:n)];
else
    y = [C'*x(n+1:end); B'*x(1:n)];
end

関数 AfunB および C を使用して、実際にスパース行列 A = [zeros(n) B; C zeros(n)] の全体を作成することなく、指定されたフラグに応じて A*x または A'*x を計算します。これは行列のスパース パターンを利用して、A*x および A'*x を計算するときにメモリを節約します。

Afun を使用して、A の最も大きい 10 個の特異値を計算します。BC および n を追加の入力として Afun に渡します。

s = svds(@(x,tflag) Afun(x,tflag,B,C,n),[1000 1000],10)
s =

  250.3248
  249.9914
   12.7627
   12.7232
   12.6988
   12.6608
   12.6166
   12.5643
   12.5419
   12.4512

A の最も大きい 10 個の特異値を直接計算して、結果と比較します。

A = [zeros(n) B; C zeros(n)];
s = svds(A,10)
s =

  250.3248
  249.9914
   12.7627
   12.7232
   12.6988
   12.6608
   12.6166
   12.5643
   12.5419
   12.4512

west0479 は、479 行 479 列の実数値のスパース行列です。この行列には、いくつかの大きい特異値と、多数の小さい特異値があります。

west0479 を読み込んで、A として格納します。

load west0479
A = west0479;

A の特異値分解を計算し、最も大きい 6 個の特異値および対応する特異ベクトルを返します。特異値の収束を確認するために、4 番目の出力引数を指定します。

[U,S,V,cflag] = svds(A);
cflag
cflag =

     0

cflag は、すべての特異値が収束したことを示します。特異値は、出力行列 S の対角成分です。

s = diag(S)
s =

   1.0e+05 *

    3.1895
    3.1725
    3.1695
    3.1685
    3.1669
    0.3038

A の完全な特異値分解を計算して、結果を確認します。A を非スパース行列に変換して、svd を使用します。

[U1,S1,V1] = svd(full(A));

対数スケールを使用して、A のすべての特異値をプロットします。

semilogy(diag(S1),'r.')
title('Singular Values of west0479')

スパース対角行列を作成して、最も大きい 6 個の特異値を計算します。

A = diag(sparse([1e4*ones(1, 8) 1e4:-1:1]));
s = svds(A)
Warning: Maximum number of iterations reached. Results may be inaccurate. 

s =

   1.0e+04 *

    1.0000
    0.9999
    0.9998
    0.9997
    0.9996
    0.9995

svds によって警告が出力されます。これは最大回数の反復を実行したが、許容誤差に収まらなかったためです。

収束の問題を解決する最も効果的な方法は、p により大きい値を使用することにより、計算に使用する Krylov 部分空間の最大サイズを増加することです。そのためには、p = 50 を指定する options 構造体を svds に渡します。

opts.p = 50;
s = svds(A,6,'largest',opts)
s =

   1.0e+04 *

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

近特異行列の最も小さい 10 個の特異値を計算します。

rng default
format shortg
B = spdiags([repelem([1; 1e-7], [198, 2]) ones(200, 1)], [0 1], 200, 200);
s1 = svds(B,10,'smallest')
Warning: Large residual norm detected. This is likely due to bad condition of
the input matrix (condition number 1.0008e+16). 

s1 =

       7.0945
       7.0945
       7.0945
       7.0945
       7.0945
       7.0945
       7.0945
       7.0945
      0.25927
   7.0888e-16

警告は、svds が適切な特異値の計算に失敗したことを示します。svds が失敗した原因は、最小と最小から 2 番目の特異値の間にギャップがあるためです。svds(...,'smallest')B の逆行列を計算する必要があり、これにより大きい数値誤差が発生しています。

比較のために、svd を使用して、正確な特異値を計算します。

s = svd(full(B));
s = s(end-9:end)
s =

      0.14196
      0.12621
      0.11045
     0.094686
     0.078914
     0.063137
     0.047356
     0.031572
     0.015787
   7.0888e-16

svds でこの計算を再現するために、B の QR 分解を行います。三角行列 R の特異値は、B と同じになります。

[Q,R,p] = qr(B,0);

R の各行のノルムをプロットします。

rownormR = sqrt(diag(R*R'));
semilogy(rownormR)
hold on;
semilogy(size(R, 1), rownormR(end), 'ro')

R の最後の要素はほぼゼロであるため、解が不安定になります。

この要素による悪い影響が解の正しい部分に及ばないように、R の最後の行を厳密にゼロに設定します。

R(end,:) = 0;

svds を使用して R の最も小さい 10 個の特異値を見つけます。結果は svd で得られたものと同等です。

sr = svds(R,10,'smallest')
sr =

      0.14196
      0.12621
      0.11045
     0.094686
     0.078914
     0.063137
     0.047356
     0.031572
     0.015787
            0

この方法を使用して B の特異ベクトルを計算するには、Q および置換ベクトル p を使用して左右の特異ベクトルを変換します。

[U,S,V] = svds(R,20,'s');
U = Q*U;
V(p,:) = V;

入力引数

すべて折りたたむ

入力行列。A は多くの場合、大規模なスパース行列です。

データ型: double
複素数のサポート: はい

計算する特異値の数。正のスカラー整数として指定します。次のいずれかの条件が満たされる場合、svds は要求された数よりも少ない特異値を返します。

  • kmin(m,n) より大きい

  • sigma = 'smallestnz'、かつ kA の非ゼロの特異値の数より大きい

k が大きすぎる場合、svds はその値を有効な k の最大値に置き換えます。

例: svds(A,2) は、A の最も大きい 2 つの特異値を返します。

特異値のタイプ。次の値のいずれかとして指定します。

オプション説明

'largest'

(既定)
最も大きい特異値

'smallest'

最も小さい特異値

'smallestnz

最も小さい非ゼロの特異値

スカラー

スカラー シフトに最も近い特異値

例: svds(A,k,'smallest') は、最も小さい k 個の特異値を出力します。

例: svds(A,k,100) は、100 に最も近い k 個の特異値を計算します。

options 構造体。次の表の 1 つ以上のフィールドを含む構造体として指定します。

オプション フィールド説明既定の設定
opts.tol

収束の許容誤差

1e-10
opts.maxit

最大反復回数

100
opts.p

Krylov 部分空間の最大サイズ

max(3*K,15)
opts.u0

左初期開始ベクトル

最大で opts.u0 または opts.v0 のいずれかを指定できます。いずれのオプションも指定しない場合、m 行 n 列の行列 A の既定は次のとおりです。

  • m < nopts.u0 = randn(m,1);

  • m >= nopts.v0 = randn(n,1);

    メモ:   svds は、プライベート乱数ストリームを使用する再現可能な方法で u0 および v0 を選択します。乱数のシードを変更しても、この randn の使用方法には影響しません。

opts.v0

右初期開始ベクトル

    メモ:   数値スカラー シフト sigma を使用する場合、svds はオプション p を無視します。

例: opts.tol = 1e-6, opts.maxit = 500 は、フィールド tol および maxit に値が設定された構造体を作成します。

データ型: struct

行列関数。関数ハンドルとして指定します。関数 Afun は、次の条件を満たさなければなりません。

  • Afun(x,'notransp') はベクトル x を受け入れて、積 A*x を返す。

  • Afun(x,'transp') はベクトル x を受け入れて、積 A'*x を返す。

    メモ:   関数ハンドルは、sigma = 'largest' (既定) の場合にのみ使用します。

例: svds(Afun,[1000 1200])

データ型: function_handle

Afun で使用される行列 A のサイズ。2 要素のサイズ ベクトル [m n] として指定します。

例: svds(Afun,[1000 1200])

出力引数

すべて折りたたむ

特異値。列ベクトルとして返されます。特異値は、非負の実数値であり、降順でリストされます。

左特異ベクトル。行列の列として返されます。Amn 列の行列で、k 個の特異値を要求した場合、U は直交列をもつ mk 列の行列になります。

特異値。対角行列として返されます。S の対角要素は非負の特異値です。Amn 列の行列で、k 個の特異値を要求した場合、Skk 列になります。

右特異ベクトル。行列の列として返されます。Amn 列の行列で、k 個の特異値を要求した場合、V は直交列をもつ nk 列の行列になります。

収束フラグ。スカラーとして返されます。値 0 は、すべての特異値が収束したことを示します。それ以外の場合、一部の特異値が収束していません。

詳細

すべて折りたたむ

ヒント

  • svds は、複数の実行にわたって再現性を確保するために、プライベート乱数ストリームを使用して既定の開始ベクトル opts.u0 および opts.v0 を生成します。svds を呼び出す前に rng を使用して乱数発生器の状態を設定しても、出力に影響しません。

  • svds を使用する方法は、小規模で密度の高い行列の特異値をいくつか求めるときに最も効率的な方法とは言えません。その問題の処理に十分なメモリがあるときには、svd(full(A)) を使用する方が迅速に処理できる場合があります。たとえば、500 行 500 列の行列の特異値を 3 個求める問題は、比較的小規模の問題であるため svd で容易に処理できます。

  • 指定された行列で svds が収束しない場合は、opts.p の値を増加することにより Krylov 部分空間のサイズを増加します。補助的な選択肢として、反復の最大回数 opts.maxit と収束の許容誤差 opts.tol を調整することも、収束の処理に役立つ場合があります。

  • 高速パフォーマンスを実現するには、k を増加すると効果が得られることがあり、特に特異値が繰り返される場合に有効です。

参照

[1] Baglama, J. and L. Reichel, “Augmented Implicitly Restarted Lanczos Bidiagonalization Methods.” SIAM Journal on Scientific Computing. Vol. 27, 2005, pp. 19–42.

[2] Larsen, R. M. “Lanczos Bidiagonalization with partial reorthogonalization.” Dept. of Computer Science, Aarhus University. DAIMI PB-357, 1998.

参考

|

R2006a より前に導入

この情報は役に立ちましたか?