Just to avoid possible confusion, let me correct a typo
(at step [2] in the example below). Apologies!
-----FW: <XFMail.111023084327.ted.harding at wlandres.net>-----
Date: Sun, 23 Oct 2011 08:43:27 +0100 (BST)
Sender: r-help-bounces at r-project.org
From: (Ted Harding) <ted.harding at wlandres.net>
To: r-help at r-project.org
Subject: Re: [R] symmetric matrix multiplication
On 23-Oct-11 07:00:07, Daniel Nordlund wrote:>> -----Original Message-----
>> From: r-help-bounces at r-project.org
>> [mailto:r-help-bounces at r-project.org]
>> On Behalf Of statfan
>> Sent: Saturday, October 22, 2011 10:45 PM
>> To: r-help at r-project.org
>> Subject: [R] symmetric matrix multiplication
>>
>> I have a symmetric matrix B (17x17), and a (17x17) square matrix A.
>> If do
>> the following matrix multiplication I SHOULD get a symmetric matrix,
>> however
>> i don't. The computation required is:
>>
>> C = t(A)%*%B%*%A
>>
>> here are some checks for symmetry
>> > (max(abs(B - t(B))))
>> [1] 0
>> > C = t(A)%*%B%*%A
>> > (max(abs(C - t(C))))
>> [1] 3.552714e-15
>>
>> Any help on the matter would be very much appreciated.
>
> Welcome to the world of floating-point calculation on
> finite precision computers. You need to read R FAQ 7.31.
> Your maximum difference is for all intents and purposes
> equal to zero.
>
> Hope this is helpful,
> Dan
>
> Daniel Nordlund
> Bothell, WA USA
In addition to Dan's comment, let me point out that you
can convert your very nearly symmetric matrix C to an
exactly (even by R's finite-precision standards) symmetric
matrix by using (C + t(C))/2. The result will differ from
the original matrix C by similar "for all intents and purposes
zero" amounts. Here is an example, using 4x4 matrices:
##[1]: The symmetric matrix B:
B <- matrix( c(
1.1, 1.2, 1.3, 1.4,
1.2, 2.2, 2.3, 2.4,
1.3, 2.3, 3.3, 3.4,
1.4, 2.4, 3.4, 4.4), byrow=TRUE, nrow=4 )
##[2]: The non-symmetric matrix B: [OOPS! Typo!!]
##[2]: The non-symmetric matrix A:
A <- matrix( c(
1.1, 1.2, 1.3, 1.4,
2.1, 2.2, 2.3, 2.4,
3.1, 3.2, 3.3, 3.4,
4.1, 4.2, 4.3, 4.4), byrow=TRUE, nrow=4 )
##[3]: An allegedly symmetric matrix C1 (constructed
## like your C):
C1 <- t(A)%*%B%*%A
##[4]: But it isn't exactly symmetric:
max(abs(C1 - t(C1)))
# [1] 5.684342e-14
##[5]: So construct an exactly symmetric version:
C2 <- (C1 + t(C1))/2
##[6]: Check that it is exactly symmetric:
max(abs(C2 - t(C2)))
# [1] 0
##[7]: And check how close it is to the original C1:
max(abs(C2 - C1))
# 1] 5.684342e-14
Hoping this helps!
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 23-Oct-11 Time: 08:43:18
------------------------------ XFMail ------------------------------
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--------------End of forwarded message-------------------------
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 23-Oct-11 Time: 08:52:38
------------------------------ XFMail ------------------------------
