Thanks for the explanation. I am not too deep into linear algebra.

T.mat, Rz.mat and Q.mat are square matrices 4 x 4.

To be more clear, here is what they look like.

h=3.39/2;

T.mat<-matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,-h,1), nrow=4) # translation of the system

alpha<-36*pi/180

cos.alpha<-cos(alpha)

minus.sin.alpha<- -1*sin(alpha)

sin.alpha<-sin(alpha)

Rz.mat<-matrix(c(cos.alpha,minus.sin.alpha,0,0,sin.alpha,cos.alpha,0,0,0,0,1,0,0,0,0,1), nrow=4)

Rx.mat<-matrix(c(1,0,0,0,0,cos.alpha, minus.sin.alpha,0,0,sin.alpha,cos.alpha,0,0,0,0,1), nrow=4)

Q.mat is a product of Rz.mat and Rx.mat with different angles.

The resultant matrix Mn.mat is used to computed the next coordinates of a helix.

center.Of.bases<-matrix(c(0.0,0.0,0.0,1),nrow=4)

Mn.mat<-qr.solve(compute.Mn.mat(T.mat,Rz.mat,Q.mat)) # This function computes the Mn.mat matrix using specific angles. Which is used to get the new coordinates.

New.center.Of.bases<-as.numeric(Mn.mat %*% center.Of.bases)

Does these lines of code tell you anything about the complexity of the problem? Let me know if I should do anything different. I was really glad hear from you since you are an expert in the area.

Cheers../Murli

-----Original Message-----

From:

[hidden email] [mailto:

[hidden email]] On Behalf Of Douglas Bates

Sent: Wednesday, July 15, 2009 7:29 AM

To: Nair, Murlidharan T

Cc:

[hidden email]
Subject: Re: [R] Matrix multiplication precision

On Wed, Jul 15, 2009 at 3:42 AM, Nair, Murlidharan T<

[hidden email]> wrote:

> Hi!!

> I am trying to multiply 5 matrices and then using the inverse of that matrix for further computation. I read about precision problems from the archives and the suggestion was to use as.numeric while computing the products.

Hmm. I'm not sure what the origins of that advice might be but I

don't think it would apply in general.

> I am still having problems with the results. Here is how I am using it

> #Mn.mat<-(T.mat %*% Rz.mat %*% Q.mat %*% Rz.mat %*% T.mat) # I was doing this in one step which I have broken down into multiple steps as below.

> Mn1.mat<-matrix(as.numeric(T.mat %*% Rz.mat),nrow=4)

> Mn2.mat<-matrix(as.numeric(Mn1.mat %*% Q.mat),nrow=4)

> Mn3.mat<-matrix(as.numeric(Mn2.mat %*% Rz.mat),nrow=4)

> Mn.mat<-matrix(as.numeric( Mn3.mat %*% T.mat),nrow=4)

I doubt very much that the as.numeric would have any effect on

precision in these cases. You are simply taking a numeric result in

its matrix form, converting it to a vector then converting the vector

back to a matrix.,

> mm <- matrix(round(rnorm(8), 3), nrow = 4)

> mm

[,1] [,2]

[1,] -0.110 2.195

[2,] 0.616 0.353

[3,] 0.589 0.970

[4,] 1.028 0.857

> as.numeric(mm)

[1] -0.110 0.616 0.589 1.028 2.195 0.353 0.970 0.857

The key in situations like this is to analyze the structure of the

computation and decide whether you can group the intermediate results.

You have written your product as

T %*% Rz %*% Q %*% Rz %*% T

What are the characteristics of T, Rz, and Q? Are they square and

symmetric, square and triangular, diagonal? Is Q orthogonal (the

factors in an orthogonal - triangular decomposition are often written

as Q and R)? Did you happen to skip a few transposes in that formula

- often formulas like that generate symmetric matrices.

And the big lesson, of course, is the first rule of numerical linear

algebra., "Just because a formula is written in terms of the inverse

of a matrix doesn't mean that is a good way to calculate the result;

in fact, it is almost inevitably the worst way of calculating the

result". You only calculate the inverse of a matrix after you have

investigated all other possible avenues for calculating the result and

found them to be fruitless.

You haven't told us what you plan to do with the inverse and that is

an important consideration., If, for example, it represents a

variance-covariance matrix, and you want to express the result as

standard deviations and correlations you would not compute the

variance-covariance matrix but stop instead at the Cholesky factor of

the inverse.

You may want to check the facilities of the Matrix package for

expressing your calculation. It is a recommended package that is

included with binary versions of R from 2.9.0 and has many ways of

expressing and exploiting structure in matrices.

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.