-
Notifications
You must be signed in to change notification settings - Fork 15
/
GramSchmidt.R
83 lines (83 loc) · 2.81 KB
/
GramSchmidt.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#' Gram-Schmidt Orthogonalization of a Matrix
#'
#' Carries out simple Gram-Schmidt orthogonalization of a matrix.
#' Treating the columns of the matrix \code{X} in the given order,
#' each successive column after the first is made orthogonal to all
#' previous columns by subtracting their projections on the current
#' column.
#'
#' @param X a matrix
#' @param normalize logical; should the resulting columns be normalized to unit length? The default is \code{TRUE}
#' @param verbose logical; if \code{TRUE}, print intermediate steps. The default is \code{FALSE}
#' @param tol the tolerance for detecting linear dependencies in the columns of a. The default is \code{sqrt(.Machine$double.eps)}
#' @param omit_zero_columns if \code{TRUE} (the default), remove linearly dependent columns from the result
#' @return A matrix of the same size as \code{X}, with orthogonal columns (but with 0 columns removed by default)
#' @author Phil Chalmers, John Fox
#' @export
#' @examples
#' (xx <- matrix(c( 1:3, 3:1, 1, 0, -2), 3, 3))
#' crossprod(xx)
#' (zz <- GramSchmidt(xx, normalize=FALSE))
#' zapsmall(crossprod(zz))
#'
#' # normalized
#' (zz <- GramSchmidt(xx))
#' zapsmall(crossprod(zz))
#'
#' # print steps
#' GramSchmidt(xx, verbose=TRUE)
#'
#' # A non-invertible matrix; hence, it is of deficient rank
#' (xx <- matrix(c( 1:3, 3:1, 1, 0, -1), 3, 3))
#' R(xx)
#' crossprod(xx)
#' # GramSchmidt finds an orthonormal basis
#' (zz <- GramSchmidt(xx))
#' zapsmall(crossprod(zz))
#'
GramSchmidt <- function (X, normalize = TRUE, verbose = FALSE,
tol=sqrt(.Machine$double.eps), omit_zero_columns=TRUE) {
B <- X
B[, 2L:ncol(B)] <- 0
if (verbose) {
cat("\nInitial matrix:\n")
printMatrix(X)
cat("\nColumn 1: z1 = x1\n")
printMatrix(B)
}
for (i in 2:ncol(X)) {
res <- X[, i]
for (j in 1:(i - 1)) res <- res - Proj(X[, i], B[, j])
if (len(res) < tol) res <- 0
B[, i] <- res
if (verbose) {
prj_char <- character(i - 1)
for (j in 1:(i - 1)) prj_char[j] <- sprintf("Proj(x%i, z%i)",
i, j)
cat(sprintf("\nColumn %i: z%i = x%i - %s\n", i, i,
i, paste0(prj_char, collapse = " - ")))
printMatrix(B)
}
}
zeros <- apply(B, 2, function(b) all(b == 0))
if (omit_zero_columns && any(zeros)){
B <- B[, !zeros]
which_zeros <- (1:ncol(X))[zeros]
message("linearly dependent column", if (length(which_zeros) > 1) "s" else "",
" omitted: ", paste(which_zeros, collapse=", "))
}
if (normalize) {
norm <- diag(1/len(B))
colnames <- colnames(B)
B <- B %*% norm
colnames(B) <- colnames
if (!omit_zero_columns && any(zeros)) B[, zeros] <- 0
if (verbose) {
cat("\nNormalized matrix: Z * inv(L) \n")
printMatrix(B)
}
}
if (verbose)
return(invisible(B))
B
}