qr

QR 分解

構文

R = qr(A)

[Q,R] = qr(A)

[Q,R,P] = qr(A)

[___] = qr(A,"econ")

[Q,R,P] = qr(A,outputForm)

[___] = qr(A,0)

[C,R] = qr(S,B)

[C,R,P] = qr(S,B)

[___] = qr(S,B,"econ")

[C,R,P] = qr(S,B,outputForm)

[___] = qr(S,B,0)

説明

R = qr(A) は、QR 分解 A = Q*R の上三角因子 R を返します。

[Q,R] = qr(A) は、m 行 n 列の行列 A に対して、A = Q*R を満たすように QR 分解を実行します。因子 R は m 行 n 列の上三角行列で、因子 Q は m 行 m 列の直交行列です。

[Q,R,P] = qr(A) はさらに、A*P = Q*R を満たす置換行列 P を返します。A が完全な場合、abs(diag(R)) が降順になるように置換行列が選択されます。

[___] = qr(A,"econ") は、前述の出力引数の任意の組み合わせを使用して、エコノミーサイズの分解を生成します。出力のサイズは、m 行 n 列の行列 A のサイズに依存します。

• m > n の場合、qr は Q の最初の n 列と R の最初の n 列のみを計算します。

• m <= n の場合、エコノミーサイズの分解は通常の分解と同じになります。

[Q,R,P] = qr(A,outputForm) は、置換情報 P を行列として返すかベクトルとして返すかを指定します。たとえば、outputForm が "vector" の場合は A(:,P) = Q*R となります。outputForm の既定値は "matrix" で、このとき A*P = Q*R になります。

[___] = qr(A,0) は qr(A,"econ","vector") と等価です。この構文の使用は推奨されません。代わりに "econ" オプションを使用してください。

[C,R] = qr(S,B) は、C = Q'*B と上三角因子 R を計算します。C と R を使用し、X = R\C として、スパース線形システム S*X = B の最小二乗解を計算できます。

[C,R,P] = qr(S,B) はさらに、R の非スパース要素が少なくなるように選択された置換行列 P を返します。C、R および P を使用し、X = P*(R\C) として、スパース線形システム S*X = B の最小二乗解を計算できます。

[___] = qr(S,B,"econ") は、前述の出力引数の任意の組み合わせを使用して、エコノミーサイズの分解を生成します。出力のサイズは、m 行 n 列のスパース行列 S のサイズに依存します。

• m > n の場合、qr は C と R の最初の n 行のみを計算します。

• m <= n の場合、エコノミーサイズの分解は通常の分解と同じになります。

[C,R,P] = qr(S,B,outputForm) は、置換情報 P を行列として返すかベクトルとして返すかを指定します。たとえば、outputForm が "vector" の場合、S*X = B の最小二乗解は X(P,:) = R\C になります。outputForm の既定値は "matrix" で、このとき S*X = B の最小二乗解は X = P*(R\C) になります。

[___] = qr(S,B,0) は qr(S,B,"econ","vector") と等価です。この構文の使用は推奨されません。代わりに "econ" オプションを使用してください。

例

Q-less QR 分解

5 行 5 列の魔方陣行列の QR 分解を求めます。1 つの出力引数を指定して、上三角因子のみを返します。

A = magic(5);
R = qr(A)

R = 5×5

  -32.4808  -26.6311  -21.3973  -23.7063  -25.8615
         0   19.8943   12.3234    1.9439    4.0856
         0         0  -24.3985  -11.6316   -3.7415
         0         0         0  -20.0982   -9.9739
         0         0         0         0  -16.0005

行列の完全な QR 分解

2 つの出力引数を指定して、魔方陣のテスト行列の完全な QR 分解を計算します。

A = magic(5);
[Q,R] = qr(A)

Q = 5×5

   -0.5234    0.5058    0.6735   -0.1215   -0.0441
   -0.7081   -0.6966   -0.0177    0.0815   -0.0800
   -0.1231    0.1367   -0.3558   -0.6307   -0.6646
   -0.3079    0.1911   -0.4122   -0.4247    0.7200
   -0.3387    0.4514   -0.4996    0.6328   -0.1774

