[18873] in s-news-athena
Re: [S] finding vectors orthogonal to a given set? (fwd)
daemon@ATHENA.MIT.EDU (Ranjan Maitra)
Thu Aug 5 11:39:48 1999
Date: Thu, 5 Aug 1999 11:35:01 -0400 (EDT)
From: Ranjan Maitra <maitra@math.umbc.edu>
To: S-news <s-news@wubios.wustl.edu>
Message-Id: <Pine.LNX.4.04.9908051134260.4654-100000@param.math.umbc.edu>
Mime-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII
Here is another suggestion I received: thanks, Rolf!
Ranjan
---------- Forwarded message ----------
Date: Thu, 5 Aug 1999 12:22:38 -0300 (ADT)
From: Rolf Turner <rolf@math.unb.ca>
To: maitra@math.umbc.edu
Subject: Re: [S] finding vectors orthogonal to a given set?
Let U be the matrix whose columns are your U_1, ..., U_k.
Assign
x <- cbind(U,diag(nrow(U))
i.e. extend your orthogonal set of vectors by a basis so that you
have a set which definitely spans R^n. Then apply Gramm-Schmidt
to the columns of x, discarding any columns which are linearly
dependent on previous columns.
The function given below will do the Gramm-Schmidt procedure for
you. Note that it is written ***naively*** --- I am no numerical
analyst --- and so it is NOT numerically sound. Perhaps someone
else will point you at a way that is numerically sound. Anyway,
this function works on toy examples.
Note that the result is orthoNORMAL, not just orthogonal.
cheers,
Rolf Turner
rolf@math.unb.ca
===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
gsl <- function(x) {
#
# Function gsl. To calculate the Gramm-Schmidt orthonormalization
# of a maximal linearly independent set of columns of a matrix x.
# Columns which are linearly dependent upon previous columns are
# discarded. This function returns the matrix of orthonormalized
# columns.
#
# If the columns of x are linearly dependent, an error is given.
#
eps <- 2*.Machine$double.eps
n <- ncol(x)
nr <- nrow(x)
rslt <- list()
k <- 0
for(j in 1:n) {
cc <- x[, j]
tmp <- if(k > 0) {
apply(matrix(unlist(lapply(rslt[1:k],
function(x, v) { sum(v * x) * x },
v = cc)), ncol = k),1,sum)
}
else 0
cc <- cc - tmp
norm <- sqrt(sum(cc * cc))
if(norm > eps) {
k <- k + 1
rslt[[k]] <- cc/norm
if(k==nr) break
}
}
matrix(unlist(rslt), ncol = k)
}
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
