0Pricing
R Academy · Lesson

Solving Linear Systems with solve()

Find solutions to Ax = b systems and compute matrix inverses.

Solving Linear Systems with solve() is a free R Academy lesson on CoddyKit — lesson 2 of 4. You can read the complete lesson below for free — then practise it hands-on in the browser with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, one of 4 lessons in the course, and your progress syncs across the web and the CoddyKit app.

Linear Systems: Ax = b

A system of linear equations can be written as Ax = b, where A is a matrix of coefficients, x is the unknown vector, and b is the right-hand side. Solving for x analytically means computing 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): The Direct Solution

solve(A, b) solves Ax = b for x. It uses LU decomposition internally, which is more numerically stable and efficient than explicitly computing A⁻¹ and then multiplying by 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): The Matrix Inverse

Calling solve(A) with only one argument returns the inverse of A: A⁻¹ such that A %*% A⁻¹ = I. Avoid using this for solving Ax=b — use solve(A,b) directly (faster and more stable).

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)

Checking Your Solution

Always verify the solution by computing A %*% x - b. Due to floating-point arithmetic, the residual won't be exactly zero, but should be near machine epsilon (~1e-15). Use norm() for a single residual magnitude.

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

Condition Number: kappa()

The condition number of A measures how sensitive the solution is to perturbations in b. A large condition number means small errors in b cause large errors in x — the system is ill-conditioned.

# 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() for Upper Triangular

backsolve(R, b) solves Rx = b where R is upper triangular using back-substitution. Much faster than general solve() for triangular systems — O(n²) vs 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() for Lower Triangular

forwardsolve(L, b) solves Lx = b where L is lower triangular using forward-substitution. Complements backsolve() and together they form the basis of LU decomposition solvers.

# 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

Multiple Right-Hand Sides

solve(A, B) where B is a matrix solves AX = B simultaneously for all columns of B. This is more efficient than calling solve(A, b) separately for each column.

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

Detecting Singular Matrices

Calling solve(A) on a singular matrix throws an error. Check det(A) or rcond(A) (reciprocal condition number) before solving. Use tryCatch() for robust code.

# 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

Least Squares with solve()

For overdetermined systems (more equations than unknowns), there's no exact solution. The least-squares solution minimizes ||Ax - b||². It solves the normal equations: 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)

Using qr.solve() for Stability

qr.solve(A, b) is numerically more stable than solve() for ill-conditioned or overdetermined systems. It uses QR decomposition instead of LU. lm() uses this internally.

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

Quick Check

Test your understanding of linear system solving in R.

Recap: Solving Linear Systems

Key takeaways: solve(A, b) solves Ax=b directly (preferred). solve(A) computes A⁻¹ (avoid for solving). backsolve(R,b) and forwardsolve(L,b) are fast for triangular systems. kappa(A) measures conditioning — large values mean unstable solutions. For overdetermined systems use qr.solve(). Always verify with 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)

Frequently asked questions

Is the “Solving Linear Systems with solve()” lesson free?

Yes — the full text of “Solving Linear Systems with solve()” is free to read here on the web, and the R Academy course includes 4 lessons in total. To practise it interactively (a built-in code editor and a 24/7 AI tutor) and unlock the rest of the R Academy course, upgrade to CoddyKit PRO.

What will I learn in “Solving Linear Systems with solve()”?

Find solutions to Ax = b systems and compute matrix inverses. You practise R Academy with hands-on code you run directly in the browser, and a 24/7 AI tutor answers your questions as you work through the lesson.

Do I need any experience to start R Academy?

No prior experience is required. R Academy on CoddyKit is structured for beginners through advanced learners; this is — lesson 2 of 4, so you can start here or from the beginning and move at your own pace.

How long does the “Solving Linear Systems with solve()” lesson take?

Most CoddyKit lessons take about 5–10 minutes. Each one is bite-sized and interactive, so you make steady progress and pick up exactly where you left off across the web and the app.

Can I write and run code in this R Academy lesson?

Yes. Every R Academy lesson includes a built-in code editor, so you write and run real code right in your browser and get instant AI feedback — no local setup required.

All lessons in this course

  1. Matrix Multiplication and Determinants
  2. Solving Linear Systems with solve()
  3. Eigenvalues and Eigenvectors
  4. SVD, QR, and Cholesky Decompositions
← Back to R Academy