R = 5×5

  -32.4808  -26.6311  -21.3973  -23.7063  -25.8615
         0   19.8943   12.3234    1.9439    4.0856
         0         0  -24.3985  -11.6316   -3.7415
         0         0         0  -20.0982   -9.9739
         0         0         0         0  -16.0005

マシンの精度内で $\mathit{A}=\mathrm{QR}$ であることを検証します。

norm(A-Q*R)

ans = 
9.5562e-15

置換された QR 分解

3 つの出力引数を指定し、QR 分解の R 因子の非ゼロ要素を減らす置換行列または置換ベクトルを返します。

スパース行列 west0479 の QR 分解を計算します。3 つの出力を指定し、$\mathrm{AP}=\mathrm{QR}$ を満たす置換行列を返します。

load west0479
A = west0479;
[Q,R,P] = qr(A);

置換行列 P について、マシンの精度内で A*P = Q*R であることを検証します。

norm(A*P-Q*R,"fro")

ans = 
3.5451e-10

次に、"vector" オプションを指定して p を置換ベクトルとして返します。

[Q,R,p] = qr(A,"vector");

置換行列 p について、マシンの精度内で A(:,p) = Q*R であることを検証します。

norm(A(:,p) - Q*R,"fro")

ans = 
3.5451e-10

分解に置換行列または置換ベクトルを使用すると、置換されていない分解に比べて、スパース入力に対する R 因子の非ゼロ要素が少なくなります。

[Q1,R1] = qr(A);
spy(R1)

spy(R)

この結果は、置換された分解により R 因子の非ゼロ要素が大幅に少なくなることを示しています。

エコノミーサイズの QR 因子を使用した線形システムの求解

係数行列のエコノミーサイズ QR 分解を使用して、線形システム $\mathrm{Ax}=\mathit{b}$ を解きます。

magic(10) の最初の 5 列を使用して、10 行 5 列の係数行列を作成します。線形方程式 $\mathrm{Ax}=\mathit{b}$ の右辺には、行列の行の和を使用します。この設定では、方程式 $\mathit{x}$ の解は 1 のベクトルになるはずです。

A = magic(10);
A = A(:,1:5)

A = 10×5

    92    99     1     8    15
    98    80     7    14    16
     4    81    88    20    22
    85    87    19    21     3
    86    93    25     2     9
    17    24    76    83    90
    23     5    82    89    91
    79     6    13    95    97
    10    12    94    96    78
    11    18   100    77    84

b = sum(A,2)

b = 10×1

   215
   215
   215
   215
   215
   290
   290
   290
   290
   290

