Menyelesaikan Sistem Linear dengan solve()

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 <- 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

# 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')

# 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)

# 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

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

# 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)  # NULL

# 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)

# 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()

# 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)