18

Rには、qr()LINPACKまたはLAPACKのいずれかを使用してQR分解を実行する関数があります（私の経験では、後者は5％高速です）。返される主なオブジェクトは、上三角行列 R (つまりR=qr[upper.tri(qr)]) に含まれる行列 "qr" です。ここまでは順調ですね。qr の下三角部分には、Q が「コンパクトな形で」含まれています。を使用して、qr 分解から Q を抽出できqr.Q()ます。の逆数を求めたいqr.Q(). つまり、私は Q と R を持っていて、それらを「qr」オブジェクトに入れたいと考えています。R は自明ですが、Q はそうではありません。目標はそれに適用することqr.solve()であり、大規模なシステムよりもはるかに高速solve()です。

4

2 に答える 2

19

序章

R は、QR 分解を計算するために、既定ではLINPACKルーチンを使用するか、指定されている場合は LAPACKルーチンを使用します。どちらのルーチンも、ハウスホルダー反射を使用して分解を計算します。mxn 行列 A は、A = QR として mxn エコノミー サイズの直交行列 (Q) と nxn 上三角行列 (R) に分解されます。ここで、Q は、t ハウスホルダー反射行列の積によって計算できます。t は小さい方です。 m-1 と n: Q = H 1 H 2 ...H t .dqrdc DGEQP3

各反射行列 H iは、長さ (m-i+1) のベクトルで表すことができます。たとえば、H 1は、コンパクト ストレージのために長さ m のベクトルを必要とします。このベクトルの 1 つを除くすべてのエントリは、入力行列の下三角の最初の列に配置されます (対角線は R 係数によって使用されます)。したがって、各リフレクションにはストレージのスカラーがもう 1 つ必要であり、これは補助ベクトル ( $qrauxR の結果で呼び出されるqr) によって提供されます。

使用されるコンパクトな表現は、LINPACK ルーチンと LAPACK ルーチンの間で異なります。

LINPACK のやり方

ハウスホルダー反射は、H i = I - v i v i T /p iとして計算されます。ここで、I は恒等行列、p iは の対応するエントリ、$qrauxv iは次のとおりです。

  • v i [1..i-1] = 0、
  • v i [i] = p i
  • v i [i+1:m] = A[i+1..m, i] (つまり、 を呼び出した後の A の下三角の列qr)

LINPACK の例

ウィキペディアの QR 分解記事の例を R で見てみましょう。

分解される行列は

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

分解を行い、結果の最も関連性の高い部分を以下に示します。

> Aqr = qr(A)
> Aqr
$qr
            [,1]         [,2] [,3]
[1,] -14.0000000  -21.0000000   14
[2,]   0.4285714 -175.0000000   70
[3,]  -0.2857143    0.1107692  -35

[snip...]

$qraux
[1]  1.857143  1.993846 35.000000

[snip...]

この分解は、2 つの Householder 反射を計算し、それらに A を掛けて R を取得することによって (隠れて) 行われました$qr。

