裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム，コンピュータ・サイエンス，統計学

Julia の小ネタ--025　QR 分解（結論）

2021年05月16日 | ブログラミング
結論から先に述べておく。

・ R の qr 分解には，デフォルトの LAPACK=FALSE と，LAPACK=TRUE がある。前者が unPivoted，後者が Pivoted。
・ 両者の結果はかなり違う。
・ Julia も両方に対応するそれぞれの関数が用意されている。デフォルトの Val(false) とVal(true) である。前者が Pivoted，後者が unPivoted。
・ R の qr 分解の結果を利用するいくつかの関数 qr.qy(), qr.qty() などがあるが，Julia では明示的に計算するだけである。
（実は，R のそれらの関数は裏で FORTRAN プログラムを利用しているのだが，単純な計算なのに読んでみると複雑怪奇）

特に指定がない限り unPiivoted QR decomposition を使えば良さそうだ。

以下の A, y について，QR 分解とそれに付随する計算を Juria と R で対比する。

A = [10 8 4
     10 3 3
      8 8 2
      3 5 9
      6 4 5]

y = [7.65
     5.88
     4.10
     2.36
     2.91]

R での Pivoted QR 分解は，LAPACK=TRUE を指定する。

R"QR2 = qr($A, LAPACK=TRUE)" # R での QR 分解の結果を QR2 に代入しておく
#=====
$qr
            [,1]       [,2]        [,3]
[1,] -17.5783958 -8.1349858 -12.1171466
[2,]   0.3626027  8.2959030   2.8239463
[3,]   0.2900821  0.1567791  -4.8166466
[4,]   0.1087808 -0.7920671  -0.3490988
[5,]   0.2175616 -0.2433863  -0.2604971

$rank
[1] 3

$qraux
[1] 1.568880 1.168779 1.681055

$pivot
[1] 1 3 2

attr(,"useLAPACK")
[1] TRUE
attr(,"class")
[1] "qr"
=====#

Julia での Pivoted QR 分解は，Val(true) を指定する。

QR = qr(A, Val(true));
QR
#=====
QRPivoted{Float64, Matrix{Float64}}
Q factor:
5×5 LinearAlgebra.QRPackedQ{Float64, Matrix{Float64}}:
 -0.56888   -0.0756797  -0.274156  -0.597487  -0.48836
 -0.56888   -0.196221    0.693239   0.36842   -0.146832
 -0.455104  -0.205194   -0.636312   0.505261   0.300988
 -0.170664   0.917519   -0.070799   0.26672   -0.229957
 -0.341328   0.268       0.185345  -0.425262   0.772315
R factor:
3×3 Matrix{Float64}:
 -17.5784  -8.13499  -12.1171
   0.0      8.2959     2.82395
   0.0      0.0       -4.81665
permutation:
3-element Vector{Int64}:
 1
 3
 2
=====#

Q，R 要素は個別に取り出せる。

QR.Q
#=====
5×5 LinearAlgebra.QRPackedQ{Float64, Matrix{Float64}}:
 -0.56888   -0.0756797  -0.274156  -0.597487  -0.48836
 -0.56888   -0.196221    0.693239   0.36842   -0.146832
 -0.455104  -0.205194   -0.636312   0.505261   0.300988
 -0.170664   0.917519   -0.070799   0.26672   -0.229957
 -0.341328   0.268       0.185345  -0.425262   0.772315
=====#

Q 要素は LinearAlgebra.QRPackedQ{Float64, Matrix{Float64}} なので，普通の行列にする必要がある場合には Matrix() を使う。

Matrix(QR.Q)
#=====
5×3 Matrix{Float64}:
 -0.56888   -0.0756797  -0.274156
 -0.56888   -0.196221    0.693239
 -0.455104  -0.205194   -0.636312
 -0.170664   0.917519   -0.070799
 -0.341328   0.268       0.185345
=====#

QR.R
#=====
3×3 Matrix{Float64}:
 -17.5784  -8.13499  -12.1171
   0.0      8.2959     2.82395
   0.0      0.0       -4.81665