A のエコノミーサイズの QR 分解を計算します。次に、x(p,:) = R\(Q'*b) で線形システム $\mathrm{QRx}=\mathit{b}$ を解きます。Q は直交であるため、この方程式は x(p,:) = R\(Q\b) と同じですが、計算が速くなります。

[Q,R,p] = qr(A,"econ","vector")

Q = 10×5

   -0.0050   -0.4775   -0.0504    0.5193    0.0399
   -0.0349   -0.5001   -0.0990   -0.1954   -0.2006
   -0.4384    0.1059   -0.4660    0.4464    0.0628
   -0.0947   -0.4151   -0.2923   -0.2542    0.5274
   -0.1246   -0.4117   -0.2812   -0.1326   -0.4130
   -0.3787    0.0209    0.2702    0.4697    0.0390
   -0.4085   -0.0017    0.2217   -0.2450   -0.2015
   -0.0648   -0.3925    0.6939    0.0669    0.1225
   -0.4683    0.0833    0.0283   -0.3038    0.5265
   -0.4982    0.0867    0.0394   -0.1822   -0.4138

R = 5×5

 -200.7112  -55.5026 -167.6040  -84.7237 -168.7997
         0 -192.1053  -40.3557 -152.4040  -39.2814
         0         0  101.3180  -89.4254   96.0172
         0         0         0   41.0248  -14.9083
         0         0         0         0   24.6386

p = 1×5

     3     1     5     2     4

x(p,:) = R\(Q'*b)

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

R の対角の片対数プロットを作成し、置換された分解によって、abs(diag(R)) が減少する R 因子が生成されることを確認します。比較のため、同じプロット内に A の特異値をプロットします。実際には、R の対角値は A の特異値と似た振る舞いをします。そのため、R の対角値は、行列 A がどの程度特異行列に近いかを判断する尺度として使用できます。

semilogy(abs(diag(R)),"-o")
hold on
semilogy(svd(A),"r-o")
legend("Diagonal of R","Singular Values of A")

スパース線形システムの求解

スパース線形システムを解き、その結果から S の列空間内にあるベクトル b の量を確認します。

10% の密度と 1 のベクトルを使用して、ランダムな 500 行 20 列のスパース行列を作成します。qr を使用して、行列を因子 R と C = Q'*b に因数分解します。

S = sprand(500,20,0.1);
b = ones(500,1);
[C,R] = qr(S,b,"econ");

結果を使用して、x = R\C で $\mathrm{Sx}=\mathit{b}$ を解きます。

x = R\C;

恒等式 ${\Vert \mathit{b}\Vert }^2={\Vert \mathrm{Sx}-\mathit{b}\Vert }^2+{\Vert \mathit{C}\Vert }^2$ について考えます。

両辺を b のノルムで除算すると、S の列空間内にある b の量を示す新しい恒等式が得られます。

$\frac{{\Vert \mathrm{Sx}-\mathit{b}\Vert }^2}{{\Vert \mathit{b}\Vert }^2}+\frac{{\Vert \mathit{C}\Vert }^2}{{\Vert \mathit{b}\Vert }^2}=1$.

1 番目の項は S の列空間内に "ない" b の量を示し、2 番目の項は S の列空間内に "ある" b の量を示します。

t1 = norm(S*x-b)^2/norm(b)^2

t1 = 
0.4000

t2 = norm(C)^2/norm(b)^2

t2 = 
0.6000

方形スパース線形システムの求解

qr を使用して、行列方程式 $\mathrm{Sx}=\mathit{B}$ を方形スパース係数行列 $\mathit{S}$ で解きます。

スパース行列 west0479 を読み込み、最初の 200 列を線形システムの方形係数行列として使用します。方程式の右辺に $\mathit{S}$ の行の和を使用します。この設定では、$\mathrm{Sx}=\mathit{B}$ の解は 1 のベクトルになります。

load west0479
S = west0479(:,1:200);
B = sum(S,2);

2 つの入力と 3 つの出力をもつ qr を使用して、$\mathrm{Sx}=\mathit{B}$ を解きます。線形システムの解は x = P*(R\C) になります。

[C,R,P] = qr(S,B);
x = P*(R\C);

マシンの精度内で $\mathrm{Sx}-\mathit{B}=0$ であることを検証します。

norm(S*x-B)

ans = 
8.3874e-11

メモ: 上三角因子 R と置換行列 P を計算するが、直交行列 Q (多くの場合、qr の呼び出しで最も計算量が多い部分) の計算は行わないというときには、B を空行列として指定します。

emptyB = zeros(size(S,1),0);
[~,R,P] = qr(S,emptyB);

入力引数

A — 入力行列

行列

入力行列。完全な行列またはスパース行列として指定します。

データ型: single | double
複素数のサポート: あり

S — 入力係数行列

行列

入力係数行列。スパース行列として指定します。2 つの入力行列がある場合、qr は線形システム S*X = B の最小二乗解を計算します。

データ型: single | double
複素数のサポート: あり

B — 右辺の行列

行列

右辺の行列。完全な行列またはスパース行列として指定します。2 つの入力行列がある場合、qr は C = Q'*B を計算します。これを使用して、線形システム S*X = B を解くことができます。

データ型: single | double
複素数のサポート: あり

outputForm — 置換出力の形状

"matrix" (既定値) | "vector"

置換出力の形状。"matrix" または "vector" として指定します。このフラグは、置換出力 P が置換行列または置換ベクトルのいずれで返されるかを制御します。このオプションを使用するには、qr に 3 つの出力引数を指定しなければなりません。

• outputForm が "vector" の場合、P は A(:,P) = Q*R を満たす置換ベクトルです。

• outputForm の既定値は "matrix" で、このとき A*P = Q*R になります。

例: [Q,R,P] = qr(A,"vector")

出力引数

Q — 直交因子

行列

直交因子。m 行 n 列の行列 A に対して A = Q*R を満たす行列として返されます。

• 完全な分解の場合、qr(A) は Q を $$Q^{H}Q=Q{Q}^H=I_m$$ を満たす m 行 m 列の直交行列として返します。

• m > n の方形行列 A の場合、エコノミーサイズの分解 qr(A,"econ") は、Q の最初の n 列と R の最初の n 行のみを計算します。Q の列は、A の列空間の正規直交基底を形成します。

マシンや MATLAB^® のリリースが異なる場合、異なる Q の列が出力されることがありますが、数値的にはいずれも正確です。Q と R の対応する行と列で符号が反転する場合がありますが、これは、反転しても式 A = Q*R の値には影響がないためです。

R — 上三角因子

行列

上三角因子。A = Q*R を満たす行列として返されます。A が非スパースであり、[Q,R,P] = qr(A) として 3 つの出力が指定された場合、R の対角は降順になります。

P — 置換情報

行列 | ベクトル

置換情報。行列またはベクトルとして返されます。P の形状は outputForm の値に依存します。また、qr は、最初の入力行列が完全とスパースのどちらであるかによっても、選択する P が満たす基準が異なります。

• 非スパース — qr は、abs(diag(R)) が降順になるように P を選択します。

• スパース — qr は、R 内の非ゼロ要素を減らすように P を選択します。

C — 線形システムの因子

行列

線形システムの因子。C = Q'*B を満たす行列として返されます。S*X = B の最小二乗解は、X = R\C です。置換出力 P が指定されている場合、解は、outputForm の値に応じて、X = P*(R\C) または X(P,:) = R\C のどちらかになります。

• outputForm が "vector" の場合、S*X = B の最小二乗解は、X(P,:) = R\C です。

• outputForm の既定値は "matrix" で、このとき S*X = B の最小二乗解は X = P*(R\C) になります。

ヒント

• 同じ係数行列を含む複数の線形システムを解くには、decomposition オブジェクトを使用します。

• 構文 [C,R] = qr(S,B) では、X = R\C の値は、S に低ランクの行列がない場合のみ S*X = B の最小二乗解になります。

参照

[1] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[2] Davis, Timothy A. “Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing Sparse QR Factorization.” ACM Transactions on Mathematical Software 38, no. 1 (November 2011): 1–22. https://doi.org/10.1145/2049662.2049670.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• コード生成では、この関数のスパース行列入力はサポートされません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意および制限:

• コード生成では、この関数のスパース行列入力はサポートされません。

• コード生成で MATLAB と異なる QR 分解が返されることがあります。方程式 Q*R = A を使用して Q 値と R 値を検証できます。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

qr 関数は GPU 配列入力をサポートしますが、次の使用上の注意および制限があります。

• スパース入力はサポートされません。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

分散配列

Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

使用上の注意および制限:

• 入力 A は非スパースでなければなりません。

詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

R2025a: 単精度のスパース入力引数の指定

入力係数行列 S は単精度として指定できます。S が single の場合、順序付けや置換ベクトルなどのインデックス付けに関連するものを除き、出力引数も single になります。

R2022b: エコノミーサイズの分解に "0" 構文は非推奨

構文 [___] = qr(A,0) と構文 [___] = qr(S,B,0) は引き続きサポートされますが、推奨されなくなりました。代わりに "econ" オプションを使用して、エコノミーサイズの分解を実行してください。

R2022a: エコノミーサイズの分解のオプション

"econ" オプションを使用して、qr でエコノミーサイズの分解を計算します。機能は、3 番目の入力が指定されない限り qr(A,0) と同じです。

R2022a: 1 出力の構文は常に上三角因子を返す

構文 R = qr(A) は、A が完全かスパースかにかかわらず、常に R を上三角行列として返します。従来は、完全な A に対し、1 出力の構文はハウスホルダー ベクトルが行列の下三角部分にある R 行列を返していました。それらのベクトルは Q 行列を部分的に定義していました。

参考

lu | chol | null | orth | qrdelete | qrinsert | qrupdate | decomposition | lsqminnorm | rank

トピック

• 分解
• 連立線形方程式