[R] symmetric matrix multiplication

[(Ted Harding)]

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:
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 ------------------------------