solve()による連立一次方程式の解法
Ax = b の連立方程式を解き、行列の逆行列を計算します。
「solve()による連立一次方程式の解法」はCoddyKit上の無料R Academyレッスンです。 これはレッスン2/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
線形システム：Ax = b
連立一次方程式は Ax = b と表せます。ここで A は係数行列、x は未知ベクトル、b は右辺ベクトルです。x を解析的に解くとは、x = A⁻¹b を計算することを意味します。
# System of equations:
# 2x + y = 5
# x + 3y = 7
# Matrix form: A %*% x = b
A <- matrix(c(2, 1,
1, 3), nrow = 2, byrow = TRUE)
b <- c(5, 7)
# What are A and b?
print(A)
print(b)
cat('We want to find x such that A %*% x = b')solve(A, b)：直接解法
solve(A, b) は Ax = b を x について解きます。内部では LU 分解を使用します。これは、A⁻¹ を明示的に計算してから b を乗算するよりも、数値的に安定で効率的です。
A <- matrix(c(2, 1,
1, 3), nrow = 2, byrow = TRUE)
b <- c(5, 7)
# Solve Ax = b
x <- solve(A, b)
print(x) # x[1] = ?, x[2] = ?
# Verify: A %*% x should equal b
residual <- A %*% x - b
print(residual) # Should be near zero
# Manual check:
# 2*(8/5) + (9/5) = 16/5 + 9/5 = 25/5 = 5 ✓
# 1*(8/5) + 3*(9/5) = 8/5 + 27/5 = 35/5 = 7 ✓solve(A)：逆行列
引数を 1 つだけ指定して solve(A) を呼び出すと、A の逆行列 A⁻¹（A %*% A⁻¹ = I を満たす行列）が返されます。Ax=b を解く場合には使用せず、直接 solve(A,b) を使ってください（より高速で安定しています）。
A <- matrix(c(4, 3,
3, 2), nrow = 2, byrow = TRUE)
# Compute inverse
A_inv <- solve(A)
print(A_inv)
# Verify: A %*% A_inv = I
A %*% A_inv # Should be identity matrix
round(A %*% A_inv, 10)
# Also A_inv %*% A = I
round(A_inv %*% A, 10)
# det(A) != 0 required for invertibility
det(A) # -1 (nonzero, so invertible)解の確認
必ず A %*% x - b を計算して解を検証してください。浮動小数点演算のため残差は厳密には 0 になりませんが、マシンイプシロン（~1e-15）に近い値になるはずです。残差の大きさを 1 つの値で確認するには norm() を使います。
A <- matrix(c(3, -1, 2,
1, 4, 0,
-2, 1, 5), nrow = 3, byrow = TRUE)
b <- c(1, 2, 3)
# Solve
x <- solve(A, b)
cat('Solution x:\n'); print(x)
# Residual check
residual <- A %*% x - b
cat('Residual vector:\n'); print(residual)
# Residual norm (should be near 0)
resid_norm <- sqrt(sum(residual^2))
cat('Residual norm:', resid_norm, '\n')
# Expected: something like 2e-16条件数：kappa()
A の条件数は、b の摂動に対して解がどの程度敏感かを表します。条件数が大きいと、b の小さな誤差が x の大きな誤差につながります。このようなシステムは悪条件です。
# Well-conditioned matrix
A_good <- matrix(c(2, 1, 1, 3), nrow = 2, byrow = TRUE)
kappa(A_good) # Small -> good
# Ill-conditioned (nearly singular) matrix
A_bad <- matrix(c(1.000, 1.001,
1.001, 1.002), nrow = 2, byrow = TRUE)
kappa(A_bad) # Very large -> bad!
# Rule of thumb: kappa > 1/machine_epsilon is trouble
.Machine$double.eps # ~2.2e-16
# For A_bad: you lose about log10(kappa) digits of precision
cat('Digits lost:', log10(kappa(A_bad)), '\n')上三角行列に対する backsolve()
backsolve(R, b) は、R が上三角行列の場合に、後退代入を使って Rx = b を解きます。三角行列のシステムでは一般的な solve() よりもはるかに高速で、計算量は O(n²) です（solve() は O(n³)）。
# Upper triangular system: Rx = b
# 2x + 3y + z = 14
# 5y + 2z = 13
# 4z = 8
R <- matrix(c(2, 3, 1,
0, 5, 2,
0, 0, 4), nrow = 3, byrow = TRUE)
b <- c(14, 13, 8)
# Solve using back-substitution
x <- backsolve(R, b)
print(x) # z=2, y=(13-4)/5=9/5, x=(14-3*9/5-2)/2
# Verify
all.equal(as.vector(R %*% x), b) # TRUE
# Compare with general solve
x_general <- solve(R, b)
all.equal(x, x_general) # TRUE (same result)下三角行列に対する forwardsolve()
forwardsolve(L, b) は、L が下三角行列の場合に、前進代入を使って Lx = b を解きます。backsolve() と相補的な関係にあり、両者を組み合わせることで LU 分解によるソルバーの基盤となります。
# Lower triangular system: Lx = b
# 3x = 6
# 2x + 4y = 10
# x + 2y + 5z = 16
L <- matrix(c(3, 0, 0,
2, 4, 0,
1, 2, 5), nrow = 3, byrow = TRUE)
b <- c(6, 10, 16)
# Solve using forward-substitution
x <- forwardsolve(L, b)
print(x) # x=2, y=(10-4)/4=1.5, z=(16-2-3)/5=2.2
# Verify
all.equal(as.vector(L %*% x), b) # TRUE
# Use case: solving L*U*x = b
# forwardsolve(L, b) -> y, then backsolve(U, y) -> x複数の右辺
B が行列である solve(A, B) は、B のすべての列について AX = B を同時に解きます。各列に対して solve(A, b) を個別に呼び出すより効率的です。
A <- matrix(c(2, 1,
1, 3), nrow = 2, byrow = TRUE)
# Solve for two right-hand sides simultaneously
B <- matrix(c(5, 7, # first system
3, 1), # second system
nrow = 2, byrow = TRUE)
# X[:,1] solves Ax = B[:,1]
# X[:,2] solves Ax = B[:,2]
X <- solve(A, B)
print(X)
# Verify both solutions
A %*% X # Should equal B
all.equal(A %*% X, B) # TRUE特異行列の検出
特異行列に対して solve(A) を呼び出すとエラーが発生します。解く前に det(A) または rcond(A)（条件数の逆数）を確認してください。堅牢なコードにするには tryCatch() を使います。
# Singular matrix (rows are linearly dependent)
S <- matrix(c(1, 2,
2, 4), nrow = 2, byrow = TRUE)
det(S) # 0 -> singular
kappa(S) # Inf
# Safe solve with tryCatch
safe_solve <- function(A, b) {
tryCatch(
solve(A, b),
error = function(e) {
cat('Matrix is singular or nearly so:\n')
cat(e$message, '\n')
return(NULL)
}
)
}
result <- safe_solve(S, c(1, 2))
print(result) # NULLsolve() による最小二乗法
過剰決定系（未知数より方程式の数が多い場合）には、厳密な解が存在しません。最小二乗解は ||Ax - b||² を最小化します。正規方程式 A'Ax = A'b を解くことで求められます。
# Overdetermined: 4 equations, 2 unknowns (y = a + b*x)
set.seed(1)
x_vals <- c(1, 2, 3, 4)
y_vals <- c(2.1, 4.0, 5.9, 8.2) # approx y = 0 + 2x
A <- cbind(1, x_vals) # Design matrix (4x2)
b <- y_vals
# Normal equations: (A'A) beta = A'b
AtA <- crossprod(A)
Atb <- crossprod(A, b)
beta_ols <- solve(AtA, Atb)
print(beta_ols) # Intercept ~0.1, slope ~2.0
# Residual sum of squares
y_hat <- A %*% beta_ols
rss <- sum((b - y_hat)^2)
cat('RSS:', rss)安定性のための qr.solve()
悪条件または過剰決定のシステムでは、qr.solve(A, b) は solve() より数値的に安定です。LU 分解の代わりに QR 分解を使用します。lm() は内部でこの方法を使用します。
# For overdetermined system, qr.solve is preferred
x_vals <- c(1, 2, 3, 4, 5)
y_vals <- c(1.9, 4.1, 6.0, 7.8, 10.1)
A <- cbind(1, x_vals)
b <- y_vals
# qr.solve handles overdetermined systems directly
beta_qr <- qr.solve(A, b)
print(beta_qr) # intercept, slope
# Equivalent to:
beta_lm <- coef(lm(y_vals ~ x_vals))
all.equal(beta_qr, beta_lm, check.names = FALSE) # TRUE
# For well-determined square systems, solve() is fine
# For overdetermined or ill-conditioned: use qr.solve()理解度チェック
R における線形システムの解法について理解度を確認しましょう。
まとめ：線形システムの解法
重要なポイント： solve(A, b) は Ax=b を直接解くため、推奨されます。solve(A) は A⁻¹ を計算するため、解法には使用しないでください。backsolve(R,b) と forwardsolve(L,b) は三角行列のシステムに適しています。kappa(A) は条件の良さを測定し、値が大きいほど解が不安定です。過剰決定系には qr.solve() を使います。必ず A %*% x - b で解を検証してください。
# Summary of solve() functions:
A <- matrix(c(3, 1, 1, 2), nrow = 2)
b <- c(9, 8)
# Solve Ax = b
x <- solve(A, b)
print(x) # c(2, 3)
# Check conditioning
kappa(A) # Small -> well-conditioned
# Verify
max(abs(A %*% x - b)) # Near zero
# For triangular systems:
R <- matrix(c(2, 3, 0, 4), nrow = 2, byrow = TRUE)
backsolve(R, c(8, 4)) # x=c(1, 1)
# For least squares (overdetermined):
# qr.solve(design_matrix, y)よくある質問
「solve()による連立一次方程式の解法」レッスンは無料ですか？
はい。「solve()による連立一次方程式の解法」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「solve()による連立一次方程式の解法」で何を学びますか？
Ax = b の連立方程式を解き、行列の逆行列を計算します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか？
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン2/4です。
「solve()による連立一次方程式の解法」レッスンにはどのくらい時間がかかりますか？
ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか？
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。
このコースのすべてのレッスン
- 行列の乗算と行列式
- solve()による連立一次方程式の解法
- 固有値と固有ベクトル
- SVD、QR、コレスキー分解