2

I have recently written code for finite element method. Since my algorithm is generating a symmetric matrix, after gathering data for each element the resulting matrix should be symmetric.

However, after I run loop over element, the resulting global matrix is not symmetric. The basic code structure is like this:

A=zeros(dof,dof)

for (each element)
    loc_A = v'*(diagonal matrix)*v
         % (v is 1xN row vector)
         % loc_A is symmetric matrix

    A(K,K) = A(K,K)+loc_A
                 % (K is Nx1 column vector)
end

Since I add symmetric matrices, the resulting matrix should also be symmetric.

However, the resulting matrix is not symmetric (I checked it with issymmetric(A)). I think this is because of round-off error. If I add A=(A+A')/2, resulting matrix is symmetric (of course...). However I don't want to do such thing. Is there any other remedy which makes A symmetric naturally without any post-processing?

rayryeng
  • 102,964
  • 22
  • 184
  • 193
Dohyun
  • 642
  • 4
  • 15
  • 4
    You are correct where the addition of symmetric matrices also generates a symmetric matrix. Your theory is right but I highly suspect it is something wrong with your code. Instead of pseudo-code, please show us the actual code snippet that you're using. My gut feeling is that you're not coding something properly. – rayryeng Jun 06 '16 at 00:46
  • The real question is: how much is the matrix non-symmetric? I mean something like `(A-A.')/norm(A)` or something. If this is very small, go ahead and just symmetrize it by hand. If it's not very small, find your bug. Also, make sure your matrices are real. – Andras Deak -- Слава Україні Jun 06 '16 at 01:01
  • 4
    `A(K,K)` doesn't give you `length(K)` elements like you think it does. It give you all permutations of `K` and `K`. You want `A(sub2ind(size(A), K,K))` instead – Suever Jun 06 '16 at 01:05
  • 1
    Perhaps this could be the issue: If you are dealing with complex numbers, the transpose (') is different from (.'). The former also conjugates values. – crazyGamer Jun 06 '16 at 01:41
  • While you have the indexing bug as noted, your diagonal matrix must be a scalar. Besides that, for other cases you might want to equate your entries not average them. Use `(triu(A,1))' + triu(A)`. Hence you have a symmetric numerical error contaminating the entries. – percusse Jun 06 '16 at 01:41
  • Why don't you want to do `A = (A+A')/2`? It's a very common and perfectly acceptable thing to do. – Sam Roberts Jun 06 '16 at 10:36
  • The difference is very small (that is why i suspected this is due to the round-off error). My code works just fine for result. But if matlab recognizes my matrix as symmetric matrix, then it will use special solver which is faster than general matrix. Since my resulting matrix is quite large, I suspected that (A+A')/2 is also lots of computation. Since I'm adding many small matrix to large matrix, I wondering is there any priori processing (to small matrix) which makes resulting matrix symmetric. – Dohyun Jun 07 '16 at 01:39

1 Answers1

1

This is due to round-off error (unless you have an obvious bug in your code). This issue is normally solved by A = (A+A')/2 (this is not an expensive operation) or you can tell the solver that your global matrix is symmetric (some solvers (i.e. MUMPS) have such an option.)

berkin
  • 33
  • 5