=====#

R 要素は普通の Matrix{Float64} である。

R と Julia の結果を比較すると，R の $qr 要素の上三角行列が Julia の .R に対応している。それ以外は似ても似つかぬというところか。

QR 分解は色々な目的で使われる。R では，qr.coef(), qr.qy(),  qr.qty(), qr.resid(), qr.fitted(), qr.solve() がある。

まず，R のそれぞれの関数の結果が，Julia でどのように計算できるかを示す。

R"qr.qy(QR2, $y)"　# QR2 は，前述の R での QR 分解の結果
#=====
           [,1]
[1,] -8.7521665
[2,] -2.2212452
[3,] -5.2286747
[4,]  3.7594407
[5,]  0.9684106
=====#

Julia では以下のようにする。

QR.Q * y
#=====
5-element Vector{Float64}:
 -8.752166455226726
 -2.221245180668157
 -5.228674673737823
  3.759440684122533
  0.968410597907452
=====#

R"qr.qty(QR2, y)"
#=====
            [,1]
[1,] -10.9589067
[2,]   0.3712013
[3,]  -0.2576622
[4,]  -0.9409527
[5,]  -1.6605395
=====#
Julia では以下のようにする。 
QR.Q' * y
#=====
5-element Vector{Float64}:
 -10.958906708515896
   0.3712012834183269
  -0.2576622285703076
  -0.9409527247136643
  -1.6605395228828508
=====#

QR.Q' * y と QR.Q \ y は同じである。

QR.Q \ y
#=====
5-element Vector{Float64}:
 -10.958906708515896
   0.3712012834183275
  -0.25766222857030807
  -0.9409527247136651
  -1.660539522882851
=====#

A * coef = y を解く。

R"qr.solve($A, $y)" # QR 分解をして方程式を解くことを明示的に要求している
# [1] 0.57427561 0.05349411 0.02653560

R"qr.coef(QR2, $y)" # qr.solve(） と同じ
# [1] 0.57427561 0.05349411 0.02653560

R"lm(y ~ A + 0)$coefficients" # 定数項なしの重回帰分析と同じことである
#         A1         A2         A3
# 0.57427561 0.05349411 0.02653560

Julia では以下のようにする。 単純だA \ y # Julia が最適と判断した方法で解を求めている
#=====
3-element Vector{Float64}:
 0.5742756089206342
 0.05349411094648444
 0.026535602880579313
=====#

qr.fitted(), qr.resid() は LAPAC=TRUE の QR 分解の結果からは求められない。

R"qr.fitted(QR2, $y)" # not supported for LAPACK QR
R"qr.resid(QR2, $y)" # not supported for LAPACK QR


############################################################

R における unPivoted QR 分解は，LAPACK=FALSE を指定する（デフォルト）。

R"QR3 = qr($A, LAPACK=FALSE)" # R での QR 分解の結果を QR3 に代入しておく
#=====
RObject{VecSxp}
$qr
            [,1]         [,2]       [,3]
[1,] -17.5783958 -12.11714664 -8.1349858
[2,]   0.5688801   5.58343597  4.1958365
[3,]   0.4551041  -0.38764215  7.1566027
[4,]   0.1706640  -0.50356818 -0.8505587
[5,]   0.3413281   0.06747076 -0.3595714

$rank
[1] 3

$qraux
[1] 1.568880 1.769156 1.383743

$pivot
[1] 1 2 3

attr(,"class")
[1] "qr"
=====#

unPivoted なので，$pivot 要素は 1, 2, 3 のように昇順となる。

Julia では，unPivoted QR のために Val(false) を指定する（デフォルト）。

QR = qr(A, Val(false))
#=====
LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}
Q factor:
5×5 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:
 -0.56888    0.198228   -0.203947  -0.64671   -0.421014
 -0.56888   -0.697277    0.181347   0.350418  -0.185742
 -0.455104   0.445145   -0.498843   0.534798   0.244691
 -0.170664   0.525132    0.755705   0.24034   -0.257403
 -0.341328  -0.0243439   0.324937  -0.339413   0.813706
