The Gram-Schmidt process is a classic algorithm in linear algebra that takes a set of linearly independent vectors and produces an orthogonal (or orthonormal) set spanning the same subspace. It underpins QR decomposition, which is used throughout statistics, signal processing, and numerical optimization. R provides multiple ways to perform this process, from writing the loop yourself to calling built-in and package functions.

How the Algorithm Works

Given vectors v1, v2, ..., vn, the Gram-Schmidt process builds orthogonal vectors u1, u2, ..., un one at a time. The first orthogonal vector is simply u1 = v1. Each subsequent vector uk is computed by subtracting from vk its projection onto every previously computed orthogonal vector. To obtain an orthonormal basis, you divide each uk by its Euclidean norm.

In plain terms: you take the next vector, remove the parts of it that point in the directions you have already covered, and what remains is the new orthogonal direction.

Manual Implementation with Loops

Here is a straightforward implementation in R using a for loop. The input is a matrix whose columns are the vectors to orthogonalize.

1gram_schmidt <- function(V) {
2  n <- ncol(V)
3  U <- matrix(0, nrow = nrow(V), ncol = n)
4  U[, 1] <- V[, 1]
5
6  for (k in 2:n) {
7    U[, k] <- V[, k]
8    for (j in 1:(k - 1)) {
9      # Subtract the projection of V[,k] onto U[,j]
10      proj <- (sum(V[, k] * U[, j]) / sum(U[, j] * U[, j])) * U[, j]
11      U[, k] <- U[, k] - proj
12    }
13  }
14
15  # Normalize each column to get an orthonormal basis
16  Q <- apply(U, 2, function(col) col / sqrt(sum(col^2)))
17  return(Q)
18}

Numeric Example

Consider three vectors in R3.

1V <- matrix(c(
2  1, 1, 0,
3  1, 0, 1,
4  0, 1, 1
5), nrow = 3, byrow = FALSE)
6
7Q <- gram_schmidt(V)
8print(round(Q, 4))

Expected output (an orthonormal matrix):

1        [,1]    [,2]    [,3]
2[1,]  0.5774  0.7071  0.4082
3[2,]  0.5774 -0.7071  0.4082
4[3,]  0.0000  0.0000 -0.8165

You can verify orthonormality by checking that Q transposed times Q is close to the identity matrix.

1print(round(t(Q) %*% Q, 10))
2#      [,1] [,2] [,3]
3# [1,]    1    0    0
4# [2,]    0    1    0
5# [3,]    0    0    1

QR Decomposition via qr()

R has a built-in qr() function that performs QR decomposition, which internally uses a numerically stable variant of Gram-Schmidt (Householder reflections). The Q matrix from QR decomposition is exactly the orthonormal basis you would get from Gram-Schmidt.

1decomp <- qr(V)
2Q_builtin <- qr.Q(decomp)
3R_matrix  <- qr.R(decomp)
4
5print(round(Q_builtin, 4))
6print(round(R_matrix, 4))

Note that the signs of columns in qr.Q() may differ from the manual implementation. Both are valid orthonormal bases; the sign convention is arbitrary.

Using the pracma Package

The pracma package provides a dedicated gramSchmidt() function that returns both the Q and R matrices directly.

1# install.packages("pracma")  # run once if not installed
2library(pracma)
3
4result <- gramSchmidt(V)
5print(round(result$Q, 4))
6print(round(result$R, 4))

This function uses the classical Gram-Schmidt algorithm, making its output more directly comparable to the manual loop above. For large matrices where numerical stability matters, qr() with Householder reflections is generally preferred.

Common Pitfalls

• Numerical instability with classical Gram-Schmidt: The textbook algorithm accumulates rounding errors when vectors are nearly parallel. For production work, use the modified Gram-Schmidt variant or the built-in qr() function, which employs Householder reflections.
• Linearly dependent input vectors: If the input vectors are not linearly independent, one of the orthogonal vectors will become a zero vector, causing a division-by-zero when normalizing. Always verify that your matrix has full column rank before applying the process.
• Sign differences between methods: qr.Q() may flip the sign of entire columns compared to a manual implementation. Both results are correct orthonormal bases. Do not assume signs will match across different functions.
• Column vs. row orientation: R stores matrices in column-major order, and the convention throughout R's linear algebra functions is that each column is a vector. Passing a matrix where vectors are arranged as rows will produce incorrect results without an error.
• Forgetting to install pracma: The gramSchmidt() function lives in the pracma package, which is not part of base R. Calling it without library(pracma) or installing the package first produces a "could not find function" error that can be confusing for beginners.

Summary

• The Gram-Schmidt process converts linearly independent vectors into an orthogonal or orthonormal basis by iteratively subtracting projections.
• A manual R implementation uses nested for loops over the matrix columns and is useful for learning the algorithm.
• The built-in qr() function provides a numerically stable QR decomposition that yields the same orthonormal Q matrix, and it should be preferred for real-world computations.
• The pracma package offers gramSchmidt() for a direct, readable call that returns both Q and R.
• Always check that input vectors are linearly independent and be aware that sign conventions may differ between methods.