Rで線形代数:極分解

Rでデータサイエンス

極分解

任意の\(n\)次元複素行列\(A\)は適当な半正定値Hermite行列\(H\)とユニタリ行列\(U\)を選んで\[A=HU\]と書ける。\(A\)が正則ならばこのような分解は一意的である。これを\(A\)の極分解という。

  • 参考引用資料: 山本哲朗(2010),『行列解析の基礎』,サイエンス社,pp.68-69

任意の正方行列は次のように対称行列と直交行列の積に分解される。それぞれを右極分解、左極分解という。\[A=QL=MQ\]ここで\[L^T=L,\,M^T=M,\,Q^TQ=QQ^T=I\]

library(dplyr)
func_matrix_decomposition_polar <- function(A) {
    # 正則行列を対象
    # rankの確認
    qr(A)$rank == unique(dim(A))
    # 右極分解
    AtA <- t(A) %*% A
    AtA
    D <- diag(sqrt(svd(AtA)$d))
    D
    U <- svd(AtA)$u
    U
    V <- svd(AtA)$v
    V
    L <- U %*% D %*% t(V)
    L
    ## 平方根の確認
    {
        {
            L %*% L
        } %>%
            round(., 10)
    } == round(AtA, 10)
    ## 対称行列の確認
    isSymmetric(L)
    # 左極分解
    AAt <- A %*% t(A)
    AAt
    D <- diag(sqrt(svd(AAt)$d))
    D
    U <- svd(AAt)$u
    U
    V <- svd(AAt)$v
    V
    M <- U %*% D %*% t(V)
    M
    ## 平方根の確認
    {
        {
            M %*% M
        } %>%
            round(., 10)
    } == round(AAt, 10)
    ## 対称行列の確認
    isSymmetric(M)
    # 直交行列 Q
    Ql <- A %*% solve(L)
    Ql
    Qr <- solve(M) %*% A
    Qr
    ## 同一の確認
    Ql == Qr
    Q <- Ql
    ## 直交の確認
    QQt <- Q %*% t(Q) %>%
        round(10)
    QtQ <- t(Q) %*% Q %>%
        round(10)
    # 分解の確認
    QL <- Q %*% L %>%
        round(10)
    MQ <- M %*% Q %>%
        round(10)
    return(list(r = qr(A)$rank, L = L, M = M, Q = Q, QQt = QQt, QtQ = QtQ, QL = QL, MQ = MQ, A = A))
}
# 例
A <- sample(0:100, 16) %>%
    matrix(nrow = 4)
func_matrix_decomposition_polar(A = A)
$r
[1] 4

$L
          [,1]     [,2]     [,3]     [,4]
[1,] 130.91353 53.12639 33.52622 34.58650
[2,]  53.12639 63.33704 27.34280 39.13281
[3,]  33.52622 27.34280 76.65692 18.19565
[4,]  34.58650 39.13281 18.19565 31.08561

$M
          [,1]     [,2]       [,3]       [,4]
[1,] 123.11441 41.29442 28.1320638 55.1017169
[2,]  41.29442 53.88216 18.6725648 52.3528268
[3,]  28.13206 18.67256 48.7096406 -0.5415076
[4,]  55.10172 52.35283 -0.5415076 76.2868867

$Q
           [,1]       [,2]        [,3]        [,4]
[1,] 0.07653764  0.7968117  0.46722571  0.37541073
[2,] 0.68667274  0.2068850  0.08996212 -0.69107594
[3,] 0.15207682 -0.5414025  0.82131017  0.09594567
[4,] 0.70674973 -0.1708011 -0.31473252  0.61014364

$QQt
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

$QtQ
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

$QL
     [,1] [,2] [,3] [,4]
[1,]   81   82   67   54
[2,]   80   25   23   12
[3,]   22    0   55    2
[4,]   94   42    6   31

$MQ
     [,1] [,2] [,3] [,4]
[1,]   81   82   67   54
[2,]   80   25   23   12
[3,]   22    0   55    2
[4,]   94   42    6   31

$A
     [,1] [,2] [,3] [,4]
[1,]   81   82   67   54
[2,]   80   25   23   12
[3,]   22    0   55    2
[4,]   94   42    6   31

参考引用資料

  1. http://www.mm.civil.tohoku.ac.jp/renzokutai/2_hizumi.pdf
  2. http://rtweb.math.kyoto-u.ac.jp/preprint/waka.pdf
  3. http://cse.naro.affrc.go.jp/takezawa/r-tips/r/20.html
  4. https://yiboyang.com/posts/2016/Apr/01/polar-decomposition-example/

最終更新

Sys.time()
[1] "2024-03-29 12:31:23 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.542'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者