R factor:
3×3 Matrix{Float64}:
 -17.5784  -12.1171   -8.13499
   0.0       5.58344   4.19584
   0.0       0.0       7.1566
=====#

R での $qr 要素の上三角行列が Julia の .R に対応している

R のいくつかの関数の結果が Julia でどのように計算できるかを示す。

R"qr.qy(QR3, $y)"　# QR3 は R での QR 分解の結果
# [1] -6.7739172 -7.4219216 -0.9351772  4.6987444  0.1448075
Julia ではこうなる。
QR.Q * y
#=====
5-element Vector{Float64}:
 -6.773917239743266
 -7.421921643467008
 -0.9351771532583504
  4.698744389504694
  0.14480746343604878
=====#

R"qr.qty(QR3, y)"
# [1] -10.9589067   0.4100200   0.1899048  -1.1146946  -1.5492706

Julia ではこうなる。
QR.Q' * y
#=====
5-element Vector{Float64}:
 -10.958906708515896
   0.4100199946407921
   0.18990476781269594
  -1.1146946162626405
  -1.5492706186064318
=====#

QR.Q' * y と QR.Q \ y は同じである。

QR.Q \ y
#=====
5-element Vector{Float64}:
 -10.958906708515896
   0.41001999464079364
   0.18990476781269755
  -1.1146946162626423
  -1.5492706186064322
=====#
回帰係数を求めるとき。
A * coef = y を解くのは，LAPACK=TRUE のときと同じである。

R"qr.solve($A, $y)" # QR 分解をして方程式を解くことを明示的に要求している
# [1] 0.57427561 0.05349411 0.02653560

R"qr.coef(QR3, $y)" # qr.solve(） と同じ
# [1] 0.57427561 0.05349411 0.02653560

R"lm($y ~ $A + 0)$coefficients" # 定数項なしの重回帰分析と同じことである
#   `#JL`$A1   `#JL`$A2   `#JL`$A3
# 0.57427561 0.05349411 0.02653560
Julia での回帰係数
A \ y # Julia が最適と判断した方法で解を求めている
#=====
3-element Vector{Float64}:
 0.5742756089206342
 0.05349411094648444
 0.026535602880579313
=====#
予測値と残差の計算
R"qr.fitted(QR3, $y)" # LAPACK=TRUE の QR 分解では求まらなかった
# [1] 6.276851 5.982845 5.075229 2.229118 3.792308

R"qr.resid(QR3, $y)" # LAPACK=TRUE の QR 分解では求まらなかった
# [1]  1.3731486 -0.1028452 -0.9752290  0.1308822 -0.8823081

それぞれ，重回帰分析 lm($y ~ $A + 0) による予測値と残差である。

R"""
res = lm($y ~ $A + 0)
print(fitted(res))
print(resid(res))
"""
# fitted(res)
#        1        2        3        4        5
# 6.276851 5.982845 5.075229 2.229118 3.792308
# resid(res)
#         1          2          3          4          5
# 1.3731486 -0.1028452 -0.9752290  0.1308822 -0.8823081

Julia での予測値
A * (A \ y)
#=====
5-element Vector{Float64}:
 6.276851388300535
 5.982845230687532
 5.075228964698108
 2.229117807419539
 3.79230811171264
=====#
Julia での残差

y .- A * (A \ y)
#=====
5-element Vector{Float64}:
  1.3731486116994658
 -0.10284523068753249
 -0.9752289646981085
  0.13088219258046108
 -0.8823081117126397
 =====#

頑張ってます。クリックお願いします。

コメント    この記事についてブログを書く
 
 
 
 
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« シーザー暗号の解読--その2 | トップ | Julia に翻訳--228　2軸グラフ »
最新の画像［もっと見る］

コメントを投稿

「ブログラミング」カテゴリの最新記事

プロフィール

プロフィール画像

バックナンバー

カレンダー

2024年10月
日 月 火 水 木 金 土
1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31

カテゴリー

雨雲の動き

  • RSS2.0