SVD、QR、コレスキー分解
svd()、qr()、chol() を使って、次元削減と因子分解を行います。
「SVD、QR、コレスキー分解」はCoddyKit上の無料R Academyレッスンです。 これはレッスン4/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
なぜ行列分解を使うのか
行列分解では、A をより単純な行列の積に分解します。これにより、隠れた構造（ランクや条件の良さ）が明らかになり、効率的な計算が可能になります。また、統計学や機械学習における PCA、回帰、最適化の基盤となります。
# Three fundamental decompositions:
# 1. SVD: A = U D V' (any matrix)
# 2. QR: A = Q R (any matrix)
# 3. Cholesky: A = L L' (symmetric positive definite)
# Each serves a purpose:
# SVD -> PCA, pseudoinverse, rank, image compression
# QR -> linear regression, Gram-Schmidt, eigenvalues
# Cholesky -> fast solve for SPD systems, simulation
A <- matrix(c(4, 3, 2,
3, 6, 1,
2, 1, 5), nrow = 3, byrow = TRUE)
cat('Matrix ready for decomposition')
print(A)SVD：svd() 関数
特異値分解（SVD）は、任意の m×n 行列を A = U D V' と分解します。U と V は直交行列で、D は減少する順に並んだ非負の特異値を対角成分に持つ対角行列です。
A <- matrix(c(1, 2, 3,
4, 5, 6), nrow = 2, byrow = TRUE)
# Compute SVD
svd_result <- svd(A)
# Components:
svd_result$d # singular values (decreasing)
svd_result$u # left singular vectors (2x2)
svd_result$v # right singular vectors (3x2)
cat('Singular values:', svd_result$d, '\n')
cat('Rank of A:', sum(svd_result$d > 1e-10), '\n')SVD による再構成
SVD から A = U %*% diag(d) %*% t(V) として A を再構成できます。低ランク近似では上位 k 個の特異値と特異ベクトルだけを残します。これは画像圧縮やノイズ除去の基盤となる方法です。
A <- matrix(c(1, 2, 3,
4, 5, 6,
7, 8, 9), nrow = 3, byrow = TRUE)
res <- svd(A)
U <- res$u; d <- res$d; V <- res$v
# Full reconstruction: A = U diag(d) V'
A_reconstructed <- U %*% diag(d) %*% t(V)
all.equal(A, A_reconstructed) # TRUE
# Rank-1 approximation (best rank-1 approx by Eckart-Young)
A_rank1 <- d[1] * U[, 1] %*% t(V[, 1])
print(round(A_rank1, 2))
# Frobenius error of rank-1 approx
sqrt(sum((A - A_rank1)^2))SVD：擬似逆行列とランク
SVD により、擬似逆行列 A⁺ = V D⁺ U' を求められます。D⁺ では、0 でない各特異値をその逆数に置き換えます。これにより、A が正方行列でない場合やフルランクでない場合でも、最小二乗の意味で Ax = b を解けます。
A <- matrix(c(1, 2, 3,
4, 5, 6), nrow = 2, byrow = TRUE)
# Pseudoinverse via SVD
svd_res <- svd(A)
tol <- 1e-10
# Invert nonzero singular values
d_inv <- ifelse(svd_res$d > tol, 1/svd_res$d, 0)
A_pinv <- svd_res$v %*% diag(d_inv) %*% t(svd_res$u)
print(A_pinv) # 3x2 pseudoinverse
# Least-squares solution: x = A+ b
b <- c(7, 8)
x_ls <- A_pinv %*% b
print(x_ls)
# Compare with MASS::ginv
# MASS::ginv(A) should give the same resultQR 分解：qr()
QR 分解では、A = QR と分解します。Q は直交行列（Q'Q = I）、R は上三角行列です。線形回帰、固有値の計算（QR アルゴリズム）、グラム・シュミット直交化などに使用されます。
A <- matrix(c(1, 2,
3, 4,
5, 6), nrow = 3, byrow = TRUE)
# Compute QR decomposition
qr_res <- qr(A)
# Extract Q and R
Q <- qr.Q(qr_res) # 3x2 orthogonal
R <- qr.R(qr_res) # 2x2 upper triangular
cat('Q is orthogonal: Q\'Q = I\n')
print(round(t(Q) %*% Q, 10)) # Identity
cat('R is upper triangular:\n')
print(round(R, 6))
# Reconstruct A = Q %*% R
all.equal(A, Q %*% R) # TRUE回帰における QR 分解
QR 分解を使って線形回帰の正規方程式を解くと、数値的に優れた結果が得られます。Q'Q = I なので、正規方程式 A'Ax = A'b は R'Q'QRx = R'Q'b となり、Rx = Q'b に簡略化できます。これは後退代入で解きます。
# Generate regression data
set.seed(7)
n <- 20
x1 <- rnorm(n); x2 <- rnorm(n)
y <- 1 + 2*x1 - 0.5*x2 + rnorm(n, sd = 0.5)
# Design matrix
X <- cbind(intercept = 1, x1 = x1, x2 = x2)
# QR solve: numerically stable
beta_qr <- qr.solve(X, y)
print(beta_qr) # ~c(1, 2, -0.5)
# Compare with lm()
beta_lm <- coef(lm(y ~ x1 + x2))
all.equal(beta_qr, beta_lm, check.names = FALSE) # TRUE
# The qr.solve internals use backsolve on R
cat('QR regression is lm()\'s default method')コレスキー分解：chol()
コレスキー分解では、対称正定値行列を A = L L' と分解します。L は下三角行列です。R の chol() は上三角因子 U を返すため、A = t(U) %*% U となります。
# Symmetric positive definite matrix
A <- matrix(c(4, 2,
2, 3), nrow = 2, byrow = TRUE)
# Cholesky decomposition: A = t(U) %*% U
U <- chol(A)
print(U)
# Note: chol() returns UPPER triangular in R
# Reconstruct
all.equal(A, t(U) %*% U) # TRUE
# Attempting Cholesky on non-PD matrix throws error
# Catch it:
result <- tryCatch(
chol(matrix(c(1, 3, 3, 2), nrow = 2)),
error = function(e) cat('Error:', e$message)
)コレスキー分解による解法
対称正定値システム（共分散行列など）では、コレスキー分解による解法は LU 分解の 2 倍高速です。逆行列を計算するには chol2inv() を使うか、コレスキー因子と backsolve/forwardsolve を使用します。
# SPD system: typical in statistics (covariance matrices)
set.seed(1)
X <- matrix(rnorm(30), nrow = 10, ncol = 3)
A <- crossprod(X) # X'X is always SPD
b <- c(1, 2, 3)
# Cholesky factor
U <- chol(A) # Upper triangular
# Solve A x = b using two triangular solves
# A = t(U) %*% U -> solve t(U) y = b, then U x = y
y <- forwardsolve(t(U), b)
x <- backsolve(U, y)
print(x)
# Verify
all.equal(as.vector(A %*% x), b)
# Or: inverse via chol2inv
A_inv <- chol2inv(U)
all.equal(x, as.vector(A_inv %*% b))PCA における SVD
SVD による PCA では、まずデータ行列 X を中心化し、中心化した X の SVD を計算します。右特異ベクトル V は主成分の方向を表し、特異値 d は説明される分散（d²/(n-1)）と関係します。
set.seed(42)
n <- 50; p <- 3
X <- matrix(rnorm(n * p), nrow = n)
# Add correlation
X[, 2] <- 0.8 * X[, 1] + 0.6 * X[, 2]
# Center
X_c <- scale(X, center = TRUE, scale = FALSE)
# SVD of centered data
svd_res <- svd(X_c)
# Variance explained by each PC
var_explained <- svd_res$d^2 / (n - 1)
prop_var <- var_explained / sum(var_explained)
cat('Proportion of variance explained:\n')
print(round(prop_var, 3))
# Compare with prcomp
pc <- prcomp(X)
all.equal(prop_var, pc$sdev^2 / sum(pc$sdev^2))SVD による行列ランク
行列のランクは、無視できない特異値の数に等しくなります。SVD はランクを判定するうえで最も数値的に信頼できる方法です。d の値のうち、しきい値（通常は max(d) * tol）を超えるものを数えるだけで求められます。
# Full-rank matrix (rank 3)
A <- matrix(c(1, 0, 0,
0, 2, 0,
0, 0, 3), nrow = 3)
svd(A)$d # c(3, 2, 1) - all nonzero -> rank 3
# Rank-deficient matrix (rank 2)
B <- matrix(c(1, 2, 3,
2, 4, 6, # row 2 = 2 * row 1
0, 1, 1), nrow = 3, byrow = TRUE)
svd(B)$d # Third value is ~0
# Rank function using SVD
matrix_rank <- function(M, tol = 1e-10) {
d <- svd(M)$d
sum(d > tol * d[1]) # relative tolerance
}
matrix_rank(B) # 2多変量正規分布のサンプリングにおけるコレスキー分解
多変量正規分布 N(μ, Σ) からサンプリングするには、まず Σ のコレスキー因子 U を計算し、標準正規乱数を X = μ + Z %*% U と変換します。MASS::mvrnorm() は内部でこの方法を使用しています。
set.seed(99)
# Covariance matrix
Sigma <- matrix(c(4, 2,
2, 3), nrow = 2)
mu <- c(1, 5)
n_samples <- 500
# Cholesky factor
U <- chol(Sigma)
# Sample standard normals
Z <- matrix(rnorm(n_samples * 2), nrow = n_samples)
# Transform: X ~ N(mu, Sigma)
X <- sweep(Z %*% U, 2, mu, FUN = '+')
# Check empirical covariance
print(round(cov(X), 2)) # Should be close to Sigma
print(round(colMeans(X), 2)) # Should be close to mu理解度チェック
R における SVD、QR 分解、コレスキー分解について理解度を確認しましょう。
まとめ：SVD、QR 分解、コレスキー分解
重要なポイント： SVD（svd()）は d、u、v を返し、PCA、ランクの計算、擬似逆行列、低ランク近似に使用されます。QR 分解（qr()）は直交行列 Q と上三角行列 R を返し、回帰や数値安定性の向上に使用されます。コレスキー分解（chol()）は A = t(U)%*%U を満たす上三角行列 U を返し、SPD システムに最も高速です。これらはそれぞれ異なる計算目的に適しており、A を正確に再構成できます。
A <- matrix(c(4, 2, 2, 3), nrow = 2)
# SVD
res <- svd(A)
all.equal(A, res$u %*% diag(res$d) %*% t(res$v)) # TRUE
# QR (on a design matrix)
X <- cbind(1, c(1,2,3,4,5))
qr_res <- qr(X)
all.equal(X, qr.Q(qr_res) %*% qr.R(qr_res)) # TRUE
# Cholesky
U <- chol(A)
all.equal(A, t(U) %*% U) # TRUE
# Use cases:
# svd() -> PCA, rank, pseudoinverse
# qr() -> regression, stability
# chol()-> SPD solve, MVN samplingよくある質問
「SVD、QR、コレスキー分解」レッスンは無料ですか？
はい。「SVD、QR、コレスキー分解」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「SVD、QR、コレスキー分解」で何を学びますか？
svd()、qr()、chol() を使って、次元削減と因子分解を行います。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか？
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン4/4です。
「SVD、QR、コレスキー分解」レッスンにはどのくらい時間がかかりますか？
ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか？
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。
このコースのすべてのレッスン
- 行列の乗算と行列式
- solve()による連立一次方程式の解法
- 固有値と固有ベクトル
- SVD、QR、コレスキー分解