> p = Aqr$qraux   # for convenience
> v1 <- matrix(c(p[1], Aqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.8571429
[2,]  0.4285714
[3,] -0.2857143

> v2 <- matrix(c(0, p[2], Aqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.9938462
[3,] 0.1107692

> I = diag(3)   # identity matrix
> H1 = I - v1 %*% t(v1)/p[1]   # I - v1*v1^T/p[1]
> H2 = I - v2 %*% t(v2)/p[2]   # I - v2*v2^T/p[2]

> Q = H1 %*% H2
> Q
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

上記で計算された Q が正しいことを確認しましょう。

> qr.Q(Aqr)
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

いいね！QR が A と等しいことも確認できます。

> R = qr.R(Aqr)   # extract R from Aqr$qr
> Q %*% R
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

LAPACKのやり方

ハウスホルダー反射は、H i = I - p i v i v i Tとして計算されます。ここで、I は恒等行列、p iは の対応するエントリ、$qrauxv iは次のとおりです。

  • v i [1..i-1] = 0、
  • v i [i] = 1
  • v i [i+1:m] = A[i+1..m, i] (つまり、 を呼び出した後の A の下三角の列qr)

R で LAPACK ルーチンを使用する場合、別のひねりがあります。列のピボットが使用されるため、分解によって別の関連する問題が解決されます。AP = QR、ここで、P は順列行列です。

LAPACKの例

このセクションでは、前と同じ例を実行します。

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> Bqr = qr(A, LAPACK=TRUE)
> Bqr
$qr
            [,1]        [,2]       [,3]
[1,] 176.2554964 -71.1694118   1.668033
[2,]  -0.7348557  35.4388886  -2.180855
[3,]  -0.1056080   0.6859203 -13.728129

[snip...]

$qraux
[1] 1.289353 1.360094 0.000000

$pivot
[1] 2 3 1

attr(,"useLAPACK")
[1] TRUE

[snip...]

$pivotフィールドに注意してください。それに戻ります。ここで、情報から Q を生成しAqrます。

> p = Bqr$qraux   # for convenience
> v1 = matrix(c(1, Bqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.0000000
[2,] -0.7348557
[3,] -0.1056080


> v2 = matrix(c(0, 1, Bqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.0000000
[3,] 0.6859203


> H1 = I - p[1]*v1 %*% t(v1)   # I - p[1]*v1*v1^T
> H2 = I - p[2]*v2 %*% t(v2)   # I - p[2]*v2*v2^T
> Q = H1 %*% H2
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

ここでも、上記で計算された Q は、R が提供する Q と一致します。

> qr.Q(Bqr)
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

最後に、QR を計算しましょう。

> R = qr.R(Bqr)
> Q %*% R
     [,1] [,2] [,3]
[1,]  -51    4   12
[2,]  167  -68    6
[3,]   24  -41   -4

違いに気づきましたか？QR は、Bqr$pivot上記の順序で列が並べ替えられた A です。

于 2011-09-17T01:49:27.740 に答える
0

OPが尋ねるのと同じ問題を調査しましたが、それは可能ではないと思います。基本的に、OP の問題は、明示的に計算された Q があれば、H1 H2 ... Ht を復元できるかどうかです。QRをゼロから計算しなければこれが可能だとは思いませんが、そのような解決策があるかどうかも知りたいです.

OPと同様の問題がありますが、別のコンテキストでは、反復アルゴリズムで列や行を追加して行列 A を変更する必要があります。最初に、QR は DGEQRF を使用して計算されるため、コンパクトな LAPACK 形式になります。行列 A が新しい行などで変更された後、既存の R の最小対角のゼロ以外の要素を消滅させ、新しい R を構築するリフレクターまたはローテーターの新しいセットをすばやく構築できますが、今では H1_old のセットがありますH2_old ... Hn_old および H1_new H2_new ... Hn_new (および同様にタウ) は、単一の QR コンパクト表現に混同することはできません。私が持っている2つの可能性は、おそらくOPには同じ2つの可能性があります。

  1. 最初に計算するときでも、更新ごとに余分なフロップを犠牲にして計算するときでも、常に Q と R を明示的に分離して維持しますが、必要なメモリは十分に制限された状態に保ちます。
  2. コンパクトな LAPACK 形式に固執しますが、新しい更新が行われるたびに、これらすべての更新リフレクターのミニ セットのリストを保持します。システムを解く時点で、大きな Q'*c を実行します。つまり、H1_u3*H2_u3*...*Hn_u3*H1_u2*H2_u2*...*Hn_u2*H1_u1*H2_u1...*Hn_u1*H1*H2* です。 ...*Hn*c ここで、ui は QR 更新番号です。これは潜在的に多くの乗算を行う必要があり、メモリを追跡する必要がありますが、間違いなく最速の方法です。

David からの長い回答は基本的に、コンパクトな QR 形式とは何かを説明していますが、明示的に計算された Q と R を入力として持つこのコンパクトな QR 形式に到達する方法は説明していません。

于 2012-03-27T14:06:09.270 に答える