[18874] in s-news-athena
Re: [S] finding vectors orthogonal to a given set? (fwd)
daemon@ATHENA.MIT.EDU (Nicholas Lange)
Thu Aug 5 11:54:00 1999
Date: Thu, 5 Aug 1999 11:48:55 -0400 (EDT)
From: Nicholas Lange <nick@mclean.harvard.edu>
To: Ranjan Maitra <maitra@math.umbc.edu>
Cc: S-news <s-news@wubios.wustl.edu>
In-Reply-To: <Pine.LNX.4.04.9908051134260.4654-100000@param.math.umbc.edu>
Message-Id: <Pine.LNX.3.96.990805114338.18646L-100000@eos.mclean.org>
Mime-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII
Hello Ranjan,
All good suggestions, and I'd add that Gram-Schimdt is of course
equivalent to performing separate linear regressions successively on the
columns, so that lsfit()or lm() could easily do the job, i.e., replace the
jth column with estimated residuals obtained from its OLS regression on
preceeding columns as j = 1,...,k.
Best regards,
Nick
-------------------------------------------------------------------------
Nicholas Lange Statistical Neuroimaging Laboratory
and Laboratory for Molecular Pharmacology
Associate Professor Mailman Research Center, McLean Hospital
Department of Psychiatry 115 Mill Street, Belmont, MA 02178
Harvard Medical School
Chief Biostatistician voice: 617-855-2139 fax: 617-855-3479
McLean Hospital e-mail: nick@mclean.harvard.edu
-------------------------------------------------------------------------
On Thu, 5 Aug 1999, Ranjan Maitra wrote:
>
> 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
>
-----------------------------------------------------------------------
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
